Skip to Main Content

Mathematical Modeling of the Vascular System

Ju Liu
Ingrid S. Lan
Alison L. Marsden

Communicated by Notices Associate Editor Reza Malek-Madani

Article cover

The Vascular Network

The cardiovascular system is a closed circuit in which the heart pumps blood from the left and right ventricles to the systemic and pulmonary circulations. The two circulations work synchronously to deliver oxygenated blood to the body and the heart itself and to route deoxygenated blood to the lungs for reoxygenation. In each circulation, blood travels away from the heart through the larger elastic arteries and smaller muscular arteries and arterioles before reaching the capillary network; within the capillaries, an exchange of oxygen, carbon dioxide, nutrients, and metabolites occurs between the blood and surrounding tissue; finally, blood returns to the heart through venules and veins. The distensibility of the vessels gives rise to the propagation of pressure and flow waves through the fluid at a finite speed. Importantly, since the diameters of capillaries and red blood cells are of the same order, the multiphase nature of blood must be considered when modeling the microcirculation. In this article, we restrict our discussion to mathematical modeling of the larger vessels with a continuum approach. Interested readers may refer to QDMV19 for a discussion on mathematical modeling of the heart.

The vascular network exhibits complex topological, geometric, physiologic, and material properties that differ greatly from subject to subject, necessitating a patient-specific approach to mathematical modeling. Furthermore, blood vessels constitute an intricate biological system that dynamically adapts to the body’s changing needs. In conjunction with the nervous system, the vascular endothelium acutely regulates the vessel tone and thus the distribution of flow to different tissues. In response to sustained changes in hemodynamic loading, arteries chronically adapt in the form of growth and remodeling (G&R), yielding changes in both the geometry and material properties. Experimental studies have, for example, demonstrated the tendency of the arterial wall to thicken in response to sustained increases in blood pressure and that of the lumen to enlarge in response to sustained increases in flow HR02.

Figure 1.

(A) Illustration of thermodynamic potentials: the internal energy , Helmholtz free energy , Gibbs free energy , and enthlapy . 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 to describe incompressible behavior. (C) The Helmholtz free energy can be transformed into the Gibbs free energy via a Legendre transformation. The resulting system has a saddle-point nature, which requires special mathematical techniques for analysis.

Proper functioning of the cardiovascular system can, however, be disrupted by either congenital or acquired diseases. These range from structural abnormalities present at birth (e.g., underdeveloped ventricles, blockages, and holes) on the congenital front, to structural changes that develop (e.g., atherosclerotic stenoses and aneurysms) on the acquired front. Despite profound advances in treatments and medical imaging, cardiovascular diseases remain the leading cause of death worldwide, and congenital heart defects remain a leading cause of birth defect-associated infant illness and death in the United States. Gold-standard interventions are still met with troubling complications. In particular, stents are frequently met with restenosis, bypass grafts experience alarmingly high rates of failure, and many diseases require a “watch and wait” treatment strategy. As alternative interventional approaches are explored, there is a pressing need to quantitatively understand and predict disease initiation and progression as well as failure modes observed in the clinic. The performance of novel medical device designs and surgical approaches must also be accurately predicted and evaluated.

It is clear that hemodynamics, vascular wall biomechanics, and cellular biochemical responses are all deeply intertwined, and comprehensive mathematical modeling of the vascular network on all scales presents tremendous research opportunities. Just as modern-day aeronautical and automotive design rely heavily on predictive mathematical modeling, we envision a future in which patient-specific modeling of cardiovascular diseases becomes a routine component of preventive care, diagnostic care, and treatment planning. In this article, we discuss recent research efforts and our perspectives on the most pressing open research directions in mathematical modeling of the vascular system.

A Unified Continuum Modeling Framework

Cardiovascular modeling requires consideration of multiphysics phenomena, including the interaction of fluids, solids, and other relevant physics, such as electrophysiology, active contraction, or the transport of reactive chemical species. Computational modeling of fluid and solid subproblems has conventionally employed dichotomous formulations for the two subdomains, creating challenges and disparities for modeling the coupled system. Unifying the formulations in a single continuum framework would indeed be ideal. Traditionally, the constitutive relations for finite elasticity have been described in terms of the Helmholtz free energy (also known as the strain energy), likely because it is a function of quantities that can be conveniently measured and controlled in laboratories, including temperature, specific volume, and the number of molecules. The Helmholtz free energy, however, ceases to be a proper potential in the incompressible limit, as the specific volume becomes constrained as a constant; accordingly, ill-conditioned matrix problems arise in numerical analyses for incompressible materials. To address this issue, a Legendre transformation can be performed on the Helmholtz free energy with respect to the specific volume (Figure 1). The resulting thermodynamic potential, namely the Gibbs free energy, can then be used to describe material behavior in both the incompressible and compressible regimes. For incompressible flows, it yields the incompressible Navier-Stokes equations; for compressible flows, the Gibbs free energy leads to the pressure primitive variable formulation, which has been deemed robust for all-speed flows. Following classical principles of mechanics, the governing equations based on the Gibbs free energy can be written in the arbitrary Lagrangian-Eulerian (ALE) description as follows:

wherein is the density, is the pressure, is the isothermal compressibility coefficient, represents the spatial coordinate with respect to a suitable reference frame that is chosen case by case and that moves with domain velocity , is the spatial coordinate in the Eulerian domain, is the velocity, is the deviatoric component of the Cauchy stress, and is the body force per unit mass. Whereas the volumetric behavior is completely described by and , the isochoric behavior is determined by . Constitutive modeling refers to the design of an appropriate form for the Gibbs free energy, which dictates the analytical forms of , , and and therefore closes the system. As will be revealed, equations 12 properly describe both viscous fluids and elastic solids, and we are actively working to extend the framework to inelastic solids as well. This unified framework bridges the gap between techniques for computational fluid dynamics (CFD) and structural dynamics, and further simplifies the analysis of fluid-structure interaction (FSI) coupled problems by presenting a unified problem requiring only a single solution strategy LM18.

