Cardiovascular fluid dynamics: a journey through our circulation

Abstract

This article presents a broad overview of the fluid mechanics of the human cardiovascular system. Beginning in the heart, we travel through the main features of our circulation to highlight important functions and diseases where fluid mechanics plays a central role. Of particular focus is the role of computational modelling in uncovering the dynamic flow phenomenon throughout our body, its association with cardiovascular disease mechanisms and progression and its importance in clinical treatment planning. We also emphasize the multiscale nature of the cardiovascular system, and associated challenges. The main aim of this review is to highlight progress and ongoing challenges in our understanding of cardiovascular haemodynamics, as well as the future outlook for translating the current state-of-the-art to widespread clinical application and improved patient outcomes.

Impact Statement

Cardiovascular disease is the leading cause of death worldwide. Fluid dynamics plays a central role in the functioning of the human cardiovascular system, as well in the progression and treatment of its associated diseases. This review highlights advances in our understanding of cardiovascular fluid dynamics, and ongoing challenges in translating engineering insights into clinical treatments for cardiovascular disease.

1. Introduction

The significance of the human cardiovascular system has intrigued philosophers and scientists for at least 3500 years. This is evidenced by the description of the metou system of vessels driven by central pulsations in the ancient Egyptian Ebers Papyrus, which is the earliest known documentation of humans’ knowledge of their circulation, dated 1550 BC (Ritner 2006). Writings on the cardiovascular system are also present in various other ancient schools of thought, including by Greek physicians of the Hippocratic Corpus, the book Nei Jing in ancient Chinese medicine and the Charaka Samhita from ancient Indian Ayurvedic medicine (Fung 1984; Patwardhan 2012; Lemonnier et al. 2017).

We now know that the adult human circulation system consists of billions of blood vessels spanning approximately 100 000 km in length, through which blood is driven by the pumping of the heart (Jones 1969). The heart pumps roughly once every second, or approximately 2 billion times over the average lifespan, through blood vessels with diameters ranging from roughly 3 cm to 5 $\mathrm {\mu }$m (Secomb 2017; United Nations 2022). The dynamics of blood flow within vessels, which occurs on the scale of seconds, has implications for cardiovascular health on the scale of years (Beere, Glagov & Zarins 1984; Giddens, Zarins & Glagov 1993). This wide range of values underscores the highly multi-scale nature of this interconnected system, which, from the fluid dynamics perspective, features phenomena that simultaneously occur at Reynolds numbers of $O(10^3)$ in the aorta and $O(10^{-3})$ in capillary networks (see table 1).

Table 1. Approximate values of the Reynolds number and Womersley number in different anatomical regions. Womersley numbers were estimated based on a pulsatile flow frequency of 1 Hz. Reported Reynolds numbers are the mean and peak values within cardiac cycles, based on the following references: Aortic flows – Ku (1997), Stein & Sabbah (1976), Stalder et al. (2011), Poelma, Watton & Ventikos (2015) and Ha et al. (2018). Coronary arteries – Kajiya et al. (1987), Perktold, Nerem & Peter (1991), He & Ku (1996) and Zhang et al. (2015). Carotid arteries – LoGerfo et al. (1981), Ku & Giddens (1983) and Lee et al. (2008). Pulmonary arteries – Sloth et al. (1994). Inferior vena cava – Joseph, Voit & Frahm (2020). Arterioles and capillaries – Secomb (2017).

In the last few decades, fluid dynamics and its effect on surrounding tissues and micro-scale biological responses have emerged as key factors in our understanding of cardiovascular function, its associated diseases and their treatments (Ku 1997). For example, diseases related to atherosclerosis, which cause the narrowing of certain arteries and are the leading cause of death globally, have been significantly correlated with regions of low and oscillatory shear stress as well as the direction and topology of the shear stress field exerted on vessel walls (Giddens et al. 1993; Arzani & Shadden 2018; Kok et al. 2019; Morbiducci et al. 2020; Hoogendoorn et al. 2021; Mazzi et al. 2021). Complex flow features, such as vortices, jets and helical flow, have been reported in blood vessels throughout our bodies and are crucial determinants of heart, valve and vascular function, as well as the development of life-threatening blood clots, tears and other diseases in vessel walls (Sforza, Putman & Cebral 2009; Morbiducci et al. 2011; Seo et al. 2014; Liu et al. 2015; Sotiropoulos, Le & Gilmanov 2016; De Nisco et al. 2019). While these interactions between blood flow and cardiovascular health are predominantly analysed in the context of Newtonian fluid dynamics, blood ceases to behave as a continuum as we progress from these inertial flows in large vessels to Stokesian flows in the microcirculation. At these scales, the dynamics of suspended cells and other constituents plays several important roles – for example, in deadly diseases such as sickle cell disease and malaria, as well as in the formation of blood clots, which can be beneficial in preventing bleeding but also dangerous when they impede flow within blood vessels (Freund 2014; Secomb 2017).

In many of these scenarios, flow-induced forces play an important role in mechanobiological sensing at the cellular level, which elicits cell signalling processes that impact long-term disease progression (Humphrey 2008). Moreover, in the treatment of several cardiovascular diseases, studies have shown improved outcomes from treatments guided by haemodynamic metrics, such as flow and pressure, rather than anatomical imaging alone (Powers 1991; Pijls et al. 2010; Taylor, Fonte & Min 2013). The recognition of the importance of haemodynamics to cardiovascular health has led to the rapid growth of image-based computational fluid dynamics (CFD) models, which simulate haemodynamics in anatomical models based on patient-specific imaging data, to analyse normal and pathological blood flow in the cardiovascular system, and correlate this haemodynamics with diseases and treatments (Taylor, Hughes & Zarins 1998; Marsden 2014; Mittal et al. 2016).

This review aims to take the reader on a journey through the intricacies of the human cardiovascular system (figure 1). Along the way, key flow phenomena of clinical and academic interest will be discussed, and the multi-scale nature of the circulatory system will be highlighted. Although much of cardiovascular health is a result of the interaction between blood flow, tissue mechanics, cellular function, genetics and other environmental and lifestyle factors, we focus here on phenomena of interest to fluid dynamics practitioners. We will also draw attention to computational modelling techniques that have advanced our understanding of cardiovascular haemodynamics and health, as well as challenges thereof.

Figure 1. An overview of the key features in the human cardiovascular system that we discuss in this review. The left panel shows a cross-section of the heart with relevant anatomical features labelled. The right panel schematically shows the order in which we will visit the different parts of the cardiovascular system.

The journey begins at the heart of the circulatory system, so to speak, in its main pumping chamber called the left ventricle.

2. Left ventricle

The left ventricle is one of the four heart chambers, and (in the normal circulation) one of two pumps that push blood on its long journey to ensure the nourishment and oxygenation of organs and tissues. It is located in the heart's lower left portion and is roughly conical in shape, with a thick muscular wall that is essential for generating the necessary force to pump oxygen-rich blood throughout the body. During one cardiac cycle, the left ventricle undergoes roughly four phases: passive filling (diastole), isovolumetric contraction, ejection (systole) and isovolumetric relaxation. In the diastolic phase, the ventricular muscle is relaxed, with oxygenated blood flowing in from the left atrium. The mitral valve, located between the left atrium and the left ventricle, is open, while the aortic valve, which connects the left ventricle and the aorta, is sealed due to the higher pressure in the aorta. As the ventricular muscles contract, the pressure increases sharply from 3–12 to 100–140 mmHg within a short time span of 200 ms. The volume of blood within the chamber remains constant since two valves are closed in the isovolumetric contraction phase. Then, as the left ventricular pressure exceeds that in the aorta, the aortic valve opens and blood is pumped into the aorta. A subsequent muscle relaxation lowers the pressure, and after isovolumetric relaxation, the mitral valve reopens, refilling the left ventricle. This cycle produces what is known as the ‘pressure-volume loop’, with the enclosed area quantifying the work done during each heart beat.

Due to the substantial change in one cardiac cycle, the flow inside the left ventricle is strongly unsteady. For example, during the rapid passive filling phase, the mitral jet, emanating from the mitral valve (diameter ${\approx }2\ \text {cm}$), generally has a peak velocity of 0.8 $\text {m s}^{-1}$ and attains a Reynolds number of 4000. For a classical jet, this would entail a turbulent dynamics. An annular vortex is formed downstream of the mitral orifice resulting from the roll-up of the shear layer emanating from the leaflets of the mitral valve. This three-dimensional vortex ring develops into a large-scale, asymmetric, clockwise rotating vortex at the centre of the left ventricle by the end of diastole (Le et al. 2012; Le & Sotiropoulos 2012; Fortini et al. 2013). This facilitates filling by reducing convective losses (Charonko et al. 2013). Its rotation also supports the ejection of blood flow toward the aortic outflow track during systole with minimal flow energy dissipation (Pedrizzetti & Domenichini 2005).

Considering this intricate haemodynamics, both in vivo imaging as well as computational modelling offer promising ways to non-invasively quantify left ventricular haemodynamics and assess both healthy and pathological behaviour. The imaging of such cardiac flows is traditionally performed using Doppler echocardiography and, more recently, using ‘four-dimensional (4-D) flow MRI’. The latter is a relatively recent development in cardiac magnetic resonance imaging (MRI) that allows the imaging of time-resolved flows and anatomies in three dimensions (Markl et al. 2012; Stankovic et al. 2014; Azarine et al. 2019; Garcia, Barker & Markl 2019). Four-dimensional flow MRI has been used extensively in the imaging of cardiac flows as well as flows within several large blood vessels to identify aberrant flow patterns and their relationship to cardiac, vascular and valvular diseases. While current clinical use remains confined to large academic centres, new imaging protocols and analysis techniques are poised to move this imaging into the clinical standard of care. However, a challenge with non-invasive in vivo flow imaging is its limited resolution, which does not allow fine-grained imaging of all flow features as well as flows in smaller blood vessels. Computational modelling is therefore a viable alternative to quantify left ventricular haemodynamics accurately and non-invasively.

