Mathematical modeling of the vascular system
MMathematical modeling of the vascular system
Ju Liu ∗ Ingrid S. Lan † Alison L. Marsden ‡ The vascular network
The cardiovascular system is a closed circuit in whichthe heart pumps blood from the left and right ventri-cles to the systemic and pulmonary circulations. Thetwo circulations work synchronously to deliver oxy-genated blood to the body and the heart itself and toroute deoxygenated blood to the lungs for reoxygena-tion. In each circulation, blood travels away from theheart through the larger elastic arteries and smallermuscular arteries and arterioles before reaching thecapillary network; within the capillaries, an exchangeof oxygen, carbon dioxide, nutrients, and metabolitesoccurs between the blood and surrounding tissue; fi-nally, blood returns to the heart through venules andveins. The distensibility of the vessels gives rise to thepropagation of pressure and flow waves through thefluid at a finite speed. Importantly, since the diame-ters of capillaries and red blood cells are of the sameorder, the multiphase nature of blood must be consid-ered when modeling the microcirculation. In this ar-ticle, we restrict our discussion to mathematical mod-eling of the larger vessels with a continuum approach.Interested readers may refer to [QDMV19] for a dis-cussion on mathematical modeling of the heart.The vascular network exhibits complex topologi-cal, geometric, physiologic, and material properties ∗ Ju Liu is an assistant professor in the Departmentof Mechanics and Aerospace Engineering at Southern Uni-versity of Science and Technology. His email address [email protected]. † Ingrid S. Lan is a graduate student in the Department ofBioengineering at Stanford University. Her email address [email protected]. ‡ Alison L. Marsden is a professor and Vera Moulton WallCenter faculty scholar in the Departments of Pediatrics (Car-diology), Bioengineering, and Institute for Computational &Mathematical Engineering at Stanford University. Her emailaddress is [email protected]. that differ greatly from subject to subject, necessi-tating a patient-specific approach to mathematicalmodeling. Furthermore, blood vessels constitute anintricate biological system that dynamically adaptsto the body’s changing needs. In conjunction withthe nervous system, the vascular endothelium acutely regulates the vessel tone and thus the distributionof flow to different tissues. In response to sustainedchanges in hemodynamic loading, arteries chronically adapt in the form of growth and remodeling (G&R),yielding changes in both the geometry and materialproperties. Experimental studies have, for example,demonstrated the tendency of the arterial wall tothicken in response to sustained increases in bloodpressure and that of the lumen to enlarge in responseto sustained increases in flow [HR02].Proper functioning of the cardiovascular systemcan, however, be disrupted by either congenital or ac-quired diseases. These range from structural abnor-malities present at birth (e.g. underdeveloped ventri-cles, blockages, and holes) on the congenital front, tostructural changes that develop (e.g. atheroscleroticstenoses and aneurysms) on the acquired front. De-spite profound advances in treatments and medicalimaging, cardiovascular diseases remain the leadingcause of death worldwide, and congenital heart de-fects remain a leading cause of birth defect-associatedinfant illness and death in the United States. Gold-standard interventions are still met with troublingcomplications. In particular, stents are frequentlymet with restenosis, bypass grafts experience alarm-ingly high rates of failure, and many diseases requirea “watch and wait” treatment strategy. As alterna-tive interventional approaches are explored, there is apressing need to quantitatively understand and pre-dict disease initiation and progression as well as fail-ure modes observed in the clinic. The performance of1 a r X i v : . [ q - b i o . T O ] F e b 𝜂, 𝑃, 𝑁= Θ𝜂 + 𝜇𝑁 𝐺 Θ, 𝑃, 𝑁 = 𝜇𝑁𝐼 𝜂, 𝐽, 𝑁= Θ𝜂 − 𝑃𝐽 + 𝜇𝑁 𝐴 Θ, 𝐽, 𝑁= −𝑃𝐽 + 𝜇𝑁 A CB
Figure 1: (A) Illustration of thermodynamic potentials: the internal energy I , Helmholtz free energy A ,Gibbs free energy G , and enthlapy H . The blue arrows indicate their relations via Legendre transformations.(B) The infinite slope of the pressure-volume curve illustrates the inability of the Helmholtz energy todescribe incompressible behavior. (C) The Helmholtz free energy can be transformed into the Gibbs freeenergy via a Legendre transformation. The resulting system has a saddle-point nature, which requires specialmathematical techniques for analysis.novel medical device designs and surgical approachesmust also be accurately predicted and evaluated.It is clear that hemodynamics, vascular wall biome-chanics, and cellular biochemical responses are alldeeply intertwined, and comprehensive mathemati-cal modeling of the vascular network on all scalespresents tremendous research opportunities. Just asmodern-day aeronautical and automotive design relyheavily on predictive mathematical modeling, we en-vision a future in which patient-specific modeling ofcardiovascular diseases becomes a routine componentof preventive care, diagnostic care, and treatmentplanning. In this article, we discuss recent researchefforts and our perspectives on the most pressingopen research directions in mathematical modelingof the vascular system. A unified continuum modeling frame-work
Cardiovascular modeling requires consideration ofmultiphysics phenomena, including the interaction offluids, solids, and other relevant physics, such as elec- trophysiology, active contraction, or the transportof reactive chemical species. Computational model-ing of fluid and solid sub-problems has convention-ally employed dichotomous formulations for the twosub-domains, creating challenges and disparities formodeling the coupled system. Unifying the formula-tions in a single continuum framework would indeedbe ideal. Traditionally, the constitutive relations forfinite elasticity have been described in terms of theHelmholtz free energy (also known as the strain en-ergy), likely because it is a function of quantitiesthat can be conveniently measured and controlledin laboratories, including temperature, specific vol-ume, the number of molecules. The Helmholtz freeenergy, however, ceases to be a proper potential inthe incompressible limit, as the specific volume be-comes constrained as a constant; accordingly, ill-conditioned matrix problems arise in numerical anal-yses for incompressible materials. To address this is-sue, a Legendre transformation can be performed onthe Helmholtz free energy with respect to the specificvolume (Figure 1). The resulting thermodynamic po-tential, namely the Gibbs free energy, can then be2sed to describe material behavior in both the incom-pressible and compressible regimes. For incompress-ible flows, it yields the incompressible Navier-Stokesequations; for compressible flows, the Gibbs free en-ergy leads to the pressure primitive variable formu-lation, which has been deemed robust for all-speedflows. Following classical principles of mechanics, thegoverning equations based on the Gibbs free energycan be written in the arbitrary Lagrangian-Eulerian(ALE) description as follows,0 = β θ ∂p∂t (cid:12)(cid:12)(cid:12)(cid:12) χ + β θ ( v − ˆ v ) · ∇ x p + ∇ x · v , (1) = ρ ∂ v ∂t (cid:12)(cid:12)(cid:12)(cid:12) χ + ρ ( v − ˆ v ) · ∇ x v − ∇ x · σ dev + ∇ x p − ρ b , (2)wherein ρ is the density, p is the pressure, β θ :=( ∂ρ/∂p ) /ρ is the isothermal compressibility coeffi-cient, χ represents the spatial coordinate with respectto a suitable reference frame that is chosen case bycase and that moves with domain velocity ˆ v , x isthe spatial coordinate in the Eulerian domain, v isthe velocity, σ dev is the deviatoric component of theCauchy stress, and b is the body force per unit mass.Whereas the volumetric behavior is completely de-scribed by ρ ( p ) and β θ ( p ), the isochoric behavior isdetermined by σ dev . Constitutive modeling refers tothe design of an appropriate form for the Gibbs freeenergy, which dictates the analytical forms of ρ ( p ), β θ ( p ), and σ dev and therefore closes the system. Aswill be revealed, equations (1)-(2) properly describeboth viscous fluids and elastic solids, and we are ac-tively working to extend the framework to inelasticsolids as well. This unified framework bridges the gapbetween techniques for computational fluid dynamics(CFD) and structural dynamics, and further simpli-fies the analysis of fluid-structure interaction (FSI)coupled problems by presenting a unified problem re-quiring only a single solution strategy [LM18]. Blood flow modeling
Vascular flow patterns are strongly influenced by con-duit geometries. Patient-specific hemodynamic mod-eling, which has opened the door to virtual treatment planning on a patient-by-patient basis, requires accu-rate geometric domains to be constructed from clin-ical image data, which are commonly acquired fromcomputed tomography angiography or magnetic res-onance imaging. Until recently, model constructionhas been a manual process requiring extensive userintervention in vessel centerline identification on aset of two-dimensional (2D) image slices, 2D imagesegmentation, and lofting of the segmentations intoa three-dimensional (3D) anatomic bounded volume.Recent advances in machine learning show promise inaccelerating the model building workflow [MWM19].Upon construction of the computational domain,CFD, which has found extensive application inweather forecasting, visual effects in digital media,and aircraft and automobile design, can be used tosimulate blood flow. The mathematical foundationof CFD is the Navier-Stokes equations, a set of non-linear partial differential equations (PDEs) describ-ing the conservation laws for viscous fluids. We notethat while blood is in fact a suspension of blood cellsin plasma that exhibits shear-thinning behavior, theNewtonian fluid model is widely accepted to be agood approximation in the large vessels. Referringback to the unified governing equations (1)-(2), in-compressibility suggests a constant ρ and β θ = 0, andthe Newtonian fluid model suggests σ dev := 2 µ ε dev ,in which µ is the dynamic viscosity, and ε dev is thedeviatoric part of the rate-of-strain tensor. Further-more, the mesh need not be moved for stationaryfluid domains. As a result, the ALE coordinates χ reduce to the Eulerian coordinates x , and ˆ v = .From a theoretical perspective, the nonlinear termin equation (2) presents challenges associated withthe existence and regularity of the solution of theNavier-Stokes equations, which remains among theClay Mathematics Institute’s unsolved MillenniumPrize Problems [Fef06]. A common approach hasthus been to approximate the Navier-Stokes equa-tions via high-fidelity numerical methods, such as fi-nite element, finite volume, or finite difference meth-ods, and extract physically meaningful informationfor analysis and prediction. Model reduction tech-niques, such as the reduced basis method for solv-ing complex parameterized PDEs in low-dimensionalspaces, have also been developed to reduce the com-3 BC D
Figure 2: CFD simulation results of an idealized medical device using the residual-based VMS formulation:Instantaneous velocity magnitude in cm/s for (A) Re = 500, (B) Re = 3500, and (C) Re = 5000. (D)Validation of these results against experimental data.putational cost.The development of CFD techniques has been amajor research thrust in applied mathematics andfluid mechanics in the past several decades. Thissubject is rich and merits a review article on its own;here we restrict comments to finite element formu-lations suitable for complex geometries and coupledproblems. In recent years, the variational multiscale(VMS) formulation has received growing attention asa powerful technique for simulating 3D fluid prob-lems. The VMS formulation is a priori consistenton all scales and can thus robustly achieve higher-order accuracy for both advection- and diffusion-dominated flows. In addition to its mathemati-cal properties, numerical studies indicate that VMSproperly captures statistical quantities relevant toflow physics, such as the energy spectrum and decayof kinetic energy in isotropic and wall-bounded turbu-lent flows [HSF18]. Like the classical stabilized meth-ods, VMS further provides a mechanism to circum-vent the Ladyzhenskaya-Babuˇska-Brezzi (LBB) con-dition associated with the divergence-free constraint,thus enabling arbitrary element types to be used formodeling incompressible flows. In contrast to
LBB-stable formulations [BBF13] which require specific el-ement types for velocity and pressure and is oftenimplementationally inconvenient, these stabilized for-mulations conveniently enable equal-order interpola-tion, largely simplifying both automatic mesh gener-ation and finite element implementation from a prac-tical standpoint. Indeed, the VMS formulation has been used in patient-specific blood flow simulationssince the late 1990s. Figure 2 shows validation of sim-ulated laminar, transitional, and fully turbulent flowsin an idealized medical device against experimentaldata, demonstrating applicability of the VMS formu-lation for various flow regimes.Over the years, the rise in high-performance com-puting and high-resolution imaging has been paral-leled by significant advances in mathematical model-ing techniques for flow simulations, and several open-source software packages are readily available. Onerepresentative example is the SimVascular project.In clinical practice, CFD has already been trans-lated into decision-making in realms including di-agnosis, risk stratification, and treatment planning.As a representative example, HeartFlow, Inc. hasachieved notable commercial success with its frac-tional flow reserve-CT technology for evaluating coro-nary lesions, which received clearance for clinical usefrom the US Food and Drug Administration (FDA) in2019. By simulating the normalized pressure gradientthrough a coronary artery blockage, HeartFlow haseliminated the need for invasive coronary angiogramsin 61% of patient cases.
Arterial wall modeling
In addition to blood flow modeling, arterial wall mod-eling is equally critical to understanding the initiationand progression of cardiovascular diseases. Unlike en-gineering materials, biological tissues are highly de-4 φ φ C i r c u m f e r en t i a l d i r e c t i on Axial directionY XZ
A B D
Figure 3: Numerical tensile testing of a 3D collagen-reinforced vascular tissue model in the axial and cir-cumferential directions. (A) Problem setting. (B) Resulting Cauchy stress distribution in an axially-testedspecimen. (C) Resulting Cauchy stress distribution in a circumferentially-tested specimen. (D) Load-displacement curves demonstrating exponential growth in the stiffness provided by collagen fibers.formable and behave as incompressible anisotropicvisco-hyperelastic materials under physiological load-ing conditions. Moreover, biological tissues have beenfound to chronically adapt to the biomechanical en-vironment through G&R. From a modeling perspec-tive, two distinct phenomena must be captured—the passive nonlinear viscoelasticity and the active
G&Rof the soft tissue.A Lagrangian framework is commonly adopted forthe solid problem, as it enables accurate prescriptionof boundary conditions. Again referring back to theunified governing equations (1)-(2), the ALE coordi-nates χ thus reduce to the Lagrangian coordinates X , and ˆ v = v . Furthermore, the displacement field u can conveniently be obtained by the kinematic re-lation d u /dt = v . Finally, the constitutive relationsfor ρ and σ dev are given in terms of the Gibbs freeenergy G as ρ = (cid:18) ∂G∂p (cid:19) − , σ dev = J − F (cid:18) ∂G∂ C (cid:19) F T = J − ∂G∂ F F T , in which the deformation gradient F , the Jacobiandeterminant J , and the right Cauchy-Green tensor C are given by F := ∂ u /∂ X + I , J := det ( F ) , C := F T F . Computational modeling of incompressible nonlin-ear elasticity in complex geometries has long beena challenge. As mentioned above, the unified frame-work enables extension of the VMS formulation tothe solid problem, thereby circumventing the LBBcondition and enabling equal-order interpolation forincompressible finite elasticity. We do, however, notethat LBB-stable formulations are equally attractive,as their energy stability can be demonstrated a prioriin the following form, ddt (cid:90) Ω X ρ (cid:107) V h (cid:107) + ρ G ich ( ˜ C h ) d Ω X = (cid:90) Ω X ρ V h · B d Ω X + (cid:90) Γ X V h · H d Γ X , in which V h and ˜ C h are discrete approximations to V and ˜ C := J − / C ; V , B , and H are the ve-locity, body force, and surface traction in the La-grangian configuration; and G ich is the isochoric partof the Gibbs free energy. This energy stability esti-mate serves as an “insurance plan” for the numericalscheme, guaranteeing reliability of the numerical re-sults.5igure 4: Illustration of the various constrained mixture configurations. A vector e on B αs is mapped to G α ( s ) e on Ω s x and F αs ( t ) e on Ω t x .In addition to incompressibility, soft tissues alsoexhibit directionally dependent, or anisotropic, ma-terial properties as a result of collagen fiber orien-tations. The directional effect of a reinforcing fibercan be characterized by a structural tensor, and wecan elegantly represent the energy potential G usinginvariants of the strain tensor and structural tensors.Figure 3 illustrates the resulting Cauchy stress dis-tributions of a 3D collagen-reinforced vascular tissuemodel subjected to numerical tensile testing in thecircumferential and axial directions. Recently, wehave additionally considered the inelastic response ofsoft tissues in our modeling of the passive mechanicalbehavior.The active G&R can be modeled by the con-strained mixture theory, which was proposed tomathematically describe the evolution of arterial wallcomposition, morphology, and thus stiffness in re-sponse to sustained alterations in hemodynamic load-ing [HR02]. The theory conceptualizes the arte-rial wall as a mixture of structurally significant con-stituents, such as collagen, elastin, and smooth mus-cle cells. These constituents continually turnover atdifferent production and removal rates, engenderingthe following set of evolution equations, in which thesuperscript α indexes the constituents, ρ α ( t ) = ρ α (0) Q α ( t ) + (cid:90) t m α ( s ) q α ( s, t ) ds. Here, ρ α is the constituent’s density, Q α ( t ) is thefraction of the constituent present at time t = 0 thatsurvives to time t , m α ( s ) is the production rate attime s ∈ [0 , t ], and q α ( s, t ) is the fraction of the con-stituent produced at time s that survives to time t .The overall density of the tissue can then be definedas ρ = (cid:80) α ρ α , and we can naturally define mass frac-tions as φ α := ρ α /ρ . These evolution equations canbe considered as an additional set of constitutive lawsthat yield source and sink terms for the mass balanceequation (1). To account for each constituent’s con-tribution to the overall mechanical behavior, the en-ergy potential is regarded as a mass-weighted averageof the constituent potentials, G = (cid:80) α φ α G α .Furthermore, the natural (or stress-free) configu-ration of each constituent is also allowed to evolveseparately, suggesting that constituents produced atdifferent times with different natural configurationscan coexist. Newly produced constituents are de-posited into the extant tissue at some prestretch, aprocess that can mathematically be described witha constituent-specific mapping ψ αs from the naturalconfiguration B αs to the tissue body Ω s x at time s .We denote its corresponding tangent map as G α ( s ).The evolution of the tissue body is described witha separate mapping ϕ t from a material point X attime 0 to a point x ∈ Ω t x at time t , such that allconstituents are kinematically constrained to deformtogether without any relative motion among the con-6 !" R C " R !" R C " R !" C " R R !" R 𝐀 ≈ 𝐀 ! + 𝐀 " +𝐀 +𝐀 % + $ &’() !" 𝐊 & Velocity magnitude: 𝒪 10 $ cm/sViscosity: 𝒪 10 %$ poiseResistance:𝒪 10 & ~𝒪 10 ’ dyn s/cm ’ Wall stiffness: 𝒪 10 ( dyn/cm $ 𝒬 Density: 𝒪 10 ) g/cm * Figure 5: Illustration of the structure of sub-matrix A in cardiovascular FSI problems. The sub-matrixinvolves contributions from the transient term ( A m ), convection term ( A c ), viscous term ( A vis ), wall stiffness( A s ), and boundary traction terms ( K l ). These terms are scaled by physical parameters that span severalorders of magnitude, reported here in the centimeter-gram-second units.stituents. We denote the tangent map of ϕ t as F ( t ),also known as the deformation gradient. In practice,computing the tissue body’s elastic stress requires ameasure of the change in geometry from B αs at time s to Ω t x at a later time t . This is achieved via a com-position of the mappings, such that a point y ∈ B αs is mapped to a point x = ϕ t ◦ ϕ − s ◦ ψ αs ( y ) ∈ Ω t x . Wecan finally represent the tangent map from B αs to Ω t x as F αs ( t ) = F ( t ) F − ( s ) G α ( s ) . Figure 4 illustrates the kinematics of the constrainedtissue mixture. To close the G&R system, we definethe deviatoric component of the Cauchy stress neededin the momentum balance equation (2) as σ dev ( t ) = J − ( t ) (cid:88) α ∂ ( φ α G α ) ∂ F ( t ) F T ( t ) , in which the energy φ α G α is expressed through a hereditary integral,( φ α G α ) ( t ) = ρ α (0) ρ ( t ) Q α ( t ) G α ( F α ( t ))+ (cid:90) t m α ( s ) ρ ( t ) q α ( s, t ) G α ( F αs ( t )) ds. Studies of the constrained mixture theory demon-strate exciting potential for predicting the structuraland morphological evolution observed in cardiovascu-lar diseases, such as abdominal aortic aneurysms andmyocardial hypertrophy. However, major challengesstill exist for extension of this theory to 3D applica-tions, as tracking the full evolutionary history mayquickly become computationally intractable.
Linear solver
Spatiotemporal discretization of nonlinear PDEs andthe Newton-Raphson method eventually boil downto repeatedly solving a linear system with millions,or even billions, degrees of freedom. In computa-tional science and engineering, solving linear systemsoften comprises the most time-consuming portion of7he analysis and thus requires careful design. The lo-cal nature of discrete differential operators yields lin-ear systems of sparse matrices that can be compactlystored in special representations, such as the com-pressed sparse row format. When used in conjunc-tion with preconditioners, iterative methods such asKrylov subspace methods present the most effectivesolution procedure for sparse matrix problems on su-percomputers. Importantly, the design of a precondi-tioner also dictates the overall efficiency, robustness,and scalability. Consider a fully implicit scheme forthe PDEs (1)-(2), which yields a system of nonlin-ear algebraic equations that can be solved with theconsistent Newton-Raphson method. The resultinglinear system is associated with the following matrix A of a 2 × A = (cid:20) A BC D (cid:21) . The sub-matrices A , B , and C respectively repre-sent the discrete convection-diffusion-reaction opera-tor, gradient operator, and divergence operator, eachwith additional numerical modeling terms from theVMS formulation. The sub-matrix D contains a massmatrix scaled by the isothermal compressibility coef-ficient and additional VMS modeling terms. Inter-ested readers may refer to [LM19] for explicit formu-las of the sub-matrices. Block factorization may beperformed on A to yield matrices of lower triangular,diagonal, and upper triangular structure, A = LDU = (cid:20) I OCA − I (cid:21) (cid:20) A OO S (cid:21) (cid:20)
I A − BO I (cid:21) , in which S := D − CA − B is the Schur complement.The design of a preconditioner for A therefore re-duces to solving smaller systems associated with A and S . This design concept is closely related to theChorin-Teman projection method and can be consid-ered an algebraic procedure wrapping the projectionmethod within an iterative solver. Whereas S is as analgebraic manifestation of the pressure Poisson equa-tion, A contains several modeling contributions, asillustrated in Figure 5 for cardiovascular FSI prob-lems. These contributions in A contain physical pa-rameters that span several orders of magnitude, fur-ther complicating the definition of S , which involves the inverse of A and is thus dense. The matrix-free technique and multigrid method represent two effec-tive strategies that can be adopted. The matrix-freetechnique takes advantage of the fact that iterativemethods do not require the algebraic definition of amatrix and instead only require the matrix’s actionon a vector. The non-sparse contributions thereforedo not need to be explicitly constructed. As an ex-ample, we outline the matrix-free definition of S inAlgorithm 1. Algorithm 1
The matrix-free algorithm for the mul-tiplication of S with a vector x . Compute the matrix-vector multiplication ˆ x ← D x . Compute the matrix-vector multiplication ¯ x ← B x . Solve for ˜ x from the linear system A ˜ x = ¯ x . Compute the matrix-vector multiplication ¯ x ← C ˜ x . return ˆ x − ¯ x .The matrix-free method defines how the two sub-matrices can be solved using the Krylov subspacemethod. Nonetheless, to further accelerate conver-gence, preconditioners can be provided by construct-ing sparse approximations to the sub-matrices. Inter-estingly, both sub-matrices have significant contribu-tions from elliptic operators. Their sparse approxi-mations can thus be regarded as discrete elliptic op-erators and be addressed effectively by the multigridmethod, which scales almost linearly with respect tothe degrees of freedom for these types of matrices[ESW14].In summary, the matrix A is solved via a LDU block decomposition that reduces the problem tosolving the smaller block matrices A and S , both ofwhich can be solved iteratively without actual alge-braic definitions. To accelerate convergence, sparseapproximations of the two sub-matrices can be con-structed to generate multigrid preconditioners with-out losing scalability on supercomputers.8 B C
Figure 6: Comparison of the WSS magnitude on a patient-specific pulmonary arterial model computed with(A) a rigid-wall simulation and (B) a FSI simulation. (C) Pressure over one cardiac cycle on the inlet (red)and two selected outlet surfaces (green and blue). Rigid-wall and FSI results are indicated by dashed andsolid lines, respectively.
Multiphysics modeling of FSI
The coupling of biofluids and biosolids is of criticalimportance in cardiovascular modeling. The Navier-Stokes equations alone are insufficient for predictingthe flow behavior, as such a model (sometimes re-ferred to as the rigid-wall model) neglects the defor-mation of the vessel wall, contractions of the heart,and motion of the heart valves, all of which give riseto a deforming fluid domain. Furthermore, the hemo-dynamic wall shear stress (WSS) and pressure aretwo important mechanical loads dictating the G&Rresponse of the vessel wall. As an illustrative ex-ample, Figure 6 compares the WSS magnitude andpressure computed from rigid-wall and FSI models.When neglecting the vessel wall compliance, hemo-dynamic quantities are consistently overpredicted.As discussed above, the unified continuum mod-eling framework enables us to describe both thefluid and solid sub-problems with the same set ofgoverning equations, albeit with different forms for σ dev , ρ ( p ), and β θ . Given the deforming fluid do-main in cardiovascular FSI problems, the fluid sub-problem requires an additional set of equations todescribe the ALE mesh motion. Most commonly, theharmonic extension algorithm or the pseudo-linear-elasticity algorithm is used to determine ˆ v in equa-tion (2). A Lagrangian description is maintained for the solid sub-problem. Importantly, the unifiedframework enables monolithic coupling of the fluidand solid sub-problems with uniform spatiotempo-ral discretization via the VMS formulation and thegeneralized- α scheme. We note that the single re-sulting linear system for the FSI problem is preciselyof the 2 × + Geometric multiscale modeling
Given the computational expense associated withsolving 3D problems, modeling the entire cardiovas-cular system with 3D models is intractable, even withmodern-day computing facilities. In addition, thespatiotemporal resolution offered by 3D models is notalways necessary for practical problems of interest. Itis indeed possible to exploit the morphology of vesselsto derive simplified models of reduced spatial resolu-tion. These 1D or 0D reduced models can either be9 egregated algorithmSolve for mesh motion N e w t o n - R a ph s o n P r o c e du r e Convergence over the entire continuum body
Exit
Calculate solid displacement 𝑢 Update fluid mesh "𝑢 ! Call reduced modelsolver for coupling 𝝈 " = 𝝈 𝝈 " = 𝝈 $’( Solid sub-domain 𝑢, 𝑝, 𝑣
Fluid sub-domain "𝑢 ! , 𝑝, 𝑣Solve 𝛽 ) ̇𝑝 + ∇ 1 𝒗 = 0𝜌 ̇𝒗 − ∇ 1 𝝈 " + ∇𝑝 − 𝜌𝒃 = 0 using VMS,Gen−α,& block PC Figure 7: Illustration of the Newton-Raphson procedure for a monolithically coupled FSI problem. The blueand magenta colors indicate procedures performed in the solid and fluid sub-domains, respectively. The soliddisplacement can be updated consistently using a segregated algorithm and subsequently used to determinethe ALE mesh motion in the fluid sub-domain.used as standalone models or as models coupled to a3D model as boundary conditions.In the 1D approach, which originated from Leon-hard Euler’s seminal work [Eul62], the 3D Navier-Stokes equations posed on a compliant axisymmetrictube are integrated over the cross-section, such thatthe spatial resolution is collapsed to a single axialdimension along the vessel, ∂A∂t + ∂Q∂z = 0 ,ρ ∂Q∂t + ρ ∂∂z (cid:18) α Q A (cid:19) + A ∂P∂z + K QA − µ ∂ Q∂z = 0 , where the three unknowns are the volumetric flowrate Q , the spatially-averaged pressure P , and thecross-sectional area A . Here, the parameters α and K depend on the assumed velocity profile over thecross-section. For a parabolic profile, for example,the momentum flux correction coefficient α = 4 / K = 8 πµ . A constitutivewall model describing the functional dependence of P on A , such as the following simple algebraic relation, is required to close the system, P ( z, t ) = P ext ( z, t ) + ψ ( A ( t, z )) ,ψ ( A ) = Eh √ π − ν √ A − √ A A , where P ext is the external pressure, A is the referencearea when P = P ext , E is the Young’s modulus, h isthe thickness, and ν is the Poisson’s ratio. Assuming A >
0, the above system is strictly hyperbolic withthe following two distinct real eigenvalues, λ , = α QA ± (cid:115) c + (cid:18) QA (cid:19) α ( α − ,c := (cid:115) Aρ dψ ( A ) dA = (cid:115) Eh √ πA ρ (1 − ν ) A . Since c (cid:29) αQ/A in hemodynamic applications, thewaves travel in two distinct directions. To properlycapture the wave propagation phenomena, the Lax-Wendroff scheme or discontinuous Galerkin methodcan be used to solve this set of nonlinear hyperbolicequations.10 VLARVLungsRA P LV A BCoronary CircuitsP LV Coronary CircuitsP LV LVLARVLungsRA P LV LVLARVLungsRA P LV A BCoronary CircuitsP LV Coronary CircuitsP LV LVLARVLungsRA P LV Figure 8: Geometric multiscale modeling of the en-tire cardiovascular system as a closed-loop system.(A) 3D-0D coupling: the aorta and coronary arter-ies, which are generated from patient-specific imagedata, comprise the 3D domain; the peripheral vascu-lature and chambers of the heart (left and right atriaand ventricles LA, RA, LV, RV) are modeled as 0Dcomponents serving as boundary conditions at the3D inlet and outlets. (B) 3D-1D-0D coupling: Theascending aorta and coronary arteries are preservedin the 3D domain, while the remainder of the aortaand its daughter vessels are replaced with a 1D modelgenerated from the vessel centerlines. The spatial resolution can be eliminated altogetherby further integrating the 1D system over a segmentof the vessel with length l and assuming a smallReynolds number such that the nonlinear convectiveterm can be neglected. With the above choice for ψ ( A ), the resulting 0D model can be derived as fol-lows, C d P dt + Q d − Q p = 0 ,L d Q dt + R Q + P d − P p = 0 , wherein R := KlA , C := 2 A √ A (1 − ν ) lEh √ π , L := ρlA , represent the viscous resistance, vessel compliance,and blood inertance, respectively. Here, the two un-knowns are the longitudinally averaged flow Q andpressure P , and the subscripts p and d denote theproximal and distal values. The above equations alsoarise in electrical circuits or hydraulic networks, andit is indeed popular to establish an electric-hydraulicanalogy, in which the current is analogous to flow,and voltage to pressure. The 0D models can beconstructed as arbitrary combinations of resistance,capacitance, and inductance elements in series andin parallel with proper matching conditions. De-spite the lack of any spatial dependence and thus thefailure to capture wave propagation phenomena, 0Dmodels comprised of compartments representing dis-tinct portions of the vasculature or chambers of theheart are often sufficient for approximating flow dy-namics in the global circulation. Furthermore, thereduction of the 3D governing equations to ordinarydifferential-algebraic equations reduces the computa-tional time by several orders of magnitude. As istrue for many modeling techniques, there is a cleartrade-off between accuracy and speed.Regardless of the choice of model, boundary con-ditions must be properly assigned. In particular,boundary conditions for a 3D FSI model must be con-sistent with the wave propagation dynamics withoutintroducing spurious reflections into the 3D domain.Reduced models of the downstream vasculature thus11epresent an extremely fitting choice. Geometric mul-tiscale modeling, or the coupling of dimensionallyheterogeneous models, further offers an efficient ap-proach for modeling the cardiovascular system, inwhich the vast majority of computing resources arereserved for 3D FSI modeling of a localized region ofinterest. Figure 8 illustrates both 3D-0D and 3D-1D-0D coupling. In the latter, a portion of the initial 3Ddomain is replaced with a 1D model generated fromthe vessel centerlines. We note that just as 3D FSIproblems require appropriate transmission conditionson the fluid-solid interface, geometric multiscale mod-eling requires mathematically and physically soundtransmission conditions. Interested readers may re-fer to [QDMV19] for an elaboration on this topic. Optimization and uncertainty quantifi-cation
Given reliable solution techniques for the forwardproblem of a multiscale multiphysics system, one canfurther invoke the mathematical concept of optimiza-tion to improve designs for medical devices and surgi-cal interventions (commonly referred to as shape op-timization) or to identify modeling parameters. Con-ceptually, optimization seeks to minimize or maxi-mize an objective function, given a set of input pa-rameters subject to certain constraints. Examplesof shape optimization include optimization of theradii, attachment angles, and attachment locations ofshunts and grafts in single-ventricle palliation, suchthat maximized pulmonary flow and an unskewedhepatic flow distribution are achieved while preserv-ing low pressures in the vena cava [YFS + InitializationSearchPollStop
Improved?Improved?Convergent?NoNoNo YesYesYesRefinemesh 𝒙 𝒌 𝒑 𝒑 𝒑 Figure 9: Flowchart for shape optimization using thesurrogate management framework.mechanobiological system of interest.Coupling cardiovascular simulations to optimiza-tion algorithms is particularly challenging, as eachobjective function evaluation requires a solution ofthe forward problem. Given the computationalcost and complexity, the application of conventionalgradient-based optimization methods, which requirenumerical determination of gradients via the finitedifference or adjoint-state method, remains cost pro-hibitive. The surrogate management framework
A B C D
Figure 10: The SimVascular pipeline (from left toright): medical image segmentation, model construc-tion, mesh generation, and simulation.12SMF) has thus emerged as a black-box derivative-free method robust for solving expensive cardiovas-cular optimization problems. The SMF is comprisedof two essential components. The SEARCH step usesa surrogate function, or an approximation of the ob-jective function with reduced computational cost, toidentify design points that are likely to improve theobjective function; the POLL step evaluates points ina positive spanning set of directions around the cur-rent best point in the parameter space (Figure 9). Ittherefore combines the efficiency of methods based onresponse surfaces with the convergence properties ofpattern search methods [Mar14]. However, we notethat to date, cardiovascular problems have been lim-ited to idealized geometries, as automation of the en-tire modeling, meshing, and simulation pipeline re-mains challenging for complex geometries that aredifficult to parameterize.Given the large number of inputs to cardiovascu-lar simulations, several sources of uncertainty exist.These pertain to clinical measurements, image-basedgeometries, material properties, and boundary con-ditions, just to name a few. As simulations are in-creasingly incorporated into the FDA approval pro-cess for medical devices and diagnostic tools, rigor-ously assessing the impact of uncertainty on simula-tion predictions must be considered a priority. Uncer-tainty quantification (UQ) is precisely the mathemat-ical field that seeks to propagate these input uncer-tainties forward, such that simulation outputs are ul-timately quantified with probability density functionsand confidence intervals. Monte Carlo (MC) sam-pling, one of the first approaches proposed for UQ, isunbiased and flexible with arbitrary and potentiallycorrelated inputs. Despite these appealing features,its slow convergence rate makes it cost prohibitive,especially considering the computational expense as-sociated with 3D simulations. Recent work has ex-tended MC to multilevel multifidelity MC, whichmakes use of different spatial resolutions (mesh size)and model fidelities (3D, 1D, 0D), to achieve reducedvariance for a fixed computational cost [FGS + The SimVascular project
To date, SimVascular [simvascular] remains theonly fully open-source software offering a completepipeline (Figure 10) from medical image segmenta-tion to 3D model construction, meshing, and patient-specific vascular FSI simulation with either ALE ora reduced linear membrane formulation [UWM + Open challenges and future directions
Having summarized the state of the art in mathemat-ical modeling of the vascular system, we make use ofthis final section to discuss open challenges worthy offurther investigation.
Multiphysics and multiscale modeling.
Giventhe dynamic interplay among hemodynamics, vascu-lar wall biomechanics, and cellular biochemical re-sponses, there is a compelling need for mathemati-cal models characterizing the impact of hemodynam-ics on mechanotransduction pathways and thus thewall composition and material properties. In addi-tion to the deformation induced by mechanical loads,cell-mediated changes involved in G&R must also bemodeled. Such a mathematical view of all scalesis critical for directly translating disease conditionsinto models with quantifiable parameters and out-puts. For example, predictive modeling of thrombusformation requires the flow-mediated transport phe-nomena to be integrated with microscale coagulationkinetics. More refined multiphysics and multiscalemathematical models are indeed needed for clinicalapplications.
Geometric parameterization.
Despite advancesin optimization algorithms for cardiovascular simula-tions, practical optimization applications with com-plex geometries are still hindered by the lack ofgeometric parameterization and manipulation tech-niques compatible with existing computer-aided de-sign (CAD) frameworks. Extending optimizationstudies beyond idealized geometries would require di-rect manipulation of meshes without introducing dis-continuities or severe mesh distortion. Recent effortsin combining CAD with computer-aided engineering(CAE) demonstrate great promise for seamless inte-gration and automation of geometric manipulationand physics-based simulations [HE10].
Verification and validation.
As mathematicalmodels for multiphysics phenomena become increas-ingly sophisticated, software implementations corre-spondingly grow in complexity. A recent study as-sessing the variability among CFD results from sev-eral research groups found rather large discrepanciesfor a relatively simple clinical test case [ea18]. Theresults from this multi-laboratory challenge signifiedthat in the face of wide adoption of medical simula-tion technology, there is a pressing need for rigorousverification of numerical techniques and equally rig-orous validation of mathematical models. A close col-laboration among applied mathematicians, engineers,experimentalists, and clinicians is essential for estab-lishing an accreditation system for mathematical andcomputational modeling of the vascular system.
Acknowledgments
The authors acknowledge Melody Dong for thepreparation of Figure 5. This work is supportedby the Guangdong-Hong Kong-Macao Joint Labo-ratory for Data-Driven Fluid Mechanics and En-gineering Applications under the award number2020B1212030001, the National Institutes of Healthgrants R01HL139796 and R01EB018302, and the Na-tional Science Foundation grant 1663671.
References [BBF13] D. Boffi, F. Brezzi, and M. Fortin,
Mixed fi-nite element methods and applications , Springer,2013.[ea18] K. Valen-Sendstad et al.,
Real-World Variabilityin the Prediction of Intracranial Aneurysm WallShear Stress: The 2015 International AneurysmCFD Challenge , Cardiovascular Engineering andTechnology (2018), 544–564.[ESW14] H.C. Elman, D.J. Silvester, and A.J. Wathen, Fi-nite elements and fast iterative solvers , Second,Oxford University Press, 2014.[Eul62] L. Euler,
Principia pro motu sanguinis per arte-rias determinando , Opera postuma (1862), 814–823.[Fef06] C.L. Fefferman,
Existence and smoothness of theNavier-Stokes equation , The Millennium PrizeProblems (2006), 67. [FGS +
20] C.M. Fleeter, G. Geraci, D.E. Schiavazzi, A.M.Kahn, and A.L. Marsden,
Multilevel and mul-tifidelity uncertainty quantification for cardio-vascular hemodynamics , Computer Methods inApplied Mechanics and Engineering (2020),113030.[HE10] T.J.R. Hughes and J.A. Evans,
Isogeomet-ric analysis , Proceedings of the internationalcongress of mathematicians 2010, 2010, pp. 299–325.[HR02] J.D. Humphrey and K.R. Rajagopal,
A con-strained mixture model for growth and remod-eling of soft tissues , Mathematical models andmethods in applied sciences (2002), 407–430.[HSF18] T.J.R. Hughes, G. Scovazzi, and L.P. Franca, Multiscale and stabilized methods , Encyclope-dia of Computational Mechanics, Second Edition(2018), 1–64.[LM18] J. Liu and A.L. Marsden,
A unified continuumand variational multiscale formulation for flu-ids, solids, and fluid-structure interaction , Com-puter Methods in Applied Mechanics and Engi-neering (2018), 549–597.[LM19] ,
A robust and efficient iterative methodfor hyper-elastodynamics with nested block pre-conditioning , Journal of Computational Physics (2019), 72–93.[Mar14] A.L. Marsden,
Optimization in cardiovascularmodeling , Annual Review of Fluid Mechanics (2014), 519–546. [MVCF +
13] M.E. Moghadam, I.E. Vignon-Clementel, R.Figliola, A.L. Marsden, and Modeling Of Con-genital Hearts Alliance (MOCHA) Investiga-tors.,
A modular numerical method for implicit0D/3D coupling in cardiovascular finite elementsimulations , Journal of Computational Physics (2013), 63–79.[MWM19] G. Maher, N. Wilson, and A.L. Marsden,
Accel-erating cardiovascular model building with con-volutional neural networks , Medical & Biologi-cal Engineering & Computing (2019), 2319–2335.[QDMV19] A. Quarteroni, L. Ded`e, A. Manzoni, and C. Ver-gara, Mathematical modelling of the human car-diovascular system: data, numerical approxima-tion, clinical applications , Vol. 33, CambridgeUniversity Press, 2019.[SAM10] S. Sankaran, C. Audet, and A.L. Marsden,
A method for stochastic constrained optimiza-tion using derivative-free surrogate patternsearch and collocation , Journal of Computa-tional Physics (2010), 4664–4682.[simvascular]
SimVascular .[UWM +
17] A. Updegrove, N.M. Wilson, J. Merkow, H. Lan,A.L. Marsden, and S.C. Shadden,
SimVascu-lar: An Open Source Pipeline for Cardiovascu-lar Simulation , Annals of Biomedical Engineer-ing (2017), 525––541.[YFS +
13] W. Yang, J.A. Feinstein, S.C. Shadden, I.E.Vignon-Clementel, and A.L. Marsden,
Opti-mization of a Y-Graft Design for Improved Hep-atic Flow Distribution in the Fontan Circula-tion , Journal of Biomechanical Engineering (2013), 011002.(2013), 011002.