Electro-fluid-mechanics of the heart

Abstract

This article presents an overview of the dynamics of the human heart and the main goal is the discussion of its fluid mechanic features. We will see, however, that the flow in the heart can not be fully described without considering its electrophysiology and elastomechanics as well as the interaction with the systemic and pulmonary circulations with which it is strongly connected. Biologically, the human heart is similar to that of all warm-blooded mammals and it satisfies the same allometric laws. Since the Paleolithic Age, however, humans have improved their living conditions, have modified the environment to satisfy their needs and, more recently, have developed advanced medical knowledge which has allowed triple the number of heartbeats with respect to other mammals. In the last century, effective diagnostic tools, reliable surgical procedures and prosthetic devices have been developed and refined leading to substantial progress in cardiology and heart surgery with routine clinical practice which nowadays cures many disorders, once lethal. Pulse duplicators have been built to reproduce the pulsatile flow and ‘blood analogues’, have been realized. Heart phantoms, can attain deformations similar to the real heart although the active contraction and the tissue anisotropy still can not be replicated. Numerical models have also become a viable alternative for cardiovascular research: they do not suffer from limitations of material properties and device technologies, thus making possible the realization of truly digital twins. Unfortunately, a high-fidelity model for the whole heart consists of a system of coupled, nonlinear partial differential equations with a number of degrees of freedom of the order of a billion and computational costs become the bottleneck. An additional challenge comes from the inherent human variability and the uncertainty of the heart parameters whose statistical assessment requires a campaign of simulations rather than a single deterministic calculation; reduced and surrogate models can be employed to alleviate the huge computational burden and all possibilities are currently being pursued. In the era of big data and artificial intelligence, cardiovascular research is also advancing by exploiting the latest technologies: equation-based augmented reality, virtual surgery and computational prediction of disease progression are just a few examples among many that will become standard practice in the forthcoming years.

1. Introduction

Most of us worry about the heart only when something goes wrong and fortunately that happens very seldom. If we only reflected on some figures, however, heart reliability would look truly exceptional: with a weight of $0.3$ kg (less than half a percent of the total body mass) and a power consumption of approximately $8$ W, it pumps blood at a rate of $5\ {\rm l}\ \min ^{-1}$ through the circulatory system, whose cumulative length exceeds $10^{5}$ km, lifelong. The heart performs this task by contracting $40$ million times per year at a rate of $70$ beats per minute (b.p.m.) although, depending on environmental conditions or intense physical exercise, it self-adapts to yield up to a six-fold flow rate, with respect to resting conditions, without us being conscious of it.

Indeed, humans have always been aware of the heart presence owing to its continuous rhythmic beat and perhaps they have also guessed its importance for life looking at the immobile hearts of hunted prey. The understanding of the exact role of this organ in body physiology, however, has taken a long time and independent rediscoveries.

Apparently, already in $2600$ BC, the Chinese Yellow Emperor, Nei Ching, in the Canon of Medicine wrote: ‘The blood current flows continuously in a circle without a beginning or end and never stops’ and that ‘all the blood is under control of the heart’ (Mackay & Mensah 2004).

Unfortunately, this brilliant intuition, incredibly ahead of its time, remained confined within China for two millennia, before Aristotle (in the fourth century BC) conceived the heart as ‘the source of all movement, since the heart links the soul with the organs of life’ which can be considered correct only philosophically.

Galen (in the second century AD) advanced the knowledge by discovering the pulmonary circulation although he considered the liver as the main circulation organ and identified the heart as the source of the ‘innate heat’ which was distributed throughout the body by blood (May 1968); this misconception was the rational for bleeding as a treatment for fever.

As centuries progressed, evidence accumulated against Aristotle's view although his writings dominated all fields of knowledge and the ipse dixit statement (he said it himself) prevented any further discussion.

The age of Renaissance revived the empirical observation and, through dissection of corpses and animals, a good understanding of the heart structure was achieved; the drawings by Leonardo da Vinci clearly show how detailed the knowledge of chambers, valves and connected veins or arteries was (see figure 1) although the notion of the heart physiology still did not differ much from Galen's ideas.

Figure 1. Drawings by Leonardo da Vinci (about 1513) of details of an ox's heart: the anatomy of valves, heart chambers and main vessels does not differ significantly from the present knowledge (adapted from Roberts & Keele 1984).

Finally, in the 17th century, William Harvey rediscovered the function of the heart when he neatly stated that ‘It has been shown by reason and experiment that blood by the beats of the ventricles flows through the lungs and heart and is pumped to the whole body … the blood in the animal body moves around in a circle continuously and … the action or function of the heart is to accomplish this pumping. This is the only reason for the motion and beat of the heart’ and also ‘The heart's one role is the transmission of the blood and its propulsion, by means of the arteries, to the extremities everywhere’ (O'Malley, Poynter & Russell 1961).

The revolutionary observations of William Harvey paved the way for a systematic study of blood circulation and the findings of capillaries and red blood cells, by Marcello Malpighi (in 1661), provided the missing link between arterial and venous blood while Richard Lower (in 1669) and Antoine Lavoisier (in 1777) explained the ${\rm CO}_2$–${\rm O}_2$ exchange in the lungs, thus completing the picture. In the same period, Galvani (1791), using dissected frogs, proved that electricity could trigger muscle contraction. The 18th century produced a number of complementary independent discoveries, although it was only in the middle of the 19th century that Purkinje (Cavero, Guillon & Holzgrefe 2017) reported the existence of a web of fibres (the nerves) within the myocardium whose role in the transmission of electrical signals was understood in the early 20th century (Tawara 1906). However, His (1949) proved that the heart chambers were electrically connected and that the heart could beat independently of the central nervous system.

In the 20th century, technology appeared in cardiology: W. Einthoven, in 1901, built the first electrocardiograph which weighed $270$ kg, occupied two rooms and required five people to operate it: he won the Noble Prize in 1924. He was a forerunner of telemedicine as in 1905, using a telephone cable, he transmitted an electrocardiogram from a hospital to his office at a distance of $1.5$ km (Mackay & Mensah 2004).

In 1931, A. Hyman developed the first cardiac pacemaker which imparted electrical impulses to the heart by transthoracic needles. The first heart defibrillation was performed in 1947 by C. Beck who, following the advice of the physiologist C.J. Wiggers, restored the normal heart rhythm by electric shocks generated by an external device.

The second half of 20th century saw the combination of technology and cardiac surgery, with the first aortic valve replacement with a mechanical prosthesis and a totally implantable pacemaker in 1960, an artificial heart in 1969 and a permanent device in 1982. Miniaturized axial pumps, employed as ventricular assist devices, were developed in 2000 and 2002.

Concurrently, diagnostic tools for cardiac imaging progressed quickly with the first echocardiography by C.H. Hertz and I. Edler in the 1950s, computerized tomography for the early diagnosis of stroke in the 1970s and magnetic resonance imaging in the early 1980s.

Nowadays, cardiology and cardiac surgery are mature disciplines whose relevance is easily understood considering the clinical, social and economic implications of cardiovascular disorders. In humans, heart-related diseases are the main cause of death worldwide being responsible for the loss of approximately eighteen million people per year (the Covid-19 pandemic has claimed, in the first two years, five million victims!). According to the World Health Organization, these disorders have been the leading cause of death, about one third of the total, globally for the last fifteen years and, although the mortality is declining, the decrease rate is slowing down. The trend for the related morbidity is even worse as, in fact, the number of cardiopath patients remains very high and, in the EU alone, over sixty million people live with a heart disease (Italianer et al. 2021).

Cardiovascular disorders (CVDs) entail dramatic changes in patients lives making it difficult to carry out normal activities, pursue hobbies or maintain social contacts. Patients of working age have their productivity and, consequently, income reduced thus deteriorating their quality of life.

In most of the countries with westernized economies, the incidence of CVDs is approximately $2\,\%$ of the population which costs $7$–$10\,\%$ of all health care expenditure. Owing to ageing of the population and the worsening of co-morbidity factors (smoking, obesity, diabetes), the incidence of CVDs is predicted to increase in the next decades yielding unsustainable health care costs before the end of this century.

Despite the explosive growth of CVD research since the 1970s, advances in clinical practice and treatments appear to have come to a halt, and finding alternative strategies has become imperative (Fiedler 2015). Novel approaches to the problems are required and the combination of engineering with cardiology is providing unprecedented opportunities both for the comprehension of the heart physiology and for the design and testing of new procedures and devices.

Sophisticated laboratory experiments have been setup to analyse the mechanics of biological tissues, the dynamics of blood in heart chambers and vessels, the fluid/structure interaction of heart valves and the electrophysiology of the myocardium. A similar effort has been devoted to study the interaction of the heart with prosthetic devices and diagnostic tools with the aim of improving the prognosis of heart impairments and their early detection.

The impressive and unrestrainable growth of computer power has further boosted cardiovascular research with computational models, virtual prototypes, digital twins and, more recently, machine learning and ‘big data’ techniques.

The main aim of this paper is to give an overview of the current knowledge and of the main issues in the field with a particular look at the fluid mechanics implications.

The paper is organized as follows. Section  2 describes the heart and circulation of vertebrates and, successively, discusses some allometric laws for the cardiovascular systems of homeothermal mammals. Section 3 illustrates human heart physiology with the main structures and its energy balance. Section 4 is devoted to blood features while § 5 describes some characteristics of cardiac tissues. Experimental and numerical heart models are illustrated in § 6 with a discussion of the main modelling issues. In § 7, we present the main features of cardiovascular flows in physiological and pathological conditions; in particular, in the latter case, we consider cardiomyopathies, prosthetic heart valves and related haemolysis and ventricular assist devices. The article is concluded in § 8 with a discussion of the future outlook and closing remarks.

2. Hearts and circulation

The heart distributes blood to tissues and organs bringing oxygen and nutrients to all the cells and removing their waste. In large organisms, this is efficiently achieved by the circulatory system, a complex network of ‘tubes’ of different cross-sections, which through advection allows the blood to travel along distances of metres within a few seconds: needless to say that without such a pumping organ, only elementary sub-millimetric organisms, with a diffusion-based distribution system, could exist.

In fact, even focussing only on animals with a circulatory system, comparing the hearts of any two of them might be difficult and even misleading on account of all possible variations. The circulation can be open, if the blood circulates freely in the tissues like in insects, or closed, when it flows within veins and arteries. The circulation is said to be incomplete when there is partial mixing between oxygenated and de-oxygenated blood or complete otherwise. Finally, the circulation is single when it crosses the heart only once in a cycle or double if the heart is crossed twice.

In figure 2, we show the circulation layouts for the vertebrate groups showing some analogies but also relevant differences.

Figure 2. Blood circulation schemes for: (a) fish; (b) amphibians and reptiles (the latter have an additional aorta); (c) mammals and birds. Arterial (oxygenated) blood is represented in red, venous blood in blue. The letter ‘A’ indicates the atrial chambers, ‘V’ the ventricular counterpart and ‘AA’ the additional aorta originating from the right side and shunting the main aorta. Only crocodilian reptiles have this structure. (d) Mechanical equivalent of the circulation of mammals.

Fish have a single-loop circulation with a heart made of one atrium and one ventricle which pumps blood to the gills to be oxygenated and from there to the rest of the body. The venous part of the circulation then collects deoxygenated blood and delivers it back to the heart atrium for a new cycle (figure 2a).

Amphibians have a more complex circulation with two loops, one for the blood oxygenation and another for its distribution to the body. To feed two circuits, the heart has two atria but only one ventricle where arterial and venous blood come in contact although they mix only partially (figure 2b). Reptiles all have a double circulation and an additional ‘aorta’ originating from the pulmonary circuit and shunting the two circulations; the heart structure is not homogeneous with some species (turtles and lizards) similar to amphibians and others (crocodiles) to mammals.

Birds and mammals also have a double circulation but it is complete and the heart has four chambers organized into two separate pumps that feed the pulmonary and systemic circulations without mixing venous and arterial blood (figure 2c).

In all cases, the heart propels the blood to arteries and veins in a precise sequence and the correct flow direction is assured by valves located in the heart and operated passively by the flow itself. The pumping action is achieved by a rhythmic, active contraction of the heart and the work done on the flow balances the viscous energy losses, mainly in the capillaries, within each cycle.

As the body mass and the size of an animal increases, so do the extension of the capillary network, the length of the circuit, total volume of blood, energy losses and therefore the size and the power of the heart.

Of course, depending on the particular vertebrate group, the body-size affects on the circulation are different and making analogies among them might be incorrect. Accordingly, to make meaningful comparisons, here we focus on homeothermic mammals with the aim of verifying how the human heart performs with respect to that of similar organisms and if its features are determined only by biology or they are also altered by external factors.

2.1. Allometric laws

Using a mechanical analogy, being the two pumps and the two exchangers (pulmonary and systemic capillary beds) crossed by the same flow rate (in series), we can think of the mammals cardiovascular system as a closed hydraulic circuit driven by a volumetric pump operating and connected to a single exchanger (figure 2d). Following Dawson (2001), the idealized model represents the heart as a cylindrical (or conical) volume of radius $a$, length $l$ and tissue thickness $h$, which undergoes a displacement $\delta$ during a periodic contraction of frequency $f$. Except for the capillaries, all the connecting vessels are lumped into a single tube of radius $r_a$ and length $L_a$. It is also assumed that all the capillaries have the same radius $r_c$ and length $L_c$ and their total number is $n_c$.

This model is very useful to derive allometric scalings although it is clearly an oversimplification of the real circulatory system. In fact, here we concentrate all the volume variation in the heart, all the fluid inertia in large vessels and all the energy losses in the capillaries. From a fluid mechanics point of view, such a system could not work since blood flow is incompressible and the volume contraction of the heart could not be accommodated by the other parts of the hydraulic circuit. Keeping in mind that, in reality, every element of the circulation contributes to different degrees to compliance, inertia and losses of the system, in this section, we will rely on the lumped parameter model to derive some interesting scaling laws and postpone the detailed description of the heart to § 3.

Let $M$ be the body mass of a mammal, the empirical evidence, supported by more than a century of experimental measurements (see Dawson (2001) and references therein), suggests the following allometric relations:

(2.1a–c)\begin{equation} a^{2}l \propto M,\quad r_a^{2} L_a \propto M,\quad n_cr_c^{2} L_c \propto M, \end{equation}

which are the volumes of the heart chambers, of connecting vessels and of the capillaries, respectively. It is worthwhile noticing that the sum of these three quantities yields the total volume of blood contained in the body that should also increase linearly with $M$: this is confirmed by direct measurements on mammals whose blood volume is $6$–$7\,\%$ of the body mass independently of the size (Schmidt-Nielsen 2004).

The fact that also the total heart mass is linearly correlated with $M$ yields

(2.2a–c)\begin{equation} a \propto M^{1/3},\quad l \propto M^{1/3},\quad h \propto M^{1/3}. \end{equation}

The blood propelled at each contraction (stroke volume) has a similar scaling $SV \sim al\delta \propto M$ which gives $\delta \propto M^{1/3}$, which is confirmed by direct measurements.

This volume is pushed, with a mean velocity $U$, into a vessel of effective radius $r_a$ and length $L_a$ with a period $T \sim f^{-1}$; balancing inertial and pressure forces yields a pressure difference

(2.3)\begin{equation} P_I \sim \frac{\rho alUL_a f^{2}}{r_a^{2}}, \end{equation}

with $\rho$ the blood density which is the same for all mammals. The quantity $P_I$, estimated as the mean arterial pressure, shows such a weak dependence on the body size that it is often considered constant (figure 3d). The values $P_I \approx 107$ mmHg for a $M=10^{-2}$ kg mouse and ${\approx }130$ mmHg for a $M=7\times 10^{2}$ kg horse (or ${\approx }150$ mmHg for a $M=4\times 10^{3}$ kg elephant) support this hypothesis although it is definitely counterintuitive since $P_I$ is thought of as overcoming the effects of gravity and of the vascular tree resistance. Clearly, the former can not be relevant as gravity acts equally on the ascending and descending flows thus balancing its action; direct confirmation of this argument comes from astronauts during space flight whose blood pressure decreases by less than $10\,\%$ despite effective gravity drops to zero (Fu et al. 2019). Concerning the friction losses, these certainly depend on the extension of the capillary network which Miller (2018) showed could be anyway perfused also at much lower mean arterial pressure.

Figure 3. (a) Heart rate (in b.p.m.) versus body mass (in kilograms) for some homeothermal mammals. The line is the power law $f \propto M^{-1/4}$. (b) Total capillary length (in metres) versus body mass (in kilograms). The line is the power law $n_cL_c \propto M^{5/6}$. (c) Radius of the pulmonary capillaries (in millimetres) versus body mass (in kilograms). The line is the power law $r_c \propto M^{1/12}$ (replotted and adapted from Dawson 2001). (d) Mean arterial pressure (in millimetres of mercury) versus body mass (in kilograms). The line is the power law $r_c \propto M^{0}$ (replotted and adapted from McMahon & Bonner 1983).

In fact, Miller (2018) suggests that the value $P_I \approx 100$ mmHg is fixed essentially by the glomerular capillary pressure at the level of kidneys that must be $P_G \approx 60$ mmHg for the ultrafiltration process to be efficient and $P_G$ is determined by largely invariant aspects of blood chemistry, plasma oncotic and capsular hydrostatic pressures which are very much the same for all mammals.

We recall now that the blood flows in a closed circuit; therefore, after a complete loop, a fluid particle must have the same (mean) pressure as at the previous cycle; this implies that the gain $P_I$ produced by the heart must be balanced by the losses in the capillaries $P_V$ (and since the former is independent of $M$, the same must hold true for the latter). However, $P_V$ can be estimated via Poiseuille flow through $n_c$ circular pipes of radius $r_c$ and length $L_c$ with the same total flow rate as the heart. This condition yields

(2.4)\begin{equation} P_I \sim P_V \sim \frac{\mu l a L_c f U}{n_c r_c^{4}}, \end{equation}

with $\mu$ the blood dynamic viscosity which is very similar for all mammals.

Furthermore, it is well established from measurements that the heart rate of mammals varies approximately with the $-1/4$ power of their mass ($f \propto M^{-1/4}$) and this is confirmed in figure 3(a). The same kind of agreement is observed for the rate of oxygen consumption which shows a $\propto M^{3/4}$ relationship; since oxygen is ultimately diffused to the tissues by the capillaries, its consumption can be related to their number and length to obtain $r_cn_c \propto M^{5/6}$ (Dawson 2001), which is again supported by the data of figure 3(b).

The above relations can be combined together to obtain $L_a \sim M^{1/4}$, $r_a \sim M^{3/8}$ and $r_c \sim M^{1/12}$, which are all confirmed by direct measurements and, to save space, we report only the last of them in figure 3(c).

The relation $f \sim M^{-1/4}$ of figure 3(a) has been further explored by Bassil, Zarzoso & Noujaim (2017) who confirmed the same scaling for the single phases of the heartbeat and explained the $-1/4$ power law by observing that the specialized electrical conduction system of the heart carries the signal through a branching network with terminal units which are size invariant for all mammals.

The weak dependence of the capillary radius $r_c$ on the body mass is instead due to the dimensions of the red blood cells whose size is approximately constant for mammals: the diameter of erythrocytes for mice is ${\approx }5.8\ \mathrm {\mu }{\rm m}$, for humans is ${\approx }7.6\ \mathrm {\mu }{\rm m}$, for horses is ${\approx } 6.0\ \mathrm {\mu }{\rm m}$ and for whales is ${\approx } 8.0\ \mathrm {\mu }{\rm m}$. There is not a unique explanation for the reason why the dimension of these cells is size independent although there are indications that this value is optimal to reduce the work of the heart to deliver oxygen to all the body cells (Schmidt-Nielsen 2004).

An interesting consequence of these allometric scalings is that the volumetric blood flow rate through the circulatory system is $\dot {Q} \approx al\delta f \propto M^{3/4}$ that, combined with $r_a \propto M^{3/8}$, yields a size independent blood velocity in the large vessels $U_a \sim \dot {Q}/r_a^{2}\propto M^{0}$. This is confirmed experimentally by Seymour et al. (2019) who found for the largest mammal arteries a slope equal to $2$ of the derivative of $\log \dot {Q}$ with respect to $\log r_a$.

However, the flow rate $\dot {Q}$ must cross also the capillary network thus yielding a velocity $U_c \sim \dot {Q}/(n_cr_c^{2}) \propto M^{-1/24}$ which is indistinguishable from a constant on account of the data scatter. This ‘constant’ velocity is presumably fixed by diffusivity of ${\rm O}_2$ and ${\rm CO}_2$ from the capillaries to the surrounding tissues which is also a size-independent process; we must admit, however, that we have found neither an experimental confirmation nor a confutation of this result.

From the above allometric scalings, it appears that the hearts of all mammals behave according to the same laws, and humans do not show any special feature. To complete the discussion, we report, in figure 4(a,b), some data analysed by Zhang & Zhang (2009) for the heart rate ($f$ in b.p.m.), life expectancy ($LE$ in years) and total number of heartbeats ($HB$) of $16$ mammals spanning seven orders of magnitude in weight, from the mouse (${\approx }3\times 10^{-2}$ kg) up to the whale (${\approx }1.7\times 10^{5}$ kg). As already mentioned, the smallest animals have the largest $f$ and, as indicated by the data, they have also the shortest lifetimes. In fact, the mouse has $f \simeq 550$ b.p.m. with a $LE\simeq 2$ years, while the whale has $f \simeq 30$ b.p.m. with $LE \simeq 30$ years. The remarkable observation, however, is that all the other mammals in between these two extrema show a consistent behaviour with small deviations with respect to the empirical correlation $LE \approx 33.5-0.0535\times f$, the only exception being humans! In fact, humans, with an average $f=70$ b.p.m. and a life expectancy of $LE\simeq 73$ years, are definitely outliers with a $LE$ three times bigger than other mammals (figure 4a).