Figure 2.

CFD simulation results of an idealized medical device using the residual-based VMS formulation: instantaneous velocity magnitude in cm/s for (A) , (B) , and (C) . (D) Validation of these results against experimental data.

Blood Flow Modeling

Vascular flow patterns are strongly influenced by conduit geometries. Patient-specific hemodynamic modeling, which has opened the door to virtual treatment planning on a patient-by-patient basis, requires accurate geometric domains to be constructed from clinical image data, which are commonly acquired from computed tomography angiography or magnetic resonance imaging. Until recently, model construction has been a manual process requiring extensive user intervention in vessel centerline identification on a set of two-dimensional (2D) image slices, 2D image segmentation, and lofting of the segmentations into a three-dimensional (3D) anatomic bounded volume. Recent advances in machine learning show promise in accelerating the model building workflow MWM19.

Upon construction of the computational domain, CFD, which has found extensive application in weather forecasting, visual effects in digital media, and aircraft and automobile design, can be used to simulate blood flow. The mathematical foundation of CFD is the Navier-Stokes equations, a set of nonlinear partial differential equations (PDEs) describing the conservation laws for viscous fluids. We note that while blood is in fact a suspension of blood cells in plasma that exhibits shear-thinning behavior, the Newtonian fluid model is widely accepted to be a good approximation in the large vessels. Referring back to the unified governing equations 12, incompressibility suggests a constant and , and the Newtonian fluid model suggests , in which is the dynamic viscosity, and is the deviatoric part of the rate-of-strain tensor. Furthermore, the mesh need not be moved for stationary fluid domains. As a result, the ALE coordinates reduce to the Eulerian coordinates , and . From a theoretical perspective, the nonlinear term in equation 2 presents challenges associated with the existence and regularity of the solution of the Navier-Stokes equations, which remains among the Clay Mathematics Institute’s unsolved Millennium Prize Problems Fef06. A common approach has thus been to approximate the Navier-Stokes equations via high-fidelity numerical methods, such as finite element, finite volume, or finite difference methods, and extract physically meaningful information for analysis and prediction. Model reduction techniques, such as the reduced basis method for solving complex parameterized PDEs in low-dimensional spaces, have also been developed to reduce the computational cost.

The development of CFD techniques has been a major research thrust in applied mathematics and fluid mechanics in the past several decades. This subject is rich and merits a review article on its own; here we restrict comments to finite element formulations suitable for complex geometries and coupled problems. In recent years, the variational multiscale (VMS) formulation has received growing attention as a powerful technique for simulating 3D fluid problems. The VMS formulation is a priori consistent on all scales and can thus robustly achieve higher-order accuracy for both advection- and diffusion-dominated flows. In addition to its mathematical properties, numerical studies indicate that VMS properly captures statistical quantities relevant to flow physics, such as the energy spectrum and decay of kinetic energy in isotropic and wall-bounded turbulent flows HSF18. Like the classical stabilized methods, VMS further provides a mechanism to circumvent the Ladyzhenskaya-Babuška-Brezzi (LBB) condition associated with the divergence-free constraint, thus enabling arbitrary element types to be used for modeling incompressible flows. In contrast to LBB-stable formulations BBF13 which require specific element types for velocity and pressure and is often implementationally inconvenient, these stabilized formulations conveniently enable equal-order interpolation, largely simplifying both automatic mesh generation and finite element implementation from a practical standpoint. Indeed, the VMS formulation has been used in patient-specific blood flow simulations since the late 1990s. Figure 2 shows validation of simulated laminar, transitional, and fully turbulent flows in an idealized medical device against experimental data, demonstrating applicability of the VMS formulation for various flow regimes.

Over the years, the rise in high-performance computing and high-resolution imaging has been paralleled by significant advances in mathematical modeling techniques for flow simulations, and several open-source software packages are readily available. One representative example is the SimVascular project. In clinical practice, CFD has already been translated into decision-making in realms including diagnosis, risk stratification, and treatment planning. As a representative example, HeartFlow, Inc. has achieved notable commercial success with its fractional flow reserve-CT technology for evaluating coronary lesions, which received clearance for clinical use from the US Food and Drug Administration (FDA) in 2019. By simulating the normalized pressure gradient through a coronary artery blockage, HeartFlow has eliminated the need for invasive coronary angiograms in 61% of patient cases.

Figure 3.

Numerical tensile testing of a 3D collagen-reinforced vascular tissue model in the axial and circumferential directions. (A) Problem setting. (B) Resulting Cauchy stress distribution in an axially-tested specimen. (C) Resulting Cauchy stress distribution in a circumferentially-tested specimen. (D) Load-displacement curves demonstrating exponential growth in the stiffness provided by collagen fibers.

Arterial Wall Modeling

In addition to blood flow modeling, arterial wall modeling is equally critical to understanding the initiation and progression of cardiovascular diseases. Unlike engineering materials, biological tissues are highly deformable and behave as incompressible anisotropic visco-hyperelastic materials under physiological loading conditions. Moreover, biological tissues have been found to chronically adapt to the biomechanical environment through G&R. From a modeling perspective, two distinct phenomena must be captured—the passive nonlinear viscoelasticity and the active G&R of the soft tissue.

A Lagrangian framework is commonly adopted for the solid problem, as it enables accurate prescription of boundary conditions. Again referring back to the unified governing equations 12, the ALE coordinates thus reduce to the Lagrangian coordinates , and . Furthermore, the displacement field can conveniently be obtained by the kinematic relation . Finally, the constitutive relations for and are given in terms of the Gibbs free energy as

in which the deformation gradient , the Jacobian determinant , and the right Cauchy-Green tensor are given by

Computational modeling of incompressible nonlinear elasticity in complex geometries has long been a challenge. As mentioned above, the unified framework enables extension of the VMS formulation to the solid problem, thereby circumventing the LBB condition and enabling equal-order interpolation for incompressible finite elasticity. We do, however, note that LBB-stable formulations are equally attractive, as their energy stability can be demonstrated a priori in the following form:

in which and are discrete approximations to and ; , , and are the velocity, body force, and surface traction in the Lagrangian configuration; and is the isochoric part of the Gibbs free energy. This energy stability estimate serves as an “insurance plan” for the numerical scheme, guaranteeing reliability of the numerical results.

Figure 4.

Illustration of the various constrained mixture configurations. A vector on is mapped to on and on .

In addition to incompressibility, soft tissues also exhibit directionally dependent, or anisotropic, material properties as a result of collagen fiber orientations. The directional effect of a reinforcing fiber can be characterized by a structural tensor, and we can elegantly represent the energy potential using invariants of the strain tensor and structural tensors. Figure 3 illustrates the resulting Cauchy stress distributions of a 3D collagen-reinforced vascular tissue model subjected to numerical tensile testing in the circumferential and axial directions. Recently, we have additionally considered the inelastic response of soft tissues in our modeling of the passive mechanical behavior.

The active G&R can be modeled by the constrained mixture theory, which was proposed to mathematically describe the evolution of arterial wall composition, morphology, and thus stiffness in response to sustained alterations in hemodynamic loading HR02. The theory conceptualizes the arterial wall as a mixture of structurally significant constituents, such as collagen, elastin, and smooth muscle cells. These constituents continually turnover at different production and removal rates, engendering the following set of evolution equations, in which the superscript indexes the constituents:

Here, is the constituent’s density, is the fraction of the constituent present at time that survives to time , is the production rate at time , and is the fraction of the constituent produced at time that survives to time . The overall density of the tissue can then be defined as , and we can naturally define mass fractions as . These evolution equations can be considered as an additional set of constitutive laws that yield source and sink terms for the mass balance equation 1. To account for each constituent’s contribution to the overall mechanical behavior, the energy potential is regarded as a mass-weighted average of the constituent potentials, .

Furthermore, the natural (or stress-free) configuration of each constituent is also allowed to evolve separately, suggesting that constituents produced at different times with different natural configurations can coexist. Newly produced constituents are deposited into the extant tissue at some prestretch, a process that can mathematically be described with a constituent-specific mapping from the natural configuration to the tissue body at time . We denote its corresponding tangent map as . The evolution of the tissue body is described with a separate mapping from a material point at time to a point at time , such that all constituents are kinematically constrained to deform together without any relative motion among the constituents. We denote the tangent map of as , also known as the deformation gradient. In practice, computing the tissue body’s elastic stress requires a measure of the change in geometry from at time to at a later time . This is achieved via a composition of the mappings, such that a point is mapped to a point . We can finally represent the tangent map from to as

Figure 4 illustrates the kinematics of the constrained tissue mixture. To close the G&R system, we define the deviatoric component of the Cauchy stress needed in the momentum balance equation 2 as

in which the energy is expressed through a hereditary integral,

Studies of the constrained mixture theory demonstrate exciting potential for predicting the structural and morphological evolution observed in cardiovascular diseases, such as abdominal aortic aneurysms and myocardial hypertrophy. However, major challenges still exist for extension of this theory to 3D applications, as tracking the full evolutionary history may quickly become computationally intractable.

Figure 5.

Illustration of the structure of submatrix in cardiovascular FSI problems. The submatrix involves contributions from the transient term (), convection term (), viscous term (), wall stiffness (), and boundary traction terms (). These terms are scaled by physical parameters that span several orders of magnitude, reported here in the centimeter-gram-second units.

Linear Solver

Spatiotemporal discretization of nonlinear PDEs and the Newton-Raphson method eventually boil down to repeatedly solving a linear system with millions, or even billions, degrees of freedom. In computational science and engineering, solving linear systems often comprises the most time-consuming portion of the analysis and thus requires careful design. The local nature of discrete differential operators yields linear systems of sparse matrices that can be compactly stored in special representations, such as the compressed sparse row format. When used in conjunction with preconditioners, iterative methods such as Krylov subspace methods present the most effective solution procedure for sparse matrix problems on supercomputers. Importantly, the design of a preconditioner also dictates the overall efficiency, robustness, and scalability. Consider a fully implicit scheme for the PDEs 12, which yields a system of nonlinear algebraic equations that can be solved with the consistent Newton-Raphson method. The resulting linear system is associated with the following matrix of a block structure:

The submatrices , , and respectively represent the discrete convection-diffusion-reaction operator, gradient operator, and divergence operator, each with additional numerical modeling terms from the VMS formulation. The submatrix contains a mass matrix scaled by the isothermal compressibility coefficient and additional VMS modeling terms. Interested readers may refer to LM19 for explicit formulas of the submatrices. Block factorization may be performed on to yield matrices of lower triangular, diagonal, and upper triangular structure,

in which is the Schur complement. The design of a preconditioner for therefore reduces to solving smaller systems associated with and . This design concept is closely related to the Chorin-Teman projection method and can be considered an algebraic procedure wrapping the projection method within an iterative solver. Whereas is as an algebraic manifestation of the pressure Poisson equation, contains several modeling contributions, as illustrated in Figure 5 for cardiovascular FSI problems. These contributions in contain physical parameters that span several orders of magnitude, further complicating the definition of , which involves the inverse of and is thus dense. The matrix-free technique and multigrid method represent two effective strategies that can be adopted. The matrix-free technique takes advantage of the fact that iterative methods do not require the algebraic definition of a matrix and instead only require the matrix’s action on a vector. The nonsparse contributions therefore do not need to be explicitly constructed. As an example, we outline the matrix-free definition of in Algorithm 1.

Figure 6.

Comparison of the WSS magnitude on a patient-specific pulmonary arterial model computed with (A) a rigid-wall simulation and (B) an 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 and solid lines, respectively.

Algorithm 1.

The matrix-free algorithm for the multiplication of with a vector .

