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 (Reference Dawson2001), 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 (Reference Dawson2001) 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 Reference Schmidt-Nielsen2004).

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. Reference Fu, Shibata, Hastings, Platts, Hamilton, Bungo, Stenger, Ribeiro, Adams-Huet and Levine2019). Concerning the friction losses, these certainly depend on the extension of the capillary network which Miller (Reference Miller2018) 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 Reference Dawson2001). (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 Reference McMahon and Bonner1983).

In fact, Miller (Reference Miller2018) 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 Reference Dawson2001), 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 (Reference Bassil, Zarzoso and Noujaim2017) 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 Reference Schmidt-Nielsen2004).

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. (Reference Seymour, Hu, Snelling and White2019) 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 (Reference Zhang and Zhang2009) 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 Reference Zhang and Zhang2009). 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.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 (Reference Katz2011) 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 Reference Piper and Preuse1993). 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!

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. Reference Amabili, Balasubramanian, Ferrari, Franchini, Giovanniello and Tubaldi2020). 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. (Reference Brindise, Busse and Vlachos2018), 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 Reference Brindise, Busse and Vlachos2018).

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. Reference Campo-Deaño, Dullens, Aarts, Pinho and Oliveira2013). Brookshier & Tarbell (Reference Brookshier and Tarbell1993) 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. (Reference Carneiro, Lima, Campos and Miranda2021) 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. Reference Yazdi, Geoghegan, Docherty, Jermy and Khanafer2018): (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 (Reference Azhari, Beyar and Sideman1999) 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 Reference Withey and Koles2008). 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 Reference Arcaute and Wicker2008). 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 Reference Moxham2003). 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. (Reference Querzoli, Fortini, Espa and Melchionna2016) 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. Reference Querzoli, Fortini, Espa and Melchionna2016).

Although not devoted to visualizations or measurements, a peculiar set-up is that described by Leopaldi et al. (Reference Leopaldi, Wrobel, Speziali, van Tuijl, Drasutiene and Chitwood2018) 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. (Reference Fortini, Querzoli, Espa and Cenedese2013) 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. (Reference Fortini, Querzoli, Espa and Cenedese2013), 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. Reference Elsinga, Scarano, Wieneke and van Oudheusden2006). 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. (Reference Poelma, Mari, Foin, Tang, Krams, Caro, Weinberg and Westerweel2011) were able to measure the three-dimensional flow in a curved pipe while Versluis et al. (Reference Versluis, Stride, Lajoinie, Dollet and Segers2020), 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. (Reference Montalba, Urbina, Sotelo, Andia, Tejos, Irarrazaval, Hurtado, Valverde and Uribe2018) 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 (Reference Del Corso, Verzicco and Viola2022).

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. Reference Clayton, Bernus, Cherry, Dierckx, Fenton, Mirabella, Panfilov, Sachse, Seemann and Zhang2011):