However, performing realistic and accurate simulations still poses considerable challenges; accurate simulations of left ventricular flow require patient-specific geometries, physiological wall motion generated by active tissue mechanics and fully coupled fluid–solid interaction with the inclusion of the mitral and aortic valves at the inlet and outlet. Significant progress has been made by using various simplifying assumptions to reduce the complexity of these models given specific research interests. For example, the flow can be simplified using lumped-parameter models (see § 7) to efficiently obtain pressure and volume curves, but these models are inherently incapable of providing the full flow field inside the chamber (Hirschvogel et al. 2017; Augustin et al. 2021). The geometry of the left ventricle can be idealized (Cheng, Oertel & Schenkel 2005; Domenichini, Pedrizzetti & Baccani 2005; Meschini et al. 2018; Dedè, Menghini & Quarteroni 2021) or patient specific (Le & Sotiropoulos 2012; Vedula et al. 2016; Hirschvogel et al. 2017). In this context, it is worth noting that the interior surface of the cardiac chambers is not smooth. It features papillary muscles, trabeculae of different sizes and false tendons, presenting a sponge-like structure extending into the chamber (see § 11). From simulations that have incorporated these detailed endocardial structures from in vivo imaging, they are thought to play a role in ensuring complete closure of the mitral valve, efficient filling and ejection of blood as well as reduced stress in the myocardium wall (Madu & D'Cruz 1997; Kulp et al. 2011; Chnafa, Mendez & Nicoud 2016; Lantz et al. 2016; Vedula et al. 2016; Sacco et al. 2018b). Another simplification in studies that have focused on ventricular filling is the absence of the mitral and aortic valves (Cheng et al. 2005), or simplified valves incorporated as resistive lumped-parameter elements (Dedè et al. 2021). More recently, 3-D valves, simulated using Arbitrary Lagrangian–Eulerian or immersed boundary methods, have been used to capture realistic flow dynamics in the healthy or pathological left ventricle (Vedula et al. 2016; Meschini et al. 2018; Meschini, Viola & Verzicco 2019). Finally, the ventricular wall motion can also be simplified as prescribed motion on moving meshes or boundaries constructed from imaging data (Domenichini et al. 2005; Chnafa, Mendez & Nicoud 2014). However, to model the contraction and relaxation of tissues accurately, fluid–solid interaction simulation of the left ventricle can also be expanded to include electrophysiology and active contraction (Vigmond et al. 2008; Quarteroni et al. 2017; Viola, Meschini & Verzicco 2020). We recommend the review of Verzicco (2022), which describes the physiology and numerical methods coupling electrophysiology, solid and fluids of the whole heart in detail.

The high Reynolds number in ventricular flow raises questions about the validity of the laminar flow assumption. Research, including both in vivo studies (Zajac et al. 2015) and numerical methods (Chnafa et al. 2014; Dedè et al. 2021), indicates that turbulence occurs during late diastole and is initiated by the breakup of large vortex structures. Turbulence models, such as large eddy simulation and Reynolds-averaged Navier–Stokes, have shown improved results in left ventricular haemodynamic analyses compared with laminar flow models and experimental data, due to their ability to capture transient or turbulent states involving small-scale phenomena that cannot be retrieved under the laminar hypothesis (Chnafa et al. 2016; Jahanzamin, Fatouraee & Nasiraei-Moghaddam 2019). In fact, a recent study by Saqr et al. (2020) demonstrated that turbulence is inherent to cardiovascular flows, not only in the left ventricle but throughout the body. However, it is important to note that fully developed turbulence is not typically observed under physiological conditions (Chnafa et al. 2014; Dedè et al. 2021). Large eddy simulation approaches are generally better suited for transitional flow regimes than Reynolds-averaged Navier–Stokes approaches, since the latter assume that turbulence is fully developed. In the context of finite-element modelling, the variational multiscale method (Hughes et al. 1998; Hughes, Mazzei & Jansen 2000; Bazilevs et al. 2007), which is a stabilization technique to control numerical oscillations associated with high Reynolds numbers, models turbulence in a similar manner to large eddy simulation and has been widely applied in both cardiac flow (Quarteroni et al. 2017; Dedè et al. 2021; Zingaro et al. 2021) and cardiovascular flow simulations (Hughes et al. 2020; Liu et al. 2020).

Left ventricular haemodynamics is significantly affected by the functioning of the mitral valve, which ensures the correct direction of blood flow from the left atrium to the left ventricle. Mitral valve disease is the most common valvular heart disorder, particularly in ageing populations, with a prevalence of more than 10 % in people over 75 years (Nkomo et al. 2006). The two main types of mitral valve disease are mitral valve stenosis and mitral valve regurgitation. Mitral valve stenosis occurs due to the narrowing or stiffening of the mitral valve, which prevents it from fully opening. As a result, flow from the left atrium to the left ventricle is obstructed and, in turn, the heart pumps less efficiently. Mitral valve stenosis causes a pressure increase in the left atrium, which induces atrial wall thickening and often atrial dilatation. The augmented atrial pressure may also produce blood congestion in the pulmonary veins (Maeder et al. 2018). From CFD simulations, the mitral jet, generated during diastole, is seen to shrink in size and strengthen when the stenosis gets more severe. As a consequence, the kinetic energy of the flow, the shear stresses acting on the tissues and the transvalvular pressure drop increase (Meschini et al. 2019). This also increases mitral valve regurgitation, where the incomplete closure of the mitral valve allows blood to flow backward into the atrium. As mitral valve regurgitation worsens, the heart must work harder to pump blood to the body, leading to a progressive increase in left ventricle volume with normalization of wall stress. In terms of the fluid dynamics, large-scale flow recirculation, smaller vortices scattered throughout the left ventricle and the presence of a regurgitant jet which promotes increased turbulence in left ventricle have been observed in patients with mitral regurgitation (Al-Wakeel et al. 2015; Bennati et al. 2023).

Common therapies for valvular diseases include valve repair or replacement with prosthetic heart valves. Present-day prosthetic heart valves can be classified into two broad categories, rigid-leaflet or mechanical valves, and valves with flexible leaflets made either of synthetic materials or biological tissue (the latter are called bioprosthetic valves). However, all known heart valve prostheses have been associated with various clinical complications due to the non-physiological flow patterns they induce upon implantation. Experimental techniques such as laser Doppler velocimetry and particle image velocimetry have been extensively employed to study prosthetic heart valve flows in vitro, showing that the transmitral vortex formation could be influenced by the saddle annulus dynamics and also by the angle of valve opening (Kheradvar & Falahatpisheh 2012). Fluid–solid interaction simulations have also revealed drastically different flow fields created by different prosthetic heart valves, and identified non-physiologic flow patterns causing blood cell damage and platelet activation (Borazjani 2013). Moreover, studies have shown that, with impaired left ventricular function, the flow disturbances introduced by the bi-leaflet mechanical valve reduce the penetration capability of the mitral jet and weaken the recirculation in the ventricular apex (Meschini et al. 2018). These effects are discussed in detailed prior reviews on heart valve fluid mechanics (Peskin 1982; Yoganathan, Chandran & Sotiropoulos 2005; Sotiropoulos et al. 2016).

Another pathology of the left ventricle is cardiomyopathy (figure 2a), which refers to a condition where the heart muscle becomes enlarged, thickened or stiff, consequently diminishing the heart's ability to pump blood and potentially leading to heart failure. The most prevalent type is dilated cardiomyopathy, which is characterized by progressive enlargement and deterioration of the left ventricle's contraction capacity. The reduced cardiac function is reflected in a weaker diastolic mitral jet which penetrates less into the wider left ventricle, and a smaller vortex in the centre of the cavity that does not support flow redirection toward the outflow, resulting in increased residual volume during pumping (Mangual et al. 2013; Eriksson et al. 2016; Riva et al. 2022).

Figure 2. (a) Average flow and coagulation patterns in the left ventricle for two patients with cardiomyopathy, from patient-specific simulations of flow coupled with biochemical coagulation reactions. Note the different flow fields and coagulation patterns, highlighting the importance of patient-specific risk. Courtesy of J.-H. Seo and R. Mittal, Johns Hopkins University. (b) Simulated flow through normal and bicuspid aortic valves. Note the asymmetric jets through bicuspid valves. The middle panel shows the more common phenotype (fusion of the left and right coronary cusp leaflets) and the right panel shows the less common phenotype (fusion of the right and non-coronary cusp leaflets) (Kaiser et al. 2022).

Left ventricular assist devices (LVADs) have become a key element for the treatment of patients suffering from terminal heart failure. These devices serve different purposes, such as a bridge to transplant (Miller et al. 2007; Aaronson et al. 2012), a bridge to recovery or as long-term therapy (Slaughter et al. 2009). The LVADs pump blood from the left ventricle to the aortic root, providing a mechanical counterpart that works in conjunction with the patient's own left ventricle. Modern LVADs employ continuous flow physiology, which means they continuously pump blood throughout the cardiac cycle. The LVADs provide significant unloading of the left ventricle, manifested by decreased end-diastolic pressure, end-diastolic volume and pressure-volume area (Uriel et al. 2016). Simultaneously, a decrease in the left ventricular volume is evident almost immediately in patients on LVAD support (McCarthy et al. 1995). Besides the left ventricle, device implantation substantially modifies the haemodynamic conditions within the aorta. Due to the high velocity introduced by the pump, the aortic valve may remain persistently closed, especially when native cardiac ejection is low. Moreover, the adverse event profile of contemporary left ventricular assist devices includes suction events (where devices unload too much so that the ventricular chamber collapses) (Vollkron et al. 2007; Griffin & Katz 2014), thrombosis within the device (Eckman & John 2012), right ventricular failure, device-related infections, stroke, dysrhythmia and aortic incompetence (Park et al. 2012; Kirklin et al. 2015). Computational fluid dynamics has been employed to characterize and optimize numerous different designs for LVAD functional performance towards improved long-term outcomes (Fraser et al. 2011). Computational fluid dynamics results identify that, with device support, the aortic root is largely a region of vascular stasis, with significantly higher shear stress on the leaflet tips in patients with aortic incompetence (Shad et al. 2021). Strong correlations between device outflow graft orientation and wall shear stress, pressure and turbulent energy dissipation in the ascending aorta are found. Based on these flow patterns, a more distal insertion has been recommended to produce at least occasional parallel flow through the aortic valve and ensure intermittent aortic valve opening (Osorio et al. 2013; Karmonik et al. 2014).

3. Aortic flow

Oxygenated blood pumped by the left ventricle to the rest of the circulation first passes through the aorta, the largest blood vessel in the human body at roughly 3 cm diameter. Continuing our Lagrangian journey, from the outlet of the left ventricle, the aorta ascends towards the thoracic region of the body, where it arches around and then descends into the abdomen. Along the way, several arteries branch off to carry blood to the heart, head, torso, kidneys and other organs. The aortic root features a valve which consists of three passive leaflets that are housed within three corresponding pouches in the aortic wall. These pouches are called aortic sinuses. The aortic valve, aortic sinuses and the curvature of the aortic arch interact with the flow emanating from the left ventricle to exhibit an array of interesting fluid mechanics that has fascinated scientists for centuries. This is famously evidenced by the illustrations of da Vinci in 1513 showing regions of recirculation around the aortic valve (Boon 2009).

Of particular interest has been the fluid mechanics of the opening and closure of the aortic valve and its relation to the geometry of the aortic sinus (Yoganathan, He & Jones 2004; Sotiropoulos et al. 2016). The three passive leaflets of the aortic valve respond to the haemodynamics within the upstream left ventricle and the downstream aorta to provide minimal flow resistance while open during ventricular contraction at systole, and to minimize backflow during ventricular filling at diastole. The opening is driven by the pressure difference between the contracting left ventricle and the aorta, i.e. the leaflets open when the pressure within the left ventricle exceeds the aortic pressure. Interestingly, the closure of the aortic valve begins before the end of this ejection phase – with an estimated 75 % of the valve closure occurring while blood continues to flow from the left ventricle to the aorta (Bellhouse & Talbot 1969). This ensures that the valve is completely closed with minimal reversed flow at the end of systole.

The hydrodynamic mechanism driving this has intrigued numerous researchers (Bellhouse & Bellhouse 1968; Peskin & Wolfe 1978), starting with the first modern study of this dynamics by Henderson & Johnson (1912). The seminal works of Bellhouse and colleagues highlighted the importance of the sinus vortex, which results from flow over the top of the aortic valves into the aortic sinuses, in initiating aortic valve closure using laboratory experiments accompanied by elegant theoretical calculations of the pressure induced by these vortices on the leaflet surfaces (Bellhouse 1969; Bellhouse & Talbot 1969). They also hypothesized the role of these vortices in maintaining the valve open position through a balance between sinus flow and the flow through the valve. While the physiological significance of the sinus vortex is still a topic of investigation (Moore & Dasi 2014, 2015), it is remarkable that the contemporary results exhibit a striking parallel with the illustrations of Da Vinci over 400 years earlier (Gharib et al. 2002).

The interaction of the flow with the valve leaflets and the aortic sinus is also important to the pathophysiology of aortic valve disease. Pathological conditions relating to aortic root diseases and congenital defects, such as bicuspid aortic valves where the valve forms with two leaflets instead of three, lead to abnormal flow patterns. In these cases, studies have reported asymmetric jets emanating from the valve opening, which impinge on the aortic wall to cause unphysiological shear stresses (Girdauskas et al. 2011; Barker et al. 2012; Verma & Siu 2014; Saikrishnan, Mirabella & Yoganathan 2015; Kaiser et al. 2022). This is shown in figure 2(b). The replacement of diseased aortic valves with mechanical and prosthetic valves also alters the haemodynamics. This is most commonly done using mechanical valves that consist of flat pivoting metal plates to mimic leaflets. The interaction of flow with the rigid leaflets gives rise to highly complex unphysiological flow within the ascending aorta, such as interacting vortex rings and shear layers emanating from the tips of the flat-plate leaflets (Dasi et al. 2009; Sotiropoulos et al. 2016). Moreover, studies have reported increased shear regions within the housing of prosthetic valves that can damage blood cells (Dasi et al. 2009). Prosthetic valves with flexible leaflets also exhibit the unphysiological flow patterns described above, albeit to a lesser degree depending on the properties of the leaflet and housing. In addition, geometrical alterations to the aortic root through surgical procedures to treat aortic valve diseases have been shown to result in secondary flows that exhibit turbulent transition and associated stress fluctuations (Stein & Sabbah 1976; Peacock 1990; Moore & Dasi 2014). The interaction of these unphysiological flow features with the valve leaflets is believed to exacerbate the stiffening or calcification of aortic valve leaflets as well as complications at the aortic wall. While the exact mechanism leading to leaflet calcification is not clear, there are strong correlations between flow-induced mechanical stresses, such as oscillatory shear stresses from flow recirculation, and the inflammation of endothelial cells on the leaflets (Ge & Sotiropoulos 2010; Gould et al. 2013). The flow through asymmetrically calcified or stenosed valves has been shown to develop into a jet that impinges on the aortic wall, where the curvature of the aortic arch introduced helical Dean vortex flows, intermittent turbulent stresses and audible sounds called murmurs (Zhu, Seo & Mittal 2018, 2019; Bailoor et al. 2021).

The rich haemodynamics associated with aortic flows has led to efforts to non-invasively quantify these flow in vivo using 4-D flow MRI in search of haemodynamic biomarkers of disease (Von Knobelsdorff-Brenkenhoff et al. 2014; Burris & Hope 2015; Garcia et al. 2019; Zimmermann et al. 2023). In this context, 4-D flow MRI is a promising method to non-invasively evaluate haemodynamic patterns that can aid in the diagnosis of valvular diseases as well as in the evaluation of prosthetics. However, the challenge with 4-D flow MRI is accurately resolving the flow dynamics and associated wall interactions, which has led to the emergence of computational methods to simulate these flows. The first method aimed at simulating the fluid dynamics of cardiac valves was introduced by Peskin and colleagues just over 50 years ago (Peskin 1972; Peskin & McQueen 1980). This technique, which is now referred to as the immersed boundary method, continues to be the mainstay of computational models for cardiac valves (Mittal & Iaccarino 2005; Griffith & Patankar 2020). Since its introduction, the immersed boundary method has been coupled with finite-element formulations of the mechanics and material properties of valve leaflets (Boffi et al. 2008; Griffith & Luo 2017; Kaiser et al. 2021). Sharp-interface immersed boundary methods have been developed, which can resolve the surface of the leaflets more accurately and therefore provide improved surface force estimates (Gilmanov & Sotiropoulos 2016; Seo et al. 2020; Bailoor et al. 2021). In the finite-element context, these immersed boundary fluid–structure interaction methods are broadly referred to as immersed finite-element or immersogeometric methods (Hsu et al. 2014, 2015; Kamensky et al. 2015).

These methods have been extensively used to study the haemodynamics of native as well as prosthetic valves in healthy and diseased states (see Votta et al. (2013) and Sun, Mao & Griffith (2019) for comprehensive reviews of applications). In fact, one such simulation-based product has recently received United States Food and Drug Administration approval to support clinical decision making for transcatheter aortic valve replacements (www.dasisim.com). However, an ongoing challenge for more widespread clinical acceptance of these simulation tools is the lack of validation with in vivo flow fields. The majority of validation studies have focused on the haemodynamics of mechanical heart valves measured in laboratory experiments (Dasi et al. 2007; Ge et al. 2008; Nobili et al. 2008; Guivier-Curien, Deplano & Bertrand 2009; Jun et al. 2014). In comparison, fewer studies have validated simulation models of heart valves with flexible leaflets (Wu et al. 2016; Sigüenza et al. 2018; Lee et al. 2020), especially with patient-specific flow fields (Emendi et al. 2021; Kaiser et al. 2023).

While a major focus of studies on aortic haemodynamics has been on the function of the aortic valve, the influence of aortic flow on the wall of the aorta has important health consequences. In particular, shear stresses induced on the aortic wall due to the jet emanating from the aortic valve as well as the resultant turbulent flow within the aorta can damage the epithelial cell lining of the aortic wall. The peak Reynolds number in the aorta can be as high as ${\sim }4000$, which is significantly higher than the transitional Reynolds number in pipes (Ku 1997). However, interesting work has shown that the nature of cardiac pulsatility, including the specific waveform of aortic flow, is close to the optimal behaviour required to combine sufficient flow rates with low wall shear stresses. This is due to the combination of a short ejection phase followed be a prolonged diastolic phase that inhibits turbulence (Scarselli et al. 2023).

Studies have also reported that 60 %–80 % of patients with bicuspid aortic valves develop dilation of the aortic wall (Verma & Siu 2014). While the exact mechanism leading to aortic dilation is unknown, altered haemodynamics associated with the impinging jet arising from bicuspid valves is thought to play a role (Girdauskas et al. 2011; Barker et al. 2012; Rodríguez-Palomares et al. 2018; Dux-Santoy et al. 2020). Furthermore, aortic dilation increases the risk of aortic dissection, which is a tear in the inner wall of the aorta that introduces a parallel flow channel within the aortic wall, referred to as the false lumen (in contrast to the native ‘true’ lumen), and poses the risk of rupture of the aortic wall (Verma & Siu 2014). Although the clinical management of aortic dissections is primarily based on anatomical imaging (Spinelli et al. 2018), studies have shown that haemodynamics plays a significant role in the prognosis (Bäumler et al. 2020). For example, while the development of clots to occlude the false lumen has been shown to improve patient outcomes, the structure of these clots can cause elevated pressure within the false lumen and increase rupture risk (Tsai et al. 2008; Berguer et al. 2015).

Another pathological condition that poses the risk of aortic wall rupture is the occurrence of dilations in the abdominal aorta, referred to as abdominal aortic aneurysms (Sakalihasan, Limet & Defawe 2005). For patients with aortic dissections and aortic aneurysms, the stratification of rupture risk depends on thrombus formation due to flow stagnation as well as stress distributions induced on the aortic wall through its interaction with the haemodynamics (Bluestein et al. 1997; Vorp et al. 2001; Francis 2002; Sho et al. 2004). A challenge in this context is the quantification of aortic wall mechanics, which is often heterogeneous and difficult to measure in vivo. Conversely, there are also haemodynamic risks associated with the constriction of the aortic wall. This narrowing of the descending aorta, which is referred to as aortic coarctation and accounts for 6 %–8 % of congenital heart defects, can lead to aortic rupture, hypertension and stroke (Singh et al. 2015). These are therefore active avenues for the future application of patient-specific computational models to aid in clinical decision making (for example, Bluestein et al. 1997; Salsac, Sparks & Lasheras 2004; Les et al. 2010; Alimohammadi et al. 2015; Chen et al. 2016; Qiao et al. 2019; Bäumler et al. 2020, among many others). Mourato et al. (2022) provide a comprehensive review of fluid–structure interaction computational models pertinent to aortic wall diseases.

4. Coronary haemodynamics

While the pumping action of the heart ensures the delivery of blood to the far reaches of our circulatory system, the coronary arteries bear the unenviable responsibility of delivering nutrition and oxygen-laden blood to the heart muscle itself. The importance of the coronary circulation is underscored by the fact that the heart has the highest oxygen consumption per gram of any organ; it extracts 70 %–80 % of the oxygen delivered, compared with 30 %–40 % utilized by skeletal muscles (Goodwill et al. 2017). It is no wonder that coronary artery disease, which inhibits the transport of blood to the myocardium due to the narrowing or occlusion of coronary arteries, is the leading cause of death worldwide (Ralapanawa & Sivakanesan 2021; Tsao et al. 2022).

The left and right coronary arteries arise on either side of the aortic sinus and carry 3 %–5 % of the volume of blood pumped by the heart into the myocardial tissue (Guyton & Hall 2006). They have diameters of ${\approx }3$ mm at their origin on the aortic sinus, and proceed through a series of bifurcations that leads to small arteries and arterioles with diameters ${\approx }100\ \mathrm {\mu }$m, and eventually capillaries with diameters ${\approx }10\ \mathrm {\mu }$m (Lee & Smith 2012). Although the peak Reynolds number in coronary arteries is typically ${<}500$, coronary flow can exhibit higher Reynolds number jets downstream of constrictions that occur due to disease. This has consequences for the progression of coronary artery disease as well as the stability of calcified plaques that cause these constrictions (Falk, Shah & Fuster 1995; Samady et al. 2011). A peculiar aspect of the coronary circulation is that, while maximum flow through most vascular networks in the body occurs during systole, a majority of the flow through coronary arteries occurs during diastole. This is due to the so-called ‘cross-talk’ between the contracting myocardial tissue and the coronary circulation from large vessels to arterioles and capillaries embedded within the myocardial tissue (Westerhof et al. 2006). The disparate length scales involved in the coronary circulation and its interaction with the beating (and moving) heart, combined with the importance of coronary flows to cardiovascular function, present a host of challenges in understanding, imaging, as well as modelling coronary artery blood flows.

Most studies on coronary artery haemodynamics have focused on the flow through the large arteries that can be clinically imaged, either invasively using coronary angiography or non-invasively using computed tomography. Efforts to measure and model coronary artery flow are motivated by several randomized clinical trials that have firmly established the importance of considering haemodynamics in clinical care. Coronary artery disease causes the narrowing and stiffening of coronary artery walls, which limits the delivery of blood to the myocardium. This has been observed to occur in regions of the vasculature that experience low and oscillating wall shear stress in laboratory experiments as well as animal models (Giddens et al. 1993). Coronary artery disease is usually treated either by inserting stents within affected arteries to relieve the narrowing, or in cases with more severe occlusion, surgically bypassing the location of narrowing in arteries by using grafts (Yusuf et al. 1994; Windecker et al. 2014). In both cases, the long-term outcomes for patients receiving these treatments has been shown to significantly improve when treatments are guided by haemodynamic measurements of pressure loss within the affected artery, rather than based on anatomical measures of the narrowing in the artery (Meijboom et al. 2008; Tonino et al. 2009; Pijls et al. 2010). There is also evidence that short as well as long-term treatment outcomes are dependent on local blood flow patterns at the site of stenting or grafting, such as flow reversal and competitive flow between the native diseased artery and bypass graft (Pagni et al. 1997; Sabik et al. 2003; Tokuda et al. 2008), as well as regions of stagnation or abnormal wall shear stress (Malek, Alper & Izumo 1999; Chatzizisis et al. 2008; Humphrey 2008).

While the major focus of studies on coronary artery blood flows has been on the flows within the coronary arteries themselves, the location of coronary arteries and their proximity to the aortic valve also raises interesting questions. Owing to their origin within the aortic sinus, researchers have long hypothesized about the influence of recirculating flow within the aortic sinus (due to interaction with the aortic valve) on coronary circulation, and vice versa. Bellhouse, Bellhouse & Reid (1968) used a laboratory model to relate flow stagnation on the ridge of the aortic sinus to high pressure that would be beneficial to driving coronary flow. They showed a significant drop in coronary flow when the coronary arteries originate outside the aortic sinus. There is also a difference in the fluid dynamics between the two sinus lobes that contain the origin of coronary arteries, or coronary sinuses, and the lobe that does not. Moore & Dasi (2015) showed that the presence of coronary arteries originating from the coronary sinuses pulls the sinus vortex closer to the walls. This leads to higher flow velocities, shear stresses and washout in the coronary sinuses compared with the non-coronary sinus. In addition, the aortic leaflets corresponding to coronary sinuses open farther into the sinus and exhibit better mechanics. These haemodynamic differences have physiological significance because low recirculation velocities lead to less washout in the non-coronary sinus, and low and oscillatory shear stresses have been linked to calcification (Chandra, Rajamannan & Sucosky 2012). Interestingly, these observations also present a possible mechanism for the observation that the non-coronary sinus is often the first to calcify – eventually leading to aortic valve calcification (Freeman & Otto 2005).

The importance of fluid mechanics in these critical clinical questions has spurred development of patient-specific computational models of coronary artery blood flows for personalized risk assessment and treatment planning. Current state-of-the-art methods use non-invasive clinical imaging to build patient-specific coronary artery anatomical models (Updegrove et al. 2017). The flow within these anatomies is simulated using either high-fidelity 3-D flow simulation methods, or reduced-order methods that evaluate the flow field either only along the centreline of the arteries or at specifically chosen points in the anatomy, such as bifurcations, by assuming Poiseuille flow throughout (see § 7). Owing to the complex geometries of vascular trees, 3-D simulations are usually performed using stabilized finite-element methods with unstructured meshes (Taylor et al. 1998; Marsden & Esmaily-Moghadam 2015). A significant challenge is that these image-based models are usually restricted to the large vessels in the coronary tree that are visible from volumetric clinical imaging. Therefore, the effect of the smaller downstream arteries and their interaction with the contracting myocardium needs to be modelled appropriately. This is usually accomplished by coupling the flow within the large arteries included in these anatomical models with lumped-parameter boundary conditions at their outlets (Kim et al. 2010). These boundary conditions are designed to account for the resistance as well as dilation of downstream blood vessels not included in the model. They also recapitulate the key features of coronary physiology in which coronary flows peak at diastole due to the ‘intramyocardial pressure’ experienced by the coronary circulation via its interaction with myocardial contractions. In addition, the compliance of the coronary artery walls as well as their movement with the contracting myocardium pose additional challenges in the modelling of coronary haemodynamics. Although the former is commonly included in models using simplified as well as high-fidelity fluid–structure interaction techniques, such as the coupled momentum method (Figueroa et al. 2006) and arbitrary Lagrangian–Eulerian (Seo et al. 2020), studies have shown that the effect of wall compliance can be small in certain situations (Torii et al. 2009; Eslami et al. 2020). However, although the movement of coronary arteries due to myocardial contraction can have significant effects on the haemodynamics, this is rarely included in simulations (Santamarina et al. 1998).

As expected, the modelling techniques and boundary conditions used to simulate coronary haemodynamics include several parameters relating to the coupling between the microcirculation and large vessels that cannot be measured on a patient-specific basis. Importantly, they determine how the flow entering the coronary artery tree is split amongst all the vessels in the tree. While there are physiological and empirical measures relating coronary artery diameters to the flow splits amongst the arteries (Murray 1926; Zhou, Kassab & Molloi 1999), these are not patient specific and do not account for vascular disease. It is not surprising that parameters affecting the amount of flow through each coronary artery have a measurable impact on our ability to predict pressure losses resulting from narrowing or aneurysms in these arteries (Wellnhofer et al. 2010; Sankaran et al. 2016; Vardhan et al. 2019). Advances in imaging methods, combined with patient-specific modelling, are now aiming to directly measure myocardial perfusion to address this need (Nieman & Balla 2020; Menon et al. 2024).

In spite of the challenges outlined above, computational modelling of coronary artery haemodynamics is a powerful tool for non-invasive clinical risk assessment and treatment planning. The computation of the pressure loss resulting from diseased coronary arteries using patient-specific computational models has successfully predicted disease severity and reduced unnecessary invasive procedures in clinical trials (Koo et al. 2011; Min et al. 2012; Taylor et al. 2013; Mathew, Gottbrecht & Salerno 2018). This has led to the approval of the first simulation-based risk assessment tool, $FFR_{{CT}}$, for clinical use by the United States Food and Drug Administration (www.heartflow.com). In addition, computational modelling also provides insights not accessible in clinical imaging. For example, the multi-directionality and topological skeleton of wall shear stress, i.e. regions of expansion and constriction in the wall shear stress field, have emerged as possible predictor of atherosclerosis (Arzani & Shadden 2018; Kok et al. 2019; Hoogendoorn et al. 2021; Mazzi et al. 2021). In addition, computational modelling has highlighted the role of coronary flow helicity in promoting healthy wall shear stress (Morbiducci et al. 2007; De Nisco et al. 2019), as well as other insights into wall shear stress patterns that are related to the progression of coronary artery disease and coronary artery bypass graft failure (Nordgaard et al. 2010; Ramachandra, Kahn & Marsden 2016; Khan et al. 2020; Candreva et al. 2022).

The importance of haemodynamics in risk assessment and treatment has also been demonstrated in other afflictions of the coronary arteries. For example, Kawasaki disease is a paediatric inflammatory condition that leads to coronary artery aneurysms in some untreated patients (McCrindle et al. 2017). A risk associated with these aneurysms is clot formation as a result of flow stagnation and recirculation. While current clinical guidelines recommend treatment based purely on aneurysm size, recent studies have demonstrated the use of computational modelling to reveal blood flow patterns within aneurysms that might provide a better indicator of patients’ risk (Sengupta et al. 2012, 2014; Grande Gutierrez et al. 2019; Menon et al. 2023). The use of computational modelling for non-invasive clinical risk assessment has also been explored in patients affected by anomalous aortic origin of coronary arteries, a congenital defect that causes the coronaries to arise from outside their normal aortic sinus location (Jiang et al. 2022). In addition to clinical risk stratification, there have also been efforts to employ patient-specific computational models of coronary haemodynamics for virtual planning of treatments and surgical procedures (Balossino et al. 2008; Sankaran et al. 2012; Morlacchi & Migliavacca 2013; Ballarin et al. 2017; Seo et al. 2021). Due to the importance of haemodynamics in the clinical diagnoses, treatment planning and disease progression of such a central part of our circulatory system, these efforts to develop personalized and predictive models of coronary artery haemodynamics, as well as developments in imaging methods to measure coronary flows, have the potential to make a significant impact on cardiovascular health.

5. Cerebral haemodynamics

Following the branching off of coronary arteries at the aortic sinuses, the next set of arteries to arise from the aorta is located downstream at the aortic arch. Of these, the most important are the left and right common carotid arteries, the latter bifurcating from the brachiocephalic artery. The carotid arteries are responsible for supplying the face, neck and, crucially, the brain. The left and right internal carotid arteries, which are branches of the left and right common carotid arteries, travel up the neck into the base of the skull, after which they bifurcate into a network of cerebral arteries that perfuse the brain. In humans, the brain accounts for only 2 %–3 % of body mass, but receives approximately 15 %–20 % of cardiac output and consumes 20 % of the available oxygen in the body under normal conditions (Claassen et al. 2021). The brain's high metabolic rate and the limited ability for anaerobic metabolism in neurons underscore the importance of maintaining healthy cerebral blood flow. The disruption of cerebral blood flow therefore has severe consequences that can lead to life-threatening strokes, permanent brain damage, as well as other maladies related to impaired removal of metabolic wastes from the brain (Tarasoff-Conway et al. 2015).

One of the most widely investigated topics relating to cerebral circulation, is the growth, rupture and treatment of cerebral aneurysms. These aneurysms typically form in the Circle of Willis, which is a traffic circle-like junction of several arteries at the base of the brain. Research has shown that the most common locations of aneurysm development coincide with regions of elevated flow and wall shear stress, usually occurring due to pathological vessel narrowing or at bifurcations between arteries (Sforza et al. 2009). Haemodynamics also plays a role in the growth and rupture of cerebral aneurysms due to the interaction between flow-induced stresses and cellular function on blood vessel walls. In fact, computational modelling has revealed that a significant proportion of cerebral aneurysms exhibit flow instabilities that lead to high-frequency oscillations even at modest Reynolds numbers of ${\sim }300$ (Khan et al. 2021). Interestingly, both low and high wall shear stress have been implicated in unfavourable aneurysm progression, with each pathological state being associated with a different vascular response (Shojima et al. 2004; Jou et al. 2008; Cebral et al. 2011; Xiang et al. 2011). While there is general agreement on the mechanism of aneurysm initiation – triggered by cell damage caused by high flow and wall shear stress states that have been observed in animal models (Stehbens 1989) – there is disagreement on the mechanisms associating abnormal wall shear stress with aneurysm progression and rupture (Xiang et al. 2014). Nevertheless, these observed links between haemodynamics and aneurysm progression, and the inability of in vivo imaging to measure critical haemodynamic stimuli, have accelerated the use of computational models to identify haemodynamic predictors of aneurysm progression (Chung & Cebral 2015). Studies have introduced several haemodynamic metrics to quantify aneurysm rupture risk, such as high oscillatory shear index and particle residence time (Gonzalez et al. 1992; Xiang et al. 2011). The risk of aneurysm rupture has also been correlated with the presence of complex vorticity patterns (Xiang et al. 2011) as well as the size, strength and location of the jet that emanates from the neck of the aneurysm and impinges on the aneurysm wall (Cebral et al. 2005).

Haemodynamics has also played a central role in the treatment of cerebral aneurysms to prevent rupture. Increasingly, clinical interventions have moved away from surgery, towards the implantation of endovascular devices that promote the formation of stabilizing clots within aneurysms (Brisman, Song & Newell 2006; Brinjikji et al. 2011). The development of such clots, either forming spontaneously or promoted by implanted devices, has been shown to reduce the incidence of rupture (Etminan & Rinkel 2016). The most commonly implanted devices are coils and flow diverters, both of which impede flow into the aneurysm and lead to flow stagnation, low wall shear stress and, consequently, clot formation (Ngoepe et al. 2018). Although the biochemistry of clot formation is currently still under investigation, computational haemodynamics models have proven to be remarkably effective at predicting the potential for clot formation and its location based on flow stagnation and residence time (Sadasivan et al. 2002; Rayz et al. 2010; De Sousa et al. 2016), as well as simulating the process of clot formation by coupling fluid and chemical transport with models of thrombogenesis (Bedekar et al. 2005; Ngoepe & Ventikos 2016).

However, broader clinical acceptance of computational models for patient-specific risk stratification and thrombosis prediction has been slow. One reason for this is the lack of clear mechanobiology that relates the observed haemodynamic patterns to pathophysiology, which makes conflicting findings regarding high and low wall shear stress difficult to interpret (Xiang et al. 2014). In the prediction of thrombosis, an interesting challenge is the wide range of time scales involved – the fluid transport occurs at the scale of seconds, clot formation occurs in minutes, while clot maturation and stabilization occurs over months (Szikora, Turányi & Marosfoi 2015). Moreover, these interacting time scales influence the properties of blood itself, which changes during the process of clotting. Cerebral haemodynamics models also contend with the familiar challenges of patient specificity in terms of model parameters and boundary conditions (Levitt et al. 2014; Steinman & Pereira 2019). An added difficulty is that the biochemical coagulation profile of each person is different, and this needs to be accounted for in patient-specific simulations of clot formation (Ngoepe et al. 2018). Finally, the lack of many large patient cohort studies has hampered their translation to clinical use (Sforza et al. 2009).

A unique aspect of the cerebral circulation that also benefits from patient-specific haemodynamic measurements and simulations is the presence of collateral arteries due to the anatomy of the Circle of Willis. Collaterals provide alternate flow paths that allow cerebral circulation to be robust to occlusions in specific parts of the vascular anatomy. However, they also complicate the clinical management of cerebral stenoses, because anatomical measurements of vessel narrowing do not correlate well with their effect on cerebral blood flow due to the presence of alternate flow paths via collaterals (Powers 1991). Furthermore, cerebral circulation exhibits a high degree of autoregulation, i.e. the blood vessels adapt to changes in upstream pressure and flow to preserve homeostatic cerebral blood flow volumes (Claassen et al. 2021). The study of the haemodynamics associated with cerebral and carotid stenosis and their treatments at the patient-specific level has therefore emerged as another area for the application of personalized imaging and computational simulations based on imaging of patient-specific cerebral artery anatomies (Birchall et al. 2006; Liang et al. 2011).

Of particular interest is the haemodynamics of the growth and rupture of atherosclerotic plaque in carotid stenosis due to the risk of dislodged plaque travelling downstream, occluding other parts of the circulation and causing stroke (Carr et al. 2013; Mukherjee et al. 2016). The location of these plaques has been shown to most often correlate with regions of low wall shear stress and flow disruption (Zarins et al. 1983; Ku et al. 1985). The development of carotid re-stenosis following treatment has also been correlated with wall shear stress patterns that exhibit high temporal variation and strong expansion/contraction through the analysis of the wall shear stress topological skeleton (Morbiducci et al. 2020). In addition, the narrowing of these arteries due to the presence of plaque leads to accelerated and highly transient downstream flows that exhibit elevated shear stress and stress gradients as well as turbulent fluctuations, all of which are thought to cause cell damage and plaque rupture (Davies et al. 1986; Gertz & Roberts 1990; Loree et al. 1991; DePaola et al. 1992; Lee et al. 2008).

Finally, a fascinating and relatively unexplored haemodynamic phenomenon associated with cerebral circulation is the transport and clearance of brain interstitial fluid, which occupies interstitial regions of the brain surrounding blood vessels and neurons and carries metabolic wastes and toxins (Shetty & Zanirati 2020). In particular, the impaired clearance of amyloid betas is a hallmark of Alzheimer's disease (Tarasoff-Conway et al. 2015). While the bulk flow of interstitial fluid has recently been identified as one of several simultaneous mechanisms for the clearance of such wastes, the mechanisms driving this flow are currently unknown. Due to the fact that interstitial fluid envelopes cerebral vasculature, studies have suggested that arterial wave propagation as a result of pulsatile blood flow are a major driver of interstitial fluid transport (Iliff et al. 2013; Ladrón-de Guevara et al. 2022). However, this has been debated due to the strength of these pulsations, with larger amplitude and lower-frequency motions caused by vasomotion proposed as an alternate mechanism (Diem et al. 2017; van Veluw et al. 2020). In spite of these disagreements, the role of circulatory phenomena in the clearance of toxins, as well as the aforementioned haemodynamic phenomena associated with aneurysms and strokes, have important consequences for brain health and are promising avenues for ongoing and future research.

6. Peripheral vascular flows

As we travel beyond the main arteries discussed so far, the circulatory system begins to branch off into the peripheral vascular system, which consists of the blood vessels outside of the heart and brain. The peripheral vascular system plays a critical role in transporting blood to and from the body's tissues, and it can be affected by various conditions. The most notable is peripheral artery disease, which affects over 8 million individuals aged 40 and older in the United States (Aday & Matsushita 2021). This disease is a slow and progressive disorder characterized by narrowing, blockage or spasms in the arteries, commonly caused by atherosclerosis – a condition marked by plaque accumulation on the arterial walls. This buildup of plaque narrows the arteries, impeding blood flow and often leading to complete occlusion. Since peripheral artery disease is a subtype of atherosclerotic disease, it is not surprising that people with peripheral artery disease have a higher risk of other atherosclerotic diseases such as coronary artery disease and stroke (Aday & Matsushita 2021).

Similar to the cerebral flow, wall shear stress has been emphasized to play an important role in peripheral artery disease: low wall shear stress may promote the build-up of plaque while high wall shear stress could increase the vulnerability of plaques (Samady et al. 2011; Casa, Deaton & Ku 2015). Simulations of local flow through thrombosis have quantified abnormal mean and oscillatory wall shear stress, which, when combined, may serve as risk factors for rupture (Xu et al. 2016; Ferrarini et al. 2021; Wang, Serracino-Inglott & Feng 2021). Additionally, simulations have highlighted the importance of patient specificity of peripheral vascular anatomies. Wood et al. (2006) investigated the differences in flow patterns, curvature and tortuosity in the superficial femoral artery between young men and women. Their MRI-based computational modelling revealed that men exhibited greater tortuosity and curvature in the superficial femoral artery, along with lower mean wall shear stress, compared with women, which might contribute to the understanding of the two- to five-fold higher incidence of peripheral vascular disease in men.

Finally, a unique challenge of modelling peripheral artery blood flow is that lower limb movement significantly affects the haemodynamics in the lower peripheral vessel. To study this, Colombo et al. (2020) mimicked hip rotation, knee flexion and complete movement of walking with idealized femoropopliteal arteries and moving-boundary models. Interestingly, they showed that leg movements induce higher time-averaged wall shear stress, with flow rate amplitude and walking period being the most influential parameters that affect wall shear stress.

7. Microvascular and capillary flows

From centimetre-scale blood vessels arising from the heart, to micron-scale capillaries perfusing the various organs in the body, the interacting scales in the circulatory system span roughly six orders of magnitude. For example, distal capillary beds, consisting of millions of micron-sized blood vessels, respond to biomechanical stimuli and metabolic demand to change the flow within upstream millimetre-sized blood vessels, which can be as much as four-fold in the coronary circulation (Wilson et al. 1990), via autoregulatory mechanisms (Goodwill et al. 2017). The flow of blood, which enters the aorta and its branches exhibiting Newtonian-like flow in the large vessels, ceases to behave as a continuum fluid as it traverses the $\mathrm {\mu }$m-sized arterioles and capillaries. This leads to a nonlinear as well as non-monotonic relationship between fluid viscosity and vessel diameter (Secomb 2017). The study of such multi-scale interaction necessitates multi-scale and multi-physics modelling approaches to tackle the prohibitive computational cost associated with simulating blood flow in these disparate scales.

The earliest physiologically motivated models of flow in the circulatory system considered the system as a so-called ‘lumped-parameter network’, drawing analogies between electric and hydraulic circuits (Westerhof et al. 1969). These highly simplified models, also known as ‘zero-dimensional models’, have basic building blocks comprising resistors (R) to model viscous resistance within blood vessels derived from Poiseuille flow, capacitors (C) to model vessel wall distensibility and inductors (L) to model fluid inertia. The simplified equations that encapsulate these physical phenomena to relate flow ($Q$) and pressure ($P$), which are analogous to electric current and voltage, respectively, of blood with density $\rho $ and viscosity $\mu $ within a blood vessel with length $l$, radius $r$, wall thickness $h$ and Young's modulus $E$ are given by

(7.1) \begin{equation} \left.\begin{array}{c} \text{Resistance: } {\rm \Delta} P = RQ, \quad \text{where } R = \dfrac{8\mu l}{{\rm \pi} r^4},\\ \text{Capacitance: } {\rm \Delta} \dot{P} = \dfrac{Q}{C}, \quad \text{where } C = \dfrac{3{\rm \pi} l r^3}{2Eh},\\ \text{Inductance: } {\rm \Delta} P = L\dot{Q}, \quad \text{where } L = \dfrac{\rho l}{{\rm \pi} r^2}. \end{array}\right\} \end{equation}