Figure 4. (a) Heart rate (in b.p.m.) versus life expectancy (in years) for some homeothermal mammals. The dotted line is the best logarithmic fit $LE \simeq 63\text {--}9.8 \ln (f)$. (b) Heartbeats in a lifetime (in units $HB=f\times LE$, the latter quantity expressed in minutes) versus heart rate (replotted from Zhang & Zhang 2009). In panels (a,b), the red circle is the human life expectancy in the Paleolithic Age. (c) Human life expectancies during the ages (data collected from https://ourworldindata.org/life-expectancy).

Similar indication comes from figure 4(b) showing the total number of heartbeats in a lifetime ($HB=f\times LE$); it appears that all mammal hearts are made to beat a given number of times $HB \approx 5\times 10^{8}$–$1.5\times 10^{9}$ which can be ‘spent’ in a short lifespan with a fast heart rate or vice versa. Once again, humans are the outliers with their ${\approx }3\times 10^{9}$ beats that almost triples the highest figure of other mammals.

Thus a legitimate question is whether humans are a biological exception or the anomaly is rather due to the improved living conditions in which the environment has been partially adapted to their own needs. To this aim, figure 4(c) reports the human life expectancy during various periods from the Paleolithic Age up to the present day; clearly man has not been always long-lived and the largest $LE$ increases have been achieved only within the last century when medical research has experienced impressive progress. If we report in figures 4(a) and 4(b) the Paleolithic life expectancy (rather then the present one), humans perform just like the others thus suggesting that our heart has evolved biologically like any mammals although our capability to modify the environment and the progress of medical science has considerably improved the overall picture.

Several comments are in order on the significance of the above discussion: animal data are unavoidably affected by large uncertainties since heart rates are often measured on anaesthetized or restrained specimens. Life expectancies vary considerably from wild to captivity and, equally important, life duration is determined by many concurrent factors (diseases, predation, fighting, food shortage) which affect different groups of similar animals to different degrees; the heart performance is just one of these factors and maybe not always be the most important.

Human lifespans are also uncertain and often even biased: those of the oldest ages have been obtained from scarce statistics, while in the Classic and Medieval ages, only selected populations (Greeks–Romans and Europeans, respectively) are represented. In addition, the averages are affected by child mortality and peerage/plebs differences: for example, the late Medieval English population had an average $LE \simeq 30$ years; however, if a member of the aristocracy survived the 21st year, the life expectancy increased up to $64$ years. Within these large scatters, the data can give only an indication although those relative to $1900$, $1950$ and $2017$ are more reliable since they are obtained as world averages.

We can close this section by concluding that the human heart is biologically similar to that of other mammals whose features follow allometric scalings derived from simple basic laws of mechanics. Nevertheless, thanks to better living conditions and the progress of clinical and surgical practice, the total number of heartbeats in a lifetime has almost tripled for humans during the last century: pushing this limit beyond the present value and extending the human life expectancy is the implicit motivation driving the research on the cardiovascular system.

3. Heart physiology and structure

The reductionist approach of the previous section was aimed only at comparing the human heart to those of other mammals and at deriving scaling laws from basic principles. That method, however, is clearly inadequate to describe the heart structure since it lumps all its complex functions into a periodic volumetric variation: in this section, we give a detailed description of the heart physiology and its architecture.

The pumping function of the heart operates synergistically with the connected vessels shown in figure 5. Here, the schematic of the system includes arrows indicating the blood direction and the names of the main elements: since the circulation forms a closed loop, any location can be assumed as a starting point.

Figure 5. Main structures and ‘piping’ of the heart: RA, right atrium; RV right ventricle; LA left atrium; LV left ventricle; Tv, tricuspid valve; Pv, pulmonary valve; Mv, mitral valve; Ao, aortic valve; Ivc, inferior vena cava; Svc, superior vena cava; Rpa, right pulmonary artery; Lpa, left pulmonary artery; Rpv, right pulmonary veins; Lpv, left pulmonary veins; Ao. aorta; Da, descending aorta; Ba, brachiocephalic artery; Ca, common artery; Sa, subclavian artery.

The superior and inferior venae cavae (Svc and Ivc) collect deoxygenated blood from the whole body and direct it to the right atrium (RA); from there, blood flows to the right ventricle (RV), initially driven by its inertia and eventually by an active RA contraction that completes the ventricle filling. A contraction of the RV then follows to increase the blood pressure therein thus causing the closing of the tricuspid valve (Tv) and the opening of the pulmonary valve (Pv). This allows the lungs to be fed by deoxygenated blood through the left and right pulmonary arteries (Lpa and Rpa).

Within the lungs, carbon dioxide/oxygen exchange takes place and oxygenated blood returns to the left atrium (LA) through the left and right pulmonary veins (Lpv, Rpv). Similarly to the right counterpart, the LA feeds the left ventricle (LV) and, when it is filled, a vigorous contraction of the latter closes the mitral valve (Mv), opens the aortic valve (Av) and squeezes oxygenated blood into the aorta (Ao), which perfuses the entire body. The head and the superior limbs are fed by the epiaortic arteries (brachiocephalic, left common and left subclavian: Ba, Ca and Sa, respectively) while all the internal organs and the inferior limbs receive blood from the descending aorta (Da).

The above main arteries branch hierarchically into smaller arteries and eventually to arterioles and capillaries which, by diffusion, distribute oxygen and nutrients to all the cells of the body. At the same time, carbon dioxide and waste products are collected by the blood in the capillaries which converge to venules, minor veins and large veins up to the inferior and superior venae cavae that convey the blood again to the right atrium for a new cycle.

The complex cardiac dynamics are the result of coordinated actions of highly specialized parts that operate synergistically: to help with the description of the main mechanisms, we report some representative views of the heart starting from the anterior (figure 6a) which shows the appearance of the heart if the chest were open by a cut of the sternum.

Figure 6. Representative views and sections of the heart: (a) anterior view; (b) vertical section of all four chambers; (c) transverse section of the ventricles; (d) vertical section of the complete left side; (e) transverse section of the fibrous skeleton with the valves. A mitral valve; B tricuspid valve; C pulmonary valve; D aortic valve. In panel (a), section planes are reported with the same colours as the frames of the other panels.

A vertical section, as in figure 6(b), reveals the two pairs of chambers of the right and left heart; an evident difference between atria and ventricles is the thin smooth walls of the former compared with the thick corrugated counterpart of the latter. The reduced thickness of the atria ($\approx$2–3 mm) is due to their function of conveying blood to the ventricles whose filling occurs for approximately $3/4$ of the volume just passively owing to the inertia of the blood stream. Only at the end of the filling process, a synchronous contraction of the atria completes this phase.

The ventricles have much thicker walls ($\approx$4–10 mm) since their contraction must close the atrioventricular valves to prevent atrial regurgitation, and increase the blood pressure to high enough values to open the semilunar valves and feed the pulmonary and systemic circulations. The high ventricular pressures during contraction require the atrioventricular valves to have fibrous strands tethering the free margins of their leaflets to prevent them from everting into the atria (figure 6b). These strands are called cordae tendineae and they are connected to protrusions of the myocardium (papillary muscles) which contract with the ventricle to enhance the pulling action. The inner surface of the ventricles has many grooves called trabeculae carnae; their exact function has not been completely addressed yet although there are indications that they enhance cardiac contractility (Munro et al. 2018).

Figure 6(c) evidences another interesting feature which is the different thickness of the left and right ventricular walls. This is easily understood considering that the right heart drives the pulmonary circulation whose pressure in the largest artery does not exceed the value $20$ mmHg ($2650$ Pa); in contrast, the systemic circulation is much more extended and the mean aortic pressure is approximately $100$ mmHg (13 300 Pa).

These pressures produce cyclic forces, up to several Newtons, on the heart valves (figure 6b,d) and on the tissues. To withstand these high forces, heart valves and tissues connect to a fibrous skeleton, embedded in the plane of figure 6(e), to create a stiff frame which supports the loads.

The muscular part, the myocardium, is a complex composite structure with many different types of cells, the myocytes, which allow the active contraction and the transmission of electrical impulses. Atrial and ventricular myocytes are specialized for contraction triggered by an activation potential which is also transmitted to neighbouring cells.

In addition to the contractile myocytes, there are other cells which have negligible contractility but are specialized as excitatory and conductive fibres. They trigger the heartbeat and assure the correct activation timing is achieved throughout the myocardium: these cells constitute the electrophysiology system, the real ‘orchestra conductor’. In physiological conditions, the heartbeat starts from the sinoatrial (SA) node, made up of self-depolarizing cells, which act as the primary pacemaker (figure 7). The electrical signal then travels, at a velocity of approximately $1\ {\rm m}\ {\rm s}^{-1}$, along the internodal pathways and the bundle of Bachmann thus depolarizing and contracting both atria within $t\approx 80$ ms. The fibrous skeleton is an electrical insulator and prevents the activating potential from moving from the atria to the ventricles; the atrioventricular (AV) node is the only passage for the signal and there the conduction velocity is only ${\leqslant }0.05\ {\rm m}\ {\rm s}^{-1}$, thus delaying the transmission of $\approx$100 ms before moving beyond. This delay is key for the filling of the ventricles which is completed by the atrial contraction while the former are still relaxed.

Figure 7. (a) Scheme of the main branches of the electrophysiologic system of the heart. ‘SA’ is the abbreviation for sino-atrial, ‘AV’ stands for atrio-ventricular. (b) Physiologic electrocardiogram tracing for a heart rate of $\approx$70 b.p.m.

Once the signal has reached the bundle of His, it speeds up to $2\ {\rm m}\ {\rm s}^{-1}$ and, following the left and right bundle branches in the septum, at $t\approx 190$ ms, the ventricle apex is depolarized. From there, the signal further accelerates along the Purkinje fibres reaching speeds of up to $4\ {\rm m}\ {\rm s}^{-1}$ and rapidly reaches the ventricular myocytes where it propagates at $1\ {\rm m}\ {\rm s}^{-1}$: by time $t \approx 230$ ms, both ventricles are fully depolarized with the contraction lasting $\approx$180 ms. This time interval is called the refractory period during which no new stimuli can trigger a new contraction.

It is worth mentioning that although the SA node is the main pacemaker, it is not the only one and other parts can start the cycle in case the first fails: the AV node has an intrinsic firing rate of $40$–$60$ b.p.m., the bundle of His $\approx$40 b.p.m. and the Purkinje fibres $\approx$20 b.p.m. In physiologic conditions, however, these latent pacemakers do not operate since the SA node ($60$–$80$ b.p.m.) has the fastest frequency and the shortest refractory period, thus setting the contraction rate for the whole heart.

The depolarization and repolarization activity of the heart generates feeble electric fields that can be detected by sensors placed on the skin to produce the typical electrocardiogram tracing (figure 7b). Each wave marks a specific event of the heart activity and cardiologists are trained to recognize pathologic behaviours looking at deviations of the tracing from the standard shape. The first P-wave of figure 7(b) is generated by the atrial depolarization and it has usually a weak intensity ($\approx$0.25 mV). The interval in between the P-wave and the QRS-complex reveals the time taken by the signal to cross the AV-node; the activity is too weak to be detected but its duration gives relevant information. The Q-wave is generated by the depolarization of the interventricular septum, the R-wave is due to the apex region which, owing to its mass, generates the most intense signal. The negative S-wave is produced by the depolarization of the left ventricle region farthest from the apex. Meanwhile, atria have repolarized but their activity is hidden by the ventricles contraction. When ventricles repolarize, they generate the T-wave and the time interval QT shortens when the heart rate increases. Sometimes the T-wave is followed by a small U-wave caused by the repolarization of the papillary muscles which are the latest to relax.

The contraction of the myocardium produces a thrust on the blood and the heart valves ensure the resulting flow is in the appropriate direction. There are two heart valves for each side of the heart (figure 6), one between atrium and ventricle and one separating ventricle and main vessel. Each valve is made of flapping, thin and tough leaflets of connective tissue lined with endocardium, the same coating as the heart chambers, and connected to the fibrous skeleton through a stiff annulus.

All valves are passive and are activated as a result of the hydrodynamic loads, mostly pressure, produced by the blood motion which, in turn, is determined by the valves configuration (fluid/structure interaction). In the right heart, the atrioventricular valve is the tricuspid, it has three leaflets whose free edges are connected to the ventricular papillary muscles by a web of cordae tendineae which prevent their prolapse into the atrium. This valve opens during ventricle diastole, as it relaxes and blood flows from the atrium, while it closes during systole when the contraction increases the ventricular blood pressure; this pressure rise causes the pulmonary (or pulmonic) valve to open and push blood into the common trunk to reach the lungs. As the ventricle contraction ceases, blood pressure therein drops and the tricuspid valve opens while the counterpressure within the pulmonary arteries seals the semilunar valve.

The left side of the heart has very similar dynamics, with the mitral valve between atrium and ventricle and the aortic valve between ventricle and aorta, the only difference being the larger pressures that these valves have to withstand which, as will be discussed in § 7.5, is the reason why they are more subjected to morbidity.

We can better appreciate the synergistic nature of the heart by looking more closely at the dynamics of these valves. In fact, the aortic valve could not operate without a suitable impedance of the downstream circulation and compliance of the aorta which stores part of the ventricle mechanical work via elastic deformation of its walls; this ensures the blood in the aorta remains pressurized during diastole and produces the counterpressure that allows the aortic valve to close. This clearly shows that the heart and aorta operate together and the pumping of the former could not be possible if the latter were replaced by a rigid pipe.

As an aside, we note that during mitral and tricuspid valve closure their leaflets approach while moving towards the atria and even after coapting, they keep retreating until the tension balances the pressure difference. This causes a physiologic backflow called ‘false regurgitation’ which is necessary for the proper sealing of the atrioventricular valves (Collia, Zovatto & Pedrizzetti 2019); the same phenomenon, however, produces also a pressure spike that, in the long term, could damage the thin atrial walls. To avoid this problem, each atrium is provided with a histologically distinct appendage which operates as a decompression chamber during ventricular systole or in the event of increased atrial pressure (Al-Saady, Obel & Camm 1999). While in the past, atrial appendages have been considered an insignificant portion of the cardiac anatomy, today, their role in the heart physiology has been recognized; they have a variety of different morphologies which have pathological relevance in case of atrial arrhythmia since clots tend to form preferentially within appendages (Garcia-Villalba et al. 2021).

It has been mentioned several times that the heart distributes oxygenated blood and nutrients to all the body to make possible the metabolic processes. However, the myocardium itself also must be perfused and this happens through the coronary circulation which is another example of the highly intertwined heart dynamics. Even though the heart accounts only for $\approx$0.5 % of the body mass, it absorbs $\approx$5 % of the cardiac output to perfuse the oxygen-eager myocardium.

The aortic root has a bulge with three lobes, the Valsalva sinuses, and from two of them, the left and right coronary arteries depart. These are $3$–$4$-mm diameter arteries that branch and taper up to a network of capillaries with the incredible density of $2500$ per ${\rm mm}^{2}$ of myocardium; the exact configuration of the various branches has a large variability among humans, even if there are some common features of the main tracts which are reported in figure 8(a).

Figure 8. (a) Scheme of the main coronary arteries and their branches. Rca, right coronary artery; Lca, left coronary artery; Lad, left anterior descending artery; Lcx, left circumflex artery; Lma, left marginal artery; Da, diagonal artery; Aia, anterior interventricular artery; Rma, right marginal artery; Pia, posterior interventricular artery; Ana, atrioventricular nodal artery. (b) Aortic (upper) and coronaric (lower) velocity profiles by means of pulse wave Doppler echo-cardiography recorded in vivo in a young healthy man with a heart rate of $\approx$70 b.p.m. (adapted from de Tullio, Pedrizzetti & Verzicco 2011b).

One problem with the heart perfusion is that the intake of the coronary arteries is in the aortic root and blood is pumped into the aorta during systole, when the vigorous myocardium contraction closes the capillary network and prevents the tissues from being perfused. However, thanks to the elastic energy accumulated in the aorta deformation, blood remains under pressure also during diastole, when the ventricles relax and perfusion is possible; this is known as the tissue pressure effect and it is confirmed by direct ultrasound measurements of the blood velocity as reported in figure 8(b) showing the out of phase fluxes. This is another example of interdependent phenomena that could not be possible if the heart did not operate in combination with other vascular elements.

Deoxygenated blood is removed from the heart muscle by several cardiac veins and most of them converge to the coronary sinus, a complex muscular vein whose outlet ends directly into the right atrium with a small (Thebesian) valve. Similarly to the coronary arteries, cardiac veins also have a variable anatomy, such as the presence of the Vieussens valve between the coronary sinus and the great cardiac vein only in $80\,\%$ of subjects and with a number of leaflets that can range from one to three.

Coronary circulation has many more sophisticated regulatory mechanisms since the heart has limited anaerobic capacity and an oxygen shortage of just a few minutes can result in permanent damage. For example, unlike skeletal muscles which use $30$–$40\,\%$ of the available oxygen, the myocardium can extract up to $80\,\%$ of it from the blood. In case an increased flow rate is needed, coronary vasculature can dilate to reduce resistance. A comprehensive discussion of all the coronary circulation details can be found in Goodwill et al. (2018).

If the electric system of the heart is the ‘orchestra conductor’, the Wiggers diagram of figure 9 is the executed ‘musical score’: in the various plots, all the relevant phases of the heart physiology are represented through the pressures in the left and right sides, the volumes of the heart chambers, the ECG signal and that of the phonocardiogram (PCG).

Figure 9. Wiggers diagram for one representative heartbeat at $\approx$70 b.p.m. The letters in the figure indicate the various phases of the cycle according to the following legenda: a, left atrium systole; b, isovolumetric left ventricle contraction; c, left ventricle contraction and early ejection phase; d late ejection phase; e isovolumetric left ventricle relaxation; f rapid diastolic filling; g slow diastolic filling; h mitral valve closes; i, aortic valve opens; j aortic valve closes; k mitral valve opens; l right atrium systole; m right ventricle isovolumetric contraction; n right ventricle contraction and early ejection phase; o late ejection phase; p right ventricle isovolumetric relaxation; q tricuspid valve closes; r pulmonary valve opens; s pulmonary valve closes; t tricuspid valve opens; u rapid diastolic filling; v slow diastolic filling. ECG is the signal of the electrocardiogram, PCG is the phonocardiogram used to monitor the heart through it sounds.

The latter is used to monitor the heart function through its sounds which are produced by the motion of the tissues and the vortical blood flow. The basic approach consists of the traditional auscultation by the stethoscope although the PCG high-fidelity recording also allows the detection of feeble murmurs, and the frequency analysis of the signal permits quantitative assessment of specific phases (Mizuno et al. 2015).

Assuming (arbitrarily) the beginning of the heart cycle in the middle of systole, figure 9 shows the blood pressures in the atria just slightly above those in the ventricles that increase their volumes because of the flowing blood. As they expand, the tension of the myocardium increases and the inertia of the stream is not enough to counter its resistance, thus an active contraction of the atria is needed to complete filling the ventricles. The SA-node then triggers an electrical impulse which quickly propagates across the atria and depolarizes them; this is evidenced by the P-wave of the ECG (figure 7b), followed by the pressure increase in the atria and the ventricular expansion to their maximum volume. Meanwhile, the activation potential has crossed the barrier of the AV-node and it can propagate along the interventricular septum (Q-wave of the ECG), the apex (R-wave) and the ventricular walls (S-wave). The myocardium depolarization activates its contraction that initially does not induce a significant volume variation (‘b’ and ‘m’ are isovolumic contractions) although the ventricular blood pressure sharply increases thus inducing the impulsive closure of mitral and tricuspid valves. The smashing of the valve leaflets, produced by the atrioventricular counterpressure, is so violent that its sound is the first we can hear in a heartbeat (the ‘lub’ of the ‘lub–dub’ onomatopoeia) by just putting our ear on a chest and this is the 1st sound of the PCG trace. As the ventricular pressure exceeds the value in the pulmonary trunk and the aorta, the pulmonary and aortic valves, separately, open and blood flows to the vessels. Ventricles volume now decreases reaching a minimum at the end of systole when the repolarization of the myocardium (T-wave of the ECG) starts the diastolic phase. Initially, the muscle relaxation is not accompanied by a volume expansion (‘e’ and ‘p’ are the isovolumic relaxations) although there is a rapid ventricular pressure drop that shuts the semilunar valves thus allowing the blood in the vessels to maintain the high pressure needed to feed the circulations. The valve closure produces the ‘dub’ of the ‘lub–dub’ sound and the second spike of the PCG trace.

As the ventricle relaxation proceeds, its pressure keeps decreasing until it drops below the atrial values and the atrioventricular valves open. The rapid filling of the ventricles starts (referred to as early wave or E-wave) and their rapid volume increases together with the intense blood recirculations sometimes producing a feeble 3rd sound in the PCG. The cycle is now complete and a new one can start.

In addition to the features just discussed, here we note that although the left and right sides of the heart operate jointly, they are not perfectly synchronous owing to different structural details and the finite conduction velocity of the electric signal propagation. For example, the waves produced by the atria contraction have the right one starting before the left since the electrical signal originates from the SA-node located on the right atrium (figure 7). Also the isovolumic contraction and relaxation have shorter durations for the right side owing to the smaller thickness of the tissue. Similar comments apply to the instants in which the atrioventricular valves open or close, when atrial and ventricular pressures cross or in which the semilunar valves operate.

In the bottom panel of figure 9, we have reported only the main events of the PCG referred to as the 1st (S1), 2nd (S2) and 3rd (S3) sounds; in fact, they consist of composite signals which are generated by several concurrent events. For example, S1 is initiated by the contraction of the ventricles (systole) and the loudest components are produced by the closure of the mitral and tricuspid valves. The former precedes by the latter $20$–$30$ ms and this produces two distinct peaks known as splitting of the first heart sound.

Also, the S2 consists of two main components which are the closure of the aortic and pulmonary valve separated by more than $20$ ms that characterizes the splitting of the second sound. There are additional components, although of small intensity, produced by haemodynamic events in the vessels immediately downstream of the semilunar valves.

Finally, there is the rare feeble 3rd sound which is generated by the vortical blood flow sweeping the compliant left ventricle walls during the passive filling phase. Features like pitch, intensity, duration and distance from S2 of this sound are used as diagnostic indicators to detect abnormal ventricular function.

Finally, among the time evolutions of the heart chamber volumes (third plot of figure 9), that of the left ventricle gives relevant information about the heart functionality: the difference between the maximum ($V_{M}$) and minimum ($V_m$) volumes during a cycle yields the quantity of blood $SV$ pumped in a beat, the stroke volume, while the ratio $EF=SV/V_{M}\times 100$ is called the ejection fraction. Values $50\,\% \leqslant EF\leqslant 70\,\%$ are considered physiological and indicate a healthy heart, while $EF < 50\,\%$ evidences impaired pumping capabilities in which blood has a residence time longer than normal in the ventricle; $EF < 35\,\%$ is considered life threatening. The product of the stroke volume times the heart rate yields the cardiac output expressed in litres per minute, $CO=HR\times SV \approx 5\ {\rm l}\ \min ^{-1}$ in healthy subjects.

3.1. Furnace or pump?

The heart, just as any other muscle, converts chemical energy into mechanical work; however, every real transformation has an efficiency $\eta =P_w/P < 1$, with $P_w$ and $P$ mechanical and chemical powers, and thus only part of the energy is used to pump blood and heat is also produced in the process.

From the volume and pressure variations during the cycle reported in figure 9, it appears that myocardium contraction/relaxation can occur in two main modes: the first (isometric) tightens/loosens the fibres without varying their length to raise/lower the pressure in the heart chambers. The second (isotonic) occurs at an approximately constant load and pressure but with decreasing/increasing fibres length and thus varying the chambers’ volume.

A common way to represent the chambers’ dynamics is to report, in a Clapeyron plane ($V$–$p$), the data during the cycle and an idealized example for the left ventricle is shown in figure 10(a); in the $V$–$p$ space, where a perfect isometric transformation is a vertical line while an isotonic one is horizontal. A heartbeat cycle is thus a rectangle and the enclosed area is the ejection work (EW) which goes from the tissues to the fluid if the loop is walked in a counterclockwise direction (vice versa otherwise). Also of interest is the triangular area to the left of the cycle which can be thought of as the potential elastic energy (PE) stored by the tissue contraction which is not directly converted into work. Indeed, owing to the elastic properties of the myocardium, the diastolic filling is not completely isotonic and the end diastolic pressure–volume relation (EDPVR) produces a preload at the beginning of the isovolumic contraction. For the same reason, the systolic emptying of the ventricle also has a similar dynamics, which is however complicated by the interaction with the arterial elasticity downstream of the aortic valve, and their interaction determines the afterload at the end of the systolic ejection.

Figure 10. (a) Schematic pressure–volume cycle for a left ventricle during a cycle: EW is the ejection work; PE is the potential energy stored in the myocardium; EDPVR is the end-diastolic-pressure–volume-relationship; ESPVR is the end-systolic-pressure–volume-relationship. (b) The same as panel (a) but for realistic data: LV left ventricle; RV right ventricle; LA left atrium; RA right atrium.

Of course, the real isovolumic phases are also not perfectly isometric and representative plots are shown in figure 10(b) for the four heart chambers. A first striking difference is the pressure ranges at which left and right ventricles operate, which are due to the different extensions of the systemic and pulmonary circulations, and this is also reflected by the different thickness of the myocardium walls. The two atria have similar cycles with a distinctive self-crossing loop entailing positive areas, given by the active contractions, and negative counterparts produced during the passive filling of the ventricles.

By adding (algebraically) the areas enclosed by all loops, a total work of $\approx$1.55 J ($\approx$1.2, $\approx$0.2, $\approx$0.075 and $\approx$0.075 J, for LV, RV, LA and RA, respectively) per heartbeat is found which, with a period $T\simeq 0.860$ s ($70$ b.p.m.), yields a mechanical work of $P_w \approx 1.8$ W.

However, a cardiac output of ${\approx }5.5\ {\rm l}\ \min ^{-1}$ and mean prevalences of $120$ mmHg and $18$ mmHg for the left and right circulations yield a power of $P_w \approx 1.7$ W, which is consistent with the previous value.

As mentioned above, this power originates from the metabolic activity of the myocardium and only a fraction of it is used to pump blood. A detailed analysis of the energy consumption during the heartbeat would require a knowledge of the relations energy–mechanical work (Fenn effect), force–velocity and energy–heat for the various regions of the heart (see Katz (2011) for a thorough discussion). These relations are not easy to determine and, even when available, they are usually obtained from measurements on isolated myocardium specimens or simplified models. However, even if they cannot be readily used to make quantitative estimates, they show that heat is produced both when the muscle contracts without shortening and when it reduces its length. Indeed, a minimum amount of heat is produced also when the muscle is completely at rest because of protein synthesis and basal metabolism.

A simple way to estimate the total heart power is by considering the myocardium oxygen consumption, estimated as ${\approx }13\ {\rm O}_2{\rm ml}\ \min ^{-1}$ every $100$ g, combined with the observation that all metabolic processes produce an average energy of ${\approx }20$ J per $1\ {\rm O}_2 {\rm ml}$ (Piper & Preuse 1993). An average adult human heart has a mass of approximately $300$ g, $20\,\%$ of which is fat and $10\,\%$ connective tissue. Assuming a net mass of ${\approx }200$ g of muscular fibres, we have an energy consumption of ${\approx }520$ J min⁻¹ which yields a power $P=8.7$ W. This figure is compatible with a coronary flow ${\approx }5\,\%$ of the total and a global heat power of $100$ W dissipated by an adult human when considering that the myocardium can extract more than $50\,\%$ of the oxygen available from the blood while skeletal muscles do not exceed 30–35 %.

We can finally compare the blood pumping power $P_w$ with the total power $P$ to compute the heart efficiency $\eta \simeq 0.2$ implying that $80\,\%$ of the chemical energy used by the heart eventually produces heat and only $20\,\%$ mechanical work.

A final comment is in order about the function of the heart: while there is no doubt that the organ physiology is aimed at pumping blood through the circulations, its heat production is definitely not negligible since it employs $4/5$ of the total power. In retrospect, the old debate about whether the heart is a pump or a furnace does not have a unique answer and we can say that even if Galen and Leonardo da Vinci did not recognize the pumping function of the heart, they were ‘$80\,\%$ correct’ by conceiving it as a heat source!

4. The blood

The synergistic dynamics of the heart applies to all the involved elements and its working fluid, the blood, is no exception. During systole, the compression work of the myocardium squeezes blood into the large arteries whose deformation stores part of the work as elastic energy. During diastole, this energy is returned to blood as kinetic energy of the flow. Concurrently, the same amount of energy is dissipated via friction as, to reach all body cells, the network of arteries, capillaries and veins has a total length of ${\approx }10^{5}$ km and a range of diameters $3$ cm–$5\ \mathrm {\mu } {\rm m}$ through which the flow is continuously sustained.

Haemodynamics is a complex subject and several comprehensive books and reviews have been devoted to it (Fung & Zweifach 1971; Cokelet 1980; Fung 1984; Baskurt et al. 2007; Nichols, O'Rourke & Vlachopoulos 2011). Here we will introduce only the main features and limit the discussion to those involved in the heart physiology.

To accomplish its function, blood has a liquid matrix in which deformable cells are suspended and its dynamics depends on the considered spatial and time scales. The liquid phase is a diluted aqueous solution with proteins and salts, the plasma, and it behaves as a Newtonian fluid of similar viscosity as pure water. The corpuscular fraction is made up of platelets (PTC), white blood cells (WBC) and red blood cells (RBC) with a volumetric fraction, the haematocrit, in the range $35\,\%$–$50\,\%$. Just for comparison, a suspension of rigid spheres at these concentrations shows an effective viscosity $10$–$100$ times that of the liquid matrix, while the whole blood, thanks to cells deformability, has an asymptotic viscosity (for shear rates ${\geqslant }50\ {\rm s}^{{-1}}$) of approximately $4$–$5$ times that of plasma, thus making possible its circulation within the limited mechanical power of the heart. A cubic millimetre (a ${\rm \mu}$l) of whole blood contains ${\approx }3\times 10^{5}$ PTC, ${\approx }7\times 10^{3}$ WBC and ${\approx }5\times 10^{6}$ RBC; therefore, the latter dominate the properties of the cellular phase.

Red blood cells are biconcave circular disks, ${\approx }8\ \mathrm {\mu }{\rm m}$ in diameter and ${\approx }2\ \mathrm {\mu }{\rm m}$ thickness, highly specialized to deliver oxygen from the lungs to the tissues and bring back carbon dioxide. For this purpose, they contain haemoglobin, which forms weak chemical bonds with ${\rm O}_2$ and ${\rm CO}_2$ molecules, within an elastic membrane with a large surface-to-volume area. RBCs are highly deformable to fit into capillaries of diameters smaller than their size, and this is accomplished by a specific cytoplasmic structure in combination with the absence of a nucleus.

Blood rheology is quite complex and highly dependent on the specific flow dynamics. In vessels of diameter larger than approximately $80\text {--}100\ \mathrm {\mu } {\rm m}$, blood can be treated as a unique continuum fluid showing a non-Newtonian shear-thinning behaviour. This property is mostly attributed to RBC deformability with the viscosity values depending exponentially on the haematocrit (Linderkamp, Stadler & Zilow 1992). In figure 11(a), we report some measurements performed in a Taylor–Couette rheometer for several haematocrit values by Brooks, Goodwin & Seaman (1970) with curves that can be well described by the empirical fit

(4.1)\begin{equation} \mu (|E|) = \mu_\infty + (\mu_0-\mu_\infty) [1+(\lambda |E|)^{2}]^{b}, \end{equation}

with $|E| = {\rm sym} (\boldsymbol {\nabla }{\boldsymbol {u}})$ the modulus of the rate-of-strain tensor and $\mu _0$, $\mu _\infty$, $\lambda$, $b$ haematocrit-dependent parameters (Siginer, De Kee & Chhabra 1999).

Figure 11. (a) Blood dynamic viscosity $\mu$ versus shear-rate $|E|$ for several values of the haematocrit (data replotted from Linderkamp et al. 1992). The thick solid lines are the $\mu (|E|)$ function given by the Carreau–Yasuda model (Siginer et al. 1999) for a haematocrit of $36\,\%$ and $48\,\%$, respectively. (b) Comparison of steady, laminar velocity profiles in a cylindrical pipe for blood at $27\,\%$ haematocrit (line) and a Newtonian fluid (dashed line) (data replotted from Goldsmith & Marlow 1979).

For each haematocrit, the blood viscosity has a plateau $\mu _0$ for vanishing shear-rates ($|E|$) then decreases and eventually attains an asymptotic value $\mu _\infty$ for large $|E|$. The initial plateau is caused by the propensity of RBCs to aggregate and form rouleaux, stacks of cells which, in turn, tend to organize into networks thus dramatically increasing blood viscosity at low shear rates. These aggregations however are reversible and they break up when subjected to relatively low mechanical loads thus yielding a viscosity reduction for increasing $|E|$.

It has long been conjectured (Charm & Kurland 1974) that rouleaux disaggregation progressed until only isolated cells were left and this resulted in the asymptotic viscosity $\mu _\infty$. However, a closer look at the radial flow distribution in a circular pipe reveals that the actual dynamics has additional important features. In fact, even if we assumed as a starting point a simple Poiseuille parabolic profile, this would give a linear rate-of-strain with the maximum at the wall and zero at the pipe axis; this non-uniform strain distribution causes rouleaux disaggregation to be more likely far from the axis thus producing also a radial viscosity gradient pointing towards the pipe axis. Moreover, RBCs deform under the flow shear leading to an accumulation at the axis (Baskurt & Meiselman 2003); haematocrit values are therefore maxima at the pipe centre and minima near the boundary with the RBC-depleted layer next to the walls which substantially contributes to the effective viscosity decrease for high shear rates.

In figure 11(b), we show the comparison of the steady, laminar velocity profile in a pipe of radius $R$ for a Newtonian fluid with the actual profile for blood at $27\,\%$ haematocrit (Goldsmith & Marlow 1979); the plug-like behaviour in the region around the axis is the result of the local high viscosity, while the RBC lean layers next to the pipe wall allow for high velocity gradients.

Even if the steady rheological blood properties might look complex, they are not fully representative of real conditions. In fact, blood flow is unsteady as the heart pumps at a given frequency and the large arteries deform as a response of the periodic pressure wave. As a result, RBCs are subjected to time-varying loads which interact with viscous and inertial properties of the cells and with aggregation/disaggregation processes to produce a viscoelastic behaviour (Cokelet 1980).

Even without considering the elasticity of the walls, the pulsatile flow alone is considerably more complex than the steady counterpart since the solution depends on the Womersley number $Wo=\sqrt {R^{2}\omega /\nu }$, which gives the relative importance of unsteady inertial and viscous forces. For $Wo < 2$, the flow is viscosity dominated and the resulting flow is parabolic with an amplitude modulated in time with pulsation $\omega$. In contrast, for $Wo > 2$, unsteady effects take over and they dominate the flow at the pipe centre while the near-wall region is still viscous; the resulting velocity profiles are then non-monotonic, with inflection points and uneven distribution of shear stresses. For the human arterial system, we have $Wo \approx 10$ in the aorta, $\approx$4 in large arteries, $\approx$1 in small arteries and ${\approx }10^{-2}$ in capillaries thus including a wide variety of different behaviours.

In the tiniest arteries and capillaries, however, still additional factors have to be considered as the vessel diameter becomes comparable with the red blood cell size. This is quantitatively evident from the plot of figure 12(a) showing the effective blood viscosity normalized by that of pure plasma as a function of the vessel diameter $D$. It can be seen that already for $D < 300\ \mathrm {\mu }{\rm m}$, the viscosity ratio starts decreasing even if only for $D \lesssim 100\ \mathrm {\mu }{\rm m}$ a steep drop is observed. The curve has a minimum for $D \approx 8\text {--}10\ \mathrm {\mu }{\rm m}$ and a new increase is found for smaller diameters. This phenomenon, known as the Fahraeus–Lindqvist effect (Fahraeus & Lindqvist 1931), can be explained by looking at the RBC configurations in blood flows within circular pipes of decreasing diameter in figure 13.

Figure 12. (a) Ratio of kinematic viscosity of the blood with that of plasma for the steady flow in a cylindrical pipe of radius $D$. (b) Ratio of the blood haematocrit ($H_b$) for a steady flow through a pipe of diameter $D$, with the reference value ($H_0$) for the blood in the feeding tank. The grey areas are the uncertainty regions covering the scatter of different measurements. (Data replotted from Baskurt et al. 2007).

Figure 13. Cartoons of RBC arrangements for the steady flow in a circular pipe: (a) $D>d$; (b) $D\approx d$; (c) $D< d$. Here, $D$ and $d$ are the tube and RBC diameters, respectively.

The cartoons show that when the vessel diameter $D$ ${\gg }d$, RBCs can arrange in multiple concentric layers and occupy a substantial volume fraction, thus behaving as described previously and producing a viscosity $3$–$4$ times that of the plasma. However, when $D$ is of the same order as $d$, the cells align in a single row thus decreasing the volume occupation and the haematocrit; this spatial arrangement prevents the cells from interacting and aggregating, yielding a viscosity close to that of the liquid matrix. Finally, in capillaries, red blood cells are deformed to fit into vessels with $D< d$; they attain the peculiar parachute shape which increases again the volume occupation and the haematocrit. In this configuration, energy is needed for the deformation of RBCs which slide along the capillary wall; the result is an effective viscosity increase for $D< d$.

This simple explanation is confirmed by measurements of the haematocrit variation with $D$ (figure 12b) indeed showing a minimum for $D \approx d$, referred to as the Fahraeus effect (Fahraeus 1929). It is worth mentioning that part of the haematocrit decrease with $D$ might also be caused by the RBC arrangement shown in figure 11(b) since vessels always branch from the side walls of larger ones and, thus, they are fed from blood layers which contain less RBCs than the bulk flow.

In all the above discussion, we have not considered the effect of WBCs on the whole blood viscosity on account of their small volume concentration. They certainly have a negligible role in large vessels, although in the microcirculation, the situation might be different, especially in pathological conditions. In fact, WBC sizes and properties depend on the type and status of the cells; they can be activated by chemical or mechanical stimuli and either migrate out of blood vessels to reach inflamed tissues (extravasation) or transiently block a microcirculation channel. Similar arguments apply to platelets whose size is approximately one-fourth of an RBC and in physiologic conditions contribute negligibly to viscosity. However, they are particularly sensitive to shear stress which activates them making possible their adhesion to damaged vessels endothelium or the formation of clots by the capture of impaired RBCs.

Before concluding this section, we wish to mention that blood also exhibits viscoelastic properties which can become relevant in pulsatile cardiovascular flows. In fact, the reversible aggregation and disaggregation of rouleaux, and the periodic deformation of RBCs under pulsatile shear, leads to cyclic storage and release of elastic energy in addition to the already mentioned shear-thinning viscosity. The former behaviour, however, depends also on flow frequency and on the multi-axial strain components and, therefore, it is usually modelled by a linear combination of modes (from $1$ to $5$), accounting for the different relaxation time scales of the fluid (Campo-Deaño et al. 2013), each one being governed by an Oldroyd-like differential equation. Every equation, in turn, depends on empirical parameters which are determined by a fitting against experimental measurements; an example of such derivation is illustrated in Campo-Deaño et al. (2013) for real blood and for blood analogues.

5. The tissues

Heart tissues are extremely composite both in terms of anatomy and properties which are anisotropic and heterogeneous; they are made of different types of specialized cells organized into unique structures aimed at performing specific functions. On account of such a complexity, a thorough description of all their features within a JFM Perspectives paper (written by a fluid dynamicist) is unrealistic. Therefore, we will discuss only those aspects needed to understand the pumping function of the heart and leave further details to specialized textbooks such as Katz (2011) or Piper & Preuse (1993).

The heart is enclosed in a non-compliant bag, the pericardium, whose external surface is in contact with the surrounding organs, mainly the lungs. The internal surface contains a small amount of fluid which reduces the friction with the epicardium, the outermost layer of the contracting heart.

The active muscular part, the myocardium, is the main and thickest layer consisting of cardiomyocytes, complex cells mechanically and electrically connected through special intercellular bridges called intercalated discs (figure 14a). There are many types of myocytes depending on the function and heart region. ‘Working myocytes’ are the contracting units, and they are mostly found in ventricles and atria respectively of bigger and smaller dimensions. A different type of myocyte are the ‘Purkinje cells’ which do not contract much but are specialized for fast conduction of electrical signals. ‘Nodal cells’ are yet another different myocyte responsible for pacemaker and conduction delay activities to initiate and coordinate the contraction of the various regions.

Figure 14. (a) Sketch of a network of myocytes with some basic structures. (b) Example of time evolution of transmembrane potential and mechanical tension for a myocyte.

The myocardium is also rich in fibroblasts which outnumber cardiomyocytes even though they account for a smaller volume fraction owing to their reduced dimensions; they are responsible for the maintenance of the extra cellular matrix, a complex network of collagen and elastin, which is the scaffold of the cardiac tissue and determines most of its passive mechanical properties.

The innermost layer, the endocardium, is a connective tissue lining heart chambers and valves; being directly in contact with blood, it is biologically similar to the endothelial cells of veins and arteries. An obvious function of the endocardium is to act as a physical barrier which protects the myocardium from the vigorous blood flow. Another important function is to dynamically sense the signals coming from the flow (mainly pressure and shear stresses) to stimulate and produce lineages of different cardiovascular cell types like valves, endothelial cells of atria and ventricles or the lining of coronary vasculature (Dye & Lincoln 2020).

The building block of the active heart muscle is the cardiac myocyte whose electrical excitation and mechanical contraction have been thoroughly investigated and can be shortly explained by referring to figure 14(b). At the basic resting state, the cell has a negative transmembrane electrical potential which prevents its contraction. When stimulated by neighbouring cells, membrane channels open and allow sodium, calcium and potassium ions to move across the intracellular and extracellular space beginning the positive-feedback mechanism of depolarization and activating the cell contraction. This ion flux (mainly ${\rm Ca}^{{2+}}$) increases the transmembrane potential to positive values, depolarizes the cell and induces contraction. The rapid depolarization is followed by a long plateau phase in which the potential slowly declines for more than $150$ ms. The repolarization phase starts immediately after with the closure of the ${\rm Ca}^{{2+}}$ channels and the attainment of the initial negative potential after approximately $250$–$300$ ms. When the depolarization process is initiated, the cell enters a refractory period during which it can not be excited again and this avoids tetany which would result in permanent contraction (incompatible with life). It is worth mentioning that the graphs of figure 14(b) are just illustrative examples since the specific timings and amplitudes of action potential and mechanical tension are dependent on the cell type and, therefore, on the region of the heart. Each myocyte depolarizes when it is reached by an electrical impulse and, in turn, it generates a new activation signal which is transmitted to the neighbouring cells through the junctional intercalated disks (cell-to-cell conduction). This implies that the electrical stimulus, and the induced muscle contraction, propagate preferentially along precise directions which make the overall tissue behaviour very anisotropic.

The first source of anisotropy is the cardiomyocyte itself, a cylinder of diameter $15\text {--}25\ \mathrm {\mu }{\rm m}$ and length $50\text {--}150\ \mathrm {\mu }{\rm m}$ which shortens its axis upon contraction; although it reduces its length by no more than $10\,\%$, the presence of different layers and the particular orientation of the fibres (figure 15) allow the ventricles to reduce their volumes by up to $70\,\%$ during systole.

Figure 15. (a) Sketch of the myocardium winding path for the atrial and ventricular walls according to the model by Torrent-Guasp et al. (2001). (b) Sketch of the resulting arrangement of the ventricles.

Locally, myocytes align the axes along a common direction and establish end-to-end junctions to form fibres. These, in turn, are organized into $4$–$6$ cells thick plane sheets (‘laminae’) separated by layers of connective tissue (Clayton et al. 2011) to reduce their coupling. These laminae have different orientations in space and across the myocardium thickness thus yielding a very anisotropic structure as hinted by figure 15. In particular, Pashakhanloo et al. (2016) report very recent measurements showing that the orientation of the fibrous bundles in the atria is homogeneous across the tissue thickness in some regions although in others, the orientation angle changes by up to $90^{\circ }$. In contrast, the ventricles are more complex having more fibre layers with an orientation which rotates continuously from the endocardium to the epicardium by up to $150^{\circ }$. This heterogeneity has implications both for the electrical conduction and the mechanical response of the tissue which considerably complicate its macroscopic description (Guan et al. 2020).

Although the individual behaviour of the cardiac myocyte is well understood and characterized, a healthy human left ventricle consists of ${\approx }5.5\times 10^{9}$ myocytes (Beltrami et al. 2001) and modelling even a single heart chamber starting from a single-cell description is practically impossible. Continuum models are thus needed in which the collective mechanical dynamics of the tissue is prescribed and the myocardium is often modelled as an orthotropic, hyperelastic material although it shows also viscoelastic features evidenced by hysteretic loops in multiple load cycles.

Most of the available information on the mechanical properties of heart muscle has been obtained by experimental tensile tests of small myocardium specimens under different conditions. As shown by Fung (1993), the curve of the tensile stress $\sigma$ (defined as the ratio of the force with the resting cross-section) versus the stretch $\lambda$ (the ratio of the deformed with the resting lengths) depends on many parameters such as temperature, pH of the environment, timings and frequencies of the load, just to mention a few. Nevertheless, in physiological conditions, a recurrent feature is a slope of the stress–stretch relationship which grows linearly with the stress $\sigma$ and can therefore be written as ${{\rm d}\sigma }/{{\rm d} \lambda } = \alpha (\sigma +\beta )$ with $\alpha$ and $\beta$ being empirical tuning parameters; upon integration, it yields $\sigma +\beta = C \times e^{\alpha \lambda }$ with the constant $C$ determined by the imposition of a stress $\sigma _0$ for a stretch $\lambda _0$. Examples of this hyperelastic nonlinear behaviour are given in figure 16 with the data for biaxial tension tests of the human myocardium, fitted with a law similar to the one above (Zhao 2018). If the myocardium specimen is small enough, there is a clear fibre- and a perpendicular cross-fibre-direction and they show different mechanical properties, as is evident in figure 16.

Figure 16. Stress versus stretch data for biaxial tension tests of human myocardium: $\circ$ mean fibre direction $\square$ cross-fibre direction. The lines are exponential-like constitutive laws of the family anisotropic continuum stored energy derived along the arguments of Fung (1993). (Data adapted from Zhao 2018).

It is worth mentioning that the specific values of the stress–strain relationship depend also on the heart region from which the specimen is extracted as the myocytes, the composition of the layers and the orientation of the fibres across the tissue thickness are quite different depending on the function of the tissue.

In addition to the inhomogeneous character of the myocardium and the natural variability among humans, other factors affecting the measurement of the mechanical properties are the preparation and conservation of the samples and the test method (uniaxial tensile, biaxial tensile, triaxial shear, constant versus time-varying loads, etc.); the interested reader is referred to Sperelakis (1989) and Holzapfel & Ogden (2009) for more detailed discussions.

6. Heart models

On account of the inherent heart complexity, cardiovascular research relies heavily on models to study blood flow and tissue dynamics in physiological and pathological conditions. In fact, models are more accessible than the real counterparts, they can be used to test hypotheses or predictions, to validate assumptions or to evaluate the performance of devices. However, for the results to be meaningful, a model has to replicate as close as possible the real system and for the heart, this might be exceedingly difficult. For this reason, hardware laboratory models usually aim at reproducing only specific parts of the heart or single details of the dynamics. Computer models are instead more flexible although at the price of huge computational costs if high-fidelity is required.

6.1. Laboratory experiments

A fundamental pillar of any meaningful experiment is its dynamic similarity with the real problem which is verified only for precise sets of model parameters. Let $q$ be a generic quantity of interest, which will depend through an unknown function $f$ on several variables; for cardiovascular flows, these can be restricted to density $\rho$, viscosity $\mu$ and velocity $U$ of the blood, a spatial scale $L$ of the problem, its periodic time scale $T$ and also some measure of the elasticity of the tissues $E$, their density $\rho _T$ and a reference (intramural) pressure difference $\Delta p$ between the inner flow and the outer environment. It is worthwhile to point out that this is just a small set of variables selected to discuss the main issues in the design of an experimental test although, depending on the specific problem, additional quantities or different sets might come up.

Another important observation is that the geometrical similarity is a fundamental prerequisite of the dynamic similitude. Thus, the laboratory model must be a faithfully scaled replica of the real system. The above length scale $L$ must be therefore representative of all the lengths which completely define the problem geometry; they can easily amount to several tens (for example, one of them is $\lambda _m$ in figure 17) and be difficult to downscale.

Figure 17. Example of experimental duplication (right) of the a specific part of the whole heart (left). The aortic root and aortic valve are considered for the evaluation of the prosthetic valve performance.

Applying the Buckingham $\varPi$-theorem to the relation $q=f(\rho, \mu, U, L, T, E, \rho _T, \Delta p)$ and using the triplet $\rho$, $U$ and $L$, to non-dimensionalize the remaining variables, we can write $\varPi _q =F[\rho U L/\mu, L/(TU), \Delta p/(\rho U^{2}), \rho _T/\rho, E/(\rho U^{2})]$, where $\varPi _q$ is the dimensionless quantity of interest and $F$ an unknown function to be determined by experiments. Here, $\rho U L/\mu =Re$ is the Reynolds number, $L/TU=St$ the Strohual number (which can be rearranged in $\sqrt {2{\rm \pi} Re\times St}=Wo$ to obtain the Womersley number) and $\Delta p/(\rho U^{2})=Eu$ the Euler number.

An exact dynamic similitude requires the non-dimensional groups for the model and the real phenomena to match perfectly: denoting by the subscript $m$ the quantities relative to the model, we can write, for example, for the Reynolds number, $\rho _m U_m L_m/\mu _m = \rho UL/\mu$, which binds the choice of some experimental parameters.

The scale factor $f_S = L_m/L$ is usually equal to unity for cardiovascular problems as $L \approx {{O}}({\rm cm})$ and there is no need to downsize the model. Instead, a major constraint comes from the blood since it is neither practical nor safe (and ethical) to use the real fluid in experiments and alternatives, with different properties, are adopted. For example, the (asymptotic) kinematic viscosity of blood is $\nu = \mu /\rho = 3.5\text {--}4.5\times 10^{-6}\ {\rm m}^{2}\ {\rm s}^{-1}$ while if water is used for experiments, it results in $\nu _m \approx 10^{-6}\ {\rm m}^{2}\ {\rm s}^{-1}$, implying that the same Reynolds number is obtained if $U_m = U(\nu _m/\nu ) \approx U/4$; velocities in the model must be therefore smaller than in the real heart. This relation can be plugged into the similitude of the Strohual number to obtain $T_m = T (U/U_m) \approx 4 T$ yielding an experiment with longer time scales; this is beneficial for experimental measurements as slow unsteady dynamics can be properly resolved in time with less effort than a fast one.

For rigid mechanical heart valves (figure 18b), the exact dynamic similitude can be obtained using the real prosthesis as a model. In fact, the inertia of the leaflets is the same as in reality and the experiment has velocities smaller by the factor $\nu /\nu _m$ and time scales longer by the inverse factor $\nu _m/\nu$; it turns out that the non-dimensional parameters formed by the ratio of unsteady ($F_U \sim \rho U L^{3}/T$) or inertial ($F_I\sim \rho U^{2} L^{2}$) fluid loads with the leaflet forces ($F_L \sim \rho _T UL^{3} /T$) are $F_L/F_I = (\rho _T/\rho )L/(UT)$ or $F_L/F_U = \rho _T/\rho$: they either depend on the product $UT$ or are independent of $U$ and $T$, thus preserving the similitude.

Figure 18. (a) Dacron tube used as aortic graft and detail of the fabric structure. (b) Experimental picture of a mechanical aortic valve in a model aortic root.

A more intuitive way to come to the same conclusion is to note that for rigid leaflets, the material elasticity $E$ does not play a role and, since they are wet on both sides by the same fluid, the same holds true for the intramural pressure $\Delta p$. In this context, the original relation reduces to $q=f(\rho, \mu, U, L, T, \rho _T)$, which in non-dimensional form reads $\varPi _q = F(Re, St, \rho _T/\rho )$, and since the densities of water and blood differ by no more than $5\,\%$, the dynamics similitude is preserved.

However, for deformable structures, $E$ is relevant and the group $E/(\rho U^{2})$ clearly indicates that, with $U$ and $U_m$ being different, the material properties of the model and real device can not be the same if the dynamic similarity has to be preserved.

This can be understood also considering that if the intramural pressure plays a role, matching the Euler number requires $\Delta p_m = \Delta p (\rho _m U_m^{2})/(\rho U^{2}) \approx \Delta p/16$, and using the same tissue elasticity with a largely reduced intramural pressure would produce unrealistic small deformations.

An important consequence of the above estimates is that a cardiovascular experiment using water as working fluid, in dynamic similitude for the flow, can not preserve the similitude also for the structural part when it maintains the same elastic properties as the real system. This is particularly unfortunate if the experiment is designed to assess the performance of an existing prosthetic device, whose material selection and construction process is the result of complex biomedical research and sophisticated industrial procedures, since a model of the same size but with different materials can not be easily manufactured.

For example, when a biological prosthetic aortic valve is to be evaluated by an experiment, it is relatively easy to place the real prosthesis in a full scale aortic root and produce a pulsatile water flow with the same Reynolds and Strohual (or Womersley) numbers as the real system (figure 17). For the flow to be visible, the aortic root has to be made using a transparent material (usually silicon rubber) whose elastic modulus $E_m$ does not match that of the aortic wall $E$ (here we are neglecting the hyperelastic and orthotropic nature of the tissues). Nevertheless, it is still possible to reproduce the correct aortic root deformations by adjusting the wall thickness $\lambda _m$ so as to have $\lambda _m E_m/\Delta p_m = \lambda E/\Delta p$ which, with the Young–Laplace equation, implies the same local radii of curvature for the model and real system.

Unfortunately, this strategy can not be used for the valve leaflets since the model is the real prosthesis ($\lambda = \lambda _m$ and $E=E_m$) that, however, is subjected to different hydrodynamic loads. Similar problems occur if the prosthesis is the aortic wall itself as it is usually a Dacron (polyethylene terephthalate; PET) vascular graft with a complex fabric structure (figure 18a) aimed at reproducing the orthotropic mechanical properties of the natural tissue (Amabili et al. 2020). It must be noted, however, that the Dacron fibre is quite rigid and, for the physiological pressure differences, the stretching is very limited; considering the negligible inertia of the graft with respect to that of the fluid therein, the flow motion turns out to be well reproduced by the experiment even if the internal loads of the structure are not representative of the real problem.

From the above discussion, it is clear that performing laboratory experiments on cardiovascular flows is very challenging because of the large number of variables and the unavoidable technical limitations. Although, in principle, dimensional analysis would allow the design of an experiment in perfect dynamic similitude, in practice, this is impossible and distorted similitudes are the rule. In the latter case, only the main parameters are matched and, for the flow analysis, experiments must have the same Reynolds and Strouhal (or Womersley) numbers as real phenonena. In this case, representative flow structures can be obtained and, with a judicious alteration of the elastic properties, the structure deformations are also in the appropriate range. In contrast, representing correctly the internal stresses of the biological tissues is exceedingly difficult, not only for the presence of the Euler number, but mainly owing to their anisotropic and nonlinear elasticity which is not reproduced by laboratory experiments. A noticeable exception is that of rigid structures, such as mechanical heart valves, since only Reynolds and Strouhal numbers matter and experiments can be performed in exact dynamic similitude.

6.1.1. Blood analogues

If the internal stresses of the tissues or the structure dynamics is the aim of the experimental test, the above distorted dynamic similitudes are not acceptable and alternative fluids, which match better the blood properties while still making possible the execution of experimental tests, must be employed. A common blood analogue is an aqueous solution of glycerol with concentrations in weight in the interval $35$–$45\,\%$ to span the physiological range of variation of the blood dynamic viscosity (figure 19a). The resulting fluid has an acceptable value of density (figure 19b) and it remains fully transparent thus allowing optical measurement techniques, such as particle image velocimetry (PIV), which are particularly useful to characterize quantitatively the flow field in terms of recirculating vortices or separation regions. For the same reason, test models are often made of transparent materials and silicone rubber (polydimethysiloxane, PDMS) is a common choice. An extra bonus of water/glycerol solutions is that the refractive index (RI) of the fluid matches that of PDMS which reduces optical distortions and allows for optical measurements. As reported in Brindise et al. (2018), fluid properties can be further tuned to those of blood by additives like urea or NaI in case of specific needs; in this case, urea is preferable to NaI as it tends to decrease the density and increase the RI thus better mimicking blood and further reducing optical interference with PDMS whose refractive index is ${\textit{RI}} \simeq 1.415$.

Figure 19. Dynamic viscosity (a), density (b) and refractive index (c) of several aqueous solutions as a function of the solute concentration (in weight). $\square$ (red) for glycerol; $\circ$ (blue) for urea; $\blacktriangle$ (green) for NaI. (Data adapted from Brindise, Busse & Vlachos 2018).

In small vessels or when the model contains a wide range of spatial scales, the non-Newtonian blood properties become relevant and the experimental fluid must also respect this feature. It has been shown that the addition to the above mixtures of quantities below $0.1\,\%$ in weight of Xanthan gum yields a fluid with a shear-thinning behaviour similar to that of figure 11. Furthermore, when subjected to pulsatile flow conditions, this fluid is able to reproduce the full viscoelastic behaviour of real blood (Campo-Deaño et al. 2013). Brookshier & Tarbell (1993) expressed by imaginary quantities the unsteady shear stress $\hat {\tau }$ and the resulting unsteady shear rate $\hat {S}$ produced in the fluid. Their ratio defines a complex viscosity $\hat {\mu }=\hat {\tau }/\hat {S}$ whose real part is the dynamic viscosity while the imaginary part is the oscillatory viscosity related to the elastic fluid behaviour; for a range of shear rates $1\text {--}1000\ {\rm s}^{-1}$ and a frequency of $2$ Hz, they showed results in good agreement with blood for physiologic haematocrit values.

At even smaller scales, such as the flow in tiny arteries or capillaries, the interactions of red blood cells become dominant and analogue suspensions must be used. Carneiro et al. (2021) report the results obtained with a suspension of PDMS microparticles with an average diameter of approximately $7\ \mathrm {\mu } {\rm m}$; using a concentration of $21\,\%$ in weight at $20\,^{\circ }{\rm C}$, they obtained similar steady, oscillatory and extensional rheological behaviour as blood at $37\,^{\circ }{\rm C}$ and $41\,\%$ haematocrit.

6.1.2. Cardiovascular phantoms

Another key ingredient of a cardiovascular model is the three-dimensional (3-D) mock-up replicating the blood domain to be investigated, possibly also with the related tissue compliance.

Apart from the actual evaluation of prosthetic devices, like valves or vascular grafts for which the real item is directly tested in the experimental apparatus, the fabrication of a cardiovascular phantom entails three main phases (Yazdi et al. 2018): (i) definition of the geometry; (ii) realization of the mould; (iii) manufacture of the phantom and extraction.

The definition of the phantom geometry is a crucial preliminary phase which heavily affects the validity and the applicability of the final results; since the same issues apply also to computational models this point is discussed with some detail. A geometry can be idealized, anatomic or patient specific and each choice has advantages and disadvantages: for example, a healthy left heart ventricle is an elongated chamber with an average end diastolic volume of $140\pm 20$ ml, an approximately axial symmetry and a diameter of $50\pm 9$ mm and a simple idealized geometry could be half of a prolate ellipsoid with two tubes on the upper flat surface to accommodate mitral and aortic valves (figure 20a). This geometry can be prescribed analytically, is extremely easy to model and it does not need particular CAD software, although it misses most of the peculiar features of a left ventricle and cannot be used to answer specific questions; for this reason, idealized geometries are employed to perform benchmarks, to assess the performance of new experimental apparatus or when the focus of the analysis is on a different region (for example, if the ventricle is used only to direct the flow through a valve or into the aorta).

Figure 20. Different models of left ventricles: idealized model (a); anatomic model (b); patient-specific model (c) and dissection of a real left ventricle human ventricle (d).

An alternative approach is to model the ventricle with a more realistic geometry using medical atlases or extracting features from databases of medical images. For example, Azhari, Beyar & Sideman (1999) defined a uniparametric helical shape descriptor to characterize quantitatively the 3-D geometry of the left ventricle; they then performed the same in vivo scan on small groups of homogeneous volunteers (ten healthy males of age $27.2\pm 3$ years, six with myocardial infarction and three with hypertrophic pathologies) and produced different averaged data for each group. The final result consisted of three anatomical left ventricle geometries containing all the physiologic anatomical details (and the pathologic ones for the second and third group) but without the tiny peculiarities related to the variability of humans. These models become more general as the size of the group from which scans are obtained increases although they are not related to any individual and can not be used to study a specific patient (see figure 20b).

Finally, there is the possibility to extract the actual ventricular geometry from the imaging (MRI or CT scan) of a specific patient to reproduce most of the features of that particular individual (figure 20c): the evident advantage is that the experimental results could be applied directly to the case of interest while the drawback is the reduced generality of the study.

A possible problem of extracting a patient-specific configuration from MRI or CT scans is that these are, in turn, the result of a complex post-processing procedure which converts two-dimensional (2-D) images into volumetric renderings. In fact, for each planar section, by thresholds on the image contrast, boundaries are identified (segmentation) and connected through interpolation (lofting) to produce a 3-D geometry. There are several open-source and licensed software which automatically differentiate the tissues and produce regular StereoLithography (STL) files ready for 3-D printing or for computations (Withey & Koles 2008). Unfortunately, each software relies on smoothing and thresholding algorithms which make the final result non-unique; any patient-specific geometry is therefore similar to the real counterpart only to a certain degree and experienced researchers are needed to tell whether or not the differences are relevant. In figure 20(d), the dissection of a human heart left ventricle is reported to compare the impressive complexity of the heart chamber with the details of the models; clearly, even the claimed patient-specific geometry is an approximation of the real physiology and this should be considered when assessing the validity of the obtained results.

Once the model geometry has been defined, a male mould can be manufactured either by turning a solid block (figure 21a) or by 3-D printing (figure 22a); the mould is the lumen of heart chambers or arteries and the phantom material is applied, usually in liquid form, to the mould surface and then solidified by curing. The transparency of the phantom depends on the mould surface smoothness, thus turned specimens are mirror finished while 3-D printed ones are grouted and painted.

Figure 21. Plexiglass of a turned, mirror finished mould for an anatomic left ventricle (a); painting of a layer of silicone rubber (b); extracted compliant phantom (c) and phantom installed in a duplicator (d).

Figure 22. 3-D printed, grouted and painted mould for an anatomic ascending aorta with the aortic root (a); phantom installed in a duplicator (b). (Mould, courtesy of G. Querzoli).

Given the geometrical complexity of the surfaces, if the phantom is rigid, lost core casting is used with a sacrificial mould which is destructively removed by melting or chemical dissolution. However, thin walled deformable phantoms can be extracted relying on their wall compliance (figure 21c) so that the (expensive) moulds can be used more than once.

When the experiment is aimed at capturing the fluid–structure interaction, wall compliance is a key factor and this is a difficult issue for the fidelity of the model. Deformable phantoms are obtained by pouring silicone rubber in the gap between a male/female mould pair or dipping a male mould in liquid rubber and curing it. In the latter case, thickness uniformity is difficult to control owing to dripping even if a continuous spinning of the mould during curing alleviates the problem (Arcaute & Wicker 2008). As an alternative, liquid rubber can be painted in thin layers around the mould (figure 21b) and quickly cured by a thermal gun; by changing the number of layers, the wall thickness can be controlled even locally although only to a certain degree.

It is worth mentioning, however, that even if deformations can be matched by adjusting the wall thickness or altering the rubber stiffness via different curing processes, the phantom still generally does not reproduce neither the hyperelastic behaviour nor the orthotropic nature of the biological tissues. In the case of myocardium, the active contraction features are also missing; therefore, if flow and structure dynamics have to be reproduced simultaneously, the latter is usually restricted to the correct kinematics.

6.1.3. Cardiovascular pulse duplicators

As discussed above, the rhythmic contraction of the heart produces a pulsatile flow at given Reynolds and Womersley (or Stokes) numbers and with specific pressure/velocity waveforms whose shape depends on the particular location on the circulatory system. In figure 23, we report these quantities sampled in the aortic root and at the beginning of the femoral artery: it can be noticed that peak pressure increases as it is measured downstream while blood velocity decreases. The former counterintuitive behaviour is due to the pulse pressure augmentation generated by the constructive interaction of forward and backward waves travelling along the compliant arteries (Moxham 2003). However, despite the fact that artery calibre decreases downstream, peak velocity decreases as well owing to the blood spilling from the aorta through the many side arteries in between aortic root and femoral arteries.

Figure 23. Waveforms of pressure (a) and velocity (b) of the blood flow in the circulatory system during a heartbeat: red solid line for quantities sampled in the aortic root; blue solid line for quantities sampled at the beginning of the femoral artery.

Pressure and velocity waveforms in cardiovascular problems depend on the heart dynamics but also on the artery network geometry and impedance which have strong variability; therefore, any experimental apparatus aiming at coping with such problems should be able to reproduce the wide variety of possibilities.

Although the actual realization of the device can change considerably, depending on the specific replicated system, the underlying general principle is quite standard and it relies on a controllable pump capable of producing a prescribed, periodic flowrate waveform to feed a closed hydraulic circuit along which the test section is connected. Since reproducing the distributed impedance of the circulatory system is impossible, it is lumped into isolated elements which can be tuned to the need. A constant head reservoir is usually added along the circuit to regulate the reference flow pressure. Bypass pipes, multiple test sections and several sensors can be combined in their main set-up to add flexibility and adjustability.

A possible layout is shown in figure 24 with the piston (P) displaced periodically by a rotating cam, driven by an electric motor, whose profile generates the desired flowrate waveform. If the walls of the ventricular chamber (VC) are stiff enough to prevent appreciable deformation, the same waveform is transferred to the inflow of the circuit (In). In a closed tank, the pressurized air (the Windkessel W) gives the compliance and the fluid therein the inertia of the circulation, while several lamination valves provide the resistance and allow a partial or total bypass of the Windkessel. The adjustable tank (T) fixes the flow pressure at point A before it returns to VC. If needed, an auxiliary test chamber (TC) can be added along the circuit to test vascular grafts or other circulation tracts.

Figure 24. Schematic of a pulse duplicator circuit with an electric motor driving a cam and a piston (P) which enforces in a ventricular chamber (VC) the desired waveform for the flow rate. The VC is fitted with a left ventricle phantom with an aortic valve which allows the flow only into the circuit (In). The flow impedance of the circuit can be adjusted through the lamination valves ${\rm {R}}_{1}$, ${\rm {R}}_{2}$ and ${\rm {R}}_{3}$ and the Windkessel W. A constant head tank (T) fixes the flow pressure at point A of the circuit. If needed, and auxiliary test chamber (TC) can be added to perform different measurements.

A duplicator can be further enriched with additional elements to reproduce more complex phenomena; in figure 25(a), we report the layout from Querzoli et al. (2016) in which the piston, driven by a computer-controlled linear motor, forces a plenum which drives a test chamber (AC) containing a deformable aortic root with the two coronary arteries. The latter cross an auxiliary chamber (CC) pneumatically connected with the plenum (PL); therefore, when its pressure increases (systole), the same pressure in CC chokes the coronaries and prevents their perfusion from the aorta (figure 25b). However, during diastole, the pressure drops in the PL and CC, the coronary tubes open and the pressurized deformable aorta can feed the coronaries (figure 25c): this set-up therefore reproduces the pressure tissue effect which allows the myocardium perfusion only during diastole when it is fully relaxed.

Figure 25. (a) Schematic of pulse duplicator circuit similar to that of figure 24 except for a piston driving a dummy plenum (PL) which feeds the aortic chamber (AC) and controls the coronary chamber (CC) in such a way to close the coronary circuit when the aorta is fed. The circuit has two constant head tanks (${\rm {T}}_{1}$ and ${\rm {T}}_{2}$) and the rest of the symbols have the same meaning as in figure 24. (b,c) Experimental measurement of the velocity magnitude and streamtraces in the aortic root and right coronary during systole and diastole, respectively. Velocity magnitude ranges from $0$ to $1.5\ {\rm m}\ {\rm s}^{-1}$ respectively for colours from black to dark red. (Courtesy of G. Querzoli and images adapted from Querzoli et al. 2016).

Although not devoted to visualizations or measurements, a peculiar set-up is that described by Leopaldi et al. (2018) and referred to as a ‘biosimulator’. A pulse duplicator with a computer-controlled piston pump is connected to the left ventricle of a porcine heart through its apex and a series of compliant tubes, Windkessels and resistances compose the preload and afterload units which feed the left atrium and receive the ventricle output. This apparatus is designed to train cardiac surgeons to perform off-pump transventricular mitral valve repairs and can be used for other types of structural heart disease or to evaluate different intracardiac devices. The biosimulator is complemented with an echocardiographic probe, placed externally on the left atrial posterior wall, and an intracardiac 2-D endoscope into the left atrium: the former provides the surgeon with the imaging available during an actual intervention while the latter gives the real situation during the simulation for comparison.

6.1.4. Measurement techniques

Apart from standard flow rate and pressure evaluations, quantitative measurements of the flowfield inside the models are performed using different methods which can be classified as optical and non-optical, with the former requiring transparent phantoms and duplicator test sections to visualize the flow in the region of interest.

The most used optical techniques are particle image velocimetry (PIV) and particle tracking velocimetry (PTV) which are very similar except that the former measures the ensemble displacement of a group or particles while the latter tracks particles individually. In both cases, the fluid is seeded by neutrally buoyant, microscopic ($10\text {--}50\ \mathrm {\mu }{\rm m}$) particles with a Stokes number small enough to behave as passive tracers and closely follow the fluid velocity (figure 26a). A section of the flow is then illuminated and recorded by a perpendicular camera so that particles (or particle clouds) can be tracked on the plane (figure 26b). For velocities of ${\approx }O(1\ {\rm m}\ {\rm s}^{-1})$, high frame rates (${>}O(500\text {--}1000$) f.p.s.) and intense light sources (usually lasers) are needed to correctly capture particle positions on each image. Furthermore, correlation algorithms on subsets of pixels across the frames and sophisticated processing techniques are employed to filter out the particles crossing the plane or other artefacts. Once the trajectories (PTV) or the displacements (PIV) have been identified, their time derivation yields the two components of the velocity field on the illuminated plane (figure 26c).

Figure 26. Example of image processing for cardiovascular flow measurements: (a) instantaneous frame capturing the position of the particles; (b) trajectory elements computed by correlating particles position across adjacent frames; (c) velocity vectors, interpolated on a regular grid, computed by time differentiation of the trajectories; (d) instantaneous position of the structure interface computed from the image in panel (a).

A bonus of this technique is that the same algorithms used to identify the particles can be employed to track the deforming boundaries (if the refractive index is not perfectly matched); in this way, not only the velocity field but also the structure kinematics can be determined (figure 26d).

A limitation of the above approach is that the plane section must cut the flow in a region where the velocity is mainly 2-D since otherwise, most of the particles would cross the illuminated plane only for a few frames and the velocity components could not be determined. If the flow has an inherent 3-D structure, more sophisticated techniques can be used, although they entail higher experimental complexity and data processing cost.

For example, the periodic character of cardiovascular flows can be exploited to reconstruct full 3-D flow fields by measuring in phase two velocity components on perpendicular planes over multiple cycles. Fortini et al. (2013) used particle-tracking algorithms to measure 2-D velocity fields over $24$ parallel planes spaced by $2.2$ mm acquired over $50$ cycles; by combining measurements on planes parallel and perpendicular to the ventricle long axis, they were able to reconstruct the 3-D velocity field inside the ventricle, as shown in figure 27. It must be noted, however, that this approach requires the cycle to cycle variations to be small enough for sections acquired from different cycles to be averaged and combined together into a single field.

Figure 27. PTV experimental reconstruction of iso-surfaces of vertical velocity ($+0.2\ {\rm m}\ {\rm s}^{-1}$, red; $-0.2\ {\rm m}\ {\rm s}^{-1}$, blue) at $0.13T$ (a); $0.18T$ (b); $0.23T$ (c); $0.29T$ (d) ($T$ is the heartbeat period). The triangulated structure represents the left ventricle wall, the large circle is the mitral ostium, the small circle is the aortic ostium. (Images adapted from Fortini et al. (2013), Courtesy of G. Querzoli).

The limitation of measuring only two velocity components on a plane can be overcome by Stereo-PIV which, using two cameras, records the same plane from two different angles and, relying on the same principles as stereoscopic vision, reconstructs the third velocity component. Issues related to perspective distortion and defocussing are handled by mapping functions which transform the acquired image pixels into real coordinate systems and they must be calibrated depending on the specific configuration.

A final possibility is tomographic PIV in which multiple cameras record from different angles multiple planar sections and the full 3-D velocity field in a volume can be reconstructed (Elsinga et al. 2006). Cost and complexity of the set-up is however an issue and also the computational need for data processing and field reconstruction might be a limiting factor. Another problem in the application of tomographic PIV to cardiovascular flow measurements is that the volume of interest is enclosed by the deformable structure whose refractive index is never perfectly matched with that of the fluid. This implies that the recordings from the various cameras have different distortion effects that have to be considered to decrease measurement errors and this further complicates data processing and its related computational cost.

Among the non-optical techniques, ultrasound imaging is the most used for cardiovascular flows since the same methods and devices are widely employed in clinical practice for diagnostics. The principle relies on the fact that sound waves have different propagation speeds in different materials and an interface between two media reflects part of the sound producing an echo. When a probe is placed on the body and emits a wave in a given direction, every tissue produces an echo which is captured by the same probe; from the delay of echoes, the depth (distance from the probe) of each tissue can be computed and an image is generated.

Blood is anechoic thus the basic method returns only the image of the solid tissues like myocardium, valves and vessels. Nevertheless, using Doppler effects, it is possible to convert echo frequency shifts into blood velocity although only as averaged values along the wave direction. Local velocity information can be obtained combining ultrasound imaging with PIV methods: in fact, gas microbubbles are opaque to sound waves and they behave as standard particles in optical methods and the former can be used to measure velocity components on the plane of the echographic image. Poelma et al. (2011) were able to measure the three-dimensional flow in a curved pipe while Versluis et al. (2020), using air microbubbles coated by biocompatible shells, were able to avoid their fast dissolution and used this method for in vivo medical imaging and diagnostic.

In some cases, even the actual apparatus for in vivo diagnostics can be used for flow measurements; Montalba et al. (2018) report of multiple MRI scans of the pulse duplicator flow in a realistic thoracic aortic phantom with the aim of assessing the variability of peak flow, mean velocity, stroke volume and wall shear stress measurements obtained from the real device under different conditions of spatial and temporal resolutions.

6.2. Computational modelling

From the discussion of the previous sections, it appears that to faithfully duplicate the dynamics of the heart, the design of a laboratory experiment entails the fulfilment of many requirements. Some of them come from the Buckingham $\varPi$-theorem, which imposes the matching of dimensionless parameters in the modelled and real phenomena, other constraints are instead due to practical considerations, like the transparency of working fluid, phantom and duplicator test section if an optical measurement technique has to be used.

In fact, any real laboratory experiment can reproduce only some cardiac features while others are only approximated or neglected. For example, using a deformable phantom in a pulse duplicator set-up, material properties and local thicknesses can be adjusted so as to replicate the correct kinematics or even the celerity of travelling waves although the orthotropic and hyperelastic mechanics of the tissues is lost. Replicating the active contraction of the tissues with the correct delays and preferential directions while maintaining the possibility of visualizing or measuring the flow is out of the question and the experimental set-up is, by necessity, incomplete.

In principle, computational models do not suffer from these limitations provided an adequate theoretical description of the system is available. For cardiovascular problems, however, computational power is the bottleneck since different complex subsystems have to be simulated including their two-way coupling. In fact, simulating the dynamics of the heart requires a formidable computational effort owing to the electrophysiology, the active contraction of the myocardium, the elastomechanics of the tissues (produced by internal and hydrodynamic loads) and the two-way coupled haemodynamics. Such simulations can be performed only nowadays thanks to the availability of efficient and inexpensive computers whose power seems ever growing. In fact, later in this paper, we will show some results of very recent numerical simulations of the whole heart which include all these effects although they are obtained relying on the latest GPU accelerated processors and massive parallel architectures which are not commodity computers yet.

Similar to experiments, computational models also start from the definition of a geometry and the same considerations as in § 6.1.2 apply here. However, standard STL files used for rapid-prototyping, with their triangle density inversely proportional to the local radii of curvature, can not be used directly for computations as their triangle distribution is usually too coarse and uneven (figure 28a) and they miss connectivity information. Another relevant difference is that the construction of a hardware phantom requires only the definition of the fluid/structure interfaces, while a computer model needs also the various tissue thicknesses (and fibre orientation) if the structure dynamics is considered. Accordingly, as a first step, a fine and homogeneous surface tessellation is performed which is also used as a starting point for the volume discretization (figure 28b).

Figure 28. (a) STL file used for 3-D printing: red mesh for the right side of the heart, green mesh for the left side. (b) the same as panel (a) but with a finer and more regular surface meshing needed for numerical modelling and with the addition of the tissue thicknesses necessary for the structure dynamics.

With the full heart geometry at hand, the computational model can be developed starting from three main building blocks: a flow solver, a structure solver and an electrophysiology module. These elements are strongly connected, physically and computationally; therefore, the solution algorithms and their coupling (the results of each component feed the others) must be carefully designed for the feasibility of the simulation, in terms of stability, accuracy and computational cost.

6.2.1. Electrophysiology

In § 5, the excitability of general cardiac myocytes and their mechanical response has been described with the aim of giving some insights into the behaviour of the cardiac tissue; here we report how the behaviour of billions of individual cells can be homogenized to derive a continuum electrophysiology model.

All cells are electrically coupled through end-to-end junction disks. As a result, they preferentially transmit an electrical signal along their axes which is the cause of conductive anisotropy. Furthermore, specialized cells are embedded within the heart tissue and they are organized into distinctive structures with different dynamics such as one-dimensional fast bundles, a fine network of fast conduction fibres, which can be modelled as a 2-D layer, and 3-D muscular tissues (figure 29). The electrically excitable cells can however be modelled by an action potential based on the difference between intracellular and extracellular potentials, as in figure 14.

Figure 29. Representation of the three-level hierarchical model of the electrophysiology system with the one-dimensional fast bundles (orange) the 2-D Purkinje network (orange) and the 3-D myocardium (dark red) (courtesy of G. Del Corso). The picture has been obtained from a numerical simulation using $2500$ linear elements for the fast bundles, ${\approx }10^{5}$ triangles for the Purkinje network and up to $1.6\times 10^{6}$ tetrahedra for the myocardium. The integration time step was kept below $5\ \mathrm {\mu }{\rm s}$. For further details, see Del Corso, Verzicco & Viola (2022).

Within these assumptions, cardiac tissues can be represented as a continuum (syncytium) characterized by two overlapping media separated by a distributed cellular membrane. Their dynamics are governed by the ion fluxes across the myocytes membrane which can be described by reaction-diffusion equations for the intracellular and extracellular potentials coupled by ionic currents. A common formulation recasts the two equations into a generalized Ohm's law and the conservation of charge and current in the form (Clayton et al. 2011):

(6.1)\begin{equation} \left.\begin{gathered} \chi \left ( C_m \frac{\partial v}{\partial t} + I^{ion}(v,{\boldsymbol{s}})+I^{s} \right ) = \boldsymbol{\nabla} \boldsymbol{\cdot} ({{\boldsymbol{M}}^{int}} \boldsymbol{\nabla} v)+\boldsymbol{\nabla} \boldsymbol{\cdot} ({{\boldsymbol{M}}^{int}} \boldsymbol{\nabla} v^{ext}), \\ 0 = \boldsymbol{\nabla} \boldsymbol{\cdot} ({{\boldsymbol{M}}^{int}} \boldsymbol{\nabla} v + ({{\boldsymbol{M}}^{int}+{\boldsymbol{M}}^{ext}}) \boldsymbol{\nabla} v^{ext}),\\ \frac{\partial{\boldsymbol{s}}}{\partial t} = F(v, {\boldsymbol{s}}). \end{gathered}\right\} \end{equation}

Here, $v= v^{ext}-v^{int}$ is directly the transmembrane (or action) potential, $v^{ext}$ its extracellular counterpart, $\chi$ a surface-to-volume ratio of the cells and $C_m$ the membrane capacitance per unit area. Additionally, ${{\boldsymbol {M}}^{int}}$ and ${{\boldsymbol {M}}^{ext}}$ are intracellular and extracellular conductivity tensors whose rank is $1$, $2$ and $3$ respectively in one-, two- and three-dimensions and they become diagonal when expressed in the fibre-axis, sheet-fibre and cross-fibre reference frame. It is worth mentioning that, as shown by Loppini et al. (2018), the conductivity tensors not only are an intrinsic property of the cardiac tissues and their fibre orientation but they depend also on the local stress, thus further coupling electrophysiology with elastomechanics and haemodynamics.

Here, $I^{s}$ is an ionic current (per unit area) stimulus which triggers the depolarization throughout the heart; therefore, in physiological conditions, it is zero everywhere except for in the sinoatrial node. Possible exceptions are regions of intense stretch which can activate additional current contributions potentially dangerous for cardiac arrhythmias and other pathologies (Loppini et al. 2018).

In physiological conditions, the heartbeat is initiated from the SA node which is the main pacemaker; however, in case of its failure, the His bundle takes over and triggers ventricular contraction, although at a slower rate.

The above model could also account for the presence of an external pacemaker whose leads complement or replace the natural triggers: in this case, additional sites for $I^{s}$ initiation should be prescribed together with suitable activation rules.

Finally, $I^{ion}(v,{\boldsymbol {s}})$ is the ionic current per unit area crossing the cells membrane which depends on the specific type of myocyte and therefore on the heart region. The third equation of (6.1) defines the various cellular models through the state vector ${\boldsymbol {s}}$, containing tens of variables governed by ordinary differential equations, which model the ionic fluxes of all the species across the membrane to determine the total current. Common cellular models are ten Tusscher & Panfilow (2006) for ventricular myocardium, Courtemanche, Ramirez & Nattel (1998) for the atria and the bundles, and Stewart et al. (2009) for the Purkinje network.

The first two eqautions of (6.1) can be reduced to a single expression assuming that the extracellular and intracellular media have the same spatial anisotropy, which implies ${{\boldsymbol {M}}^{ext}} = \lambda {{\boldsymbol {M}}^{int}}$ with $\lambda$ a scalar quantity. The model, referred to as the ‘monodomain model’, then simplifies to a single equation:

(6.2)\begin{equation} \chi \left ( C_m \frac{\partial v}{\partial t} + I^{ion}(v,{\boldsymbol{s}})+I^{s} \right ) = \boldsymbol{\nabla} \boldsymbol{\cdot} ({{\boldsymbol{M}}} \boldsymbol{\nabla} v), \end{equation}

with ${{\boldsymbol {M}}}={{\boldsymbol {M}}^{int}} \lambda /(\lambda +1)$, which still needs to be completed by cellular models.

It is worth mentioning that cellular models describe phenomena occurring at microscopic scales whose dynamics are therefore much faster than those modelled at a continuum level by partial differential equations; unfortunately, they are strongly coupled, thus the time integration of the electrophysiology model is very stiff and penalized by the fast cellular dynamics. However, cellular models are made of ordinary differential equations and most of them are linear thus making possible their explicit analytical solution within each time step; this approach was first proposed by Rush & Larsen (1978) showing enhanced stability properties.

Once the electrophysiology model, either (6.1) or (6.2), is solved numerically over the whole heart, the 3-D distribution of the transmembrane potential is obtained at every instant; using the relations suggested by Nash & Panfilov (2004) or Niederer & Smith (2008) (see figure 14), the internal active tension along the local fibres direction is then determined from which myocardium contraction is obtained.

6.2.2. Elastomechanics of tissues

Tissue dynamics can by modelled by the Lagrangian description of displacements $\boldsymbol {w}$ with respect to a reference configuration:

(6.3a,b)\begin{equation} \rho_s \frac{{\rm d}^{2} \boldsymbol{w}}{{\rm d}t^{2}} = \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{P}} + \rho_s \boldsymbol{f}_s, \quad {\boldsymbol{PF}}^{T}={\boldsymbol{FP}}^{T}, \end{equation}