The matrix-free method defines how the two submatrices can be solved using the Krylov subspace method. Nonetheless, to further accelerate convergence, preconditioners can be provided by constructing sparse approximations to the submatrices. Interestingly, both submatrices have significant contributions from elliptic operators. Their sparse approximations can thus be regarded as discrete elliptic operators and be addressed effectively by the multigrid method, which scales almost linearly with respect to the degrees of freedom for these types of matrices ESW14.

In summary, the matrix is solved via an block decomposition that reduces the problem to solving the smaller block matrices and , both of which can be solved iteratively without actual algebraic definitions. To accelerate convergence, sparse approximations of the two submatrices can be constructed to generate multigrid preconditioners without losing scalability on supercomputers.

Multiphysics Modeling of FSI

The coupling of biofluids and biosolids is of critical importance in cardiovascular modeling. The Navier-Stokes equations alone are insufficient for predicting the flow behavior, as such a model (sometimes referred to as the rigid-wall model) neglects the deformation of the vessel wall, contractions of the heart, and motion of the heart valves, all of which give rise to a deforming fluid domain. Furthermore, the hemodynamic wall shear stress (WSS) and pressure are two important mechanical loads dictating the G&R response of the vessel wall. As an illustrative example, Figure 6 compares the WSS magnitude and pressure computed from rigid-wall and FSI models. When neglecting the vessel wall compliance, hemodynamic quantities are consistently overpredicted.

Figure 7.

Illustration of the Newton-Raphson procedure for a monolithically coupled FSI problem. The blue and magenta colors indicate procedures performed in the solid and fluid subdomains, respectively. The solid displacement can be updated consistently using a segregated algorithm and subsequently used to determine the ALE mesh motion in the fluid subdomain.

As discussed above, the unified continuum modeling framework enables us to describe both the fluid and solid subproblems with the same set of governing equations, albeit with different forms for , , and . Given the deforming fluid domain in cardiovascular FSI problems, the fluid subproblem requires an additional set of equations to describe the ALE mesh motion. Most commonly, the harmonic extension algorithm or the pseudo-linear-elasticity algorithm is used to determine in equation 2. A Lagrangian description is maintained for the solid subproblem. Importantly, the unified framework enables monolithic coupling of the fluid and solid subproblems with uniform spatiotemporal discretization via the VMS formulation and the generalized- scheme. We note that the single resulting linear system for the FSI problem is precisely of the block structure discussed above LM18. As a result, the data structures and solution methods for the FSI problem can be constructed in a manner practically identical to that of CFD (Figure 7). To enable flexible coupling between dimensionally heterogeneous models (see the next section), a modular approach is adopted within each Newton-Raphson iteration to communicate information between the 3D and reduced models in a fashion similar to the Gauss-Seidel method MVCF13.

Geometric Multiscale Modeling

Given the computational expense associated with solving 3D problems, modeling the entire cardiovascular system with 3D models is intractable, even with modern-day computing facilities. In addition, the spatiotemporal resolution offered by 3D models is not always necessary for practical problems of interest. It is indeed possible to exploit the morphology of vessels to derive simplified models of reduced spatial resolution. These 1D or 0D reduced models can either be used as standalone models or as models coupled to a 3D model as boundary conditions.

In the 1D approach, which originated from Leonhard Euler’s seminal work Eul62, the 3D Navier-Stokes equations posed on a compliant axisymmetric tube are integrated over the cross-section, such that the spatial resolution is collapsed to a single axial dimension along the vessel,

where the three unknowns are the volumetric flow rate , the spatially-averaged pressure , and the cross-sectional area . Here, the parameters and depend on the assumed velocity profile over the cross-section. For a parabolic profile, for example, the momentum flux correction coefficient , and the friction parameter . A constitutive wall model describing the functional dependence of on , such as the following simple algebraic relation, is required to close the system:

where is the external pressure, is the reference area when , is the Young’s modulus, is the thickness, and is the Poisson’s ratio. Assuming , the above system is strictly hyperbolic with the following two distinct real eigenvalues:

Since in hemodynamic applications, the waves travel in two distinct directions. To properly capture the wave propagation phenomena, the Lax-Wendroff scheme or discontinuous Galerkin method can be used to solve this set of nonlinear hyperbolic equations.

Figure 8.

Geometric multiscale modeling of the entire cardiovascular system as a closed-loop system. (A) 3D-0D coupling: the aorta and coronary arteries, which are generated from patient-specific image data, comprise the 3D domain; the peripheral vasculature and chambers of the heart (left and right atria and ventricles LA, RA, LV, RV) are modeled as 0D components serving as boundary conditions at the 3D inlet and outlets. (B) 3D-1D-0D coupling: the ascending aorta and coronary arteries are preserved in the 3D domain, while the remainder of the aorta and its daughter vessels are replaced with a 1D model generated from the vessel centerlines.

The spatial resolution can be eliminated altogether by further integrating the 1D system over a segment of the vessel with length and assuming a small Reynolds number such that the nonlinear convective term can be neglected. With the above choice for , the resulting 0D model can be derived as follows:

wherein

represent the viscous resistance, vessel compliance, and blood inertance, respectively. Here, the two unknowns are the longitudinally averaged flow and pressure , and the subscripts and denote the proximal and distal values. The above equations also arise in electrical circuits or hydraulic networks, and it is indeed popular to establish an electric-hydraulic analogy, in which the current is analogous to flow, and voltage to pressure. The 0D models can be constructed as arbitrary combinations of resistance, capacitance, and inductance elements in series and in parallel with proper matching conditions. Despite the lack of any spatial dependence and thus the failure to capture wave propagation phenomena, 0D models comprised of compartments representing distinct portions of the vasculature or chambers of the heart are often sufficient for approximating flow dynamics in the global circulation. Furthermore, the reduction of the 3D governing equations to ordinary differential-algebraic equations reduces the computational time by several orders of magnitude. As is true for many modeling techniques, there is a clear trade-off between accuracy and speed.