Lumped-parameter models use suitable arrangements of the above resistors, capacitors and inductors to provide macroscopic descriptions of blood flow rates and pressures at specific locations in hydraulic circuits via a set of ordinary differential equations satisfying global mass conservation but only partial momentum conservation (Pfaller et al. 2022). Figure 3 shows two examples of such zero-dimensional models. Although the simplicity of these models does not conventionally allow accurate modelling of complex flow topologies, such as vessel bifurcations and curvature, several studies have proposed improvements and corrections to account for flow nonlinearities while preserving the simplicity of these models (Mynard & Valen-Sendstad 2015; Chnafa et al. 2017; Mirramezani & Shadden 2020). The minimal computational expense of zero-dimensional models has led to their wide use in low-order modelling of multi-scale circulatory phenomena by effectively ‘lumping’ together the resistance and capacitance of entire micro-vascular/capillary beds into single electric circuit-like components (Vignon-Clementel et al. 2006; Kim et al. 2010).

Figure 3. Examples of multiscale computational models. (a) A model of Fontan surgery that allows the quantification of local as well as whole-body haemodynamics (Schwarz et al. 2021). (b) A coronary flow model coupling 3-D flow in the aorta and large coronary vessels, 0-D flow in smaller myocardial vessels (shown in blue) and systemic circulation and Darcy flow in the myocardial microvasculature (Menon et al. 2024).