where $\rho _s$ is the tissue density, ${\boldsymbol {P}}$ is the first Piola–Kirkhhoff stress tensor and $\boldsymbol {f}_s$ the vector of external volume forces. Note that the equation is in Lagrangian form; therefore, the time derivative of $\boldsymbol {w}$ is material (or total). The first Piola–Kirkhhoff stress tensor ${\boldsymbol {P}}={\boldsymbol {FS}}$ is related to the strain tensor ${\boldsymbol {E}} =[ \boldsymbol {\nabla } {\boldsymbol {w}} + (\boldsymbol {\nabla } {\boldsymbol {w}})^\textrm {T} + (\boldsymbol {\nabla } {\boldsymbol {w}})^\textrm {T} \boldsymbol {\nabla } {\boldsymbol {w}}]/2$ through a suitable strain energy density $\phi ({\boldsymbol {E}})$, which may also contain an additive term for the active contraction of fibres, when present (active stress formulation). Here, ${\boldsymbol {S}} = \partial \phi /\partial {\boldsymbol {E}}$ and ${\boldsymbol {F}} = {\boldsymbol {I}} - \boldsymbol {\nabla } \boldsymbol {w}$ (Gurtin 1981). For hyperelastic, orthotropic materials, Fung-type strain energy densities are often adopted in the form $\phi ({\boldsymbol {E}})= C(e^{Q}-1)$ being $Q= \sum _{i,j} b_{ij} E_{ii}E_{jj}$ or, for a bi-axial strain, $Q=\alpha E_{ff}^{2}+\beta E_{cc}^{2}$, with $E_{ff}$ and $E_{cc}$ the strain tensor components in the fibre and cross-fibre directions and $C$, $\alpha$, $\beta$ adjustable parameters (Fung 1993).