(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. (Reference Loppini, Gizzi, Ruiz-Baier, Cherubini, Fenton and Filippi2018), 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. Reference Loppini, Gizzi, Ruiz-Baier, Cherubini, Fenton and Filippi2018).

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 (Reference ten Tusscher and Panfilow2006) for ventricular myocardium, Courtemanche, Ramirez & Nattel (Reference Courtemanche, Ramirez and Nattel1998) for the atria and the bundles, and Stewart et al. (Reference Stewart, Aslanidi, Noble, Noble, Boyett and Zhang2009) 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 (Reference Rush and Larsen1978) 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 (Reference Nash and Panfilov2004) or Niederer & Smith (Reference Niederer and Smith2008) (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 Reference Gurtin1981). 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 Reference Fung1993).

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. (Reference Ambrosi, Arioli, Nobile and Quarteroni2011) 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 Reference Fedosov2010), 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 Reference de Tullio and Pascazio2016; Viola, Meschini & Verzicco Reference Viola, Meschini and Verzicco2020). 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. Reference Nealen, Muller, Keiser, Boxerman and Carlson2006) 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 Reference van Gelder1998).

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 Reference De Vita, de Tullio and Verzicco2016). 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. Reference Campo-Deaño, Dullens, Aarts, Pinho and Oliveira2013).

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. Reference Bavo, Roccatello, Iannacone, Degroote, Vierendeels and Segers2016).

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 (Reference Peskin1972). 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 Reference Verzicco and Querzoli2021).

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 (Reference Griffith and Patankar2020) entirely devoted to the application of IB methods to biological flows.

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. Reference Vignon-Clementel, Figueroa, Jansen and Taylor2006; Esmaily Moghadam et al. Reference Esmaily Moghadam, Bazilevs, Hsia, Vignon-Clementel and Marsden2011).

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 Reference Marsden2014).

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. Reference Viola, Spandan, Meschini, Romero, Fatica, de Tullio and Verzicco2022). 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. Reference de Tullio, Afferrante, Demelio, Pascazio and Verzicco2011a). (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. Reference Meschini, de Tullio, Querzoli and Verzicco2018). (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 Reference Viola, Meschini and Verzicco2021).

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. (Reference de Tullio, Cristallo, Balaras and Verzicco2009) 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. Reference Celotto, Zovatto, Collia and Pedrizzetti2019) or even by a prescribed kinematics (Ge & Sotiropulos Reference Ge and Sotiropulos2007). 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 (Reference Sotiropulos, Le and Gilmanov2016) and Mittal et al. (Reference Mittal, Seo, Vedula, Choi, Liu, Huang, Jain, Younes, Abraham and George2016).

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. Reference Le, Elbaz, van der Geest and Sotiropulos2019; Collia et al. Reference Collia, Zovatto, Tonti and Pedrizzetti2021). 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 Reference Westerhof, Lankhaar and Westerhof2008): 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. Reference Pfaller, Cruz Varona, Lang, Bertoglio and Wall2020). 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 (Reference Moulton, Hong and Secomb2017) 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.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. (Reference Viola, Spandan, Meschini, Romero, Fatica, de Tullio and Verzicco2022).

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. (Reference Harfi, Seo, Yasir, Welsh, Meyer, Abraham, Gerge and Mittal2017) 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 Reference Mitter, Shah and Thomas2017).

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 Reference Meschini, Mittal and Verzicco2021.)

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. Reference Kräuter, Reiter, Reiter, Nizhnikava, Masana, Schmidt, Fuchsjäger, Stollberger and Reiter2020).

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. Reference Sotiropulos, Le and Gilmanov2016). 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. Reference Azarine, Garcon, Stansal, Canepa, Angelopoulos, Silvera, Sidi, Marteau and Zins2019), 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. Reference Russo, Banuth, Nader and Dreyfuss2020), 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 Reference Gaine, Naeije and Peacock2014). 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. Reference Gaine, Naeije and Peacock2014).

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. Reference Meschini, de Tullio, Querzoli and Verzicco2018). 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 (Reference Nitti, De Cillis and de Tullio2022) 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. (Reference Chung, Karamanoglu and Kovacs2004) 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 Reference Chung, Karamanoglu and Kovacs2004).

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 (Reference Weintraub, Semsarian and Macdonald2017) 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. (Reference Maron, Ommen, Semsarian, Spirito, Olivotto and Maron2014) 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 Reference Gregor and Curila2015) 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 Reference Sacks, Marryman and Schmidt2009; Sotiropulos et al. Reference Sotiropulos, Le and Gilmanov2016).

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. Reference Grigioni, Morbiducci, D'Avenio, Di Benedetto and Del Gaudio2005) 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. Reference de Tullio, Cristallo, Balaras and Verzicco2009). 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. Reference Behbahani, Behr, Hormes, Steinseifer, Arora, Coronado and Pasquali2009). 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. Reference de Tullio, Cristallo, Balaras and Verzicco2009.)

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. Reference de Tullio, Cristallo, Balaras and Verzicco2009.)

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. Reference Mitlyng, Chandrashekhar, Furne and Levitt2006). A similar effect is found for platelets: Hugh & Gebton (Reference Hugh and Gebton1970) 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. (Reference Goldstone, Chiu, Baiocchi, Lingala, Patrick, Fischbein and Woo2017), 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. (Reference Cetinkaya, Poggenpohl, Bramlage, Hein, Doss, Bramlage, Schonburg and Richter2019) 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. Reference Meschini, de Tullio, Querzoli and Verzicco2018); (c) the same as panel (b) but experimental PIV measurements. (Courtesy of G. Querzoli, adapted from Querzoli, Fortini & Cenedese Reference Querzoli, Fortini and Cenedese2010.)

Indeed, Meschini et al. (Reference Meschini, de Tullio, Querzoli and Verzicco2018) 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. (Reference Querzoli, Fortini and Cenedese2010) and later confirmed by independent experiments and numerical simulations by Vu et al. (Reference Vu, Rossini, Montes, Campos, Moon, Martinez-Lezcapi, Bermejo, Del Alamo and May-Newman2019).

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. Reference Viola, Jermin, Warnock, Querzoli and Verzicco2019). 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. Reference Schaffer, Arnaoutakis, Allen, Weiss, Patel, Russel, Shah and Conte2011).

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. Reference Viola, Jermin, Warnock, Querzoli and Verzicco2019.)

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 Reference Del Corso, Verzicco and Viola2020), cardiac myocytes (Mirams, Niederer & Clayton Reference Mirams, Niederer and Clayton2020) or the left ventricle (Cai et al. Reference Cai, Ren, Wang, Xie, Zhu and Gao2021). 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. (Reference Strocchi2020), 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 (Reference Ghanem, Higdon and Owhadi2017). 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 Reference Kalnay2002).

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 (Reference Arzani and Dawson2021) 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 Reference Tanishita and Yamamoto2016). 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. Reference Li, Gao, Mangion, Berry and Luo2009); 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.