Regardless of the choice of model, boundary conditions must be properly assigned. In particular, boundary conditions for a 3D FSI model must be consistent with the wave propagation dynamics without introducing spurious reflections into the 3D domain. Reduced models of the downstream vasculature thus represent an extremely fitting choice. Geometric multiscale modeling, or the coupling of dimensionally heterogeneous models, further offers an efficient approach for modeling the cardiovascular system, in which the vast majority of computing resources are reserved for 3D FSI modeling of a localized region of interest. Figure 8 illustrates both 3D-0D and 3D-1D-0D coupling. In the latter, a portion of the initial 3D domain is replaced with a 1D model generated from the vessel centerlines. We note that just as 3D FSI problems require appropriate transmission conditions on the fluid-solid interface, geometric multiscale modeling requires mathematically and physically sound transmission conditions. Interested readers may refer to QDMV19 for an elaboration on this topic.

Optimization and Uncertainty Quantification

Given reliable solution techniques for the forward problem of a multiscale multiphysics system, one can further invoke the mathematical concept of optimization to improve designs for medical devices and surgical interventions (commonly referred to as shape optimization) or to identify modeling parameters. Conceptually, optimization seeks to minimize or maximize an objective function, given a set of input parameters subject to certain constraints. Examples of shape optimization include optimization of the radii, attachment angles, and attachment locations of shunts and grafts in single-ventricle palliation, such that maximized pulmonary flow and an unskewed hepatic flow distribution are achieved while preserving low pressures in the vena cava YFS13. On the other hand, optimization can also be used as a tool for data assimilation, which seeks to identify modeling parameters that best reproduce experimental or clinical data. For example, external tissue support parameters have previously been identified by minimizing the discrepancy between model and image contours of the vessel lumen, and G&R modeling parameters have been identified to achieve homeostatic stress states experimentally observed in the vessel. In all cases, appropriate objective functions must be determined to measure performance, a process that requires deep understanding of the clinical disease or mechanobiological system of interest.

Figure 9.

Flowchart for shape optimization using the surrogate management framework.

Coupling cardiovascular simulations to optimization algorithms is particularly challenging, as each objective function evaluation requires a solution of the forward problem. Given the computational cost and complexity, the application of conventional gradient-based optimization methods, which require numerical determination of gradients via the finite difference or adjoint-state method, remains cost prohibitive. The surrogate management framework (SMF) has thus emerged as a black-box derivative-free method robust for solving expensive cardiovascular optimization problems. The SMF is comprised of two essential components. The SEARCH step uses a surrogate function, or an approximation of the objective function with reduced computational cost, to identify design points that are likely to improve the objective function; the POLL step evaluates points in a positive spanning set of directions around the current best point in the parameter space (Figure 9). It therefore combines the efficiency of methods based on response surfaces with the convergence properties of pattern search methods Mar14. However, we note that to date, cardiovascular problems have been limited to idealized geometries, as automation of the entire modeling, meshing, and simulation pipeline remains challenging for complex geometries that are difficult to parameterize.

Given the large number of inputs to cardiovascular simulations, several sources of uncertainty exist. These pertain to clinical measurements, image-based geometries, material properties, and boundary conditions, just to name a few. As simulations are increasingly incorporated into the FDA approval process for medical devices and diagnostic tools, rigorously assessing the impact of uncertainty on simulation predictions must be considered a priority. Uncertainty quantification (UQ) is precisely the mathematical field that seeks to propagate these input uncertainties forward, such that simulation outputs are ultimately quantified with probability density functions and confidence intervals. Monte Carlo (MC) sampling, one of the first approaches proposed for UQ, is unbiased and flexible with arbitrary and potentially correlated inputs. Despite these appealing features, its slow convergence rate makes it cost prohibitive, especially considering the computational expense associated with 3D simulations. Recent work has extended MC to multilevel multifidelity MC, which makes use of different spatial resolutions (mesh size) and model fidelities (3D, 1D, 0D), to achieve reduced variance for a fixed computational cost FGS20. As an alternative to MC, polynomial chaos expansions describe smooth stochastic responses using interpolating polynomials that are mutually orthogonal with respect to the probability measure of the random inputs and can achieve up to exponential convergence rates with respect to the polynomial order. Significant research has also focused on adaptive methods appropriate for applications involving discontinuous stochastic responses. We conclude this section by emphasizing that optimization and UQ are two important mathematical tools that, when combined, enable identification of optimal designs and parameters that are robust to fluctuations in modeling choices and surgical implementation SAM10.

Figure 10.

The SimVascular pipeline (from left to right): medical image segmentation, model construction, mesh generation, and simulation.

Figure 11.

Geographical distribution of SimVascular users.

The SimVascular Project

To date, SimVascular [simvascular] remains the only fully open-source software offering a complete pipeline (Figure 10) from medical image segmentation to 3D model construction, meshing, and patient-specific vascular FSI simulation with either ALE or a reduced linear membrane formulation UWM17. It was originally developed in Charles Taylor’s lab at Stanford University and released in 2007. At the time, however, the integration of several licensed commercial components hindered new user adoption and prevented complete open source release. In 2013, the senior author of this article accordingly launched a joint revitalization effort with collaborators, aiming to integrate open-source alternatives and improve computational techniques for all stages of the pipeline. Nowadays, SimVascular has attracted over 4500 domestic and international users (Figure 11). The development team continues to host conference workshops, Youtube tutorials, and a user forum to foster further advances in cardiovascular research. It has also proven to be an effective educational tool in graduate engineering courses, elucidating fundamental principles in CFD, FSI, and human physiology. As an active software project with regularly forthcoming enhancements, it represents the state of the art in cardiovascular simulation. Current research and development efforts span finite element method development for 3D G&R, automated vessel path identification, model manipulation, optimization, UQ, and even virtual reality for visualizing simulation results in real time. Importantly, the motivation underlying these efforts is the ultimate vision of virtual patient-specific treatment planning on a clinically relevant time frame, in which physicians can make informed decisions about the optimal intervention directly from the clinic.