Rather than including the active contraction in ${\boldsymbol {P}}$, a possible alternative is to write the deformation gradient as the product ${\boldsymbol {F}} = {\boldsymbol {F}}_e{\boldsymbol {F}}_a$ (Kröner–Lee decomposition), with ${\boldsymbol {F}}_e$ the elastic part and ${\boldsymbol {F}}_a$ the virtual distortion due to the active contraction which is prescribed by a constitutive relation: this is referred to as active strain formulation and the reader is referred to Ambrosi et al. (2011) for a more extensive description.

The above elastomechanics model can be integrated numerically by finite-element methods, although the large displacements form of ${\boldsymbol {E}}$, the time-dependent dynamics and a number of elements of the order of hundreds of thousands to discretize the whole heart entail excessive computational costs.

A possible emerging alternative for the modelling of tissue dynamics is based on the interaction potential approach (Fedosov 2010), in which all tissues are described by triangles or tetrahedra and the mass is distributed among the element nodes which are connected by nonlinear springs (figure 30a). All internal loads and constraints, like incompressibility or surface preservation, are expressed by a combination of potentials $\varPhi$ whose gradient yields the internal forces (de Tullio & Pascazio 2016; Viola, Meschini & Verzicco 2020). Indicating by ${\boldsymbol {x}}_i(t)$ the instantaneous position of the $i$th node and by $m_i$ its mass, second Newton's law of motion reads