Another class of simplified vascular flow models that has enabled the modelling in large, multi-scale circulatory networks at reasonable computational expense are so-called ‘one-dimensional’ flow models (Formaggia et al. 1999, 2001). The description of flow within blood vessels is obtained by averaging the 3-D Navier–Stokes equations across the vessel cross-section and projecting this on the 1-D curve tracing the centre of blood vessels ($z$). This yields partial differential equations describing the flow rate, $Q(t,z)$, and pressure, $P(t,z)$, within vessels of cross-sectional area $S(t,z)$

(7.2)\begin{equation} \left.\begin{array}{c} \displaystyle \dfrac{\partial Q}{\partial t} + \dfrac{4}{3}\dfrac{\partial }{\partial z}\left(\dfrac{Q^2}{S}\right) + \dfrac{S}{\rho}\dfrac{\partial p}{\partial z} = Sf - \dfrac{8\mu}{\rho}\dfrac{Q}{S} + \dfrac{\mu}{\rho} \dfrac{\partial^2 Q}{\partial z^2},\\ \displaystyle \dfrac{\partial S}{\partial t} + \dfrac{\partial Q}{\partial z} = 0. \end{array}\right\} \end{equation}

The above equations are closed using constitutive models for the wall deformation that relate $P(t,z)$ to $S(t,z)$. These models provide higher-fidelity descriptions of the fluid dynamics than lumped-parameter models, particularly in vasculature with more complicated topology, while still remaining significantly less computationally expensive than full 3-D flow simulations.