Open Challenges and Future Directions

Having summarized the state of the art in mathematical modeling of the vascular system, we make use of this final section to discuss open challenges worthy of further investigation.

Multiphysics and multiscale modeling

Given the dynamic interplay among hemodynamics, vascular wall biomechanics, and cellular biochemical responses, there is a compelling need for mathematical models characterizing the impact of hemodynamics on mechanotransduction pathways and thus the wall composition and material properties. In addition to the deformation induced by mechanical loads, cell-mediated changes involved in G&R must also be modeled. Such a mathematical view of all scales is critical for directly translating disease conditions into models with quantifiable parameters and outputs. For example, predictive modeling of thrombus formation requires the flow-mediated transport phenomena to be integrated with microscale coagulation kinetics. More refined multiphysics and multiscale mathematical models are indeed needed for clinical applications.

Geometric parameterization

Despite advances in optimization algorithms for cardiovascular simulations, practical optimization applications with complex geometries are still hindered by the lack of geometric parameterization and manipulation techniques compatible with existing computer-aided design (CAD) frameworks. Extending optimization studies beyond idealized geometries would require direct manipulation of meshes without introducing discontinuities or severe mesh distortion. Recent efforts in combining CAD with computer-aided engineering (CAE) demonstrate great promise for seamless integration and automation of geometric manipulation and physics-based simulations HE10.

Verification and validation

As mathematical models for multiphysics phenomena become increasingly sophisticated, software implementations correspondingly grow in complexity. A recent study assessing the variability among CFD results from several research groups found rather large discrepancies for a relatively simple clinical test case ea18. The results from this multilaboratory challenge signified that in the face of wide adoption of medical simulation technology, there is a pressing need for rigorous verification of numerical techniques and equally rigorous validation of mathematical models. A close collaboration among applied mathematicians, engineers, experimentalists, and clinicians is essential for establishing an accreditation system for mathematical and computational modeling of the vascular system.

Acknowledgments

The authors acknowledge Melody Dong for the preparation of Figure 5. This work is supported by the National Institutes of Health grants R01HL139796 and R01EB018302, the National Science Foundation grant 1663671, and the Guangdong-Hong Kong-Macao Joint Laboratory for Data-Driven Fluid Mechanics and Engineering Applications under the award number 2020B1212030001.

References