(6.4)\begin{equation} m_i \ddot{\boldsymbol{x}}_i +\gamma_i \dot {\boldsymbol{x}}_i = {\boldsymbol{f}}_i^{I}+{\boldsymbol{f}}_i^{A}+{\boldsymbol{f}}_i^{E}, \quad {\rm with} \ {\boldsymbol{f}}_i^{I} ={-}\boldsymbol{\nabla} \varPhi({\boldsymbol{x}}_i). \end{equation}

In the above equation, ${\boldsymbol {f}}_i^{A}$ is the force produced by the active contraction of the fibres and ${\boldsymbol {f}}_i^{E}$ the external and haemodynamic forces. The term $\gamma _i \dot {\boldsymbol {x}}_i$ gives an internal friction which, being dissipative, can not be expressed in potential form.

Figure 30. (a) Sketch of a triangulated discretization of the tissue with fibres: the active contraction force ${\boldsymbol {f}}_i^{A}$ is aligned with the local fibre direction, the internal force ${\boldsymbol {f}}_i^{I}$ is given by the vector sum of all the partial forces (${\boldsymbol {f}}_i^{I'},\ldots$) of the edges concurring to the $i$th node of mass $m_i$. Here, ${\boldsymbol {f}}_i^{E}$ accounts for the external forces. (b) Structure sample discretized by random triangulation to recover statistical homogeneity; (c) the same as panel (b) but with regular triangulation which produces unphysical anisotropic directions along the aligned edges.

Figure 31. Sketch of possible spatial discretizations for a left atrium/mitral valve/left ventricle assembly: (a) body conformal mesh with the mitral valve open during diastole; (b) the same as panel (a) but with the mitral valve closed during systole; (c) the same as panel (a) but for a regular mesh to be used with immersed boundary methods.

This approach is already a standard in computer graphics (Nealen et al. 2006) and it is particularly attractive in the present context on account of its reduced computational cost compared with the standard finite element techniques. However, this is not a real continuum approach since structures are modelled as a discrete network of concentrated masses connected by (hyper)elastic elements. In this mesoscale description, the continuum features of the structure are recovered only at scales larger than the elements (figure 30b) and care must be taken not to introduce artificial anisotropy by a too regular discretization as in figure 30(c) (van Gelder 1998).

6.2.3. Haemodynamics

The mathematical model for blood dynamics is well assessed (but also the most computationally demanding) since it is given by the well-known Navier–Stokes equations which for an incompressible, viscous fluid read:

(6.5a,b)\begin{equation} \frac{\partial {\boldsymbol{u}}}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} ({\boldsymbol{uu}}) ={-}\frac{\boldsymbol{\nabla} p}{\rho} + \frac{1}{\rho}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\tau} + {\boldsymbol{f}},\quad \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{u}} = 0, \end{equation}

with ${\boldsymbol {u}}$, $p$ and $\rho$ the local velocity, pressure and density, $\boldsymbol{\tau}$ the viscous stress tensor and ${\boldsymbol {f}}$ the specific volume forces. The constitutive relation $\boldsymbol{\tau}=2 \mu \textrm {sym}(\boldsymbol {\nabla }{\boldsymbol {u}})$ is commonly used and the dynamic viscosity $\mu$ can be either constant, if blood is considered a Newtonian fluid, or a function of $\boldsymbol {\nabla } {\boldsymbol {u}}$, as in (4.1), for a shear-thinning fluid. In a relatively large system like the heart, the Newtonian fluid model is adequate and reliable although in some niche applications, for example haemolysis, the more realistic shear-thinning relation is preferable (De Vita, de Tullio & Verzicco 2016). If the viscoelastic properties of blood have to be reproduced, a separate set of differential equations must be solved to determine the time dependent viscosity field (Campo-Deaño et al. 2013).

Equation (6.5a,b) is solved with the boundary condition ${\boldsymbol {u}} = \dot {\boldsymbol {w}}$ (or alternatively ${\boldsymbol {u}} = \dot {\boldsymbol {x}}$ if the model (6.4) is used) applied at the fluid/tissue interface and this couples fluid and structure dynamics. In turn, once (6.5a,b) is solved, pressure and velocity fields are obtained within the domain of interest and from them, the hydrodynamic loads on a surface $S_i$ with local outward normal ${\boldsymbol {n}}$ are computed through

(6.6)\begin{equation} {\boldsymbol{f}}_i^{E} = \int_{S_i} ({-}p {\boldsymbol{n}}+\boldsymbol{\tau} \boldsymbol{\cdot} {\boldsymbol{n}})\,{\rm d}S. \end{equation}

These loads are then used to determine the dynamics of the tissues.

It is important to note that the evaluation of (6.6) and the integration of (6.5a,b) can only be performed numerically; therefore, a preliminary discretization of the fluid volume is required. Using standard techniques, with the mesh conforming to the boundaries, could be extremely difficult considering the complex geometrical configuration and its time dependence that would require grid morphing or remeshing. Even more troublesome is the treatment of heart valves since they are located in between different volumes which are connected when the former are open and separated when closed. This topological change of the fluid volume results in a singularity of body conformal meshes which render the equation integration terminally expensive and eventually impossible (Bavo et al. 2016).

A viable alternative is to decouple the system geometry from the volume discretization using non-body-conformal meshes and to replace the physical boundary conditions by space- and time-dependent volume forcings ${\boldsymbol {f}}$ mimicking the effects of the tissues on the flow. These schemes are referred to as immersed boundary (IB) methods and, for cardiovascular flows, they were originally proposed by Peskin (1972). IB methods allow the handling of complex geometry flows within the ease and efficiency of structured regular grids and they do not need the generation of body-fitted meshes that, for moving and deforming boundaries, entail an unbearable computational overhead. Furthermore, as the physical boundaries move across the computational nodes, the procedure is insensitive to the topological changes produced by valves opening and closing. However, the disadvantage of the IB methods is that as the valve leaflets approach each other, the fluid layer in between shrinks and the resolution of the underlying fixed mesh becomes insufficient. This problem is alleviated by using lubrication solutions and contact models whenever the fluid layer in between two immersed boundaries becomes thinner than $O(3\text{--} 4)$ mesh elements (Verzicco & Querzoli 2021).

Another related problem is that the absence of a boundary conformal mesh prevents the possibility of wall grid refinement which is necessary to properly resolve the boundary layer dynamics. To some extent, the missing grid stretching can be compensated by a finer regular grid but only up to a certain Reynolds number.

In fact, as $Re$ increases, boundary layers at the walls and flow scales in the bulk both decrease in size and, as the former shrink at a faster rate, eventually, they become the bottleneck. This is a serious limitation for the application of IB to high- $Re$ flows such as those encountered in automotive or aerospace problems; cardiovascular flows, however, evolve at Reynolds numbers of the order of a few thousands ($\approx$3000 for the left ventricle and $\approx$6000 for the aortic root), thus a brute force refinement of the regular grid is still possible within an affordable computational cost. As an aside, we note also that the flow scales of cardiovascular flows are not as small as the Reynolds number would suggest since the peak flow rate is sustained only for a small fraction of the cardiac cycle (see figure 23b) and turbulence can not fully develop up to the smallest allowed scales: this makes the IB methods an even more computationally efficient choice. For a more detailed discussion of this point, the reader is referred to the recent review by Griffith & Patankar (2020) entirely devoted to the application of IB methods to biological flows.