These 0-D and 1-D descriptions of the circulation have been widely used in multi-scale flow simulations of extensive vascular networks (Westerhof et al. 1969; Avolio 1980; Olufsen 1999). They are most often used in conjunction with high-fidelity 3-D models of blood flow in a multi-scale manner, where 3-D blood flow simulations are performed in limited vascular segments of interest and 0-D or 1-D models approximate downstream/upstream circulatory behaviour that is not included in 3-D blood flow simulations (Migliavacca et al. 2006; Vignon-Clementel et al. 2006; Kim et al. 2010; Sankaran et al. 2016; Grande Gutierrez et al. 2019; Grande Gutiérrez, Sinno & Diamond 2022; Menon et al. 2024). This is shown in figure 3. Such multi-scale models have allowed researchers to investigate the relationship between widely disparate circulatory scales, such as the effect of exercise on changes in cardiac function and pulmonary circulation (Kung et al. 2014), haemodynamic changes in distal/proximal circulation to compensate for local abnormalities within specific vascular regions (Balossino et al. 2009; Schwarz et al. 2021) and the effects of local surgical interventions (e.g. shunt sizing) on global physiologic parameters such as oxygen delivery and cardiac workload in patients with congenital heart disease (Migliavacca et al. 2006; Moghadam et al. 2012; Esmaily-Moghadam, Hsia & Marsden 2015).