[BBF13]
Daniele Boffi, Franco Brezzi, and Michel Fortin, Mixed finite element methods and applications, Springer Series in Computational Mathematics, vol. 44, Springer, Heidelberg, 2013, DOI 10.1007/978-3-642-36519-5. MR3097958Show rawAMSref\bib{Boffi2013}{book}{ author={Boffi, Daniele}, author={Brezzi, Franco}, author={Fortin, Michel}, title={Mixed finite element methods and applications}, series={Springer Series in Computational Mathematics}, volume={44}, publisher={Springer, Heidelberg}, date={2013}, pages={xiv+685}, isbn={978-3-642-36518-8}, isbn={978-3-642-36519-5}, review={\MR {3097958}}, doi={10.1007/978-3-642-36519-5}, } Close amsref.
[ESW14]
Howard C. Elman, David J. Silvester, and Andrew J. Wathen, Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics, 2nd ed., Numerical Mathematics and Scientific Computation, Oxford University Press, Oxford, 2014, DOI 10.1093/acprof:oso/9780199678792.001.0001. MR3235759Show rawAMSref\bib{Elman2014}{book}{ author={Elman, Howard C.}, author={Silvester, David J.}, author={Wathen, Andrew J.}, title={Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics}, series={Numerical Mathematics and Scientific Computation}, edition={2}, publisher={Oxford University Press, Oxford}, date={2014}, pages={xiv+479}, isbn={978-0-19-967880-8}, review={\MR {3235759}}, doi={10.1093/acprof:oso/9780199678792.001.0001}, } Close amsref.
[ea18]
K. Valen-Sendstad et al., Real-world variability in the prediction of intracranial aneurysm wall shear stress: The 2015 International Aneurysm CFD Challenge, Cardiovascular Engineering and Technology 9 (2018), 544–564.Show rawAMSref\bib{Sendstad2018}{article}{ author={et~al., K. Valen-Sendstad}, title={Real-{w}orld {v}ariability in the {p}rediction of {i}ntracranial {a}neurysm {w}all {s}hear {s}tress: {T}he 2015 {I}nternational {A}neurysm {CFD} {C}hallenge}, date={2018}, journal={Cardiovascular Engineering and Technology}, volume={9}, pages={544\ndash 564}, } Close amsref.
[Eul62]
L. Euler, Principia pro motu sanguinis per arterias determinando, Opera postuma, 1862, pp. 814–823.Show rawAMSref\bib{euler1862principia}{article}{ label={Eul62}, author={Euler, L.}, title={Principia pro motu sanguinis per arterias determinando}, journal={Opera postuma, 1862}, pages={pp. 814\ndash 823}, } Close amsref.
[Fef06]
Charles L. Fefferman, Existence and smoothness of the Navier-Stokes equation, The millennium prize problems, Clay Math. Inst., Cambridge, MA, 2006, pp. 57–67. MR2238274Show rawAMSref\bib{fefferman2006}{article}{ author={Fefferman, Charles L.}, title={Existence and smoothness of the Navier-Stokes equation}, conference={ title={The millennium prize problems}, }, book={ publisher={Clay Math. Inst., Cambridge, MA}, }, date={2006}, pages={57--67}, review={\MR {2238274}}, } Close amsref.
[FGS20]
Casey M. Fleeter, Gianluca Geraci, Daniele E. Schiavazzi, Andrew M. Kahn, and Alison L. Marsden, Multilevel and multifidelity uncertainty quantification for cardiovascular hemodynamics, Comput. Methods Appl. Mech. Engrg. 365 (2020), 113030, 37, DOI 10.1016/j.cma.2020.113030. MR4089535Show rawAMSref\bib{Fleeter2020}{article}{ author={Fleeter, Casey M.}, author={Geraci, Gianluca}, author={Schiavazzi, Daniele E.}, author={Kahn, Andrew M.}, author={Marsden, Alison L.}, title={Multilevel and multifidelity uncertainty quantification for cardiovascular hemodynamics}, journal={Comput. Methods Appl. Mech. Engrg.}, volume={365}, date={2020}, pages={113030, 37}, issn={0045-7825}, review={\MR {4089535}}, doi={10.1016/j.cma.2020.113030}, } Close amsref.
[HE10]
T. J. R. Hughes and J. A. Evans, Isogeometric analysis, Proceedings of the International Congress of Mathematicians 2010, World Scientific, 2010, pp. 299–325.Show rawAMSref\bib{Hughes2010}{inproceedings}{ label={HE10}, author={Hughes, T. J. R.}, author={Evans, J. A.}, title={Isogeometric analysis}, date={2010}, booktitle={Proceedings of the International Congress of Mathematicians 2010, World Scientific,}, pages={299\ndash 325}, } Close amsref.
[HSF18]
T. J. R. Hughes, G. Scovazzi, and L. P. Franca, Multiscale and stabilized methods, 2nd ed., Encyclopedia of Computational Mechanics, 2018, pp. 1–64.Show rawAMSref\bib{hughes2018}{article}{ label={HSF18}, author={Hughes, T. J. R.}, author={Scovazzi, G.}, author={Franca, L. P.}, title={Multiscale and stabilized methods}, journal={2nd ed., Encyclopedia of Computational Mechanics, 2018, }, pages={pp. 1\ndash 64}, } Close amsref.
[HR02]
J. D. Humphrey and K. R. Rajagopal, A constrained mixture model for growth and remodeling of soft tissues, Math. Models Methods Appl. Sci. 12 (2002), no. 3, 407–430, DOI 10.1142/S0218202502001714. MR1894260Show rawAMSref\bib{Humphrey2002}{article}{ author={Humphrey, J. D.}, author={Rajagopal, K. R.}, title={A constrained mixture model for growth and remodeling of soft tissues}, journal={Math. Models Methods Appl. Sci.}, volume={12}, date={2002}, number={3}, pages={407--430}, issn={0218-2025}, review={\MR {1894260}}, doi={10.1142/S0218202502001714}, } Close amsref.
[LM18]
Ju Liu and Alison L. Marsden, A unified continuum and variational multiscale formulation for fluids, solids, and fluid-structure interaction, Comput. Methods Appl. Mech. Engrg. 337 (2018), 549–597, DOI 10.1016/j.cma.2018.03.045. MR3801792Show rawAMSref\bib{Liu2018}{article}{ author={Liu, Ju}, author={Marsden, Alison L.}, title={A unified continuum and variational multiscale formulation for fluids, solids, and fluid-structure interaction}, journal={Comput. Methods Appl. Mech. Engrg.}, volume={337}, date={2018}, pages={549--597}, issn={0045-7825}, review={\MR {3801792}}, doi={10.1016/j.cma.2018.03.045}, } Close amsref.
[LM19]
Ju Liu and Alison L. Marsden, A robust and efficient iterative method for hyper-elastodynamics with nested block preconditioning, J. Comput. Phys. 383 (2019), 72–93, DOI 10.1016/j.jcp.2019.01.019. MR3908729Show rawAMSref\bib{Liu2019}{article}{ author={Liu, Ju}, author={Marsden, Alison L.}, title={A robust and efficient iterative method for hyper-elastodynamics with nested block preconditioning}, journal={J. Comput. Phys.}, volume={383}, date={2019}, pages={72--93}, issn={0021-9991}, review={\MR {3908729}}, doi={10.1016/j.jcp.2019.01.019}, } Close amsref.
[MWM19]
G. Maher, N. Wilson, and A. L. Marsden, Accelerating cardiovascular model building with convolutional neural networks, Medical & Biological Engineering & Computing 57 (2019), 2319–2335.Show rawAMSref\bib{Maher2019}{article}{ author={Maher, G.}, author={Wilson, N.}, author={Marsden, A. L.}, title={Accelerating cardiovascular model building with convolutional neural networks}, date={2019}, journal={Medical \& Biological Engineering \& Computing}, volume={57}, pages={2319\ndash 2335}, } Close amsref.
[Mar14]
Alison L. Marsden, Optimization in cardiovascular modeling, Annual review of fluid mechanics. Vol. 46, Annu. Rev. Fluid Mech., vol. 46, Annual Reviews, Palo Alto, CA, 2014, pp. 519–546. MR3205456Show rawAMSref\bib{Marsden2014}{article}{ label={Mar14}, author={Marsden, Alison L.}, title={Optimization in cardiovascular modeling}, conference={ title={Annual review of fluid mechanics. Vol. 46}, }, book={ series={Annu. Rev. Fluid Mech.}, volume={46}, publisher={Annual Reviews, Palo Alto, CA}, }, date={2014}, pages={519--546}, review={\MR {3205456}}, } Close amsref.
[MVCF13]
Mahdi Esmaily Moghadam, Irene E. Vignon-Clementel, Richard Figliola, Alison L. Marsden, and Modeling of Congenital Hearts Alliance (MOCHA) Investigators, A modular numerical method for implicit 0D/3D coupling in cardiovascular finite element simulations, J. Comput. Phys. 244 (2013), 63–79, DOI 10.1016/j.jcp.2012.07.035. MR3064206Show rawAMSref\bib{Moghadam2013}{article}{ label={MVCF$^+$13}, author={Esmaily Moghadam, Mahdi}, author={Vignon-Clementel, Irene E.}, author={Figliola, Richard}, author={Marsden, Alison L.}, author={Modeling of Congenital Hearts Alliance (MOCHA) Investigators}, title={A modular numerical method for implicit 0D/3D coupling in cardiovascular finite element simulations}, journal={J. Comput. Phys.}, volume={244}, date={2013}, pages={63--79}, issn={0021-9991}, review={\MR {3064206}}, doi={10.1016/j.jcp.2012.07.035}, } Close amsref.
[QDMV19]
Alfio Quarteroni, Luca Dede’, Andrea Manzoni, and Christian Vergara, Mathematical modelling of the human cardiovascular system: Data, numerical approximation, clinical applications, Cambridge Monographs on Applied and Computational Mathematics, vol. 33, Cambridge University Press, Cambridge, 2019, DOI 10.1017/9781108616096. MR3967732Show rawAMSref\bib{quarteroni2019mathematical}{book}{ author={Quarteroni, Alfio}, author={Dede', Luca}, author={Manzoni, Andrea}, author={Vergara, Christian}, title={Mathematical modelling of the human cardiovascular system}, series={Cambridge Monographs on Applied and Computational Mathematics}, volume={33}, subtitle={Data, numerical approximation, clinical applications}, publisher={Cambridge University Press, Cambridge}, date={2019}, pages={x+280}, isbn={978-1-108-48039-0}, review={\MR {3967732}}, doi={10.1017/9781108616096}, } Close amsref.
[SAM10]
S. Sankaran, C. Audet, and A. L. Marsden, A method for stochastic constrained optimization using derivative-free surrogate pattern search and collocation, J. Comput. Phys. 229 (2010), 4664–4682.Show rawAMSref\bib{Sankaran2010}{article}{ author={Sankaran, S.}, author={Audet, C.}, author={Marsden, A. L.}, title={A method for stochastic constrained optimization using derivative-free surrogate pattern search and collocation}, date={2010}, journal={J. Comput. Phys.}, volume={229}, pages={4664\ndash 4682}, } Close amsref.
[simvascular]
SimVascular, www.simvascular.org.Show rawAMSref\bib{simvascular}{misc}{ label={simvascular}, title={Sim{V}ascular, \url {www.simvascular.org}}, how={\url {www.simvascular.org}}, } Close amsref.
[UWM17]
A. Updegrove, N. M. Wilson, J. Merkow, H. Lan, A. L. Marsden, and S. C. Shadden, SimVascular: An open source pipeline for cardiovascular simulation, Ann. Biomed. Engrg. 45 (2017), 525–541.Show rawAMSref\bib{Updegrove2017}{article}{ author={Updegrove, A.}, author={Wilson, N. M.}, author={Merkow, J.}, author={Lan, H.}, author={Marsden, A. L.}, author={Shadden, S. C.}, title={Sim{V}ascular: {A}n {o}pen {s}ource {p}ipeline for {c}ardiovascular {s}imulation}, date={2017}, journal={Ann. Biomed. Engrg.}, volume={45}, pages={525\ndash 541}, } Close amsref.
[YFS13]
W. Yang, J. A. Feinstein, S. C. Shadden, I. E. Vignon-Clementel, and A. L. Marsden, Optimization of a Y-graft design for improved hepatic flow distribution in the Fontan circulation, J. Biomech. Engrg. 135 (2013), 011002.Show rawAMSref\bib{Yang2013}{article}{ author={Yang, W.}, author={Feinstein, J. A.}, author={Shadden, S. C.}, author={Vignon-Clementel, I. E.}, author={Marsden, A. L.}, title={Optimization of a {Y}-{g}raft {d}esign for {i}mproved {h}epatic {f}low {d}istribution in the {F}ontan {c}irculation}, date={2013}, journal={J. Biomech. Engrg.}, volume={135}, pages={011002}, } Close amsref.

Ju Liu is an assistant professor in the Department of Mechanics and Aerospace Engineering at Southern University of Science and Technology. His email address is liuj36@sustech.edu.cn.

Ingrid S. Lan is a graduate student in the Department of Bioengineering at Stanford University. Her email address is ingridl@stanford.edu.

Alison L. Marsden is a professor and Vera Moulton Wall Center faculty scholar in the Departments of Pediatrics (Cardiology) and Bioengineering and the Institute for Computational & Mathematical Engineering at Stanford University. Her email address is amarsden@stanford.edu.

Article DOI: 10.1090/noti2278

Credits

The opening image is courtesy of Ju Liu.

Figures 1, 4, 7, and 9–11 are courtesy of the authors.

Figure 2 is reprinted from Computer Methods in Applied Mechanics and Engineering, vol. 367, Ju Liu, Weiguang Yang, Melody Dong, Alison L. Marsden, “The nested block preconditioning technique for the incompressible Navier-Stokes equations with emphasis on hemodynamic simulations,” 2020, with permission from Elsevier.

Figure 3 is courtesy of John Wiley & Sons, Ltd., ©2019. Reprinted, with permission, from International Journal for Numerical Methods in Engineering, vol. 120, J. Liu, A. L. Marsden, and Z. Tao, “An energy-stable mixed formulation for isogeometric analysis of incompressible hyperelastodynamics,” 937–963. https://doi.org/10.1002/nme.6165.

Figure 5 is courtesy of Melody Dong.

Figure 6 is reprinted from Mechanics Research Communications, vol. 107, Ju Liu, Weiguang Yang, Ingrid S. Lan, Alison L. Marsden, “Fluid-structure interaction modeling of blood flow in the pulmonary arteries using the unified continuum and variational multiscale formulation,” 2020, with permission from Elsevier.

Figure 8A is reprinted from Computers & Fluids, vol. 142, Justin S. Tran, Daniele E. Schiavazzi, Abhay B. Ramachandra, Andrew M. Kahn, Alison L. Marsden, “Automated tuning for parameter identification and uncertainty quantification in multi-scale coronary simulations,” 2017, with permission from Elsevier.

Photo of Ju Liu is courtesy of Ju Liu.

Photo of Ingrid S. Lan is courtesy of Victor Tieu.

Photo of Alison L. Marsden is courtesy of Toni Bird.