6.2.4. Fluid–structure–electrophysiology interaction

In principle, the double time integration of (6.4) or (6.3a,b) yields the deformation field and the configuration of all heart tissues whose wet surface velocity can be used to impose the boundary conditions in (6.5a,b) for the haemodynamics while the volume, with its hierarchy of one- and two-dimensional structures, provides the domain of the electrophysiology equations. However, evolving in time, the structure entails the knowledge of the fluid loads whose behaviour depends, in turn, on the structure configuration itself. This ‘chicken and egg’ problem is known as a fluid/structure interaction (FSI) which, for cardiovascular problems, is further complicated by the addition of the electrophysiology (FSEI).

Given the strongly interconnected dynamics of the three systems, it is not obvious whether they must be solved simultaneously or sequentially and, in the latter case, from where to start.

In the first case, all models are solved together as a single large dynamical system through an iterative procedure (strong coupling), while in the second, the output of each model is used as the input for the successive one (loose coupling) in a given order. Strong coupling is stable and robust although it is computationally expensive since it requires several iterations among the solvers (de Tullio & Pascazio 2016) for every time step ($\Delta t$). However, loose coupling needs only one model solution per $\Delta t$ although it is prone to become unstable and poses strong limitations on the time step size.

Which of the two approaches is computationally more efficient depends on the specific problem and its physical parameters; Borazjani, Ge & Sotiropulos (2008) have shown that, when added mass effects are prevalent, loose coupling FSI becomes unstable and strongly penalizes the integration time step. In fact, in a heart model, valve dynamics is dominated by added mass effects as the intrinsic mass of valve leaflets is negligible with respect to that of the displaced fluid and this would suggest that strong coupling FSI is needed. However, Viola et al. (2020) have shown that, for the heart, the maximum $\Delta t$ is anyway limited by the fast electrophysiology dynamics and the high elastic frequencies of the stiffened myocardium during contraction. As a consequence, within these limitations, a loose coupling FSI maintains the integration stability and provides computational time savings of the order of $50$–$70\,\%$ while giving the same results as strong coupling FSI.

In contrast, the model of de Tullio et al. (2009) considered the dynamics of a bi-leaflet mechanical valve in a rigid aortic root and only the fluid/structure interaction time scales were present; according to the above arguments, strong coupling FSI was used for the time integration with a number of iterations between fluid and structure solvers of the order of $4$ each step.

6.2.5. Boundary conditions

As noted in § 2, human circulation is closed thus only boundary conditions at the blood/tissues interfaces are needed to determine the haemodynamics of the heart. However, modelling the whole circulation with the same detail as the heart is practically impossible, and therefore the discretization is truncated at some level of the main vessels and the missing part is replaced by additional models (figure 32). These usually take the form of boundary conditions mimicking resistance, inertia and compliance of the eliminated circulation circuits and their choice is of paramount importance for the quality of the results (Vignon-Clementel et al. 2006; Esmaily Moghadam et al. 2011).

Figure 32. Schematic representation of the total circulation (a) and of the heart with the necessary boundary conditions (b).

In the simulation of an open-loop model, like a tract of aorta or even a single side of the heart, inflow and outflow boundary conditions must be provided separately. The former are usually easier and are imposed as time dependent velocity profiles (plug flow, parabolic and Womersley solutions for rigid or deformable vessels) over all inlet sections. Outflow conditions are more delicate and they usually result in a set of ordinary differential equations, if inertia, resistance and compliance are lumped in a zero-dimensional model, or partial differential equations if the same parameters are distributed on a one-dimensional domain. Either way, the resulting relations can seldom be integrated analytically and they need numerical solutions sometimes implying iterations with the main model (Marsden 2014).

In closed-loop models, inflow and outflow are connected (figure 32) and the correct pressure/velocity conditions come from suitable coupling relations. This is usually achieved by lumped parameter networks whose structure and features can be tuned to obtain the desired inlet and outlet waveforms. In fact, these are just virtual Windkessels, analogous to those introduced in § 6.1.3, which are implemented by analytical and numerical relations and therefore can replicate functions much more complex than the hardware counterparts.

A possible problem is that boundary conditions coming from low-dimensional virtual Windkessels miss some information needed by the 3-D main model like the velocity distribution over the inflow sections or even at the outflow sections in regions of backflow; also in this case, iterative coupling procedures are needed which further increase the computational load.

6.2.6. Reduced models

Running a high-fidelity heart model accounting for all the features described in §§ 6.2.1–6.2.5 is a formidable computational task which can not be accomplished without high-performance parallel computing, possibly accelerated with the latest GPU technologies (Viola et al. 2022). Such extensive undertakings, however, are necessary only in exceptional cases while less demanding or ‘reduced’ models can be used to address specific problems or to perform routine activities. The complexity can be reduced in many different ways, depending on the question to be answered and the available computational resources.

A first simplification can be obtained by considering only a subset of the complete model, thus reducing both the geometrical complexity and the volume of the computational domain: common examples are a single heart chamber or the initial portion of the ascending aorta with natural or prosthetic valves, as shown in figure 33. In some cases, the set-ups can be further simplified by neglecting the deformability of some parts and replacing them by rigid structures; of course in each case, appropriate boundary conditions and fluid/structure interaction procedures must be selected to cope with the specific system dynamics.

Figure 33. Examples of reduced models: (a) deformable aortic root with coronary arteries and bi-leaflet mechanical aortic valve. The contours are the pressure field during diastole on the symmetry plane of the valve, showing the perfusion of the coronaries (adapted from de Tullio et al. 2011a). (b) Deformable left ventricle and mitral valve anchored to a rigid structure. The contours are the velocity magnitude with vectors on the symmetry plane, during the diastolic early wave (adapted from Meschini et al. 2018). (c) Complete deformable left heart with mitral and aortic valves and with the aorta up to the descending tract. The contours are the velocity magnitude with vectors on the symmetry plane, during the diastolic early wave (adapted from Viola, Meschini & Verzicco 2021).

A somewhat subtle configuration is that of a rigid aortic root with a mechanical aortic valve whose geometry is completely determined by the angular position of the rigid leaflets. This supposedly simple set-up does not need a structure solver since the balance of angular momentum about the hinge regulates angular velocity and position of the leaflets which, in turn, determine the blood domain configuration. However, an incompressible flow in a fully rigid domain does not have any degree of compliance resulting in an extremely stiff numerical problem in which, at every time step, flow and structure must be in perfect dynamic equilibrium for the numerical model not to diverge. The vanishing mass of the rotating leaflets, dominated by added mass effects, exacerbates the stability problems especially in the fast dynamic phases. This configuration was analysed by de Tullio et al. (2009) who reported, during the valve closure, leaflets angular velocities up to $60\ \textrm {rad}\ \textrm {s}^{-1}$. The accurate and stable integration of this fast dynamics implied time steps as small as ${\approx }3\ \mathrm {\mu }\textrm {s}$ combined with a strong-coupling FSI algorithm which required approximately $4$ iterations between fluid and solid solvers per step.

Apart from the specific set-up, heart valves are always a bottleneck for digital cardiovascular models as the flow through them has the highest velocities, their leaflets, both native and artificial, have negligible thickness and inertia and yet, when closed, they support the largest pressure differences. For these reasons, when the focus of the model is not the valve FSI or the detailed features of the flow through them, they are replaced by simpler models reproducing only the main features.

The easiest approximation is the check valve or ‘diode model’ in which a time dependent resistive volume is positioned within the valve ostium such that its resistance is vanishing when the valve is open and tends to infinity when closed. This is implemented by adding an extra term to (6.5a,b) in the form ${\boldsymbol {f}}' = -{\boldsymbol {u}} \delta ({\boldsymbol {x}})/K(t)$ with $\delta ({\boldsymbol {x}})$ giving the volumetric localization of the forcing and $K(t)$ a Darcy-like porosity, which modulates the time opening ($K(t) \rightarrow \infty$) and closing ($K(t) \rightarrow 0$). This model is computationally inexpensive although it relies on user-defined parameters for the timing and sharpness of valve effects which need additional tuning.

Furthermore, ‘diode models’ completely miss the peculiar flow patterns produced by the valve leaflets which sometimes are needed to make quantitative analysis of the downstream phenomena; an alternative procedure consists of having the actual valve geometry but driving its dynamics either by an integral momentum balance (Celotto et al. 2019) or even by a prescribed kinematics (Ge & Sotiropulos 2007). These methods yield a more realistic flow across the valve although they tend to produce spurious vorticity owing to unavoidable mismatches between the assigned leaflet position and its correct counterpart coming from local FSI equilibrium. More details on the specific topic of heart valves modelling can be found in the comprehensive reviews by Sotiropulos, Le & Gilmanov (2016) and Mittal et al. (2016).

Prescribing the kinematics of the tissues and using time-dependent configurations as boundary conditions for the flow is a strategy that can be used not only for the valves but also for the endocardium thus avoiding the integration of all the tissue dynamics, including its electrophysiology. The advantage of such a method is the possibility to perform realistic simulations at a fraction of the computational cost of the complete model; however, the results are determined by the prescribed configuration and therefore should be viewed more as a constrained flow configuration rather than as a real predictive result.

Nowadays boundary kinematics is often extracted from high-resolution four-dimensional medical imaging like MRI or CT scans thus making this technique particularly suitable for patient-specific analyses (Le et al. 2019; Collia et al. 2021). However, even the most advanced imaging devices can not give more than $\approx$30 images per heartbeat at spatial resolutions of approximately $3\times 3\times 3\ \textrm {mm}^{3}$ voxels and, therefore, space–time interpolations are needed to feed the numerical simulations that are performed on finer meshes and advanced in time with smaller time steps.

A further model reduction is more or less consciously operated when the haemodynamics is integrated on a mesh coarser than required by the flow Reynolds number. In fact, in this context, the size of the smallest flow scales is difficult to estimate as classical formulae for Kolmogorov scale ($\eta _K$) and boundary layer thickness are valid for statistically steady conditions and rigid boundaries at rest, respectively. None of these hypotheses apply to cardiovascular flows that are pulsatile, evolve within moving, deforming boundaries and are driven by an alternation of favourable stabilizing pressure gradients and unfavourable counterparts. Considering that in each heart district, the peak blood velocity is attained only for a small fraction of the heartbeat (figure 23) and assuming that the energy does not cascade up to the smallest scales, we can conjecture that the minimum spatial resolution is a multiple (say $3$–$4$) of $\eta _K$: with the physiological heart dimensions and flow parameters, we get $\approx$0.2 mm for the left ventricle during diastole and $\approx$0.1 mm for the aortic root during systole.

Discretizations with such fine meshes are seldom used and resulting haemodynamics is either stabilized by ad hoc numerical procedures (upwind discretizations, numerical viscosity, filtering) or by turbulence modelling (LES or even RANS) that are also developed in different contexts and within hypotheses not fulfilled by cardiac flows.

An even more extreme strategy to simplify the cardiac system can be the replacement of the Navier–Stokes equations with zero-dimensional Windkessels which provide the pressure exerted on the endothelial walls of the heart chambers (Westerhof, Lankhaar & Westerhof 2008): in this way, a coupled system of nonlinear partial differential equations is reduced to a handful of ordinary differential equations that can be integrated within negligible times.

Further complexity reduction can be obtained by parametrizations of the solid structure, for example, extracting from full-order-model solutions the modes of the tissue dynamics and constructing with them an observation-based model (Pfaller et al. 2020). In fact, using proper-orthogonal decomposition (POD), the modes are ordered according to their energy content and a few modes are enough to capture the relevant kinematics of the structure.

A complementary possibility is to consider only a subset of the heart chambers and describe them by simplified parametric shapes with a reduced number of degrees of freedom. An example is the model of Moulton, Hong & Secomb (2017) who considered a left ventricle of oblate spheroid shape with a myocardium of non-uniform thickness constrained by volume conservation. Aortic and mitral valves were modelled by resistive elements and the circulations by lumped-parameter Windkessels; by properly calibrating the model parameters, realistic distributions of myocardial stresses are obtained within very limited computational costs.

7. Cardiovascular haemodynamics

7.1. Physiological conditions

The periodic beating of the heart produces a pulsatile blood flow whose correct direction, from atria to ventricles and then to the main arteries, is assured by the atrioventricular and semilunar valves. These are passive structures which move mainly owing to fluid inertia and pressure differences produced by the active contraction of the myocardium and the energy stored by the elastic stretching of deformable tissues.

Figure 34. Flow structure in the heart chambers during diastole (a,c) and during systole (b,d); panels (a,b) report a planar section through the left heart while panels (c,d) through the right heart. The colour contours represent velocity magnitude (courtesy of F. Viola). The pictures have been obtained from a numerical simulation of the complete heart (electrophysiology, elastomechanics and haemodynamics), on a mesh of approximately $2\times 10^{8}$ nodes and an integration time step of $5\ \mathrm {\mu } \textrm {s}$. The numerical method is that described in Viola et al. (2022).

Figures 34(a) and 34(b) report snapshots, obtained from numerical simulations (CFD), of the blood flow in the left heart during two representative phases of the cycle which can be better understood with the help of the top panel of figure 9. The flow of figure 34(a) occurs in the region f of the Wiggers diagram and it shows the rapid passive filling of the left ventricle during early diastole (E-wave) in which both LV and RV chambers are relaxed and blood mainly flows because of its kinetic energy acquired in the pulmonary veins: this is confirmed by figure 35(d) showing CFD colour contours of pressure in the upper left heart region with velocity vectors overlaid. Mitral leaflets are open and the emerging fast mitral jet drives an intense ventricular recirculation which is believed to be beneficial for the cyclic sweeping and ‘washing’ of the corrugated endocardium (figure 20d) and to prevent the formation of stagnant flow regions at the ventricle apex which are prone to produce clots. The last argument is made quantitative by the E-wave propagation index (EPI) defined as the ratio of the mitral propagation length ($L_{mj}$) and the ventricle major axis length ($L_{lv}$) at maximum expansion: $L_{mj}$ can be estimated as $\int _0^{T_E} V_M(t)\,\textrm {d}t$, with $T_E$ the duration of the E-wave and $V_M(t)$ the jet velocity at the mitral leaflet tips which, together with $L_{lv}$, can both be evaluated by simple ultrasound Doppler imaging. Harfi et al. (2017) found that healthy left ventricles are characterized by $EPI \geqslant 1$ confirming that the mitral jet has to have enough momentum to penetrate up to ventricle apex.

Figure 35. Colour contours of instantaneous pressure overlaid with velocity vectors for the flow through the mitral and aortic valves during various phases of the cycle: (a) isovolumic systole, the mitral valve is closing and the aortic valve is still sealed; (b) isotonic systole, the mitral valve is sealed and the aortic valve open; (c) incipient diastole, the mitral valve is opening and the aortic valve is closed; (d) early diastole (E-wave), the mitral valve is fully open and the aortic valve is sealed. The pressure range has been adjusted in each panel to evidence the pressure variation across the valves. (Courtesy of F. Viola.) The data are obtained from the same numerical simulation as described in figure 34.

The large ventricular vortex is subsequently reinforced by a second mitral jet produced by the active atrial contraction (A-wave) which completes the filling of the ventricle to its maximum capacity.

The strength of E- and A-waves is measured through the peak velocity magnitude of the produced mitral jets and their ratio ($E/A$) is assumed as a diagnostic indicator of the overall cardiac health; in physiologic conditions, it results in $1 \leqslant E/A \leqslant 2$ while higher values indicate the ‘diastolic disfunction’ which is the inability of the atrium to contract or of the ventricle to properly relax in both cases leading to improper filling of the latter chamber (Mitter, Shah & Thomas 2017).

Meanwhile, the aortic valve is sealed by the counterpressure between the blood in the tensioned aorta and that in the relaxed left ventricle (figure 35c,d); this avoids aortic regurgitation and allows the feeding of systemic circulation also during diastole.

Figure 34(b) shows the same side of the heart during early systole, region c of the Wiggers diagram, with an opposite dynamics with respect to figure 34(a). In fact, the vigorous active contraction of the ventricular myocardium causes the blood pressure to rise sharply beyond the aortic value thus causing the opening of the aortic valve and driving an intense jet into the ascending aorta (see also figure 35b). During systole, blood flows into the aorta at a faster rate than it is distributed to the systemic circulation and the excess mass is stored in the aorta itself which slightly increases its volume (approximately $5\,\%$) and stretches the walls thus converting kinetic energy and pressure of the stream into elastic potential energy which is gradually returned to the flow during diastole.

Thanks to the different electrophysiology timings, while the left ventricle squeezes, the atrium is fully relaxed and the counterpressure between the two chambers closes the mitral valve, again avoiding regurgitation (figure 35a,b).

Figure 9 evidences that this part of the cycle produces, between left atrium and ventricle, the largest pressure difference of all blood circulation; in physiological conditions, it can exceed $110$ mmHg (${\approx }$14 600 Pa) which, on an ostium of diameter ${\approx }2.5$ cm, generates a force of $\approx$7 N to be countered by the sealed mitral leaflets.

A common belief is that part of this load is sustained by the pulling action of cordae tendineae which connect the free edges of the valve leaflets to the papillary muscles; however, a closer look at the system dynamics shows that it is not the case. In fact, during systole, the entire left ventricle contracts and its major axis shortens by up to $25\,\%$ thus decreasing the distance between leaflet edges and papillary muscles. Since the cordae tendineae do not contract, their pulling action ceases in the middle of systole potentially allowing mitral valve prolapse. However, we have seen that papillary muscles are reached directly by the Purkinje fibres whose conduction velocity of the activation potential (${\approx } 4\ \textrm {m}\ \textrm {s}^{-1}$) is faster than that of the myocardium muscle (${\approx }1\ \textrm {m}\ \textrm {s}^{-1}$); this implies that the former are fully contracted before the ventricle starts shortening (figure 36b) and the tethering of the cordae tendineae is exerted while the mitral valve is closing. Once the valve is sealed, the curvature of its coapted leaflets balances the pressure loads through membrane tension, as suggested by the Young–Laplace equation. At the leaflets’ root, the load is then transferred to the stiff fibrous mitral annulus and the pulling action of the cordae is no longer needed (figure 36c). This mechanism resembles the construction of an arch or a dome which needs the support of a scaffold only until the keystone is in place. Direct confirmation of the above dynamics comes from the experience of heart surgeons who refer that once the mitral valve is closed, the tension of the cordae tendineae is indeed absent (R. De Paulis Personal Communication).

Figure 36. Dynamics of a mitral valve in a model left ventricle with cordae tendineae and papillary muscles: (a) diastole; (b) early systole when the papillary muscles contract before the ventricle; (c) late systole. The two mitral valve leaflets are represented with different colours to show the coaptation during closure. Note that the leaflets’ dynamics is computed as part of the solution and not prescribed, therefore once they have coapted (panel b) they can still drift normally or slide tangentially (panel c) depending on the loads. The bold arrows indicate the direction of the blood stream to and from the ventricle. (Figure adapted from Meschini, Mittal & Verzicco 2021.)

Once more, this example clearly illustrates that the various heart systems are extremely coupled and synergistic; therefore, considering only a single aspect can give an incomplete or misleading picture.

Despite similar dynamic mechanisms, mitral and aortic valves are quite different, both in terms of geometry and operating conditions. The former is made up of an anterior (AL, towards the interventricular septum) and a posterior leaflet (PL), each one having three main segments; the valve is asymmetric, the leaflets have opposite antero-lateral curvatures, and different mobility and length. This peculiar structure produces a jet whose axis points in different directions during diastole: initially, the jet axis points to the posterior ventricle wall (figure 35c) then the jet rotates (figure 35d) and, at the peak of the E-wave, it is directed towards the apex. The initial development of the E-wave jet results in the formation of a toroidal vortex ring whose interaction with the posterior ventricular wall yields the aforementioned large-scale recirculations (see also figure 27): there are attempts to connect the mitral vortex ring features to the efficiency of the ventricle although there are no unique views and the matter is still debated (Kräuter et al. 2020).

The aortic valve, with untethered leaflets and threefold symmetry, is comparatively simpler than the mitral one although its dynamics is still very involved. The valve opens at the beginning of systole, when ventricular pressure exceeds the aortic value, and the resulting flow is a high-speed jet into the ascending aorta (figures 34b and 35b). The vorticity layers, detaching from the leaflets’ tips, separate the jet core from the region of the aortic root sinuses although the shear stresses at the interface drive a stationary toroidal structure referred to as the aortic sinus vortex or ASV. The physiological function of this structure is debated although there are indications that it helps to drive blood towards the coronary ostia for their perfusion and produces a smoother valve closure during diastole (Sotiropulos et al. 2016). The shear layer instability and its FSI with the valve induces intense leaflet fluttering which modulates, by high-frequency perturbations, the jet and triggers its transition to turbulence at the end of systole when an intense adverse pressure gradient reverts the flow and closes the valve.

In addition to the small-scale fluctuations, a large-scale swirling flow is also formed in the aorta (Azarine et al. 2019), which maintains the coherence of the stream with a constant shear stress level in the aortic arch and descending tract; the latter reduces the accumulation of low density lipoproteins and prevents the formations of plaques.

The presence of a precise shear stress range is a common condition to all the cardiovascular endothelium; every region is characterized by physiological values sensed as an indicator of regular heart functioning. In contrast, all variations become signals of abnormal conditions to which the tissue reacts with long-term remodelling (Russo et al. 2020), an irreversible condition which can trigger a negative feedback loop leading to heart failure; the adverse dynamics can originate not only by degenerative diseases or ischaemic episodes but also by surgical procedures or implantation of prostheses, all of which are events capable of modifying the native flow pattern.

The right heart is topologically similar to the left counterpart although, supporting only $15\,\%$ of the total workload, it shows significant differences. The right ventricle has a strongly 3-D structure, it is ‘wrapped’ around the anterior side of the left ventricle and is made of a lower muscular part, the sinus, and an upper outflow tract, the conus, which connects with the common trunk. The systolic contraction is complex and influenced by the left part: the highest contribution comes from the shortening of the long axis, then from the movement of the free wall towards the septum and finally from the bulging of the latter into the right ventricle cavity (Gaine, Naeije & Peacock 2014). The resulting flow is highly 3-D and a single planar section gives only an approximate view of its structure. Figure 34(c) shows a snapshot of the flow during early diastole when the pressure is in equilibrium throughout the right heart (figure 9 regions u and v): a cross-section of the superior vena cava shows the feeding of the right atrium with blood flowing inertially, through the tricuspid valve, to the ventricle. A similar flow is present also in the inferior vena cava; however, this flow is not visible since it is not captured by the section plane. In this phase, the pressure in the pulmonary arteries exceeds that of the ventricle and this keeps the pulmonary valve closed. In contrast, when systole starts (figure 9 regions m and n), the active isovolumic contraction increases pressure in the ventricle and the first effect is to close the tricuspid valve to avoid regurgitation into the atrium. As the contraction progresses, the pressure keeps increasing until the common trunk value is exceeded and the pulmonary valve opens (figure 34d) to supply blood to the lungs.

The tricuspid and pulmonary valves are quite different with the latter similar to the aortic valve while the former has three main leaflets with the free margins connected by an intricate web of cordae tendineae to papillary muscles. The right valves are less loaded and damaged than those of the left side because the low impedance pulmonary circulation works on pressure differences of only $20$–$25$ mmHg ($\approx$3000 Pa). For this reason, any right heart impairment, compensated by the stronger left side, goes initially unnoticed and allows the disease to progress and worsen because of slow remodelling processes. At later stages, when symptoms appear, for example, in the form of pulmonary hypertension, the disease, either caused by valvular, ventricular or arterial problems, is usually irreversible and with negative prognosis. The right heart has long been neglected as it was viewed as a trivial circuit bringing blood to the lungs and therefore its haemodynamics was not interesting: nowadays, the understanding has totally changed and the pulmonary circulation is considered a key synergistic component of the correct heart pumping function (Gaine et al. 2014).