Although their minimal computational expense has allowed the study of physiological interactions amongst different scales, the primary challenge associated with such multi-scale models is the relationship between their parametrization and the associated physiology they aim to capture. For example, in models consisting of 3-D descriptions of vascular haemodynamics coupled with reduced-order representations of the beating heart, it is not immediately obvious how the parameters of the heart model, usually consisting of resistors, time-varying capacitors and valves, relate to patients’ physiological cardiac volume, ejection fraction and other clinically relevant metrics. The need to associate these low-order flow models with clinical measurements for personalized modelling has led to development of parameter estimation and automated tuning schemes, which ideally should account for uncertainties (Tran et al. 2017; Fleeter et al. 2020; Seo et al. 2020; Menon et al. 2024).

While the above low-order modelling strategies can be used to provide macroscopic descriptions of flow and pressure through lumped-parameter representations of microvasculature, they are computationally intractable when modelling flow that is resolved at the scale of each blood vessel in dense capillary networks within tissues. This is due to the complexity and computational expense of modelling millions of interconnected vessels, as well as the network-like structure of capillary beds where simplified Poiseuille flow, without accounting for viscous effects at junctions between blood vessels, may not be valid. A spatially averaged/coarse-grained description of flow is therefore required in such microvascular networks. In fact, dimensional analysis as well as microfluidic measurements have shown a linear relationship between the pressure gradient and spatially averaged flow within microvasculature (Fung 1997; Stauber et al. 2017) – akin to Darcy's law governing flow through porous materials. Recent work has therefore led to the development of computational models for microvascular/capillary tissue perfusion that represent dense capillary networks and interstitial regions of tissues as porous materials.

