1. Introduction
Most of us worry about the heart only when something goes wrong and fortunately that happens very seldom. If we only reflected on some figures, however, heart reliability would look truly exceptional: with a weight of $0.3$ kg (less than half a percent of the total body mass) and a power consumption of approximately $8$ W, it pumps blood at a rate of $5\ {\rm l}\ \min ^{-1}$ through the circulatory system, whose cumulative length exceeds $10^{5}$ km, lifelong. The heart performs this task by contracting $40$ million times per year at a rate of $70$ beats per minute (b.p.m.) although, depending on environmental conditions or intense physical exercise, it self-adapts to yield up to a six-fold flow rate, with respect to resting conditions, without us being conscious of it.
Indeed, humans have always been aware of the heart presence owing to its continuous rhythmic beat and perhaps they have also guessed its importance for life looking at the immobile hearts of hunted prey. The understanding of the exact role of this organ in body physiology, however, has taken a long time and independent rediscoveries.
Apparently, already in $2600$ BC, the Chinese Yellow Emperor, Nei Ching, in the Canon of Medicine wrote: ‘The blood current flows continuously in a circle without a beginning or end and never stops’ and that ‘all the blood is under control of the heart’ (Mackay & Mensah Reference Mackay and Mensah2004).
Unfortunately, this brilliant intuition, incredibly ahead of its time, remained confined within China for two millennia, before Aristotle (in the fourth century BC) conceived the heart as ‘the source of all movement, since the heart links the soul with the organs of life’ which can be considered correct only philosophically.
Galen (in the second century AD) advanced the knowledge by discovering the pulmonary circulation although he considered the liver as the main circulation organ and identified the heart as the source of the ‘innate heat’ which was distributed throughout the body by blood (May Reference May1968); this misconception was the rational for bleeding as a treatment for fever.
As centuries progressed, evidence accumulated against Aristotle's view although his writings dominated all fields of knowledge and the ipse dixit statement (he said it himself) prevented any further discussion.
The age of Renaissance revived the empirical observation and, through dissection of corpses and animals, a good understanding of the heart structure was achieved; the drawings by Leonardo da Vinci clearly show how detailed the knowledge of chambers, valves and connected veins or arteries was (see figure 1) although the notion of the heart physiology still did not differ much from Galen's ideas.
Finally, in the 17th century, William Harvey rediscovered the function of the heart when he neatly stated that ‘It has been shown by reason and experiment that blood by the beats of the ventricles flows through the lungs and heart and is pumped to the whole body … the blood in the animal body moves around in a circle continuously and … the action or function of the heart is to accomplish this pumping. This is the only reason for the motion and beat of the heart’ and also ‘The heart's one role is the transmission of the blood and its propulsion, by means of the arteries, to the extremities everywhere’ (O'Malley, Poynter & Russell Reference O'Malley, Poynter and Russell1961).
The revolutionary observations of William Harvey paved the way for a systematic study of blood circulation and the findings of capillaries and red blood cells, by Marcello Malpighi (in 1661), provided the missing link between arterial and venous blood while Richard Lower (in 1669) and Antoine Lavoisier (in 1777) explained the ${\rm CO}_2$–${\rm O}_2$ exchange in the lungs, thus completing the picture. In the same period, Galvani (Reference Galvani1791), using dissected frogs, proved that electricity could trigger muscle contraction. The 18th century produced a number of complementary independent discoveries, although it was only in the middle of the 19th century that Purkinje (Cavero, Guillon & Holzgrefe Reference Cavero, Guillon and Holzgrefe2017) reported the existence of a web of fibres (the nerves) within the myocardium whose role in the transmission of electrical signals was understood in the early 20th century (Tawara Reference Tawara1906). However, His (Reference His1949) proved that the heart chambers were electrically connected and that the heart could beat independently of the central nervous system.
In the 20th century, technology appeared in cardiology: W. Einthoven, in 1901, built the first electrocardiograph which weighed $270$ kg, occupied two rooms and required five people to operate it: he won the Noble Prize in 1924. He was a forerunner of telemedicine as in 1905, using a telephone cable, he transmitted an electrocardiogram from a hospital to his office at a distance of $1.5$ km (Mackay & Mensah Reference Mackay and Mensah2004).
In 1931, A. Hyman developed the first cardiac pacemaker which imparted electrical impulses to the heart by transthoracic needles. The first heart defibrillation was performed in 1947 by C. Beck who, following the advice of the physiologist C.J. Wiggers, restored the normal heart rhythm by electric shocks generated by an external device.
The second half of 20th century saw the combination of technology and cardiac surgery, with the first aortic valve replacement with a mechanical prosthesis and a totally implantable pacemaker in 1960, an artificial heart in 1969 and a permanent device in 1982. Miniaturized axial pumps, employed as ventricular assist devices, were developed in 2000 and 2002.
Concurrently, diagnostic tools for cardiac imaging progressed quickly with the first echocardiography by C.H. Hertz and I. Edler in the 1950s, computerized tomography for the early diagnosis of stroke in the 1970s and magnetic resonance imaging in the early 1980s.
Nowadays, cardiology and cardiac surgery are mature disciplines whose relevance is easily understood considering the clinical, social and economic implications of cardiovascular disorders. In humans, heart-related diseases are the main cause of death worldwide being responsible for the loss of approximately eighteen million people per year (the Covid-19 pandemic has claimed, in the first two years, five million victims!). According to the World Health Organization, these disorders have been the leading cause of death, about one third of the total, globally for the last fifteen years and, although the mortality is declining, the decrease rate is slowing down. The trend for the related morbidity is even worse as, in fact, the number of cardiopath patients remains very high and, in the EU alone, over sixty million people live with a heart disease (Italianer et al. Reference Italianer2021).
Cardiovascular disorders (CVDs) entail dramatic changes in patients lives making it difficult to carry out normal activities, pursue hobbies or maintain social contacts. Patients of working age have their productivity and, consequently, income reduced thus deteriorating their quality of life.
In most of the countries with westernized economies, the incidence of CVDs is approximately $2\,\%$ of the population which costs $7$–$10\,\%$ of all health care expenditure. Owing to ageing of the population and the worsening of co-morbidity factors (smoking, obesity, diabetes), the incidence of CVDs is predicted to increase in the next decades yielding unsustainable health care costs before the end of this century.
Despite the explosive growth of CVD research since the 1970s, advances in clinical practice and treatments appear to have come to a halt, and finding alternative strategies has become imperative (Fiedler Reference Fiedler2015). Novel approaches to the problems are required and the combination of engineering with cardiology is providing unprecedented opportunities both for the comprehension of the heart physiology and for the design and testing of new procedures and devices.
Sophisticated laboratory experiments have been setup to analyse the mechanics of biological tissues, the dynamics of blood in heart chambers and vessels, the fluid/structure interaction of heart valves and the electrophysiology of the myocardium. A similar effort has been devoted to study the interaction of the heart with prosthetic devices and diagnostic tools with the aim of improving the prognosis of heart impairments and their early detection.
The impressive and unrestrainable growth of computer power has further boosted cardiovascular research with computational models, virtual prototypes, digital twins and, more recently, machine learning and ‘big data’ techniques.
The main aim of this paper is to give an overview of the current knowledge and of the main issues in the field with a particular look at the fluid mechanics implications.
The paper is organized as follows. Section 2 describes the heart and circulation of vertebrates and, successively, discusses some allometric laws for the cardiovascular systems of homeothermal mammals. Section 3 illustrates human heart physiology with the main structures and its energy balance. Section 4 is devoted to blood features while § 5 describes some characteristics of cardiac tissues. Experimental and numerical heart models are illustrated in § 6 with a discussion of the main modelling issues. In § 7, we present the main features of cardiovascular flows in physiological and pathological conditions; in particular, in the latter case, we consider cardiomyopathies, prosthetic heart valves and related haemolysis and ventricular assist devices. The article is concluded in § 8 with a discussion of the future outlook and closing remarks.
2. Hearts and circulation
The heart distributes blood to tissues and organs bringing oxygen and nutrients to all the cells and removing their waste. In large organisms, this is efficiently achieved by the circulatory system, a complex network of ‘tubes’ of different cross-sections, which through advection allows the blood to travel along distances of metres within a few seconds: needless to say that without such a pumping organ, only elementary sub-millimetric organisms, with a diffusion-based distribution system, could exist.
In fact, even focussing only on animals with a circulatory system, comparing the hearts of any two of them might be difficult and even misleading on account of all possible variations. The circulation can be open, if the blood circulates freely in the tissues like in insects, or closed, when it flows within veins and arteries. The circulation is said to be incomplete when there is partial mixing between oxygenated and de-oxygenated blood or complete otherwise. Finally, the circulation is single when it crosses the heart only once in a cycle or double if the heart is crossed twice.
In figure 2, we show the circulation layouts for the vertebrate groups showing some analogies but also relevant differences.
Fish have a single-loop circulation with a heart made of one atrium and one ventricle which pumps blood to the gills to be oxygenated and from there to the rest of the body. The venous part of the circulation then collects deoxygenated blood and delivers it back to the heart atrium for a new cycle (figure 2a).
Amphibians have a more complex circulation with two loops, one for the blood oxygenation and another for its distribution to the body. To feed two circuits, the heart has two atria but only one ventricle where arterial and venous blood come in contact although they mix only partially (figure 2b). Reptiles all have a double circulation and an additional ‘aorta’ originating from the pulmonary circuit and shunting the two circulations; the heart structure is not homogeneous with some species (turtles and lizards) similar to amphibians and others (crocodiles) to mammals.
Birds and mammals also have a double circulation but it is complete and the heart has four chambers organized into two separate pumps that feed the pulmonary and systemic circulations without mixing venous and arterial blood (figure 2c).
In all cases, the heart propels the blood to arteries and veins in a precise sequence and the correct flow direction is assured by valves located in the heart and operated passively by the flow itself. The pumping action is achieved by a rhythmic, active contraction of the heart and the work done on the flow balances the viscous energy losses, mainly in the capillaries, within each cycle.
As the body mass and the size of an animal increases, so do the extension of the capillary network, the length of the circuit, total volume of blood, energy losses and therefore the size and the power of the heart.
Of course, depending on the particular vertebrate group, the body-size affects on the circulation are different and making analogies among them might be incorrect. Accordingly, to make meaningful comparisons, here we focus on homeothermic mammals with the aim of verifying how the human heart performs with respect to that of similar organisms and if its features are determined only by biology or they are also altered by external factors.
2.1. Allometric laws
Using a mechanical analogy, being the two pumps and the two exchangers (pulmonary and systemic capillary beds) crossed by the same flow rate (in series), we can think of the mammals cardiovascular system as a closed hydraulic circuit driven by a volumetric pump operating and connected to a single exchanger (figure 2d). Following Dawson (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:
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
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
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.
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
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).
Similar indication comes from figure 4(b) showing the total number of heartbeats in a lifetime ($HB=f\times LE$); it appears that all mammal hearts are made to beat a given number of times $HB \approx 5\times 10^{8}$–$1.5\times 10^{9}$ which can be ‘spent’ in a short lifespan with a fast heart rate or vice versa. Once again, humans are the outliers with their ${\approx }3\times 10^{9}$ beats that almost triples the highest figure of other mammals.
Thus a legitimate question is whether humans are a biological exception or the anomaly is rather due to the improved living conditions in which the environment has been partially adapted to their own needs. To this aim, figure 4(c) reports the human life expectancy during various periods from the Paleolithic Age up to the present day; clearly man has not been always long-lived and the largest $LE$ increases have been achieved only within the last century when medical research has experienced impressive progress. If we report in figures 4(a) and 4(b) the Paleolithic life expectancy (rather then the present one), humans perform just like the others thus suggesting that our heart has evolved biologically like any mammals although our capability to modify the environment and the progress of medical science has considerably improved the overall picture.
Several comments are in order on the significance of the above discussion: animal data are unavoidably affected by large uncertainties since heart rates are often measured on anaesthetized or restrained specimens. Life expectancies vary considerably from wild to captivity and, equally important, life duration is determined by many concurrent factors (diseases, predation, fighting, food shortage) which affect different groups of similar animals to different degrees; the heart performance is just one of these factors and maybe not always be the most important.
Human lifespans are also uncertain and often even biased: those of the oldest ages have been obtained from scarce statistics, while in the Classic and Medieval ages, only selected populations (Greeks–Romans and Europeans, respectively) are represented. In addition, the averages are affected by child mortality and peerage/plebs differences: for example, the late Medieval English population had an average $LE \simeq 30$ years; however, if a member of the aristocracy survived the 21st year, the life expectancy increased up to $64$ years. Within these large scatters, the data can give only an indication although those relative to $1900$, $1950$ and $2017$ are more reliable since they are obtained as world averages.
We can close this section by concluding that the human heart is biologically similar to that of other mammals whose features follow allometric scalings derived from simple basic laws of mechanics. Nevertheless, thanks to better living conditions and the progress of clinical and surgical practice, the total number of heartbeats in a lifetime has almost tripled for humans during the last century: pushing this limit beyond the present value and extending the human life expectancy is the implicit motivation driving the research on the cardiovascular system.
3. Heart physiology and structure
The reductionist approach of the previous section was aimed only at comparing the human heart to those of other mammals and at deriving scaling laws from basic principles. That method, however, is clearly inadequate to describe the heart structure since it lumps all its complex functions into a periodic volumetric variation: in this section, we give a detailed description of the heart physiology and its architecture.
The pumping function of the heart operates synergistically with the connected vessels shown in figure 5. Here, the schematic of the system includes arrows indicating the blood direction and the names of the main elements: since the circulation forms a closed loop, any location can be assumed as a starting point.
The superior and inferior venae cavae (Svc and Ivc) collect deoxygenated blood from the whole body and direct it to the right atrium (RA); from there, blood flows to the right ventricle (RV), initially driven by its inertia and eventually by an active RA contraction that completes the ventricle filling. A contraction of the RV then follows to increase the blood pressure therein thus causing the closing of the tricuspid valve (Tv) and the opening of the pulmonary valve (Pv). This allows the lungs to be fed by deoxygenated blood through the left and right pulmonary arteries (Lpa and Rpa).
Within the lungs, carbon dioxide/oxygen exchange takes place and oxygenated blood returns to the left atrium (LA) through the left and right pulmonary veins (Lpv, Rpv). Similarly to the right counterpart, the LA feeds the left ventricle (LV) and, when it is filled, a vigorous contraction of the latter closes the mitral valve (Mv), opens the aortic valve (Av) and squeezes oxygenated blood into the aorta (Ao), which perfuses the entire body. The head and the superior limbs are fed by the epiaortic arteries (brachiocephalic, left common and left subclavian: Ba, Ca and Sa, respectively) while all the internal organs and the inferior limbs receive blood from the descending aorta (Da).
The above main arteries branch hierarchically into smaller arteries and eventually to arterioles and capillaries which, by diffusion, distribute oxygen and nutrients to all the cells of the body. At the same time, carbon dioxide and waste products are collected by the blood in the capillaries which converge to venules, minor veins and large veins up to the inferior and superior venae cavae that convey the blood again to the right atrium for a new cycle.
The complex cardiac dynamics are the result of coordinated actions of highly specialized parts that operate synergistically: to help with the description of the main mechanisms, we report some representative views of the heart starting from the anterior (figure 6a) which shows the appearance of the heart if the chest were open by a cut of the sternum.
A vertical section, as in figure 6(b), reveals the two pairs of chambers of the right and left heart; an evident difference between atria and ventricles is the thin smooth walls of the former compared with the thick corrugated counterpart of the latter. The reduced thickness of the atria ($\approx$2–3 mm) is due to their function of conveying blood to the ventricles whose filling occurs for approximately $3/4$ of the volume just passively owing to the inertia of the blood stream. Only at the end of the filling process, a synchronous contraction of the atria completes this phase.
The ventricles have much thicker walls ($\approx$4–10 mm) since their contraction must close the atrioventricular valves to prevent atrial regurgitation, and increase the blood pressure to high enough values to open the semilunar valves and feed the pulmonary and systemic circulations. The high ventricular pressures during contraction require the atrioventricular valves to have fibrous strands tethering the free margins of their leaflets to prevent them from everting into the atria (figure 6b). These strands are called cordae tendineae and they are connected to protrusions of the myocardium (papillary muscles) which contract with the ventricle to enhance the pulling action. The inner surface of the ventricles has many grooves called trabeculae carnae; their exact function has not been completely addressed yet although there are indications that they enhance cardiac contractility (Munro et al. Reference Munro, Shen, Ward, Ruygrok, Crossman and Soeller2018).
Figure 6(c) evidences another interesting feature which is the different thickness of the left and right ventricular walls. This is easily understood considering that the right heart drives the pulmonary circulation whose pressure in the largest artery does not exceed the value $20$ mmHg ($2650$ Pa); in contrast, the systemic circulation is much more extended and the mean aortic pressure is approximately $100$ mmHg (13 300 Pa).
These pressures produce cyclic forces, up to several Newtons, on the heart valves (figure 6b,d) and on the tissues. To withstand these high forces, heart valves and tissues connect to a fibrous skeleton, embedded in the plane of figure 6(e), to create a stiff frame which supports the loads.
The muscular part, the myocardium, is a complex composite structure with many different types of cells, the myocytes, which allow the active contraction and the transmission of electrical impulses. Atrial and ventricular myocytes are specialized for contraction triggered by an activation potential which is also transmitted to neighbouring cells.
In addition to the contractile myocytes, there are other cells which have negligible contractility but are specialized as excitatory and conductive fibres. They trigger the heartbeat and assure the correct activation timing is achieved throughout the myocardium: these cells constitute the electrophysiology system, the real ‘orchestra conductor’. In physiological conditions, the heartbeat starts from the sinoatrial (SA) node, made up of self-depolarizing cells, which act as the primary pacemaker (figure 7). The electrical signal then travels, at a velocity of approximately $1\ {\rm m}\ {\rm s}^{-1}$, along the internodal pathways and the bundle of Bachmann thus depolarizing and contracting both atria within $t\approx 80$ ms. The fibrous skeleton is an electrical insulator and prevents the activating potential from moving from the atria to the ventricles; the atrioventricular (AV) node is the only passage for the signal and there the conduction velocity is only ${\leqslant }0.05\ {\rm m}\ {\rm s}^{-1}$, thus delaying the transmission of $\approx$100 ms before moving beyond. This delay is key for the filling of the ventricles which is completed by the atrial contraction while the former are still relaxed.
Once the signal has reached the bundle of His, it speeds up to $2\ {\rm m}\ {\rm s}^{-1}$ and, following the left and right bundle branches in the septum, at $t\approx 190$ ms, the ventricle apex is depolarized. From there, the signal further accelerates along the Purkinje fibres reaching speeds of up to $4\ {\rm m}\ {\rm s}^{-1}$ and rapidly reaches the ventricular myocytes where it propagates at $1\ {\rm m}\ {\rm s}^{-1}$: by time $t \approx 230$ ms, both ventricles are fully depolarized with the contraction lasting $\approx$180 ms. This time interval is called the refractory period during which no new stimuli can trigger a new contraction.
It is worth mentioning that although the SA node is the main pacemaker, it is not the only one and other parts can start the cycle in case the first fails: the AV node has an intrinsic firing rate of $40$–$60$ b.p.m., the bundle of His $\approx$40 b.p.m. and the Purkinje fibres $\approx$20 b.p.m. In physiologic conditions, however, these latent pacemakers do not operate since the SA node ($60$–$80$ b.p.m.) has the fastest frequency and the shortest refractory period, thus setting the contraction rate for the whole heart.
The depolarization and repolarization activity of the heart generates feeble electric fields that can be detected by sensors placed on the skin to produce the typical electrocardiogram tracing (figure 7b). Each wave marks a specific event of the heart activity and cardiologists are trained to recognize pathologic behaviours looking at deviations of the tracing from the standard shape. The first P-wave of figure 7(b) is generated by the atrial depolarization and it has usually a weak intensity ($\approx$0.25 mV). The interval in between the P-wave and the QRS-complex reveals the time taken by the signal to cross the AV-node; the activity is too weak to be detected but its duration gives relevant information. The Q-wave is generated by the depolarization of the interventricular septum, the R-wave is due to the apex region which, owing to its mass, generates the most intense signal. The negative S-wave is produced by the depolarization of the left ventricle region farthest from the apex. Meanwhile, atria have repolarized but their activity is hidden by the ventricles contraction. When ventricles repolarize, they generate the T-wave and the time interval QT shortens when the heart rate increases. Sometimes the T-wave is followed by a small U-wave caused by the repolarization of the papillary muscles which are the latest to relax.
The contraction of the myocardium produces a thrust on the blood and the heart valves ensure the resulting flow is in the appropriate direction. There are two heart valves for each side of the heart (figure 6), one between atrium and ventricle and one separating ventricle and main vessel. Each valve is made of flapping, thin and tough leaflets of connective tissue lined with endocardium, the same coating as the heart chambers, and connected to the fibrous skeleton through a stiff annulus.
All valves are passive and are activated as a result of the hydrodynamic loads, mostly pressure, produced by the blood motion which, in turn, is determined by the valves configuration (fluid/structure interaction). In the right heart, the atrioventricular valve is the tricuspid, it has three leaflets whose free edges are connected to the ventricular papillary muscles by a web of cordae tendineae which prevent their prolapse into the atrium. This valve opens during ventricle diastole, as it relaxes and blood flows from the atrium, while it closes during systole when the contraction increases the ventricular blood pressure; this pressure rise causes the pulmonary (or pulmonic) valve to open and push blood into the common trunk to reach the lungs. As the ventricle contraction ceases, blood pressure therein drops and the tricuspid valve opens while the counterpressure within the pulmonary arteries seals the semilunar valve.
The left side of the heart has very similar dynamics, with the mitral valve between atrium and ventricle and the aortic valve between ventricle and aorta, the only difference being the larger pressures that these valves have to withstand which, as will be discussed in § 7.5, is the reason why they are more subjected to morbidity.
We can better appreciate the synergistic nature of the heart by looking more closely at the dynamics of these valves. In fact, the aortic valve could not operate without a suitable impedance of the downstream circulation and compliance of the aorta which stores part of the ventricle mechanical work via elastic deformation of its walls; this ensures the blood in the aorta remains pressurized during diastole and produces the counterpressure that allows the aortic valve to close. This clearly shows that the heart and aorta operate together and the pumping of the former could not be possible if the latter were replaced by a rigid pipe.
As an aside, we note that during mitral and tricuspid valve closure their leaflets approach while moving towards the atria and even after coapting, they keep retreating until the tension balances the pressure difference. This causes a physiologic backflow called ‘false regurgitation’ which is necessary for the proper sealing of the atrioventricular valves (Collia, Zovatto & Pedrizzetti Reference Collia, Zovatto and Pedrizzetti2019); the same phenomenon, however, produces also a pressure spike that, in the long term, could damage the thin atrial walls. To avoid this problem, each atrium is provided with a histologically distinct appendage which operates as a decompression chamber during ventricular systole or in the event of increased atrial pressure (Al-Saady, Obel & Camm Reference Al-Saady, Obel and Camm1999). While in the past, atrial appendages have been considered an insignificant portion of the cardiac anatomy, today, their role in the heart physiology has been recognized; they have a variety of different morphologies which have pathological relevance in case of atrial arrhythmia since clots tend to form preferentially within appendages (Garcia-Villalba et al. Reference Garcia-Villalba2021).
It has been mentioned several times that the heart distributes oxygenated blood and nutrients to all the body to make possible the metabolic processes. However, the myocardium itself also must be perfused and this happens through the coronary circulation which is another example of the highly intertwined heart dynamics. Even though the heart accounts only for $\approx$0.5 % of the body mass, it absorbs $\approx$5 % of the cardiac output to perfuse the oxygen-eager myocardium.
The aortic root has a bulge with three lobes, the Valsalva sinuses, and from two of them, the left and right coronary arteries depart. These are $3$–$4$-mm diameter arteries that branch and taper up to a network of capillaries with the incredible density of $2500$ per ${\rm mm}^{2}$ of myocardium; the exact configuration of the various branches has a large variability among humans, even if there are some common features of the main tracts which are reported in figure 8(a).
One problem with the heart perfusion is that the intake of the coronary arteries is in the aortic root and blood is pumped into the aorta during systole, when the vigorous myocardium contraction closes the capillary network and prevents the tissues from being perfused. However, thanks to the elastic energy accumulated in the aorta deformation, blood remains under pressure also during diastole, when the ventricles relax and perfusion is possible; this is known as the tissue pressure effect and it is confirmed by direct ultrasound measurements of the blood velocity as reported in figure 8(b) showing the out of phase fluxes. This is another example of interdependent phenomena that could not be possible if the heart did not operate in combination with other vascular elements.
Deoxygenated blood is removed from the heart muscle by several cardiac veins and most of them converge to the coronary sinus, a complex muscular vein whose outlet ends directly into the right atrium with a small (Thebesian) valve. Similarly to the coronary arteries, cardiac veins also have a variable anatomy, such as the presence of the Vieussens valve between the coronary sinus and the great cardiac vein only in $80\,\%$ of subjects and with a number of leaflets that can range from one to three.
Coronary circulation has many more sophisticated regulatory mechanisms since the heart has limited anaerobic capacity and an oxygen shortage of just a few minutes can result in permanent damage. For example, unlike skeletal muscles which use $30$–$40\,\%$ of the available oxygen, the myocardium can extract up to $80\,\%$ of it from the blood. In case an increased flow rate is needed, coronary vasculature can dilate to reduce resistance. A comprehensive discussion of all the coronary circulation details can be found in Goodwill et al. (Reference Goodwill, Dick, Kiel and Tune2018).
If the electric system of the heart is the ‘orchestra conductor’, the Wiggers diagram of figure 9 is the executed ‘musical score’: in the various plots, all the relevant phases of the heart physiology are represented through the pressures in the left and right sides, the volumes of the heart chambers, the ECG signal and that of the phonocardiogram (PCG).
The latter is used to monitor the heart function through its sounds which are produced by the motion of the tissues and the vortical blood flow. The basic approach consists of the traditional auscultation by the stethoscope although the PCG high-fidelity recording also allows the detection of feeble murmurs, and the frequency analysis of the signal permits quantitative assessment of specific phases (Mizuno et al. Reference Mizuno, Koichiro, Shirai and Shiina2015).
Assuming (arbitrarily) the beginning of the heart cycle in the middle of systole, figure 9 shows the blood pressures in the atria just slightly above those in the ventricles that increase their volumes because of the flowing blood. As they expand, the tension of the myocardium increases and the inertia of the stream is not enough to counter its resistance, thus an active contraction of the atria is needed to complete filling the ventricles. The SA-node then triggers an electrical impulse which quickly propagates across the atria and depolarizes them; this is evidenced by the P-wave of the ECG (figure 7b), followed by the pressure increase in the atria and the ventricular expansion to their maximum volume. Meanwhile, the activation potential has crossed the barrier of the AV-node and it can propagate along the interventricular septum (Q-wave of the ECG), the apex (R-wave) and the ventricular walls (S-wave). The myocardium depolarization activates its contraction that initially does not induce a significant volume variation (‘b’ and ‘m’ are isovolumic contractions) although the ventricular blood pressure sharply increases thus inducing the impulsive closure of mitral and tricuspid valves. The smashing of the valve leaflets, produced by the atrioventricular counterpressure, is so violent that its sound is the first we can hear in a heartbeat (the ‘lub’ of the ‘lub–dub’ onomatopoeia) by just putting our ear on a chest and this is the 1st sound of the PCG trace. As the ventricular pressure exceeds the value in the pulmonary trunk and the aorta, the pulmonary and aortic valves, separately, open and blood flows to the vessels. Ventricles volume now decreases reaching a minimum at the end of systole when the repolarization of the myocardium (T-wave of the ECG) starts the diastolic phase. Initially, the muscle relaxation is not accompanied by a volume expansion (‘e’ and ‘p’ are the isovolumic relaxations) although there is a rapid ventricular pressure drop that shuts the semilunar valves thus allowing the blood in the vessels to maintain the high pressure needed to feed the circulations. The valve closure produces the ‘dub’ of the ‘lub–dub’ sound and the second spike of the PCG trace.
As the ventricle relaxation proceeds, its pressure keeps decreasing until it drops below the atrial values and the atrioventricular valves open. The rapid filling of the ventricles starts (referred to as early wave or E-wave) and their rapid volume increases together with the intense blood recirculations sometimes producing a feeble 3rd sound in the PCG. The cycle is now complete and a new one can start.
In addition to the features just discussed, here we note that although the left and right sides of the heart operate jointly, they are not perfectly synchronous owing to different structural details and the finite conduction velocity of the electric signal propagation. For example, the waves produced by the atria contraction have the right one starting before the left since the electrical signal originates from the SA-node located on the right atrium (figure 7). Also the isovolumic contraction and relaxation have shorter durations for the right side owing to the smaller thickness of the tissue. Similar comments apply to the instants in which the atrioventricular valves open or close, when atrial and ventricular pressures cross or in which the semilunar valves operate.
In the bottom panel of figure 9, we have reported only the main events of the PCG referred to as the 1st (S1), 2nd (S2) and 3rd (S3) sounds; in fact, they consist of composite signals which are generated by several concurrent events. For example, S1 is initiated by the contraction of the ventricles (systole) and the loudest components are produced by the closure of the mitral and tricuspid valves. The former precedes by the latter $20$–$30$ ms and this produces two distinct peaks known as splitting of the first heart sound.
Also, the S2 consists of two main components which are the closure of the aortic and pulmonary valve separated by more than $20$ ms that characterizes the splitting of the second sound. There are additional components, although of small intensity, produced by haemodynamic events in the vessels immediately downstream of the semilunar valves.
Finally, there is the rare feeble 3rd sound which is generated by the vortical blood flow sweeping the compliant left ventricle walls during the passive filling phase. Features like pitch, intensity, duration and distance from S2 of this sound are used as diagnostic indicators to detect abnormal ventricular function.
Finally, among the time evolutions of the heart chamber volumes (third plot of figure 9), that of the left ventricle gives relevant information about the heart functionality: the difference between the maximum ($V_{M}$) and minimum ($V_m$) volumes during a cycle yields the quantity of blood $SV$ pumped in a beat, the stroke volume, while the ratio $EF=SV/V_{M}\times 100$ is called the ejection fraction. Values $50\,\% \leqslant EF\leqslant 70\,\%$ are considered physiological and indicate a healthy heart, while $EF < 50\,\%$ evidences impaired pumping capabilities in which blood has a residence time longer than normal in the ventricle; $EF < 35\,\%$ is considered life threatening. The product of the stroke volume times the heart rate yields the cardiac output expressed in litres per minute, $CO=HR\times SV \approx 5\ {\rm l}\ \min ^{-1}$ in healthy subjects.
3.1. Furnace or pump?
The heart, just as any other muscle, converts chemical energy into mechanical work; however, every real transformation has an efficiency $\eta =P_w/P < 1$, with $P_w$ and $P$ mechanical and chemical powers, and thus only part of the energy is used to pump blood and heat is also produced in the process.
From the volume and pressure variations during the cycle reported in figure 9, it appears that myocardium contraction/relaxation can occur in two main modes: the first (isometric) tightens/loosens the fibres without varying their length to raise/lower the pressure in the heart chambers. The second (isotonic) occurs at an approximately constant load and pressure but with decreasing/increasing fibres length and thus varying the chambers’ volume.
A common way to represent the chambers’ dynamics is to report, in a Clapeyron plane ($V$–$p$), the data during the cycle and an idealized example for the left ventricle is shown in figure 10(a); in the $V$–$p$ space, where a perfect isometric transformation is a vertical line while an isotonic one is horizontal. A heartbeat cycle is thus a rectangle and the enclosed area is the ejection work (EW) which goes from the tissues to the fluid if the loop is walked in a counterclockwise direction (vice versa otherwise). Also of interest is the triangular area to the left of the cycle which can be thought of as the potential elastic energy (PE) stored by the tissue contraction which is not directly converted into work. Indeed, owing to the elastic properties of the myocardium, the diastolic filling is not completely isotonic and the end diastolic pressure–volume relation (EDPVR) produces a preload at the beginning of the isovolumic contraction. For the same reason, the systolic emptying of the ventricle also has a similar dynamics, which is however complicated by the interaction with the arterial elasticity downstream of the aortic valve, and their interaction determines the afterload at the end of the systolic ejection.
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−1 which yields a power $P=8.7$ W. This figure is compatible with a coronary flow ${\approx }5\,\%$ of the total and a global heat power of $100$ W dissipated by an adult human when considering that the myocardium can extract more than $50\,\%$ of the oxygen available from the blood while skeletal muscles do not exceed 30–35 %.
We can finally compare the blood pumping power $P_w$ with the total power $P$ to compute the heart efficiency $\eta \simeq 0.2$ implying that $80\,\%$ of the chemical energy used by the heart eventually produces heat and only $20\,\%$ mechanical work.
A final comment is in order about the function of the heart: while there is no doubt that the organ physiology is aimed at pumping blood through the circulations, its heat production is definitely not negligible since it employs $4/5$ of the total power. In retrospect, the old debate about whether the heart is a pump or a furnace does not have a unique answer and we can say that even if Galen and Leonardo da Vinci did not recognize the pumping function of the heart, they were ‘$80\,\%$ correct’ by conceiving it as a heat source!
4. The blood
The synergistic dynamics of the heart applies to all the involved elements and its working fluid, the blood, is no exception. During systole, the compression work of the myocardium squeezes blood into the large arteries whose deformation stores part of the work as elastic energy. During diastole, this energy is returned to blood as kinetic energy of the flow. Concurrently, the same amount of energy is dissipated via friction as, to reach all body cells, the network of arteries, capillaries and veins has a total length of ${\approx }10^{5}$ km and a range of diameters $3$ cm–$5\ \mathrm {\mu } {\rm m}$ through which the flow is continuously sustained.
Haemodynamics is a complex subject and several comprehensive books and reviews have been devoted to it (Fung & Zweifach Reference Fung and Zweifach1971; Cokelet Reference Cokelet1980; Fung Reference Fung1984; Baskurt et al. Reference Baskurt, Hardeman, Rampling and Meiselman2007; Nichols, O'Rourke & Vlachopoulos Reference Nichols, O'Rourke and Vlachopoulos2011). Here we will introduce only the main features and limit the discussion to those involved in the heart physiology.
To accomplish its function, blood has a liquid matrix in which deformable cells are suspended and its dynamics depends on the considered spatial and time scales. The liquid phase is a diluted aqueous solution with proteins and salts, the plasma, and it behaves as a Newtonian fluid of similar viscosity as pure water. The corpuscular fraction is made up of platelets (PTC), white blood cells (WBC) and red blood cells (RBC) with a volumetric fraction, the haematocrit, in the range $35\,\%$–$50\,\%$. Just for comparison, a suspension of rigid spheres at these concentrations shows an effective viscosity $10$–$100$ times that of the liquid matrix, while the whole blood, thanks to cells deformability, has an asymptotic viscosity (for shear rates ${\geqslant }50\ {\rm s}^{{-1}}$) of approximately $4$–$5$ times that of plasma, thus making possible its circulation within the limited mechanical power of the heart. A cubic millimetre (a ${\rm \mu}$l) of whole blood contains ${\approx }3\times 10^{5}$ PTC, ${\approx }7\times 10^{3}$ WBC and ${\approx }5\times 10^{6}$ RBC; therefore, the latter dominate the properties of the cellular phase.
Red blood cells are biconcave circular disks, ${\approx }8\ \mathrm {\mu }{\rm m}$ in diameter and ${\approx }2\ \mathrm {\mu }{\rm m}$ thickness, highly specialized to deliver oxygen from the lungs to the tissues and bring back carbon dioxide. For this purpose, they contain haemoglobin, which forms weak chemical bonds with ${\rm O}_2$ and ${\rm CO}_2$ molecules, within an elastic membrane with a large surface-to-volume area. RBCs are highly deformable to fit into capillaries of diameters smaller than their size, and this is accomplished by a specific cytoplasmic structure in combination with the absence of a nucleus.
Blood rheology is quite complex and highly dependent on the specific flow dynamics. In vessels of diameter larger than approximately $80\text {--}100\ \mathrm {\mu } {\rm m}$, blood can be treated as a unique continuum fluid showing a non-Newtonian shear-thinning behaviour. This property is mostly attributed to RBC deformability with the viscosity values depending exponentially on the haematocrit (Linderkamp, Stadler & Zilow Reference Linderkamp, Stadler and Zilow1992). In figure 11(a), we report some measurements performed in a Taylor–Couette rheometer for several haematocrit values by Brooks, Goodwin & Seaman (Reference Brooks, Goodwin and Seaman1970) with curves that can be well described by the empirical fit
with $|E| = {\rm sym} (\boldsymbol {\nabla }{\boldsymbol {u}})$ the modulus of the rate-of-strain tensor and $\mu _0$, $\mu _\infty$, $\lambda$, $b$ haematocrit-dependent parameters (Siginer, De Kee & Chhabra Reference Siginer, De Kee and Chhabra1999).
For each haematocrit, the blood viscosity has a plateau $\mu _0$ for vanishing shear-rates ($|E|$) then decreases and eventually attains an asymptotic value $\mu _\infty$ for large $|E|$. The initial plateau is caused by the propensity of RBCs to aggregate and form rouleaux, stacks of cells which, in turn, tend to organize into networks thus dramatically increasing blood viscosity at low shear rates. These aggregations however are reversible and they break up when subjected to relatively low mechanical loads thus yielding a viscosity reduction for increasing $|E|$.
It has long been conjectured (Charm & Kurland Reference Charm and Kurland1974) that rouleaux disaggregation progressed until only isolated cells were left and this resulted in the asymptotic viscosity $\mu _\infty$. However, a closer look at the radial flow distribution in a circular pipe reveals that the actual dynamics has additional important features. In fact, even if we assumed as a starting point a simple Poiseuille parabolic profile, this would give a linear rate-of-strain with the maximum at the wall and zero at the pipe axis; this non-uniform strain distribution causes rouleaux disaggregation to be more likely far from the axis thus producing also a radial viscosity gradient pointing towards the pipe axis. Moreover, RBCs deform under the flow shear leading to an accumulation at the axis (Baskurt & Meiselman Reference Baskurt and Meiselman2003); haematocrit values are therefore maxima at the pipe centre and minima near the boundary with the RBC-depleted layer next to the walls which substantially contributes to the effective viscosity decrease for high shear rates.
In figure 11(b), we show the comparison of the steady, laminar velocity profile in a pipe of radius $R$ for a Newtonian fluid with the actual profile for blood at $27\,\%$ haematocrit (Goldsmith & Marlow Reference Goldsmith and Marlow1979); the plug-like behaviour in the region around the axis is the result of the local high viscosity, while the RBC lean layers next to the pipe wall allow for high velocity gradients.
Even if the steady rheological blood properties might look complex, they are not fully representative of real conditions. In fact, blood flow is unsteady as the heart pumps at a given frequency and the large arteries deform as a response of the periodic pressure wave. As a result, RBCs are subjected to time-varying loads which interact with viscous and inertial properties of the cells and with aggregation/disaggregation processes to produce a viscoelastic behaviour (Cokelet Reference Cokelet1980).
Even without considering the elasticity of the walls, the pulsatile flow alone is considerably more complex than the steady counterpart since the solution depends on the Womersley number $Wo=\sqrt {R^{2}\omega /\nu }$, which gives the relative importance of unsteady inertial and viscous forces. For $Wo < 2$, the flow is viscosity dominated and the resulting flow is parabolic with an amplitude modulated in time with pulsation $\omega$. In contrast, for $Wo > 2$, unsteady effects take over and they dominate the flow at the pipe centre while the near-wall region is still viscous; the resulting velocity profiles are then non-monotonic, with inflection points and uneven distribution of shear stresses. For the human arterial system, we have $Wo \approx 10$ in the aorta, $\approx$4 in large arteries, $\approx$1 in small arteries and ${\approx }10^{-2}$ in capillaries thus including a wide variety of different behaviours.
In the tiniest arteries and capillaries, however, still additional factors have to be considered as the vessel diameter becomes comparable with the red blood cell size. This is quantitatively evident from the plot of figure 12(a) showing the effective blood viscosity normalized by that of pure plasma as a function of the vessel diameter $D$. It can be seen that already for $D < 300\ \mathrm {\mu }{\rm m}$, the viscosity ratio starts decreasing even if only for $D \lesssim 100\ \mathrm {\mu }{\rm m}$ a steep drop is observed. The curve has a minimum for $D \approx 8\text {--}10\ \mathrm {\mu }{\rm m}$ and a new increase is found for smaller diameters. This phenomenon, known as the Fahraeus–Lindqvist effect (Fahraeus & Lindqvist Reference Fahraeus and Lindqvist1931), can be explained by looking at the RBC configurations in blood flows within circular pipes of decreasing diameter in figure 13.
The cartoons show that when the vessel diameter $D$ ${\gg }d$, RBCs can arrange in multiple concentric layers and occupy a substantial volume fraction, thus behaving as described previously and producing a viscosity $3$–$4$ times that of the plasma. However, when $D$ is of the same order as $d$, the cells align in a single row thus decreasing the volume occupation and the haematocrit; this spatial arrangement prevents the cells from interacting and aggregating, yielding a viscosity close to that of the liquid matrix. Finally, in capillaries, red blood cells are deformed to fit into vessels with $D< d$; they attain the peculiar parachute shape which increases again the volume occupation and the haematocrit. In this configuration, energy is needed for the deformation of RBCs which slide along the capillary wall; the result is an effective viscosity increase for $D< d$.
This simple explanation is confirmed by measurements of the haematocrit variation with $D$ (figure 12b) indeed showing a minimum for $D \approx d$, referred to as the Fahraeus effect (Fahraeus Reference Fahraeus1929). It is worth mentioning that part of the haematocrit decrease with $D$ might also be caused by the RBC arrangement shown in figure 11(b) since vessels always branch from the side walls of larger ones and, thus, they are fed from blood layers which contain less RBCs than the bulk flow.
In all the above discussion, we have not considered the effect of WBCs on the whole blood viscosity on account of their small volume concentration. They certainly have a negligible role in large vessels, although in the microcirculation, the situation might be different, especially in pathological conditions. In fact, WBC sizes and properties depend on the type and status of the cells; they can be activated by chemical or mechanical stimuli and either migrate out of blood vessels to reach inflamed tissues (extravasation) or transiently block a microcirculation channel. Similar arguments apply to platelets whose size is approximately one-fourth of an RBC and in physiologic conditions contribute negligibly to viscosity. However, they are particularly sensitive to shear stress which activates them making possible their adhesion to damaged vessels endothelium or the formation of clots by the capture of impaired RBCs.
Before concluding this section, we wish to mention that blood also exhibits viscoelastic properties which can become relevant in pulsatile cardiovascular flows. In fact, the reversible aggregation and disaggregation of rouleaux, and the periodic deformation of RBCs under pulsatile shear, leads to cyclic storage and release of elastic energy in addition to the already mentioned shear-thinning viscosity. The former behaviour, however, depends also on flow frequency and on the multi-axial strain components and, therefore, it is usually modelled by a linear combination of modes (from $1$ to $5$), accounting for the different relaxation time scales of the fluid (Campo-Deaño et al. Reference Campo-Deaño, Dullens, Aarts, Pinho and Oliveira2013), each one being governed by an Oldroyd-like differential equation. Every equation, in turn, depends on empirical parameters which are determined by a fitting against experimental measurements; an example of such derivation is illustrated in Campo-Deaño et al. (Reference Campo-Deaño, Dullens, Aarts, Pinho and Oliveira2013) for real blood and for blood analogues.
5. The tissues
Heart tissues are extremely composite both in terms of anatomy and properties which are anisotropic and heterogeneous; they are made of different types of specialized cells organized into unique structures aimed at performing specific functions. On account of such a complexity, a thorough description of all their features within a JFM Perspectives paper (written by a fluid dynamicist) is unrealistic. Therefore, we will discuss only those aspects needed to understand the pumping function of the heart and leave further details to specialized textbooks such as Katz (Reference Katz2011) or Piper & Preuse (Reference Piper and Preuse1993).
The heart is enclosed in a non-compliant bag, the pericardium, whose external surface is in contact with the surrounding organs, mainly the lungs. The internal surface contains a small amount of fluid which reduces the friction with the epicardium, the outermost layer of the contracting heart.
The active muscular part, the myocardium, is the main and thickest layer consisting of cardiomyocytes, complex cells mechanically and electrically connected through special intercellular bridges called intercalated discs (figure 14a). There are many types of myocytes depending on the function and heart region. ‘Working myocytes’ are the contracting units, and they are mostly found in ventricles and atria respectively of bigger and smaller dimensions. A different type of myocyte are the ‘Purkinje cells’ which do not contract much but are specialized for fast conduction of electrical signals. ‘Nodal cells’ are yet another different myocyte responsible for pacemaker and conduction delay activities to initiate and coordinate the contraction of the various regions.
The myocardium is also rich in fibroblasts which outnumber cardiomyocytes even though they account for a smaller volume fraction owing to their reduced dimensions; they are responsible for the maintenance of the extra cellular matrix, a complex network of collagen and elastin, which is the scaffold of the cardiac tissue and determines most of its passive mechanical properties.
The innermost layer, the endocardium, is a connective tissue lining heart chambers and valves; being directly in contact with blood, it is biologically similar to the endothelial cells of veins and arteries. An obvious function of the endocardium is to act as a physical barrier which protects the myocardium from the vigorous blood flow. Another important function is to dynamically sense the signals coming from the flow (mainly pressure and shear stresses) to stimulate and produce lineages of different cardiovascular cell types like valves, endothelial cells of atria and ventricles or the lining of coronary vasculature (Dye & Lincoln Reference Dye and Lincoln2020).
The building block of the active heart muscle is the cardiac myocyte whose electrical excitation and mechanical contraction have been thoroughly investigated and can be shortly explained by referring to figure 14(b). At the basic resting state, the cell has a negative transmembrane electrical potential which prevents its contraction. When stimulated by neighbouring cells, membrane channels open and allow sodium, calcium and potassium ions to move across the intracellular and extracellular space beginning the positive-feedback mechanism of depolarization and activating the cell contraction. This ion flux (mainly ${\rm Ca}^{{2+}}$) increases the transmembrane potential to positive values, depolarizes the cell and induces contraction. The rapid depolarization is followed by a long plateau phase in which the potential slowly declines for more than $150$ ms. The repolarization phase starts immediately after with the closure of the ${\rm Ca}^{{2+}}$ channels and the attainment of the initial negative potential after approximately $250$–$300$ ms. When the depolarization process is initiated, the cell enters a refractory period during which it can not be excited again and this avoids tetany which would result in permanent contraction (incompatible with life). It is worth mentioning that the graphs of figure 14(b) are just illustrative examples since the specific timings and amplitudes of action potential and mechanical tension are dependent on the cell type and, therefore, on the region of the heart. Each myocyte depolarizes when it is reached by an electrical impulse and, in turn, it generates a new activation signal which is transmitted to the neighbouring cells through the junctional intercalated disks (cell-to-cell conduction). This implies that the electrical stimulus, and the induced muscle contraction, propagate preferentially along precise directions which make the overall tissue behaviour very anisotropic.
The first source of anisotropy is the cardiomyocyte itself, a cylinder of diameter $15\text {--}25\ \mathrm {\mu }{\rm m}$ and length $50\text {--}150\ \mathrm {\mu }{\rm m}$ which shortens its axis upon contraction; although it reduces its length by no more than $10\,\%$, the presence of different layers and the particular orientation of the fibres (figure 15) allow the ventricles to reduce their volumes by up to $70\,\%$ during systole.
Locally, myocytes align the axes along a common direction and establish end-to-end junctions to form fibres. These, in turn, are organized into $4$–$6$ cells thick plane sheets (‘laminae’) separated by layers of connective tissue (Clayton et al. Reference Clayton, Bernus, Cherry, Dierckx, Fenton, Mirabella, Panfilov, Sachse, Seemann and Zhang2011) to reduce their coupling. These laminae have different orientations in space and across the myocardium thickness thus yielding a very anisotropic structure as hinted by figure 15. In particular, Pashakhanloo et al. (Reference Pashakhanloo, Herzka, Ashikaga, Mori, Gai, Bluemke, Trayanova and McVeigh2016) report very recent measurements showing that the orientation of the fibrous bundles in the atria is homogeneous across the tissue thickness in some regions although in others, the orientation angle changes by up to $90^{\circ }$. In contrast, the ventricles are more complex having more fibre layers with an orientation which rotates continuously from the endocardium to the epicardium by up to $150^{\circ }$. This heterogeneity has implications both for the electrical conduction and the mechanical response of the tissue which considerably complicate its macroscopic description (Guan et al. Reference Guan, Yao, Luo and Gao2020).
Although the individual behaviour of the cardiac myocyte is well understood and characterized, a healthy human left ventricle consists of ${\approx }5.5\times 10^{9}$ myocytes (Beltrami et al. Reference Beltrami2001) and modelling even a single heart chamber starting from a single-cell description is practically impossible. Continuum models are thus needed in which the collective mechanical dynamics of the tissue is prescribed and the myocardium is often modelled as an orthotropic, hyperelastic material although it shows also viscoelastic features evidenced by hysteretic loops in multiple load cycles.
Most of the available information on the mechanical properties of heart muscle has been obtained by experimental tensile tests of small myocardium specimens under different conditions. As shown by Fung (Reference Fung1993), the curve of the tensile stress $\sigma$ (defined as the ratio of the force with the resting cross-section) versus the stretch $\lambda$ (the ratio of the deformed with the resting lengths) depends on many parameters such as temperature, pH of the environment, timings and frequencies of the load, just to mention a few. Nevertheless, in physiological conditions, a recurrent feature is a slope of the stress–stretch relationship which grows linearly with the stress $\sigma$ and can therefore be written as ${{\rm d}\sigma }/{{\rm d} \lambda } = \alpha (\sigma +\beta )$ with $\alpha$ and $\beta$ being empirical tuning parameters; upon integration, it yields $\sigma +\beta = C \times e^{\alpha \lambda }$ with the constant $C$ determined by the imposition of a stress $\sigma _0$ for a stretch $\lambda _0$. Examples of this hyperelastic nonlinear behaviour are given in figure 16 with the data for biaxial tension tests of the human myocardium, fitted with a law similar to the one above (Zhao Reference Zhao2018). If the myocardium specimen is small enough, there is a clear fibre- and a perpendicular cross-fibre-direction and they show different mechanical properties, as is evident in figure 16.
It is worth mentioning that the specific values of the stress–strain relationship depend also on the heart region from which the specimen is extracted as the myocytes, the composition of the layers and the orientation of the fibres across the tissue thickness are quite different depending on the function of the tissue.
In addition to the inhomogeneous character of the myocardium and the natural variability among humans, other factors affecting the measurement of the mechanical properties are the preparation and conservation of the samples and the test method (uniaxial tensile, biaxial tensile, triaxial shear, constant versus time-varying loads, etc.); the interested reader is referred to Sperelakis (Reference Sperelakis1989) and Holzapfel & Ogden (Reference Holzapfel and Ogden2009) for more detailed discussions.
6. Heart models
On account of the inherent heart complexity, cardiovascular research relies heavily on models to study blood flow and tissue dynamics in physiological and pathological conditions. In fact, models are more accessible than the real counterparts, they can be used to test hypotheses or predictions, to validate assumptions or to evaluate the performance of devices. However, for the results to be meaningful, a model has to replicate as close as possible the real system and for the heart, this might be exceedingly difficult. For this reason, hardware laboratory models usually aim at reproducing only specific parts of the heart or single details of the dynamics. Computer models are instead more flexible although at the price of huge computational costs if high-fidelity is required.
6.1. Laboratory experiments
A fundamental pillar of any meaningful experiment is its dynamic similarity with the real problem which is verified only for precise sets of model parameters. Let $q$ be a generic quantity of interest, which will depend through an unknown function $f$ on several variables; for cardiovascular flows, these can be restricted to density $\rho$, viscosity $\mu$ and velocity $U$ of the blood, a spatial scale $L$ of the problem, its periodic time scale $T$ and also some measure of the elasticity of the tissues $E$, their density $\rho _T$ and a reference (intramural) pressure difference $\Delta p$ between the inner flow and the outer environment. It is worthwhile to point out that this is just a small set of variables selected to discuss the main issues in the design of an experimental test although, depending on the specific problem, additional quantities or different sets might come up.
Another important observation is that the geometrical similarity is a fundamental prerequisite of the dynamic similitude. Thus, the laboratory model must be a faithfully scaled replica of the real system. The above length scale $L$ must be therefore representative of all the lengths which completely define the problem geometry; they can easily amount to several tens (for example, one of them is $\lambda _m$ in figure 17) and be difficult to downscale.
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.
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$.
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).
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.
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.
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.
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.
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).
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.
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).
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.
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):
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:
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:
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
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.
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:
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
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.4. Fluid–structure–electrophysiology interaction
In principle, the double time integration of (6.4) or (6.3a,b) yields the deformation field and the configuration of all heart tissues whose wet surface velocity can be used to impose the boundary conditions in (6.5a,b) for the haemodynamics while the volume, with its hierarchy of one- and two-dimensional structures, provides the domain of the electrophysiology equations. However, evolving in time, the structure entails the knowledge of the fluid loads whose behaviour depends, in turn, on the structure configuration itself. This ‘chicken and egg’ problem is known as a fluid/structure interaction (FSI) which, for cardiovascular problems, is further complicated by the addition of the electrophysiology (FSEI).
Given the strongly interconnected dynamics of the three systems, it is not obvious whether they must be solved simultaneously or sequentially and, in the latter case, from where to start.
In the first case, all models are solved together as a single large dynamical system through an iterative procedure (strong coupling), while in the second, the output of each model is used as the input for the successive one (loose coupling) in a given order. Strong coupling is stable and robust although it is computationally expensive since it requires several iterations among the solvers (de Tullio & Pascazio Reference de Tullio and Pascazio2016) for every time step ($\Delta t$). However, loose coupling needs only one model solution per $\Delta t$ although it is prone to become unstable and poses strong limitations on the time step size.
Which of the two approaches is computationally more efficient depends on the specific problem and its physical parameters; Borazjani, Ge & Sotiropulos (Reference Borazjani, Ge and Sotiropulos2008) have shown that, when added mass effects are prevalent, loose coupling FSI becomes unstable and strongly penalizes the integration time step. In fact, in a heart model, valve dynamics is dominated by added mass effects as the intrinsic mass of valve leaflets is negligible with respect to that of the displaced fluid and this would suggest that strong coupling FSI is needed. However, Viola et al. (Reference Viola, Meschini and Verzicco2020) have shown that, for the heart, the maximum $\Delta t$ is anyway limited by the fast electrophysiology dynamics and the high elastic frequencies of the stiffened myocardium during contraction. As a consequence, within these limitations, a loose coupling FSI maintains the integration stability and provides computational time savings of the order of $50$–$70\,\%$ while giving the same results as strong coupling FSI.
In contrast, the model of de Tullio et al. (Reference de Tullio, Cristallo, Balaras and Verzicco2009) considered the dynamics of a bi-leaflet mechanical valve in a rigid aortic root and only the fluid/structure interaction time scales were present; according to the above arguments, strong coupling FSI was used for the time integration with a number of iterations between fluid and structure solvers of the order of $4$ each step.
6.2.5. Boundary conditions
As noted in § 2, human circulation is closed thus only boundary conditions at the blood/tissues interfaces are needed to determine the haemodynamics of the heart. However, modelling the whole circulation with the same detail as the heart is practically impossible, and therefore the discretization is truncated at some level of the main vessels and the missing part is replaced by additional models (figure 32). These usually take the form of boundary conditions mimicking resistance, inertia and compliance of the eliminated circulation circuits and their choice is of paramount importance for the quality of the results (Vignon-Clementel et al. Reference Vignon-Clementel, Figueroa, Jansen and Taylor2006; Esmaily Moghadam et al. Reference Esmaily Moghadam, Bazilevs, Hsia, Vignon-Clementel and Marsden2011).
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.
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. Cardiovascular haemodynamics
7.1. Physiological conditions
The periodic beating of the heart produces a pulsatile blood flow whose correct direction, from atria to ventricles and then to the main arteries, is assured by the atrioventricular and semilunar valves. These are passive structures which move mainly owing to fluid inertia and pressure differences produced by the active contraction of the myocardium and the energy stored by the elastic stretching of deformable tissues.
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.
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).
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.