7.2. Flow scales and dynamics

From the above discussion, it is clear that classifying the heart haemodynamics into canonical flow behaviours is pointless since the blood motion shows bits of many classical problems but does not fit any of them.

The main point to be considered is the nature of the forcing which is periodic in time with short and concentrated active phases (figure 9). In fact, the steep pressure gradients driving blood motion are applied within two intervals of $\approx$100 ms in a cycle of $\approx$860 ms and the resulting pulsatile flow is strongly unsteady. For example, the mitral jet, with a peak velocity of ${\approx }0.8\ \textrm {m}\ \textrm {s}^{-1}$ and an ostium diameter of $\approx$2 cm, attains a Reynolds number of $\approx$4000, which for a classical jet would entail turbulent dynamics. In this case, however, the flow is accelerated by a favourable pressure gradient for $\approx$50 ms and then strongly decelerated, by the adverse counterpart, in the subsequent $\approx$50 ms. Accordingly, the initial smooth, accelerating flow becomes transitional during the deceleration when the destabilizing pressure gradient produces an explosive growth of small flow scales (Meschini et al. 2018). Nevertheless, this transition is also limited in time and followed by a long quiescent phase of $\approx$300–400 ms during which the flow settles down, small scales decay and only large-scale structures survive. In fact, the flow is neither laminar nor turbulent showing a transitional behaviour restricted to the short deceleration phase. When averaged over a complete period, the mean Reynolds number is $\approx$400 with approximately $3/4$ of the period in which the flow is not forced and its energy decays.

A very similar dynamics for the aortic root was recently reported by Nitti, De Cillis & de Tullio (2022) who investigated the flow downstream of several valves by direct numerical simulation. During peak systole, the highest Reynolds number of the cardiovascular system was attained ($Re\simeq 6500$) and yet the energy spectra were not characterized by the small-scale content of statistically steady turbulence at the same $Re$. The reason for this behaviour was identified in the stabilizing effect of the favourable pressure gradient which accelerated the flow from rest to the peak velocity. In contrast, during deceleration, the adverse pressure gradient destabilized the flow and triggered the forward energy cascade which was evidenced by the appearance of a short inertial range in the spectra. The duration of the forward cascade, however, was too short for the energy to reach the smallest Kolmogorov scale which, for those flow parameters, was estimated as $\eta _K \approx 40\ \mathrm {\mu } \textrm {m}$.

A similar comment applies to the flow sweeping the tissues, which has often been analysed by making analogies with boundary layers, pipe flows and jets impinging on walls. As an example, if we considered the flow in the aortic root at $Re\simeq 6500$ as the flow in a circular pipe, we should conclude that the flow has a sustained turbulence with a log-law for the mean velocity profile at the wall. Once again, we note that the dynamics is very different since the flow is pulsatile, the peak Reynolds number is sustained for less than $1/5$ of the period and for $2/3$ of it, the inflow is zero. A quick look at figures 34(a) and 34(b) confirms this scenario, with recirculations in the Valsalva sinuses and a massive separation on the convex side of the arch during systole (figure 34b), and a stagnant flow during diastole (figure 34a). Thus, in addition to the strongly unsteady outer flow, the dynamics is affected by the moving/deforming boundaries and by curvature effects. Moreover, the peculiar blood rheology (figure 11b) makes it difficult to compute wall shear stresses from standard Newtonian fluid relations. The same task becomes pointless in the ventricles on account of the actual endothelium morphology (figure 20d) which resembles more a compliant network of fibres than a solid interface.

7.3. Heart rate effect

For every heartbeat, the ventricles pump into the circulations a given blood volume (stroke volume, $SV$) which, multiplied by the heart rate ($HR$), yields the volumetric flow rate or cardiac output $CO=SV\times HR$. For a healthy adult of mass $\approx$70 kg, the mean resting stroke volume is $SV \simeq 70$ ml and the heart rate $HR \simeq 70$ b.p.m. which, in the appropriate units, yield $CO \simeq 5\ \textrm {l}\ \min ^{-1}$. However, depending on external conditions or physical exertion, $CO$ and $HR$ might both increase, the former up to $5$–$6$ times and the latter $2.5$–$3$ times thus entailing a $SV$ dependence on $HR$. Indeed, a heartbeat involves electrophysiology, active contraction, passive elasticity of tissues and haemodynamics, each one with different time scales whose interplay determines a nonlinear behaviour of $CO$. The cyclic filling/emptying of a ventricle can be discussed from figure 37(a) which is nothing but the time derivative of the left ventricle volume in the third panel of figure 9 (with a different time origin). The initial filling of the ventricle (E-wave) is inertial and yields approximately $70\,\%$ of the blood volume in a time $Et$. The remaining $30\,\%$ is given by the active atrial contraction (A-wave) whose duration is $At$ and is separated from the initial phase by a time interval $Dt$ referred to as diastasis: the sum of $Et+Dt+At$ gives the diastole duration. Chung et al. (2004) have analysed, for a cohort of young healthy volunteers, the dependence of the time intervals on $HR$ as a consequence of physical exercise and figure 37(b) shows that the period reduction is mostly obtained by a decrease of the diastasis which, for $HR \geqslant 80$ b.p.m., becomes negative. Indeed, in this regime, E-wave and A-wave partially overlap as confirmed by the curve for the flow rate at $HR=100$ b.p.m. of figure 37(a) (red curve). It is worth noticing that the stroke volume ($SV$) also changes showing an increase up to $HR \approx 130$ b.p.m. followed by a saturation up to $HR \approx 160$ b.p.m. and then a decrease for even higher frequencies. The initial $SV$ increase is due to the combination of a larger telediastolic volume and a stronger active contraction; their beneficial effect, however, is countered by the diastasis shortening until the total duration of diastole is not enough to allow for the complete filling of the ventricle thus reducing $SV$. Another factor contributing to the deterioration of $SV$ is the growing blood peak velocities, both through mitral and aortic valve (figure 37d), which increases the transvalvular pressure drops and decreases the pumping efficiency. Nevertheless, as long as $SV$ decreases less than linearly with $HR$, the cardiac output can still increase with the heart rate and this is possible up to $\approx$180 b.p.m. Beyond this frequency, $CO$ also diminishes and the only monotonically increasing quantity is the power consumption.

Figure 37. (a) Blood flow rate versus time for a healthy left ventricle; the positive part of the curve is the inflow through the mitral valve, while the negative side is the outflow through the aortic valve. Thick line for $HR=60$ b.p.m., red solid line for $HR=100$ b.p.m. (b) Duration of representative time intervals (indicated in panel a) for the black curve at $HR=60$ b.p.m.) versus $HR$: $Et$ (blue) duration of E-wave; $At$ (green) duration of A-wave; $Dt$ duration of diastasis; $St$ (red) duration of systole. (c) Stroke volume ($SV$), cardiac output ($CO$) and power consumption ($P$) versus $HR$. (d) Mean blood velocity through valvular ostia for E-wave $VE$ (blue), A-wave $VA$ (green) and systole $VS$ (red). (Some of the data are adapted and replotted from Chung, Karamanoglu & Kovacs 2004).

It is important to point out that although the described phenomena and the related trends are general, the specific values are age- and condition-dependent: in fact, trained athletes reach higher values of $SV$ which saturates only beyond $HR \geqslant 200$ b.p.m. However, the maximum attainable $HR$ decreases with age and the upper limit (in b.p.m.) is given by the empirical formula $HR_{max} = 220 -A$, with $A$ the age expressed in years.

7.4. Cardiomyopathies

A relevant class of disorders affecting the cardiac muscle is referred to as cardiomyopathies and they concern the mechanical structure or the electrical function of the myocardium, both with negative effects of its pumping efficiency.

The most common is dilated cardiomyopathy, in which the left ventricle progressively enlarges and deteriorates its contracting capability: the progression of the disease is caused by myocardium remodelling with an increase of end-systolic and end-diastolic volumes. The ventricle tends to assume a spherical shape (figure 38a) and to reduce its diastolic relaxation as well as the related rapid ventricle filling; the stroke volume declines together with the ejection fraction and the cardiac output, thus increasing the risk of congestive heart failure. The impaired ventricular filling also causes an increase of end-diastolic pressure whose effect is to further increase the ventricle volume. The picture is made even worse by the cordae tendineae pulling action of an expanding ventricle on the mitral leaflets which can not properly seal during systole: this produces significant mitral regurgitation and further reduction of the cardiac output.

Figure 38. Echographic images of the heart left ventricle: (a) systolic phase of a ventricle affected by dilated cardiomyopathy; (b) the same as panel (a) but for a healthy ventricle reported for comparison. (c–f) diastolic and systolic phases for a ventricle affected by hypertrophic cardiomyopathy. LVOT indicates the left-ventricle-outflow-tract and SAM the systolic-anterior-motion. In all panels, the endocardium is evidenced by a green line while in panels (c–f), the mitral leaflets are evidenced by a red line. The thick white arrows of panels (c,e) indicate the main blood direction.

Weintraub, Semsarian & Macdonald (2017) relates the causes of the disease and the above mechanical phenomena either to genetic factors or to acquired causes such as myocarditis and prolonged exposure to alcohol and drugs. As the dilated cardiomyopathy degenerates, it tends to also involve the right ventricle with similar dynamics until the heart can not produce compensatory changes and fails because of congestion; the survival rate at $5$ years from the onset of the disease is $50\,\%$ and the major risk factor is an ejection fraction ${\leqslant }25\,\%$.

The second most common cardiomyopathy is the hypertrophic one that occurs in one out of five hundred people: it is a genetic disease which manifests as an abnormal thickening of the interventricular septum and produces a bulge on the left side (figure 38c), usually without impairing the pumping capability. In some cases, however, the asymmetric septum growth triggers a series of additional phenomena which can obstruct the left ventricular outflow tract (LVOT) and cause sudden cardiac death.

Maron et al. (2014) report that most of the patients affected by hypertrophic cardiomyopathy (HCM) present elongated mitral leaflets which have larger mobility within the ventricle. During systole, blood is squeezed into the aorta accelerating through the LVOT (figure 38e) which is restricted by the thick septum and the long anterior mitral leaflet. This produces abnormal blood accelerations and low pressures which pull the leaflets towards the septum and produce their systolic anterior motion or SAM (figure 38f); the severity of the disorder is measured by the intraventricular subvalvular pressure drop $\Delta p_{spd}$, defined as the maximum pressure difference across the LVOT.

The combination of severe septum thickening, elongated leaflets and strong ejection fractions ($EF > 60\,\%$) produces high values of $\Delta p_{spd}$ which, when exceeding the threshold of $30$ mmHg, causes the outflow tract obstruction that can result in sudden cardiac death.

A common therapy for HCM is to use beta-blockers which reduce the heart rate and the systolic pressure (Gregor & Curila 2015) thus decreasing $\Delta p_{spd}$ and the systolic anterior motion of the mitral leaflets. Another possibility is to induce, by local alcohol injection, controlled ischaemic events that cause necrosis and reduce the septal volume. Finally, surgery can be executed either by removing the bulge of the septum (septal myectomy) or by reducing the excessive length of the mitral leaflets (mitral plication) and both procedures aim at decreasing the $\Delta p_{spd}$ and the induced SAM.

A quite uncommon disease is the restrictive cardiomyopathy characterized by left ventricle stiffness and reduced filling. Unlike the previous two disorders, in this case, there are no morphological changes of the ventricles but rather haemodynamic issues. Also the primary causes of the disease are not identified and, only in some cases, it can be associated with metabolic diseases or radiotherapy effects.

Finally, arrhythmogenic cardiomyopathy occurs in approximately one out of five thousand people, affecting either the left or the right ventricle. It is characterized by electric instability of the myocardium that can degenerate in ventricular fibrillation and induce sudden cardiac death.

7.5. Prosthetic valves

Assuming an average human life expectancy of $\approx$73 years (figure 4) and a physiological rate of $70$ b.p.m., every heart performs slightly less than $3$ billion beats, which is the same number of times heart valves open and close, their leaflets are bent and stretched and their structure is cyclically loaded; similar to the rest of the heart, only a delicate balance between mechanical stresses and biochemical processes triggered by the leaflets’endothelial cells allows such an extraordinary resistance to fatigue.

Indeed, inflammatory, infectious events or age-related damage accumulation can alter this equilibrium and lead to degenerative valve pathologies that can be roughly divided into stenosis and insufficiency: in the first case, the valve does not open completely and offers a reduced cross-section to blood thus causing increased pressure losses and abnormal flow velocities. Insufficiency instead occurs when a valve does not close properly allowing back flow or regurgitation which impairs the pumping capability of the heart. These pathologies, either resulting from slow tissue remodelling or accumulation of calcification processes, usually appear in mild forms with moderate symptoms although they eventually evolve to more severe grades which require pharmacological treatment or surgery.

However, congenital valve diseases are malformations developed before birth affecting size, shape, number of leaflets or the way they are connected to the valvular ring: the resulting disorders are again stenosis and/or insufficiency, but they can now show up in severe forms from birth.

Two main parameters are used to quantify valve efficiency: the effective orifice area ($EOA$) and the regurgitant fraction ($RF$). The first is defined as the ratio of the mean forward flow rate ($\bar {Q}$) through the valve and the mean pressure drop ($\overline {\Delta p}$); unlike the pure geometrical orifice area ($GOA$), $EOA$ accounts also for the velocity and pressure time modulations and it can be estimated by the Gorlin relation $EOA=\bar {Q}/(51.6\sqrt {\overline {\Delta p}})$, with $\bar {Q}$ being expressed in $\textrm {ml}\ \textrm {s}^{-1}$, $\overline {\Delta p}$ in mmHg and $EOA$ in $\textrm {cm}^{2}$.

$RF$ is defined as the backflowing blood volume normalized by the forward counterpart, is expressed in percent and it is used to quantify the degree of insufficiency of a valve.

Physiologic values of both $RF$ and $EOA$ depend on the specific heart valve, the age and the height/weight of the patient, usually expressed by the body surface area ($BSA$ $\textrm {m}^{2}$); clinicians have compiled and refined tables and nomograms which give normal and pathologic ranges to reveal cases that need therapy.

The most common treatments of valvular diseases are surgical remodelling and replacement by the implantation of a prosthesis and in both cases, the aim is to restore $EOA$ and $RF$ as close as possible to their physiological values. Given the very different workload of the left and right heart sides, valve surgery is performed almost exclusively on mitral and aortic valves (Sacks, Marryman & Schmidt 2009; Sotiropulos et al. 2016).

From a technical point of view, designing and building a device with the same function as a heart valve is relatively easy even if biocompatibility of materials, durability requirements, reduced pressure losses and haemodynamic considerations make a prosthetic valve a very advanced piece of engineering which has considerably evolved since the first attempts. In figure 39, some past and present models are reported and an evident difference is that those in the upper row are made of rigid parts (mechanical valves) while the lower row shows valves with flexible leaflets (biological valves).

Figure 39. Past and present prosthetic heart valve models: (a) caged ball; (b) tilting disc; (c) bi-leaflet; (d) biological tri-leaflet; (e) biological stented (or sutureless); ( f) biological for transcatheter implantation. The picture from panel (f) refers to the ‘Medtronic CoreValve Evolut PRO valve’ taken from Medtronic's web site.

Apart from isolated attempts to use valves explanted from corpses, the first prosthetic valves were the ‘caged ball’ model (figure 39a) developed in the 1950s and initially implanted in the descending aorta (heterotropic position) to treat aortic valve insufficiency. Only since the 1960s, when extracorporeal circulation became available, enabling open heart surgery to be performed, have caged ball valves been implanted in their natural annular position. The working principle is based on the force experienced by a body in relative motion with a fluid: in systole, the valve opens as the forward stream drags the ball away from the ring and the cage prevents it from becoming lost in the circulation. Vice versa, in diastole, the backflow pulls the sphere back towards the ring and the valve seals. This device is extremely simple and reliable although the presence of a blunt body in the aortic lumen greatly reduces $EOA$, produces relevant pressure losses and a strongly turbulent wake whose intense stresses alter the equilibrium with the surrounding tissues and the integrity of the blood itself.

A first improvement came from the ‘tilting disc’ model developed in the 1970s (figure 39b) in which the framed disc can rotate by $55^{\circ }$–$70^{\circ }$ when dragged by the fluid, thus increasing the $EOA$ and decreasing the flow perturbation: this relatively simple design change reduced the transvalvular pressure drop from approximately $24$ mmHg, for a caged ball valve of $21$ mm diameter in the aortic position, to $14$ mmHg, for a tilting disc valve in similar conditions.

Bi-leaflet valves (figure 39c) behave even better as the two semicircular leaflets can rotate about the hinges beyond $80^{\circ }$ and their streamlined attitude further reduces flow perturbations; nowadays, this model is the only one used among mechanical valves and it can have a pressure loss as low as $10$–$11$ mmHg if properly implanted.

Regardless of the specific model, for all mechanical devices, the presence of rigid obstacles within the valve ostium reduces their haemodynamic performance with respect to biologic models which, having flexible leaflets, closely reproduce the dynamics of the native systems. This is particularly evident from figure 40, where an instantaneous snapshot of the peak systolic flow through biological and mechanical valves in the aortic position are shown: the inclined leaflets in the middle of the high-speed stream produce thin shear layers which induce high turbulence levels in the downstream region. The characteristic triple-jet flow of the bi-leaflet valve is considerably different from that of the biological valve whose larger and compact ostium generates a smoother coherent flow which maintains its regularity farther downstream. Figure 40(a) also shows the forcing of the shear layer on the recirculation in the Valsalva sinus which, instead, is weaker for a biological valve as the flexible leaflet shields it from the main flow.

Figure 40. Snapshots from numerical simulations of the flow through prosthetic aortic valve models during early systole: (a) mechanical bi-leaflet valve; (b) biological tri-leaflet valve. Colour contours range from $+1\ \textrm {m}\ \textrm {s}^{-1}$ to $-0.3\ \textrm {m}\ \textrm {s}^{-1}$ respectively from magenta to blue. The lower panels report the cross-section of the ostium with the valve in the fully open position. (Courtesy of M.D. de Tullio.)

Considering that the native valves have transvalvular pressure drops (TPDs) of the order of a few mmHg (see figure 9), biological valves, with their $TPD <10$ mmHg and negligible flow perturbations, would be the obvious choice if they had not had durability issues; in fact, after the initial use of porcine valves, the leaflets are now manufactured from bovine pericardium which, after a few years, tends to calcify and to degrade owing to fatigue accumulation. Moreover, the ‘single digit’ pressure drop of a newly implanted biological valve increases in time and, typically after $10$–$15$ years, a further valve replacement is necessary. The latter procedure has been simplified during the last decade since, unlike the model in figure 39(d), stented valves such as those shown in figure 39(e, f) can be crimped into a small section catheter, introduced through a femoral or brachial artery, and expanded inside the old valve without requiring new open heart surgery. This is referred to as TAVI (trans-catheter aortic valve implantation) or a valve-in-valve procedure which is becoming very popular especially among elderly individuals or patients who are non-eligible for open heart surgery.

Obviously, a valve expanded inside an existing one entails an EOA reduction and the procedure can not be repeated. Therefore, if the patient has a life expectation exceeding $20$–$25$ years, a bi-leaflet mechanical valve can be preferable owing to its lifetime durability. Unfortunately, this advantage comes at the price of increased stress levels in the blood whose immediate adverse effect is haemolysis.

Indeed, the thin shear layers and strong velocity fluctuations produced by the mechanical valve both contribute to the generation of abnormal stresses on the red blood cells whose viscoelastic membranes can be permanently impaired and eventually torn with the consequent release of haemoglobin in the plasma. In fact, it has been shown that high levels of mechanical stress act in combination with exposure times to cause a cumulative cell damage which can be quantified by the haemoglobin index $HI$. This quantity is defined as the ratio (in percent) of the concentration variation of the plasma-free haemoglobin (haemoglobin in plasma) after a stress $\tau$ has been applied for a time $t$, $dPHb$ (in mg (100 ml⁻¹)), with the initial total haemoglobin concentration $Hb$ (in mg (100 ml⁻¹)). Experiments, performed with different devices, have resulted in a considerable scatter of the results although there is consensus about expressing the index in the form: $HI = 100\times dPHb/Hb = C t^{\alpha } \tau ^{\beta }$, with $C$, $\alpha$ and $\beta$ coefficients to be determined experimentally. It must be noted that this expression, obtained for a constant, single-component stress, can not be applied directly to unsteady, complex flows like in figure 40 and some amendments must be introduced. The first is to consider a time-dependent stress whose instantaneous value $\tau$ acts for an infinitesimal time $\textrm {d} t$ and produces a damage increment $\textrm {d} HI$, obtained by differentiation $\textrm {d} HI = C \alpha t^{\alpha -1} \tau ^{\beta } \,\textrm {d} t$ which, after some manipulation (Grigioni et al. 2005) and integrated in time, yields:

(7.1)\begin{equation} HI(t) = C \alpha \int_{t_0}^{t} \left ( \int_{t_0}^{\xi} \tau (\xi)^{\beta/\alpha} \,{\rm d} \xi +HI(t_0) \right )^{\alpha-1} \tau(\zeta)^{\beta/\alpha} \,{\rm d}\zeta, \end{equation}

with $HI(t_0)$ the initial damage.

Another necessary change is to reduce the local multi-component stress tensor $\boldsymbol{\tau} =\{\tau _{ij}\}$ of a complex 3-D flow to the scalar quantity $\tau$ needed by (7.1). Simple possibilities are to use the tensor second invariant $\tau = \sqrt {\boldsymbol{\tau} : \boldsymbol{\tau}/2}$ or the von Mises criterion $\tau = \tau _{eq} = \sqrt {3\tau _{ij}\tau _{ij}/2}$ obtained by matching the deformation work of a fluid element in a plain shear flow with the same work for a full 3-D deformation (de Tullio et al. 2009). A more complex alternative is to model the deformation of red blood cells produced by the fluid local strain field and derive the scalar stress from the instantaneous cell configuration (strain-based model): in this context, a cell is approximated by a triaxial ellipsoid and a tensor equation is evolved in time to obtain its semi-axes (Behbahani et al. 2009). The model accounts for the cell membrane tension and the different viscosity of outer and inner fluids; however, the model contains a number of parameters which must be tuned depending on the specific blood properties.

Another key point for the application of (7.1) is that $\tau (t)$ is the stress experienced at each instant by the cell (or by the fluid particle) during its motion and therefore it is a Lagrangian quantity. In fact, high stress levels do not necessarily imply appreciable blood damage as long as they are localized in a region crossed in a short time while, in contrast, moderate stresses can produce cumulative damage if exerted for a long time. This implies that (7.1) must be computed along Lagrangian paths and the evaluation of haemolysis effects for a given device entails statistical averaging over many trajectories to sample all representative cases.

Figure 41 reports for the flow in figure 40(a) three trajectories coloured by the quantity $\tau _{eq}$ and showing very different stress levels depending on the region crossed. Figure 42 yields more quantitative information with the time evolution of the stress and the induced haemolysis index on each trajectory; the most relevant events are clearly the sweeping of solid boundaries which combine high stress levels with reduced convection velocities.