These models have been applied to study capillary flows, oxygen or nutrient exchange, as well as the transport of therapeutics in various organs such as in alveoli within the lungs (Erbertseder et al. 2012; Zurita & Hurtado 2022), within tumours (Secomb et al. 2004; Stylianopoulos & Jain 2013; d'Esposito et al. 2018), in the brain (Sweeney, Walker-Samuel & Shipley 2018; Kim et al. 2023) and in the heart (Chapelle et al. 2010; Cookson et al. 2012; Michler et al. 2013; Di Gregorio et al. 2021, 2022; Papamanolis et al. 2021; Kim et al. 2023). While the development of tissue perfusion models is relatively nascent compared with the 3-D and lumped-parameter models discussed above, in the context of flow through cardiac tissue, porous medium flow models have been coupled with high-fidelity and reduced-order models of flow in upstream coronary arteries (Hyde et al. 2014; Di Gregorio et al. 2021; Papamanolis et al. 2021; Kim et al. 2023; Menon et al. 2024) as well as the elasto-mechanics describing cardiac contraction and its effect on the pressure driving coronary flow within the cardiac tissue (Chapelle et al. 2010; Cookson et al. 2012; Di Gregorio et al. 2021). Figure 3(b) shows an example of a Darcy myocardial flow model coupled to upstream 3-D flow in coronary arteries and 0-D distal circulation models. However, as with the low-order micro-circulation models discussed above, a challenge with these approaches is the correspondence between modelling parameters, such as Darcy permeability and pressure, with quantities that can be measured in vivo using current clinical techniques. Recent advances in imaging methods hold promise for more accurate prescription of these parameters.

Delving even deeper into smallest scales of the circulatory system, i.e. the cellular level, fluid dynamics continues to play a prominent role. Typical Reynolds numbers in capillaries are $O(10^{-3})$, where Stokes flow assumptions are valid (Secomb 2017). However, the challenge at these scales is that the continuum description of blood breaks down. It is more accurately characterized as a complex fluid or suspension of predominantly red blood cells, white blood cells and platelets within plasma. These constituents exhibit a complex dynamics, such as migration away from vessel walls (creating a cell-free layer), intricate orientations and motions, as well as non-uniform concentrations within vessels as well as between different branching vessels (Chien et al. 1967a,b; Chien 1970; Fedosov, Noguchi & Gompper 2014; Takeishi et al. 2019). This occurs in response to various factors including their concentration, which varies with scale from ${\sim }50\,\%$ in most of the body to ${\sim }20\,\%$ in small vessels with $30\mathrm {\mu }$m diameters, flow-induced shear forces within vessels and their different material properties (Zhao, Shaqfeh & Narsimhan 2012; Qi & Shaqfeh 2017). The cell-scale dynamics within small vessels has several implications, including for clot formation, drug delivery and, most importantly, blood viscosity (Chien 1970). Interestingly, the dynamics and deformability of blood's constituents lead to the effective viscosity of blood being only 4–5 times that of plasma, compared with an equivalent suspension of rigid spheres which would be 10–100 times as viscous as its solvent (Verzicco 2022). Shear-induced cell migration and deformation within vessels also leads to blood acting as a shear-thinning fluid and produces a non-monotonic trend in blood viscosity with respect to vessel diameter (Pries, Neuhaus & Gaehtgens 1992). Effective blood viscosity first decreases with increasing diameter due to the cell-free near vessel walls. This is referred to as the Fåhræus–Lindqvist effect (Fåhræus & Lindqvist 1931). It reaches a minimum in vessels with diameters of the order of red blood cell diameters, where red blood cells line up in a single-file-like motion, and then increases in vessels with smaller diameters due to the deformation of red blood cells. A review of these effects is available in Secomb (2017), and efforts to simulate them are discussed in Freund (2014).

8. Venous flow

After perfusing tissues throughout the body, blood begins its return journey to the heart via the venous system. The venous circulation consists of veins that carry blood depleted of oxygen and nutrients, and stowed with metabolic wastes, from deep within organs, tissues and cells to the great veins that return to the heart. Compared with arteries, veins operate under lower pressure with low pulsatility, and therefore have thinner walls. For the large veins close to the heart, such as the vena cava, there is some mild pulsatility due to right atrial contraction and respiratory effects.

A unique feature of venous haemodynamics is the presence of one-way valves within most veins in the body, adding complexity to venous flow simulation. These valves are essential for maintaining unidirectional venous flow, especially given the low pressure under which they operate. Their dysfunction, which impairs the ability of the venous network to return blood to the heart, is a primary cause of chronic venous inefficiency and venous hypertension (Bergan et al. 2006). Although these valves operate in very low-speed flows, the formation of an axial and helical jet, the presence of flow recirculation within the sinus, and the pulsatility of flow resulting from the opening and closure of valves have been observed and are hypothesized to play roles in healthy venous function (Hamer, Malone & Silver 1981; Lurie et al. 2003; Lurie & Kistner 2013). In this context, another interesting feature is that when multiple valves are present at venous junctions, the valve openings are consistently observed to be at an angle to each other and this angle correlates with the distance between the valves (Lurie & Kistner 2012). This asymmetry naturally gives rise to helical flows that have been shown to improve efficiency and reduce flow stagnation – thus preventing thrombosis/clotting (Lurie & Kistner 2013; Chen et al. 2018).

The formation of clots within veins – or deep venous thrombosis – presents a significant health concern, with an annual incidence of approximately 1 per 1000 adults in the United States (Cushman 2007). Physical immobility is one of the most significant risk factor for deep venous thrombosis because the absence of muscular activity leads to a loss of oscillatory blood flow around valve areas in veins, thereby promoting thrombus or clot formation (Welsh et al. 2019). Subsequently, these clots can travel to the lungs and get lodged in the pulmonary arteries, leading to a potentially fatal pulmonary embolism. In its prevention, alongside anticoagulation therapy, placing a filter in the inferior vena cava serves as an alternative strategy. This mechanical device is designed to filter thrombi carried by the blood flow within the vein. Using in vitro experiments and numerical simulation, the design of these filters has been characterized and optimized with a focus on maximizing thrombus capture efficiency, minimizing areas of low wall shear stress, and ensuring minimal disruption to blood flow (Singer, Wang & Diachin 2010; Feng, Li & Feng 2022). We recommend the recent review by Watson et al. (2024) for a discussion of modelling efforts in this domain to understand flow-mediated mechanisms driving deep vein thrombosis.

Beyond the pulmonary circulation, abnormalities in venous flow can have far-reaching health effects because the entire circulatory system is intricately linked with other bodily fluid systems. For instance, global multi-scale models of human circulation with its interactions with central nervous system fluids have been developed to explore how venous irregularities might influence the development and progression of neurological diseases (Müller & Toro 2014; Toro et al. 2022).

A theme of central focus, with regard to the haemodynamics of the venous system, is palliative approaches for patients with single ventricle physiology. This is a group of congenital heart defects in which one ventricle (either the left or right) is too underdeveloped to function normally. Typically, these patients undergo a series of three surgeries between birth and 4 years of age to reconfigure their systemic and pulmonary circulations sequentially so that they are driven by the single functioning ventricle. These surgeries are primarily centred around the great veins that return blood to the heart, with the aim of maintaining passive flow in these and the pulmonary vessels due to the absence of the second functioning ventricle. While these surgeries are life saving, the associated complications have spurred the use of CFD for improving surgical designs virtually. In these applications, 3-D computational models are often integrated with lumped-parameter networks, thus offering insights into the local haemodynamics while also capturing the interactions with the broader circulatory system (Corsini et al. 2011; Kung et al. 2020; Schwarz et al. 2021).

During the stage one, Norwood surgery, the pulmonary circulation is connected to the aorta with a shunt to provide blood flow to the lungs. Early computational studies proposed and examined several variations in the position and size of the shunt, and found the larger shunts result in poorer oxygen delivery while leading to increased pulmonary flow (Migliavacca et al. 2001; Hsia et al. 2011). The effect of various degrees of pulmonary banding in the hybrid procedure has been simulated to determine the optimal systemic to pulmonary flow ratio and oxygen delivery (Corsini et al. 2011). Moreover, modelling provides valuable quantitative insight into the optimization of the neo-aortic reconstruction. With patient specific geometries, a smooth aortic arch and a larger anastomosis have been shown to result in higher energy efficiency in the blood flow (Qian et al. 2010; Itatani et al. 2012).

In the second stage, the Glenn procedure, the superior vena cava is directly connected to the pulmonary arteries. By allowing passive blood flow to the lungs from the upper body, it reduces the workload on the single ventricle. The connection can be achieved through either the hemi-Fontan or bidirectional Glenn (Douglas et al. 1999; Kogon et al. 2008), but questions remain as to which results in the best haemodynamic results. Patient-specific simulation initially indicated that the Glenn connection is haemodynamically more efficient compared with hemi-Fontan, due to its direct end-to-side anastomosis (Pekkan et al. 2009). Despite what appears to be significant local power loss differences, follow up studies found no notable distinctions in haemodynamic or physiological parameters between the two models (Kung et al. 2020).

After the third surgery, known as the Fontan procedure, the total cavopulmonary connection is established with inferior and superior vena cava connected to the pulmonary arteries. This distinct circulation causes elevated central venous pressure and reduced pulmonary pulsatility, which are associated with morbidities such as liver fibrosis, exercise intolerance and heart failure (Goldberg et al. 2011). Besides assessing the cavopulmonary connection at rest, computational frameworks have also offered valuable insights into the impacts of respiration and exercise in Fontan patients (Marsden et al. 2007; Kung et al. 2014). They reveal that the elevated resistance at connection and nonlinearly increased energy loss significantly restrict the Fontan circulation's capacity to augment cardiac output during exercise, demonstrating a limited ability to accommodate exercise conditions (Whitehead et al. 2007; Sundareswaran et al. 2008). The connection was then proposed to be offset to prevent the head-on collision of flow (De Leval et al. 1996), but this design introduces uneven hepatic flow distribution and therefore uneven growth in pulmonary vascular tree. Subsequently, the bifurcated Y-graft was introduced and improved energy loss and hepatic flow distribution compared with the offset design in computational studies that were later validated clinically (Marsden et al. 2009; Yang et al. 2012; Martin et al. 2015; Trusty et al. 2019, 2020). Figure 4(a) shows examples of flow fields resulting from these offset and Y-graft configurations. It is worth mentioning that, although power loss is a major focus of these computational studies, Fontan haemodynamics from 100 patients showed that Fontan power loss varies from patient to patient, and elevated levels are correlated with lower systemic flow and cardiac index (Haggerty et al. 2014). As their predictive power has grown, computational simulations continue to be used in investigations of new palliative techniques, such as models of tissue-engineered vascular grafts in the Fontan pathway (Schwarz et al. 2021).

Figure 4. (a) Haemodynamics associated with standard offset grafts and novel Y-grafts in the Fontan procedure to treat single ventricle physiology (Marsden et al. 2010; Yang et al. 2015). The highly disturbed flow and the split between left and right pulmonary arteries are ongoing challenges. (b) Patient-specific models for treating pulmonary stenoses. Comparing exclusively proximal stenting vs more extensive stenting showed improved pressure alleviation with extensive stenting (Lan et al. 2022).

9. Right heart

After its long journey from large blood vessels through the depths of tissues and eventually through the venous system, de-oxygenated blood returns to the right side of the heart, which consists of the right atrium and the right ventricle. The superior and inferior vena cava first direct de-oxygenated blood into the right atrium. From there, it flows into the right ventricle, which is the secondary pumping chamber of the heart, responsible for driving flow to the lungs to be re-oxygenated. Blood moving from the right atrium to the right ventricle flows through the tricuspid valve that separates these chambers, and is initially driven by its own inertia. Eventually, an active right atrial contraction completes the ventricular filling. Subsequently, the right ventricle contracts to increase the blood pressure therein, thus closing the tricuspid valve and opening the pulmonary valve to pump blood into the lungs. The left and right sides of the heart, working in parallel, operate synchronously to ensure equal volumes of blood flow, maintaining equilibrium in the circulatory system. The atrial and ventricular septums separate the right and left atria and ventricles, respectively.

Despite similarities to the left heart, the right atrium possess unique anatomical and functional traits. The right atrium connects to the right ventricle through the tricuspid valve, while the left atrium links to the left ventricle via the bicuspid mitral valve. Both atria feature trabeculated appendages that may enhance stagnant flow and promote clot formation (or thrombosis). However, the right atrial appendage has a larger junction than its left counterpart, which likely accounts for reduced thrombosis in the right atrium appendage compared with the left (see § 11) (Lang et al. 2022). Flow patterns within the right atrium are significantly influenced by the inflows from the superior and inferior vena cava. Studies have observed that these incoming flows engage in a circular motion rather than colliding directly, which results in predominantly rotational energy that tends to form a single vortex in the right atrium (Kilner et al. 2000). The kinetic energy is better preserved compared with configurations with multiple vortices or less prevalent helical flow, and conversely, atria without this single vortex flow tend to rely more on energy from atrial contractions (Dewhurst et al. 2020).

Often termed the ‘forgotten ventricle’, the right ventricle also exhibits distinct characteristics from the left ventricle. It is a thin-walled crescent-shaped chamber that wraps around the left ventricle and pumps against the low-impedance pulmonary circuit at a lower pressure of approximately 5–20 mmHg (667–2650 Pa). As such, its mass is roughly 1/6th that of the left ventricle, and its stroke work is approximately a quarter that of the left ventricle. Interestingly, although both right and left ventricular contractions are predominantly driven by shortening of the long-axis, left ventricular contractions significantly affect right ventricular function. The bulging of the ventricular septum due to high pressure in left ventricle reduces the right ventricle volume, and about 20 % to 40 % of the right ventricular systolic pressure and volume outflow result from left ventricular contraction (Santamore & Dell'Italia 1998; Gaine et al. 2014). The right ventricle is also more afterload sensitive. When afterload rises, the right ventricle experiences a significant surge in energy use and a decrease in efficiency compared with the left ventricle, causing the right ventricle's stroke volume to decrease more rapidly. This sensitivity could explain why systemic arterial hypertension is better tolerated than pulmonary arterial hypertension (Bravo et al. 2022).

The flow field in the right ventricle differs from the left ventricle largely due to its different geometry. A number of MRI and CFD studies have shown that the general flow path through the right ventricle is dominated by the smooth curvature of the chamber with a less acute change in direction between the inflow and the outflow tract (Kilner et al. 2000; Fredriksson et al. 2011). A ring-like vortex forms behind the tricuspid valve at the early diastole, which interacts with the narrower boundaries of the right ventricle, particularly on the septal side that is flatter, quickly breaking into a weakly turbulent flow pattern (Kilner et al. 2000; Markl, Kilner & Ebbers 2011). Thus, the flow in the right ventricle features a higher dissipative state, whereas the left ventricle preserves more circulation and kinetic energy throughout diastole. Additionally, its shape helps streamline the flow, promoting the formation of a helical flow towards the pulmonary orifice (Mangual, Domenichini & Pedrizzetti 2012). During systolic contraction, the residual vorticity extends along the outflow tract, adding a slightly helical pattern to the largely irrotational velocity field, and assists in mixing of the blood (Collia et al. 2021). A number of CFD simulations have also characterized the effect of endocardial structures, such as trabeculae and papillary muscles, on disrupting the main vortex ring and dissipating energy in intra-ventricular flows (Sacco et al. 2018a). However, as with the left ventricle, this is usually ignored in both simulations and imaging due to the challenge of simultaneously resolving both the macro-scale flow features and the geometric details of the endocardial surface.

Flow modelling techniques have also been expanded to describe changes linked to pathological conditions, serving as markers for impaired right ventricular function. For example, pulmonary arterial hypertension patients, who have increased pressure and diminished ejection fraction, exhibit a significant increase in right ventricle kinetic energy work density and pronounced pulmonary artery energy loss across the cardiac cycle, which correlates with right ventricle ejection fraction (Han et al. 2015). Moreover, haemodynamic analysis has revealed an increase in both right atrium and right ventricle vorticity at late diastole in pulmonary artery hypertension patients, suggesting its potential interplay between right ventricle fluid dynamics and tissue mechanics leading to pulmonary artery hypertension (Fenster et al. 2015).

10. Pulmonary circulation

The pulmonary circulation carries oxygen-depleted blood from the right side of the heart to the lungs for oxygen enrichment, and then channels the oxygen-saturated blood back to the left side of the heart. Starting from the right ventricular outflow tract, blood is pumped into the main pulmonary artery, which splits into the left and right pulmonary arteries. As these arteries penetrate the lungs, they further subdivide into smaller vessels, following a pattern analogous to the bronchial tree in the respiratory system. At their finest level, these vessels sprawl out across the alveoli, the tiny air sacs where the critical process of gas exchange occurs. Here, carbon dioxide is expelled from the blood and replaced with fresh oxygen. This newly oxygenated blood is then collected by the pulmonary veins and directed back to the left atrium of the heart, ready to be pumped to the rest of the body.

Owing to the central role of blood oxygenation and its subsequent transport, the focus of pulmonary haemodynamics studies has generally revolved around two important pathologies – pulmonary arterial hypertension and pulmonary stenosis. Pulmonary arterial hypertension is a chronic disease characterized by progressive elevation in pulmonary artery pressure due to increased pulmonary vascular resistance. The progression of pulmonary artery hypertension is driven by the remodelling of vascular tissue in response to elevated haemodynamic pressure. The wall thickness and vascular compliance decrease during this remodelling, which leads to increased right ventricle afterload and ventricle–vascular decoupling (Kheyfets et al. 2013). A primary clinical diagnostic indicator of pulmonary arterial hypertension is the pulmonary vascular resistance. It is calculated as the ratio of the pressure drop through the pulmonary circulation to pulmonary blood flow and signifies the structural changes at the pulmonary arteriole level. Interestingly, this widely used clinical metric does not reflect the fact that blood flow is pulsatile with reflective waves. Therefore, several studies have proposed pulmonary impedance to better predict patient outcomes by taking into account not only resistance but also elasticity and inertia to characterize the steady and oscillatory behaviours of blood flow and dynamic interactions between the heart and the lungs, especially under varying physiological conditions (Hunter et al. 2008; Naeije 2012; Kheyfets et al. 2015).

Pulmonary stenosis, or narrowing of the pulmonary vessel lumen, is another significant cause of pulmonary flow disruption. Peripheral pulmonary artery stenosis is commonly associated with congenital heart diseases, including Williams syndrome, Alagille syndrome and tetralogy of Fallot. Severe peripheral pulmonary artery stenosis results in right ventricular hypertension, right ventricular hypertrophy and eventually right ventricular failure. Additionally, if pulmonary blood flow is not promptly restored, it can lead to irreversible damage to the pulmonary artery (Mainwaring & Hanley 2016). Interventional procedures alleviate the stenosis with catheter-based angioplasty or surgical homograft patches. Nonetheless, such surgeries are challenging since they require long hours of cardiopulmonary bypass and address only the most proximal pulmonary arteries. Moreover, post-operative complications include persistent right ventricle hypertension, tears in the vessel wall, aneurysm formation, in-stent restenosis and even death (Cunningham et al. 2013; Felmly et al. 2023).

In light of these debilitating pathologies, both CFD and 4-D flow MRI have been used to investigate blood flow patterns that provide insights into the physiological mechanisms behind disease progression and clinical outcomes. For example, MRI data have revealed that pulmonary artery hypertension patients exhibit reduced flow, peak velocities, wall shear stress, helicity, vorticity and strain (Tang et al. 2012; Barker et al. 2015; Schäfer et al. 2017; Sieren et al. 2019). Increased recirculation and large vortex formation along the inferior curvature of the main pulmonary artery are also reported in pulmonary artery hypertension patients, which can serve as marker for pulmonary artery hypertension risk stratification. This formation of large-scale vortices stems from vessel dilation and diminished right ventricular function, creating shear layers with dramatically different velocities resulting in flow instabilities. This also explains the decreased wall shear stress at the proximal artery wall, which leads to endothelial dysfunction and maladaptive mechanobiological responses (Reiter et al. 2008, 2013).

While the resolution limitations of MRI hinder its effectiveness in distal branches where disease is thought to initiate, computational simulations can provide an in-depth assessment of haemodynamics in these farther reaches of the pulmonary arterial tree. Automatic algorithms have been developed to segment the pulmonary vasculature below 1 mm resolution (Doel, Gavaghan & Grau 2015). With more detailed anatomical models, blood flow simulations have reported decreased wall shear stress in the proximal arterial tree in severe pulmonary artery hypertension, as well as increased oscillatory shear index and relative residence time in the distal arteries, which indicate more stagnant blood flow (Tang et al. 2012; Kheyfets et al. 2015; Yang et al. 2019; Pillalamarri et al. 2021). Computational models have allowed the exploration of pulmonary artery hypertension linked to congenital heart defects, and found increased wall shear stress with increasing ventricular septal defect size in paediatric pulmonary artery trees (Dong et al. 2021). Moreover, recent work has attempted to couple models of pulmonary haemodynamics with vascular growth and remodelling to simulate mechanobiological metrics of interest and potential treatment strategies for pulmonary artery hypertension patients (Szafron et al. 2023). Image-based CFD can also serve as a virtual platform to evaluate interventions on a patient-specific basis to optimize surgical procedures. Yang, Feinstein & Vignon-Clementel (2016) and Yang et al. (2018) performed virtual pulmonary artery reconstruction for patients with Alagille syndrome and achieved accurate predictions of postoperative pulmonary artery flows via adaptive outflow boundary conditions. Recently, Lan et al. (2022) performed virtual intervention in Williams and Alagille syndrome patients and showed that stenting only proximal lesions is often insufficient to reduce overall vascular resistance and thus right ventricular pressures (figure 4b).

A challenge in such computational models is that of prescribing the boundary conditions required to accurately simulate the haemodynamics both before and after clinical interventions. For example, outflow boundary conditions substantially impact the flow distribution in the distal vessels, which is a significant predictor of pulmonary hypertension severity (Kheyfets et al. 2015). Models also need to consider flow autoregulation in distal pulmonary vasculature to correctly predict haemodynamics post-treatments (Yang et al. 2016, 2018). Another challenge is assessing blood flow in distal pulmonary arteries with diameters between 50 and 500 $\mathrm {\mu }$m, which are beyond typical modelling resolution. To address this, 3-D models have been combined with 0-D morphometric trees that consist of resistances in parallel and series based on existing pulmonary artery morphometric data of vessel segment diameters, lengths and branching patterns (Huang et al. 1996). Results from models with morphometric trees have uncovered increased wall shear stress in the distal small vessels compared with controls (Yang et al. 2019), as well as elevated pulmonary vascular resistance and hypertension due to rarefaction in the pulmonary circulation (Olufsen et al. 2012). They have also reported decreased shear stress along with increased cyclic stretch due to pulmonary vascular resistance at the microvascular level (Bartolo et al. 2022). While it is challenging to accurately recapitulate the wide spatial scales and physiological responses at these scales in computational models, continued progress is important for clinical translation, especially due to the difficulty in accurately resolving these scales using clinical imaging.

11. Left atrium

At the end of our tour, we return to the left side of the heart where oxygen-rich blood returning from the lungs via the pulmonary veins collects in the left atrium. From here, blood is directed back into the left ventricle, via contraction of the left atrium, to restart its journey through the circulatory system.

Haemodynamic studies have revealed intricate flow patterns within the left atrium. The inflow from pulmonary veins are asymmetric, with left pulmonary veins forming a rotational loop while right pulmonary veins are more aligned with the wall (Fyrenius et al. 2001). Such asymmetric filling generates asymmetric atrial flow pathways and tends to preserve blood momentum during the cardiac cycle (Kilner et al. 2000). During systole and mid-diastole, discrete left atrial vortices have been observed both from clinical (Fyrenius et al. 2001; Park et al. 2013) and computational data (Vedula et al. 2015; Zingaro et al. 2021). The formation of atrial vortices is conjectured to produce a better ‘washout’ effect in the left atrium and avoid intra-atrial flow stasis (Fyrenius et al. 2001).

In terms of the haemodynamics, the most interesting feature of the left atrium is the left atrial appendage. Anatomically, the left atrium has a smooth-walled interior, except for the left atrial appendage protruding from its upper left corner. The left atrial appendage is a small and ear-shaped pouch, featuring a ridged interior due to the presence of pectinate muscles. It is more distensible than the rest of the left atrium and is hypothesized to act as a decompression chamber during left ventricular systole and periods of high atrial pressure (Al-Saady, Obel & Camm 1999). However, its structure has a propensity to cause increased flow stasis – which can lead to the formation of blood clots and the risk of stroke. This is particularly important in the context of atrial fibrillation, which is the most prevalent and clinically important cardiac arrhythmia. In patients with this condition, the heart's electrical signals are chaotic, leading to the inefficient atrial pumping. This can result in an irregular and frequently rapid heart rate, along with the risk of the blood pooling and potential clotting. If untreated, atrial fibrillation increases the risk of heart-related complications, such as stroke or heart failure. Notably, over 90 % of left atrial thrombi (or blood clots) originate in the left atrial appendage (Al-Saady et al. 1999). Flow patterns within the left atrial appendage have long been suggested as potential predictors for these events (Agmon et al. 1999; Kamp et al. 1999).

Therefore, the integration of computational methods with medical imaging has been used to provide a comprehensive view of the haemodynamics within the left atrial appendage. Decreased mean and peak velocities in left atrial appendage have been reported for atrial fibrillation patients, thus making the washout field incomplete and more prone to fluid stagnation. A large vortex has been observed in the centre of the left atrium for atrial fibrillation patients, instead of the multiple discrete vortices in healthy left atrium (Park et al. 2013). Four typical left atrial appendage morphologies and their impacts on the haemodynamics have been investigated using CFD, showing the close relation between morphology and thrombosis, in agreement with clinical reports. By modelling flow in distinct left atrial appendage morphologies, simulation results have underscored the critical importance of patient-specific left atrial appendage geometric characteristics in influencing thromboembolic risks (Bosi et al. 2018; Masci et al. 2019). While various CFD models have been employed for the left atrium to address potential shifts from transitional to turbulent flows, it is noteworthy that, in atrial fibrillation contexts, models without turbulence remain viable (Dueñas-Pamplona et al. 2021). For more details about computational modelling of the left atrium and the left atrial appendage, we refer to the review by Valvez et al. (2023).

12. Future outlook

As we conclude our journey through the human circulatory system, we hope to have impressed upon the reader the central role that fluid dynamics plays in a wide range of circulatory phenomena. The decades of research that have uncovered these insights have been primarily motivated by one goal – to use physics and engineering insights in advancing clinical care. To that end, several ongoing challenges remain – especially in the context of widespread clinical translation of computational models.

Chief amongst these is the pressing need for more comprehensive verification and validation of realistic patient-specific computational models and their associated insights in large-scale clinical studies (Steinman et al. 2012). This is a theme that was discussed several times in preceding sections of this review, which highlights its importance in all areas of the circulation. While there have been efforts to validate computational models with in vitro flows (Emendi et al. 2021; Kaiser et al. 2023; Zimmermann et al. 2023), there is an urgent need for in vivo validation using patient-specific clinical data.

A related ongoing challenge is the need to develop methods for more personalized computational models that are informed by a variety of patient-specific clinical measurements – not just anatomical imaging, as is currently the norm. This requires the combination of accurate models with data-driven methods to tune model parameters based on routine as well as novel clinical measurements/imaging (Tran et al. 2017; Menon et al. 2024). The translation of such models to clinical use also critically depends on the use of rigorous uncertainty quantification techniques (Fleeter et al. 2020).

We should also highlight emerging tools to process and interpret the fluid dynamics of cardiovascular function. This includes methods that enable new and higher quality insights from clinical imaging, such as super-resolution and de-noising (Fathi et al. 2020). In this context, there have also been several recent efforts to enhance the utility of low-resolution clinical imaging by augmenting such measurements with physics-based and data-enabled techniques to extract the clinically relevant haemodynamics that is otherwise under-resolved (Bakhshinejad et al. 2017; Habibi et al. 2021; Zhang et al. 2022a,b). Finally, new analysis techniques focused on topological features of the flow field itself, including Lagrangian coherent structures and the topology of the wall shear stress field, have recently been shown to provide promising clinical indicators for clinical risk and disease progression (Arzani et al. 2017; Arzani & Shadden 2018; Morbiducci et al. 2020; Mazzi et al. 2021).

As important as fluid dynamics is to cardiovascular health, future treatments will hinge on more comprehensive predictive models that incorporate new findings from a range of different disciplines, including mechanobiology, genetics, and cell signalling (Schwarz et al. 2023a). This has led to the development of multi-physics models, such as whole-heart models, which combine haemodynamics with the spatio-temporal behaviour of the heart at molecular, cellular and tissue scale – from ionic activation in electrophysiology to the tissue mechanics leading to cardiac contraction – to provide a more detailed description of cardiac function (Sugiura et al. 2022; Verzicco 2022; Bucelli et al. 2023). Models combining haemodynamics with mechanobiology are also being developed to predict vascular growth and remodelling, which occur due to the interplay between haemodynamics and biomechanical stimuli experienced by vessel walls, in realistic patient-specific models (Schwarz et al. 2023b). Finally, while computational modelling has provided many of the crucial insights discussed in this review, such modelling is often not possible in time frames that allow their widespread clinical uptake. To address this, machine learning methods are being developed to automatically construct patient-specific anatomical models as well as extract new insights from routine imaging (Fathi et al. 2020; Kong, Wilson & Shadden 2021; Arzani et al. 2022). There is also a need for real-time cardiovascular simulations, through reduced-order modelling and machine learning methods (Pegolotti et al. 2023).

Efforts to address these challenges will yield validated, accurate and robust digital twins that will accelerate the promise of personalized and non-invasive in silico medicine.