Figure 41. Lagrangian trajectories of three ideal particles ($\textrm {d}\kern0.06em {\boldsymbol {x}}\, \textrm {d} t= {\boldsymbol {u}}$) flowing through the aortic root with a bi-leaflet mechanical valve. The trajectories are coloured by the local and instantaneous value of $\tau _{eq}$. The three panels refer to particles released from different initial positions. (Adapted from de Tullio et al. 2009.)

Figure 42. Time evolution of $\tau _{eq}$ and $HI$ for the trajectories within the mechanical valve of figure 41. Solid line for figure 41(a); dashed line for figure 41(b) and dotted line for figure 41(c). (Adapted from de Tullio et al. 2009.)

Mechanical valves, however, produce most of the adverse effects when they are closed since the counterpressure they support forces high speed leakage through unavoidable gaps forming at the contact between leaflets and between leaflets and housing (figure 43). Another source of leakage is the leaflet hinges where a gap of the order of $100\ \mathrm {\mu } \textrm {m}$ is intentionally created by design to generate effluent chimneys washing, with every heartbeat, the hinges and preventing stagnant blood from producing clots.

Figure 43. (a) Flow leakage through a closed bi-leaflet aortic valve during diastole, obtained by numerical simulation. Contours are velocity magnitude and colours range from $+1\ \textrm {m}\ \textrm {s}^{-1}$ to $-0.3\ \textrm {m}\ \textrm {s}^{-1}$ respectively from magenta to cyan. (The related pressure field is that reported in figure 33a.) (b) Schematic of the valve and of the hinge mechanism, where the cavity in the housing is reported in red while the spherical ear of the leaflets is indicated by a black circle. The green lines and arrows show the possible leakage of the leaflet when it is closed. The ear of the leaflet rolls without sliding into the butterfly-like cavity of the housing to avoid friction wear.

Figure 44 reports a comparison of haemolysis indices computed for the valves of figure 40 clearly showing that blood damage generated by the biological prosthesis is one order of magnitude smaller than that caused by the mechanical one. Furthermore, the time evolution of $HI$ evidences a first increase at peak systole caused by the intense velocity fluctuations and an even larger increment at early diastole during valve closure: in both cases, the mechanical valve shows worse haemodynamic behaviour.

Figure 44. Time evolution of the haemolysis index during a heartbeat for a mechanical and a biological valve in aortic position as in figure 40. The dashed curve is the aortic flow rate waveform to evidence the phases of $HI$ growth. Here, $HI$ has been computed by releasing $4\times 10^{5}$ particles per cycle over four heartbeats and phase averaging the data. (Courtesy of M.D. de Tullio.)

The result of figure 44 finds clinical confirmation in the decreased red blood cell half-life time in patients with mechanical heart valves: in fact, in physiological conditions, the RBC half-life is approximately $120$ days, after which RBCs are eliminated by the the spleen and liver and their iron recycled, while the presence of mechanical valves reduces this time to $99$ days (Mitlyng et al. 2006). A similar effect is found for platelets: Hugh & Gebton (1970) report a half-life time in physiological conditions of approximately $7$ days which reduces to $\approx$5.5 days when a prosthetic mitral valve is implanted. One of the most serious adverse effect of mechanical haemolysis is that abnormal stresses not only impair RBCs but also produce platelet activation which tend to capture the circulating damaged cells before they can be sequestrated by the spleen and liver. These conditions can lead to clot formation in the arterial circulation, increasing the risk of ischaemiae which can be lethal in case of heart infarctions and strokes. To prevent these life-threatening events, lifelong anticoagulation therapy is necessary, and this induces high level of haemorrhage risk.

Prosthetic valves are used also for mitral valve replacement although all the issues encountered for aortic valves are exacerbated by the non-circular ostium, the asymmetric native valve anatomy, the leaflets connection by cordae tendineae to papillary muscles, and its effect on the ventricular recirculation. Furthermore, the main access to the site is still the left ventricle apex which adds myocardial injury: all these factors concur to make mitral valve replacement less successful than the aortic procedure. In fact, the poorer outcome of mitral valve replacement with respect to the aortic one is confirmed by Goldstone et al. (2017), who performed an assessment of the fifteen-year survival after surgery for a large cohort of patients. This study separated the causes of mortality showing that, even after considering bleeding and reoperation, strokes were the most likely deadly events thus pointing at altered ventricular haemodynamics with increased probability of clot formation.

In contrast to aortic valve replacement, the adverse results do not depend on the valve model; in fact, Cetinkaya et al. (2019) evaluated the long-term outcome of mitral valve replacement using biological and mechanical valves, and found that, although mechanical valves are associated with an increased risk of bleeding caused by the anticoagulant therapy, there are no significant survival differences for the two valves after a follow up of ten years.

Figure 45. Numerical simulation of the flow in a model left ventricle during early diastole (E-wave) for an ejection fraction $EF=40\,\%$: (a) flow with native mitral valve; (b) flow with mechanical bi-leaflet valve in mitral position, panels (a,b) adapted from Meschini et al. 2018); (c) the same as panel (b) but experimental PIV measurements. (Courtesy of G. Querzoli, adapted from Querzoli, Fortini & Cenedese 2010.)

Indeed, Meschini et al. (2018) investigated the left ventricle flow after a mitral valve replacement and they found that none of the prosthetic valves preserved the beneficial variable direction mitral jet, caused by the native leaflets asymmetry, and all artificial valves give smaller EPI values thus yielding mitral jets with a reduced penetration towards the ventricle apex (figure 45). This is particularly evident for impaired ventricles whose ejection fraction is smaller than physiological values ($50$–$70\,\%$): in these cases, the pressure losses produced by the prosthetic valve further reduce the mitral jet momentum which, in turn, drives a weaker ventricular vortex. The extreme case is given by a mitral mechanical valve in an impaired ventricle ($EF=40\,\%$) whose weak jet is further split into three thinner jets by the open leaflets and its momentum is transferred to small-scale fluctuations; the resulting EPI is only $0.69$, much smaller than the unity threshold needed to guarantee a proper ventricle washing and a stagnant blood region at the ventricle apex is formed. Similar results were found by laboratory measurements by Querzoli et al. (2010) and later confirmed by independent experiments and numerical simulations by Vu et al. (2019).

7.6. Ventricular assist devices

The shortage of organs available for transplantation and the issues related to foreign-body rejection make the realization of a full prosthetic heart a key goal for bioengineering. However, this goal is still far from being reached. In principle, replacing the mechanical function of the heart is simple since it is just a volumetric pump with limited flow rate, prevalence and power. Nevertheless, handling blood in a living body and reproducing the self-adaptive heart capabilities entail many subtleties that, to date, have prevented the existing devices from becoming a functional artificial heart.

In fact, the current models are used only as a bridge therapy to transplantation as they need bulky external driving units which must be operated in hospitals or heavy battery packs of limited duration: additional limiting factors are the weight and size of the implanted unit, the cannulae tunnelled through the chest to convey power which are often sources of infections and haemolysis effects on the blood. Presently, the survival of implanted patients ranges from a few days to a few hundred days with a strongly reduced life quality.

Moreover, the design of ventricular assist devices (VADs) has been more successful. These devices can be implanted to patients with some remaining heart function. VADs use has gradually changed from being used exclusively as a bridge to transplantation to a destination therapy for elderly patients or a temporary support to ease the healing after major cardiac surgery or myocardial infarction. Several configurations are available to assist the right/left ventricle (RVAD/LVAD) or both (BiVAD), with LVAD representing the majority of cases.

The working principle is to reduce the load of a ventricle by a transapical implantation of a pump which shunts the flow directly to a downstream artery and restores the cardiac output to physiological values; in figure 46, a drawing is shown for a LVAD, implanted in an impaired dilated left ventricle, which drives the flow into the ascending aorta. The residual ventricular ejection fraction ($EF_r$, typically ${<}30\,\%$) multiplied by the heart rate yields the remaining cardiac output $CO_r$ which is complemented by the LVAD with a flow rate $CO_{VAD}$ that can be modulated by changing the bypass ratio. The total cardiac output $CO = CO_r+CO_{VAD}$ can be adjusted to the target value and even modified to stimulate the natural heart function or to provide additional support.

Figure 46. (a) Drawing of the implanted left ventricular assist device (LVAD) with the shunt in the ascending aorta; (b) section of a dilated left ventricle with an apical LVAD implantation.

A major problem is the use of continuous $HR=60$ bpmaxial pumps which are preferred because of their compact and lightweight design and haemodynamic properties, although they provide a steady $CO_{VAD}$: even if the total $CO$ attains physiological values, the pulsatile character of the flow is lost and the resulting haemodynamics changes, in the medium term, also affect other organs. The ventricular flow is clearly the most affected as the physiological recirculation is weakened and additional stagnant flow regions are produced. Figure 47 reports PIV experimental measurements of the flow at peak diastole for two ventricles without and with LVAD indeed confirming that, despite the higher $CO$ of the implanted case, the presence of VADs not only reduces the circulation intensity but also deteriorates its spatial coherence (Viola et al. 2019). As mentioned before, these dynamical conditions increase the risk of thrombus formation which is countered by a anticoagulant therapy in turn increasing the haemorrhagic risk: indeed, bleeding is the major short-term issue for VAD implantation and it accounts for $60\,\%$ of adverse events (Schaffer et al. 2011).

Figure 47. Instantaneous snapshots of the flow at peak diastole in a left ventricle with $CO=2.4\ \textrm {l}\ \min ^{-1}$ (a) and for a LVAD implanted ventricle with $CO>7\ \textrm {l}\ \min ^{-1}$ (b). Velocity vectors obtained by PIV measurements and colour contours of out-of-plane vorticity, red for clockwise and blue for anticlockwise: (c) for the configuration in panel (a) and (d) for the configuration in panel (b). The length $d_v$ indicates the size of the vortex from the compact vortex patch. (Adapted from Viola et al. 2019.)

8. Outlook and closing remarks

In this article, we have tried to give a comprehensive overview of cardiovascular mechanics although, given the complexity of the subject and the involvement of different disciplines, many topics have been just mentioned and others overlooked. In particular, we have focussed on those aspects of the electro-fluid-mechanics of the heart required to understand its complex haemodynamics.

Furthermore, given the pedagogical scope of JFM Perspectives papers, we have presented a systematic discussion of consolidated cardiovascular knowledge rather then an illustration of the latest results. Nevertheless, even within this context, there are emerging fields which are evolving very rapidly and it is very likely that they will be the future trends of cardiovascular research; in the following, a short overview of the most promising developments will be given together with a discussion of possible future applications.

8.1. Uncertainty quantification, virtual trials and surrogate models

Typical digital heart models consist of interdependent systems of partial differential equations enforcing conservation and balance laws constrained by additional semi-empirical relations. The equations are typically nonlinear and evolve in a complex domain whose shape and topology is part of the sought results. For these models, analytical solutions are out of the question and numerical integration is required. These models are then completed by physical and numerical parameters which usually can not be obtained by direct measurements and therefore must be inferred from medical imaging or indirect evidence. Even the definition of a cardiovascular geometry entails approximations since standard models have to cope with strong biological variability while a patient-specific configuration, which is constructed for a single individual and is one-of-a-kind, is unavoidably affected by measurement errors and artefacts of reconstruction algorithms.

Concurrently, the continuous growth of available computational power has increased the complexity and sophistication of cardiovascular models to such a level that their predicting capabilities are now limited more by the uncertainty of input parameters than by the fidelity of the biophysical description. Quantifying the propagation of the various uncertainties through the model and assessing their effect on the quantities of interest, therefore, has become an urgent need to determine the confidence on the results and to decide which elements need refinement.

Uncertainty quantification (UQ) is a cross-disciplinary approach applied in those fields where a single deterministic realization is not representative of the system dynamics: weather forecasting is an example in which tens of concurrent simulations are run using randomly perturbed initial conditions to account for their incomplete knowledge and measurement errors. These simulation campaigns produce a set of weather maps from which the likelihood of some events can be extracted via statistical processing.

Given the human variability and the epistemic uncertainty of most of the model parameters, a similar approach is key in cardiovascular flows to assess the reliability of predictions that are used to guide choices of clinicians and surgeons. Indeed there are already some attempts to apply UQ to cardiac dynamics although they are still limited to some specific aspects of the system such as the electrophysiologic system (Del Corso, Verzicco & Viola 2020), cardiac myocytes (Mirams, Niederer & Clayton 2020) or the left ventricle (Cai et al. 2021). The high computational cost of full heart simulations has prevented its analysis from an UQ perspective; the trend, however, is in that direction given the continuously growing computational power available. This is confirmed by Strocchi et al. (2020), who have compiled a cohort of $24$ geometries of four-chambers hearts from patients affected by heart failure and have made the database publicly available to stimulate UQ electromechanics studies.

Needless to say, the determination of heart geometries is only the starting point of a virtual cohort, and needs to be completed by elastomechanic properties of the tissues, electrophysiologic parameters, blood properties and impedances of pulmonary and systemic circulations to be used as boundary conditions: the UQ analysis of all these quantities entails a huge computational burden which needs careful planning of the cases to be simulated.

This undertaking, however, would definitely be worthwhile and the results extremely rewarding. In fact, running cardiovascular simulations on a digital cohort of hearts is nothing but a virtual clinical trial which, however, does not entail all the ethical and economical implications of the counterpart performed on real patients. Furthermore, in a virtual model, all relevant parameters can be isolated and varied in clean and controllable conditions, and variables can be modified individually or in small groups. In contrast, it is very unlikely to enrol real patients who fit exactly the needs. This is even more true for rare pathologies or for the extreme manifestations (the distribution tails) of more common diseases.

The high computational cost of complex multiphysics heart models and the need of trained researchers to run them constitutes an insurmountable barrier to the use of digital twins in clinical practice. Furthermore, they rely on massive parallel supercomputers, have a time-to-solution of the order of days and produce terabytes of data which need extensive post-processing to extract relevant information; this is totally incompatible with clinical decisions that must be taken, at most, within hours and obtained by commodity computers.

For these reasons, efforts are being made to produce simpler mathematical models capable of replicating some features of the high-fidelity counterparts without requiring exorbitant computational resources, lengthy elaboration times and skilled technicians. These are referred to as surrogate models and they encompass a variety of techniques all aimed at drastically decreasing the number of degrees of freedom and the mathematical complexity.

The basic idea is to consider the high-fidelity model as a black box transfer function ${\mathcal {M}}$ which maps the uncertain vector of input parameters ${\boldsymbol {x}}$ to the uncertain vector of the quantities of interest ${\boldsymbol {y}}={\mathcal {M}}{\boldsymbol {x}}$. Owing to input uncertainties, characterizing this transfer function entails a large number of simulations which is usually unaffordable; therefore, ${\mathcal {M}}$ is surrogated by another function $\hat {{\mathcal {M}}}$ which is constrained by some properties of ${\mathcal {M}}$ and calibrated by a representative set of high-fidelity model results. A common choice is the polynomial chaos expansion (PCE) in which $\hat {{\mathcal {M}}}$ is expressed as a series of orthogonal Legendre or Hermite polynomials depending on the statistical character of the input parameters. The coefficients of the series are then computed through the calibration step which requires the predictions of $\hat {{\mathcal {M}}}$ to ‘match’ (or to minimize the error) with respect to ${\mathcal {M}}$ and in practical computations, the number of coefficients is limited by truncating the surrogate function to a maximum order $p$.

Additional details and alternative techniques can be found in Ghanem, Higdon & Owhadi (2017). Here, we wish to point out that the crucial point of this technique is the calibration of $\hat {{\mathcal {M}}}$ which is nothing but a fitting which therefore depends on the reference dataset coming from clinical observations or from the high-fidelity model.

In retrospect, these surrogate models can be considered as the modern version of nomograms, multiple-entry tables on which medical doctors have relied for more than a century and that are now being replaced by digital analogues calibrated from high-quality databases.

8.2. Data assimilation and data-driven simulations

Cardiovascular models, experiments and clinical observations are complementary in nature and should go hand in hand, possibly exploiting their synergistic combination. In fact, mathematical models are based on approximate physical descriptions and are fed by uncertain parameters, experiments typically have low resolutions and noisy results, while clinical data are sparse and unevenly distributed. In this context, hybrid approaches aimed at combining information from different sources would be extremely helpful and recent applications have moved towards this direction.

Indeed, in other disciplines, these are already consolidated techniques and, once again, we mention weather forecasts in which the set of hundreds of millions to billion variables of the computational model has to be initialized by at most tens of thousands of measurements and observations collected from meteorological stations, ground-based remote sensing instruments or aeroplane and ship detections. Early attempts made using simple interpolations yielded awkward initial fields producing unphysical acoustic and gravity waves. A major improvement came from a preliminary model integration with the solution constrained at those points where measurements were available; this was the beginning of data assimilation which is nowadays performed continuously in time to generate updated initial conditions and to prevent drifts of the numerical solution (Kalnay 2002).

In simpler words, data assimilation consists of combining recent observations with contemporary model predictions to obtain a consistent actual picture of the system which can be used to initiate new simulations for the time ahead. Data assimilation can also be thought of as a way of keeping a model ‘on track’ by constantly constraining it with updated observations.

In cardiovascular dynamics, the assimilable information might come from medical imaging or other diagnostics usually with limited spatio-temporal resolution; a computational model can then be run with the solution constrained by the available data. There is a variety of possible constraints and they can be the enforcement of boundary kinematics on a portion or all the domain but also volume forcings on the flow to penalize deviations from the available velocity measurements: Arzani & Dawson (2021) describes a large number of techniques, like principal component analysis, Kalman filtering or dynamic mode decomposition, together with several examples for their application.

In figure 48, a conceptual data assimilation, from two-dimensional echo-Doppler measurements of the aortic root and upper left ventricle to a three-dimensional model of the left heart, is provided showing that by this approach, it is possible to ‘augment’ the diagnostic capability of simple tools. In fact, echographic imaging is simple and innocuous but the two-dimensional information must be complemented with the knowledge and experience of the medical doctor; however, the same data could be assimilated by a three-dimensional computational model and the numerical solution of the physics-based differential equations used to complete the picture, therefore providing a sort of augmented reality driven by traditional tools.

Figure 48. Conceptual view of data assimilation from two-dimensional echo-Doppler measurements to constrain a three-dimensional computational model.

8.3. Mechanotransduction models for long-term dynamics

The myocardium innermost layer, the endocardium, as well as that of arteries and veins, the endothelium, are constantly exposed to haemodynamic loads consisting of pressure and shear stresses. The latter are converted by the endothelial cells in biochemical signals which regulate tissue permeability, dilatation or restriction, adhesiveness to circulating leukocytes and platelets, in addition to changes of the cells’ structure and tissue remodelling: the associated time scales range from seconds to months and some processes are the result of a cumulative effect over time.

Studies at the level of individual cells have identified specialized structures, such as glycocalyxes or membrane pores, which transform stresses into biochemical processes through protein alterations and associations (Tanishita & Yamamoto 2016). Owing to the intractable high number of cells, however, when a tissue or an organ is considered, the above mechanisms have to be modelled at larger spatial scales; these are continuum mechanotransduction models that relate the sensed mechanical stimulus to the final effects on the tissue.

A simple example could be an incipient valve stenosis or minor plaque accumulation causing a local blood acceleration and, say, a $10\,\%$ increase of the local wall shear stress: this variation is unlikely to produce appreciable immediate effects on cardiac function, which can temporarily adapt to much larger pressure and flow-rate variations. Nevertheless, the long-term reaction of endothelial cells to the altered stress can induce a progressive tissue thickening, and even alter fibre alignment, thus triggering a negative feedback loop.

These processes, however, develop on time scales of the order of months or even years, while accurate computational and experimental models are limited to tens of heartbeats and therefore to evolution times of the order of minutes. This huge disparity of time scales can not be bridged by brute force computation and alternative strategies have to be found. A possibility to overcome these limitations is to couple high-fidelity simulations and experiments with continuum mechanotransduction models so that the former provide a detailed spatio-temporal distribution of the wall shear stresses over all the endothelial tissue and the latter determine the long-term evolution of the altered condition. Of course, while tissue remodelling occurs, the flow and the wall shear stresses change; therefore, the two models should be iteratively executed up to the attainment of an equilibrium configuration or to the impairment of the system physiology.

Continuum mechanotransduction models can be very useful for anticipating the long-term outcome of a surgical procedure, the evolution of abnormal haemodynamics or the progression of a disease diagnosed at the initial stage (Li et al. 2009); these innovative combinations of clinical observations and computational science are the basis of precision medicine which takes into account individual variability to predict more accurately treatments and prevention strategies.

8.4. Closure

In this JFM Perspectives paper, the dynamics of the human heart has been presented with a particular focus on its haemodynamic aspects. Still, the highly synergistic nature of the heart pumping function has required also the description of the elastodynamics of the tissues and of the electrophysiology. These systems are strongly interconnected and all cooperate to ensure organ function: the final result is a reliable and efficient volumetric pump that, using only $8$ W of power, ensures a flow-rate of $5\ \textrm {l}\ \min ^{-1}$ along an intricate circuit of cumulative length $10^{5}$ km, lifelong. Indeed, the mechanical efficiency of the pump is not extraordinary since only $2$ W are used as pumping power, with the remaining $6$ W producing precious heat which constantly warms up the blood distributed throughout the body via the circulations.

Biologically, the human heart is similar to that of all warm-blooded mammals and complies to the same allometric laws; nevertheless, improved living conditions, the ability to modify the environment and recent advances in medical research have almost tripled its lifetime in terms of heartbeats compared to the ‘natural’ situation.

Reproducing the heart dynamics, both by computational models or laboratory experiments, is a formidable challenge since each system is highly sophisticated and two-way coupled to the others: the electrical system, by the propagation of an activation potential and delaying barriers, determines the coordinated active contraction and relaxation of the myocardium which is an orthotropic, inhomogeneous and hyper- (visco)-elastic tissue. The time-dependent volume change of the heart chambers produces a flow from the atria to the ventricles and from there to the circulations, the backflow being prevented by the passive opening and closing of four heart valves. Finally, the blood is also functional to the system since it is a dense suspension with shear-thinning behaviour thus reducing its effective viscosity in high shear-rate regions and, therefore, limiting the power needed for a given flow rate.

Laboratory experiments usually focus on specific details of the heart, such as one chamber possibly with one or two valves or the coupling between one chamber and one of the main vessels: the technical complexity of pulse duplicators and optical accessibility of the system prevent more complete models from being assembled. Blood analogues reproduce well the physiological rheology while heart phantoms made of homogeneous, isotropic and elastic materials (usually transparent silicone rubber) match only the mean deformation and fail to reproduce the other details.

Numerical models, in contrast, do not suffer from these limitations and, in principle, can account for all physiological and pathological details of tissues, blood and electrophysiology to realize a ‘digital twin’ of the heart. This undertaking, however, entails huge computational complexity and costs that are only becoming possible nowadays thanks to the latest massive parallel computers and GPU accelerated architectures.

To support the medical community with useful results, models have to account for human variability and uncertainty of input parameters: this implies that single deterministic runs are of little use and extensive simulation campaigns have to be performed in the form of virtual clinical trials, thus further exacerbating the computational burden.

However, digital twins with their simulation campaigns are too complex and computationally demanding to be used in clinical practice which, therefore, has to be provided with simpler, more manageable tools executed on commodity computers: these include reduced and surrogate models that can be considered nothing but a modern version of the traditional nomograms which have supported clinical decision making for more than a century.