Patient-Specific Inverse Modeling of In Vivo Cardiovascular Mechanics with Medical Image-Derived Kinematics as Input Data: Concepts, Methods, and Applications

1. Introduction

The primary role of numerical modeling in cardiovascular biomechanics has been to predict the performance of medical devices and to estimate physiological and mechanical cues acting on tissues, such as pressure and flow-driven stresses. Given the vast experimental evidence of mechanical factors producing effects on cellular differentiation, signaling, communication, and function [1,2,3,4,5], in silico experiments have explored the role of mechanical stimuli on normal and pathological tissue growth and remodeling [6]. From a clinical research standpoint, the development of patient-specific biomechanical models could provide more accurate and detailed data leading to a better understanding of the onset and progression of cardiovascular disease [7]. In addition, computational modeling has also been proposed as a supporting tool for medical practice on a patient-specific basis, which could provide non-invasive assessments of tissue properties, structure, and mechanical loads as physiologically meaningful risk stratification factors. Such patient-specific analyses have the potential to bring immense benefits to clinical practice by supporting diagnosis, treatment planning, and predictions of the outcome of surgical procedures with minimum associated costs and risk to the patients [8].

However, for biomechanical models to provide low-risk patient-specific solutions, personalized non-invasive clinical studies must be readily available to quantify regional cardiovascular function. Current medical imaging technology, namely echocardiography and magnetic resonance imaging (MRI), offers not only anatomical information but also high-resolution kinematics data of tissue motion and blood flow [9,10,11,12]. Kinematic-derived quantities, such as peak and average strain on the myocardium and aortic walls, have shown a good correlation with clinical risk markers [13]. Nevertheless, kinematic information alone cannot provide insights about mechanical forces, stresses, and tissue material properties, which are necessary for a full understanding of healthy and pathophysiological phenomena [14].

The inverse method, or data assimilation method, is an approach that allows solving classic mechanics problems “backwards”; that is, retrieving material properties and dynamic information (stress and forces) using measured kinematic information and loading boundary conditions as input [15]. Several research groups have coupled computational-mechanics tools with medical imaging technology to retrieve relevant biomechanical and hemodynamic markers from normal and pathological human tissues and organs [16], including diverse cardiovascular components [17]. Inverse modeling approaches have the potential to become a valuable tool for the non-invasive assessment of patient-specific cardiovascular health by providing quantitative physiological metrics that cannot be directly measured in vivo but may be derived entirely from clinical evaluations and the application of biomechanical principles.

The relevance of patient-specific modeling and its potential impact on the future of personalized and predictive health care has been acknowledged by several funding agencies. In 2003, the Interagency Modeling and Analysis Group was formed from the collaboration of nine institutes of the National Institutes of Health (NIH) and three directorates of the National Science Foundation (NSF). This group released its first funding opportunity in 2004 under the title “Interagency Opportunities in Multiscale Modeling in Biomedical, Biological, and Behavioral Systems Solicitation”, which has been regularly reissued since then, and led to the creation of the Multiscale Modeling Consortium which includes over 100 projects on multiscale modeling of biological systems. The European Union initiated the “Structuring the Europhysiome” Consortium in 2006, which led to the Virtual Physiological Human project, an ongoing initiative that aims to bring together policymakers, regulatory agencies, funding bodies, industry, and research organizations towards the development of integrated computer models of the mechanical, physical, and biochemical functions of the living human body [7]. These initiatives have motivated the integration of multidisciplinary teams of biologists, physicians, and engineers who are faced with the challenge of bringing together field-specific nomenclatures, techniques, and analytical approaches.

Inverse modeling of biomechanical systems requires the confluence of state-of-the-art techniques from several disciplines including clinical care, medical imaging, simulation engineering, data analysis, and computer science (Figure 1). Inverse methods have been developed on a highly specific basis and tailored to the available clinical data, tissue/pathology of interest, and available resources; thus, inverse modeling developments only share a general data processing pipeline, while differing on the clinical data source, imaging technique, and modeling approach. Therefore, the task of designing an inverse method pipeline requires a comprehensive understanding of the process at all stages, for which familiarity with fundamental concepts and terminology is a prerequisite. The latter can be challenging due to the multidisciplinary nature not only of the method itself but also of the clinically relevant phenomena to model.

Figure 1. Timeline of microprocessor speed as a measure of computation capability, and relevant landmarks on the fields of biomechanics theory, medical imaging and simulation that make possible modern patient-specific image-based inverse modeling of the cardiovascular system. Acronyms: CT, computerized tomography; DENSE, displacement encoding with stimulated echoes; FEM, finite element method; FVM, finite volume method; IE, inverse elastostatics; FSI, fluid–structure interactions; PC, phase contrast; MRI, magnetic resonance imaging.

Inverse modeling analyses can also be applied to in vitro experimental setups. The main advantage of this approach is that input data is not limited by the available clinical tests, mechanical loads can be precisely controlled, and kinematics can be measured with high-resolution instruments. Furthermore, the target of the inverse method can be defined not only in terms of kinematic information but also in terms of controlled mechanical loads and stress measurements [18]. Moreover, the outputs of inverse modeling can be validated with controlled experimental parameters. Inverse analyses of in vitro setups have been applied to explanted animal and human tissue, and to engineered tissue constructs. Notably, inverse modeling has been applied to resolve mechanics at a cellular level. The traction force microscopy (TFM) technique was introduced by Butler et al., to estimate the force that adherent cells exert on their surroundings by solving the traction field in a hydrogel of known properties, cultured with cells, and with embedded beads as visual markers [19]. By tracking the bead displacements through microscopy, and setting known boundary conditions, the traction field is resolved by an exact solution assuming a semi-infinite medium. Further development was introduced by Tambe et al., with the monolayer stress microscopy (MSM) technique, which allowed the inverse estimation of stress fields across monolayer cellular constructs under static and dynamic conditions by inducing controlled displacements of the boundary under a motorized microscope [20]. These, and other similar techniques, have been used to explore the response of cardiovascular cells (endothelial cells, cardiomyocytes, smooth muscle cells, etc.), cellular layers, and engineered tissue to mechanical stimuli in terms of cell proliferation, migration, expression, and synthesis of extracellular matrix components [21,22,23]. The detailed and accurate results that can be retrieved from inverse analyses of in vitro experiments can provide valuable insights into cardiovascular mechanobiology. These insights contribute to the understanding of how macroscopic biomechanical factors affect the healthy or pathological growth and remodeling of cardiovascular and engineered tissues. However, the replication of in vivo physiological conditions in vitro is cumbersome, and the results of in vitro experimentation can be challenging to extrapolate to patient-specific situations. As a result, the clinical application of inverse analyses of in vitro experiments remains limited.

This article aims first to serve as a referential document for concepts and methods from all involved disciplines on patient-specific in vivo inverse modeling; and secondly, to highlight the potential clinical application of patient-specific inverse modeling in the cardiovascular research field. Specifically, we review the fundamentals of cardiovascular tissue and blood biomechanics, modeling and simulation, and medical imaging, as they relate to the inverse modeling approach and its applications in cardiovascular medicine. In Section 2, we review the application of the principles of classical continuum mechanics to the study of blood and tissue motion, with special emphasis on the constitutive equations that have been proposed to describe the mechanical behavior of cardiovascular tissue and blood. Section 3 briefly summarizes the fundamentals and main features of the finite element method (FEM) and finite volume method (FVM), as the most popular formulations for the numerical solution of biomechanical models. Next, we review the general definition of inverse problems and the available alternatives to solve inverse mechanics problems in Section 4. In Section 5, we review the working principles and main features of ultrasound (US), magnetic resonance imaging (MRI), and computational tomography (CT) imaging, giving special attention to the available techniques for the resolution of tissue and blood kinematics. Finally, in Section 6, we present a comprehensive review of the applications of imaging-based inverse modeling approaches to patient-specific human cardiovascular mechanics, including the resolution of the unloaded configuration and the estimation of tissue properties and stresses. Reviewed applications include healthy and diseased heart valves, cardiac and arterial walls, and hemodynamics of large arteries. To highlight the potential application of the inverse-modeling approach in cardiovascular medicine, we focus herein mostly on developments made in human studies, with a few mentions of relevant and pioneering studies in animals.

2. Governing Principles of Biomechanics

Modern biomechanics consists of the formulation of governing equations describing balances of mass, linear and angular momentum, and energy to biological systems and physiological processes. The human body maintains a uniform and stable temperature through homeostatic thermoregulation. Thus, contributions due to temperature change in the internal energy of the material, heat fluxes, and heat supply are typically negligible to the energy balance, which is in turn reduced to the balance between deformational energy and stress power (thermodynamic work). In classical continuum mechanics of purely mechanical processes, the balance of angular momentum directly translates to the symmetry of the stress tensor, and therefore, the relevant governing equations for most cardiovascular mechanics applications consist only of the balances of linear momentum and mass. However, most biological systems are open, and continuously interact with their surroundings, and thus the conservation principles must be handled carefully, especially with respect to tissue growth and atrophy within relevant timescales.

Given that the resolution of most in vivo medical imaging is on the scale of millimeters, only phenomena occurring at the tissue level can be directly associated with these measurements. The assumption of material continuity is reasonable for the formulation of the governing principles at this scale, leaving any additional considerations dealing with the extracellular and intracellular micro-environments to be included ad hoc with additional modeling formulations and constitutive equations.

To apply these principles, it is necessary to relate the stress tensor to kinematic measures, which is in essence the description of the mechanical behavior of the material under study. This information is provided by a constitutive equation; these can be either phenomenological equations “arbitrarily” formulated to reproduce experimental observations, or analytical expressions inspired by theoretical interactions of the material constituents at the micro or molecular scale. The selection of adequate models to describe the phenomena of interest is key to the success of any engineering analysis. The selected model must be complex enough to describe the most salient observable features at the scale of interest, while ideally being simple enough to provide a rational interpretation of its parameters and results and render a computationally tractable numerical problem. After fitting these parametrized models to experimental data, the constitutive equation can provide an additional understanding of the underlying mechanisms associated with the mechanical response of the material. Models of increased complexity usually require a larger number of parameters to be fitted, and overparametrized models can lead to solution multiplicity which obscures its interpretation and validity. For the sake of generality, the constitutive equation must also be independent of the frame of reference, comply with the second principle of thermodynamics, and yield amenable mathematical treatment and systems of equations that are solvable [24].

After formulating the constitutive relation as a function of specific unknowns (e.g., displacement or velocity fields), the resulting system of governing equations is then particularized into specific problems by the definition of temporal and spatial domains of interest and the imposition of appropriate boundary and initial conditions. To obtain a unique solution, it is necessary to constrain the problem by assigning a first-order boundary condition on at least one of the boundaries (e.g., by prescribing a known displacement or velocity). On the rest of the boundaries, higher-order boundary conditions can be applied to impose distributed forces such as a known pressure. Once solved, the result of this forward formulation is the transient spatial distribution of displacement or velocity throughout the domain as a result of the specified loads and material properties. From this kinematic information, strain and strain-rate distributions can be numerically derived, and the stress distributions retrieved via the constitutive equation.

In the following subsections, we present a review of the main features of the mechanical behavior and composition of cardiovascular tissues and blood, as well as the constitutive equations that have been developed and applied to model such behaviors. Then, we review the algorithms implemented to model fluid–structure interactions and their importance in cardiovascular mechanics simulations. Finally, we briefly discuss the modeling of biological tissue growth and remodeling by the application of the constrained mixture theory.

2.1. Structural Mechanics of Cardiovascular Tissue

The study of soft biological tissues under the framework of finite elasticity was initiated by Y.C. Fung and others in the late 1960s, setting the basis of modern biomechanics [25,26]. As in classical solid mechanics, the mechanical analyses of cardiovascular tissues are usually performed with a Lagrangian formulation of the governing principles (Figure 2). Applied forces are imposed as boundary conditions. For blood vessels, these are prescribed as transmural pressure differences that often assume a traction-free condition on the adventitial surface. More recently, however, growing attention to the role of perivascular and pericardial support and tethering has promoted the inclusion of restrictions to the displacement of the outer surface of the heart and vasculature [14,27,28]. In addition, more complex formulations of cardiovascular tissue mechanics which departs from classical elastic solids have been proposed to account for complex microstructural compositions, the inclusion of pre-stress/strains, chemically activated muscular tone, and viscous energy dissipation.

Figure 2. Representation of the modeling process of structural mechanics of cardiovascular tissue and fluid mechanics of the blood flow with a continuum mechanics approach. Structural mechanics of cardiovascular tissue are usually analyzed with a Lagrangian formulation that follows the deformation of a given portion of the tissue. Blood flow mechanics is usually analyzed with an Eulerian formulation, that is, analyzing the mass and energy balances on a fixed volume of interest through which the fluid flows.

Cardiovascular tissues are comprised of multiple layers of cells and extracellular matrix (ECM) components. The ECM is a network of macromolecules that is continuously synthesized and degraded by active cells and functionally provides them with structural and biochemical support. Typically, collagen, elastin, and fibrillin are regarded as the main structural constituents responsible for the macroscopic mechanical behavior of cardiovascular tissues [29]. Healthy cardiovascular tissue retains residual stress even when unloaded (i.e., the tissue is pre-strained relative to a reference state of zero transmural pressure). Circumferential and axial pre-stress/strain in vascular conduits have been widely established by measuring how much these tissues recoil to an open configuration when excised and cut transversely and longitudinally to relieve the residual stress [30,31]. It has been hypothesized that pre-strain plays a relevant role in balancing higher stresses on the luminal surface of blood vessels and promoting a homogenized transmural stress distribution and homeostatic equilibrium of the vascular tissue [32]. Notably, residual strain is heterogeneous and has been shown to vary with patient age and health, likely as a consequence of heterogeneous growth and remodeling and/or damage.

Additionally, cardiovascular tissue is muscular in nature and actively contracts/distends. Thus, its mechanical behavior is affected by the activation of actin-myosin sliding filaments, which depends on ion-based chemical signaling and determines the muscular tone. In the myocardium, striated muscle activation is responsible for cardiac contraction. In large arteries, contraction of smooth muscle cells regulates downstream vascular resistance, blood flow, and propagation of the pressure-pulse wave along the vascular tree.

Of note, cardiovascular tissue also exhibits viscoelastic behavior, which has been established with stress relaxation, creep, and strain-rate experiments. It has been argued that viscous energy dissipation of healthy tissue, functioning at a regular physiological rate (~1 Hz), is negligible compared to stored strain energy [33]. Nevertheless, viscoelasticity may play a critical role under pathological conditions where the deformation rate is increased, such as in atrial fibrillation, or when dealing with highly dissipative structures such as lipid pools in atherosclerotic plaques [34,35]. Despite the relevance of viscoelastic properties to pathological conditions, standardized testing protocols have yet to be developed for the exploration of its relation to disease onset and progression [36]. Notably, if viscous dissipation and inertial effects are neglected, all temporal terms in the governing equations are canceled, rendering the problem a quasi-static process (which is the most common approach applied to vascular wall mechanics).

2.1.1. Constitutive Equations

Passive Properties

Linearized elasticity falls short in describing the complex behavior of biological tissue, not only because the mechanical response is highly non-linear but also the material undergoes finite motions and deformations. In the case of non-linear behavior, it is usually convenient to employ the formulation of hyperelasticity and express the constitutive equations as the relation of a scalar stored energy density function to the deformation gradient tensor or the strain tensor (while other definitions of stretch or strain tensors are also possible and common). The scalar energy density function represents the amount of deformational energy stored per unit volume and is defined in such a way that the stress tensors can be obtained from their derivatives with respect to the strain or stretch tensor.

The passive behavior of cardiovascular tissue is characterized by an increasing resistance to deformation with strain. This behavior is represented by an increasing slope in the stress versus strain/stretch and strain energy versus strain/stretch curves, with a close to zero slope at zero strain and rapidly increasing at physiological ranges. This behavior has been attributed to the structural characteristics of the ECM components. It has been suggested that the increasing resistance to deformation with strain is owed to the progressive engagement of wavy bundles of elastin and collagen fibers to support the mechanical loads (Figure 3).

Figure 3. Representative biaxial stress-stretch behavior of healthy cardiovascular tissue. The stress and slope increase with the stretch/strain in any direction. This suggests that the cardiovascular tissue stiffens with stretch/strain which is hypothesized to be a consequence of the progressive engagement of ECM components to resist further deformation. This behavior is modeled by exponential functions of the deformation tensor components and/or invariants. Mark symbols in the figures show experimental biaxial test data of the human thoracic aorta for (a,b) young patients (20 to 35 years of age), and for (c,d) older patients (57 to 71 years of age). Solid lines represent the best-fit approximation with a four-fiber family constitutive equation. Reprinted/adapted with permission from Ref. [37], 2014, Elsevier.

Hyperelastic isotropic material models, such as the Neo–Hookean and the Mooney–Rivlin constitutive equations, can accurately describe the behavior of amorphous bodies such as lipid pools in atherosclerotic formations, and fit some portions of the pressure–volume relation of blood vessels. The symmetry of these material models allows its representation as to the linear combination of the deformation gradient invariants weighted by the material properties. This formulation allows the determination of unique material properties for a given mechanical behavior. However, these models fail to reproduce the characteristic highly nonlinear and anisotropic behavior of cardiovascular tissue (Figure 2 and Figure 3). This was originally addressed by the use of phenomenological equations, e.g., the Fung orthotropic exponential model being one of the most commonly used. Its success relies on its relative simplicity, widespread numerical implementation, and accuracy in the prediction of stress-strain curves [24,26]. Guccione et al. proposed modifications to Fung’s orthotropic model based on myofiber structure and orientation to tailor myocardial tissue behavior [38]. These phenomenological equations usually consist of two terms contributing to the strain energy function. First, the contribution of volume changes is written as a function of the determinant (third invariant) of the deformation gradient tensor. The second term is the deviatoric contribution to the strain energy function, defined to be proportional to an exponential function of the components of the strain tensor. The proportional constant sets the scale of the material stiffness, and the function of the strain tensor components defines the material anisotropy. Part of the success of Fung-like equations in describing the stiffening of cardiovascular tissue with strain relies on the exponential functional form to quantify the effect of strain increments on deformational energy. However, the fitted constants of phenomenological equations lack physical interpretation which is desirable for studies aiming to relate material properties with pathological conditions [39,40].

In the past decades, many microstructure-inspired constitutive models have been proposed to specifically suit cardiovascular tissue behavior [39,41]. Fiber-family models have been particularly successful in reproducing the anisotropic behavior of the vascular wall while keeping physiological meaning to some of the fitting constants, the Holzapfel–Ogden model and its many variants being the most popular for cardiovascular tissue. These models assume that families of 1D fibers, each with specific mechanical behavior, orientation distribution, and volume fraction, are embedded within an isotropic continuum matrix. The isotropic component of the strain energy function is usually defined as a function proportional to the first invariant of the deformation gradient tensor. The contribution of fiber families is a weighted sum of exponential functions, the weighting factors being the fiber family material parameters. The exponential functions are defined depending on deformation tensor invariants and the relative orientation of fibers to strain, such that the contribution of a fiber family is maximized if the strain deformation occurs in the direction of the fibers [37]. Again, exponential functions are employed to mimic the stiffening effect of strain on cardiovascular tissue (Figure 3), this time being directly attributed to the ECM fiber components. Many improvements to these models have been proposed to account for different coupling effects such as inextensibility of fibers or cross-linking of the fiber ensembles [42,43,44,45].

Active Properties

Adequate modeling of active contraction is key for an accurate description of cardiovascular function in general, being particularly critical for modeling the heart. Conceptually, there are two possible approaches, the active stress models are the most common and they assume the stress tensor can be decomposed as the sum of a passive and active component, while active strain models assume a product decomposition of the deformation gradient.

Most active contraction models used in the inverse analysis of ventricular mechanics through continuum mechanics are simplifications of more complex bio-chemo-mechanical models such as the work of Hunter et al. [46]. The latter proposes a four-state variable model which includes the passive elasticity of myocardial tissue, the binding of calcium ions (Ca2+) to troponin C and its release, tropomyosin movement kinetics, the myofiber length, and the kinetics of cross-bridge tension build-up under perturbation of myofilament length. In practice, the detailed information required for the evaluation of this model is out of reach, so simplified models assume the active stress acts mostly lengthways the direction of myofibers, with a magnitude that is proportional to the fiber length and the activation status. The activation status is often expressed as a time-dependent spatially heterogeneous function ranging from 0 to 1 [47]. Active strain function will impose the relative shortening lengthways of myofiber directions as a function of location and time along the cardiac cycle.

Typically, the activation state is assumed to be instantaneously homogeneous within the region of study, however, it is known that cellular activation propagates as an electrical wave and the excitation-contraction coupling poses a complex electromechanical problem [48]. This propagation has been modeled macroscopically as a reaction-diffusion problem of the electric potential through the intracellular and extracellular domains, thus known as the bidomain model. The monodomain model is a simplification that assumes the same propagation anisotropy for both domains. Some interesting research has been developed to apply inverse modeling to fit the parameters of mono and bidomain equations using patient-specific electrocardiography [49,50]. These, however, fall out of the scope of this review, and the interested reader is encouraged to study the abundant literature on the inverse problem of electrocardiography [51,52].

The application of active contraction models requires the specification of the local myofiber orientations. Patient-specific myofibers orientation can be resolved via diffusion tensor MR imaging (DT MRI) [53]. Being a relatively novel technique, these scans are rarely available from medical records of cardiovascular-disease patients, although their relevance in the biomedical field is a growing topic of discussion [54]. Therefore, myofiber orientation is either assumed a priori with simplified models or obtained through diffeomorphic transformations with the employment of precomputed cardiac atlases [55]. Bayer et al. proposed a Laplace–Dirichlet rule-based algorithm for assigning myofiber orientation to computational heart models that showed good agreement with DT MRI measurements. This algorithm consists of the resolution of the Laplace equation on the simulation domain with appropriate Dirichlet boundary conditions constrained by the following rules: the longitudinal fiber direction is parallel to the endocardial and epicardial surfaces, the longitudinal fibers rotate clockwise throughout the ventricular wall from a positive helical angle at the endocardium to a negative helical angle at the epicardium (both imposed by the user), fibers in the papillary muscles and trabeculae are assumed parallel to the long axis of these structures, the transverse fiber direction is perpendicular to longitudinal fibers, and fiber orientation in the septum is continuous with the ventricular walls [55]. Similarly, Potse et al. proposed a rule-based algorithm to define myofiber orientations assuming that longitudinal fibers are orthogonal to the local vector pointing to the shortest path between endocardium and pericardium, with a clockwise varying helical angle [56]. Rijcken et al. derived an equation for longitudinal and transverse myofiber orientation by solving an optimization problem, which maximized the ejection while maintaining fiber strain as homogeneous as possible on idealized geometries [57].

2.2. Fluid Mechanics of Blood Flow

For the study of fluids, it is more practical to implement an Eulerian formulation of the governing equations. This formulation is obtained by applying the Reynolds transport theorem to the equations for mass and momentum balance. Thus, this formulation solves the relation between flow driving forces, flow velocity, and deformation rates (Figure 2). A simplifying assumption applicable to biological systems is the incompressibility of the fluids, as most of them are either liquids or gases moving at subsonic velocities. Additionally, it is convenient to decompose the stress tensor into a spherical tensor representing the hydrostatic pressure and a deviatoric stress tensor. With this decomposition, constitutive equations can be designed to specifically relate the deviatoric stress components to the viscous dissipation of momentum.

Blood flow is generally assumed to be laminar throughout the circulatory system. The main arguments for this assumption are the pulsatile nature of the flow, the reduced dimensions of the vessels, and relatively low velocities, each contributing to the viscous effects overcoming the inertial forces and preventing turbulent random motion. However, it has been argued that transition to turbulent flows could be achieved locally in stenotic arteries. The use of a laminar model to study those cases could lead to an underestimation of wall shear stress, and stress oscillation [58,59]. Unlike most conventional engineering flows, blood flow is pulsatile and contained by compliant conduits of complex geometry. Since the 1950s, Womersley [60], McDonald [61], Taylor [62], Pedley [63], and others, developed analytical and experimental studies of pulsatile flow in mammals, identifying the most relevant parameters and features of this type of flow, thus setting the bases for modern hemodynamics.

Besides the intricacies of pulsatile flow in distensible conduits, the blood itself is a complex fluid. Blood consists of a suspension of cells in an aqueous solution of proteins and minerals called plasma. Plasma occupies approximately 55% of the blood volume, the rest being mainly occupied by red blood cells, white blood cells, and platelets. The rheological behavior of blood depends on how its constituents interact with each other and with the vessel walls, in consequence, this behavior is non-linear and highly dependent on the volumetric composition of blood, the flow conditions, and vessel dimensions. Modeling the complex interactions of blood constituents is a challenging statistical mechanics problem [64]. Some researchers have shown that cell aggregation and disaggregation are relevant to accurately describing blood rheology, especially in capillary flows where the cell size is comparable to the vessel diameter. Multiscale approaches have been successful in coupling the behavior of single cells as elastic entities with the transport equations of fluid flow, which are relevant for the study of clotting, aggregation, and platelet activation [65,66]. These approaches are computationally expensive, making them impractical for the study of large vessels.

Constitutive Equations

In the study of large and medium vessels (the ones that can be feasibly resolved with standard medical imaging), blood is often assumed to be a single-phase continuum. This approximation is reasonable given the relatively small size of cell aggregates compared to the vessel dimensions (and thickness of the boundary layer), and the relative relevance of inertia on flow motion [67]. For these cases, phenomenological constitutive equations describing the macroscopic behavior of flow are often applied. The linear Newtonian fluid is the simplest and most commonly employed model, providing reasonable results in vessels with diameters down to 200 μm [67]. Constitutive equations, such as the Casson, Herschel–Bulkley, and Carreu–Yasuda, incorporate the shear-thinning effect on apparent viscosity by introducing a yield shear stress term [68]. To account for the effect of the volumetric share of cell suspension, recent works have included the hematocrit as an independent variable for the estimation of the effective viscosity [69].

2.3. Fluid-Structure Interactions (FSI)

Mechanics of the vascular wall and hemodynamics have been mostly studied as isolated problems; however, the function of the cardiovascular system is the result of complex interactions between blood, the actively contractile cardiac tissue, and the compliant vascular walls. The interaction of fluids and solids can conceptually be achieved by coupling the boundary conditions on the interface between the solid and the fluid, such that the field of displacements, velocities, and stresses are continuous and derivable at all points in a monolithic fully coupled approach. This, however, poses many implementation difficulties for complex 3D domains that can only be solved numerically. In addition, the typically large deformations of the cardiovascular walls cannot be handled by linearized methods used in conventional engineering applications.

The immersed boundary method, introduced by Peskin, was originally developed for the study of flow around heart valves and was rapidly adopted for many other applications [70]. In this approach, the Eulerian variables of fluid dynamics describing the surrounding flow are defined on a fixed computational grid, while the Lagrangian variables, accounting for the deformation of the tissue structures, are defined in a curvilinear computational grid that can be displaced with no conforming constraints in respect to the Eulerian grid. The moving solid boundary interacts with the fixed fluid domain by means of elastic body forces which are modulated by Dirac delta-like functions [71,72]. The fictitious domain method is a generalization of the immersed boundary method, which solves the coupling of the Lagrangian and Eulerian domains by the use of Lagrange multipliers instead of the concept of body forces [73]. This method is computationally less demanding as it does not require fitting the interface boundary at the cost of impaired accuracy near the interface.

Similarly, Figueroa et al. proposed the coupled momentum method, consisting of changing the non-slip condition on the fluid boundary to a traction condition, which is strongly coupled to the degrees of freedom of modified thin-membrane elements. This allows the formulation of the solid equations on the same Eulerian frame as in the fluid equations. In consequence, the fluid–solid interface mesh remains fixed, while the boundary nodes will have nonzero velocities [74]. Another well-established method for FSI simulation is the arbitrary Lagrangian–Eulerian (ALE) algorithm, which allows the arbitrary convective motion of the computational nodes of the discretization grid with respect to a fixed reference frame. Typically, the nodes on the fluid-solid interface are treated with a Lagrangian formulation. To deal with large or heterogeneous deformations of the interface, several implementations include the re-discretization of the computational domain to avoid the influence of ill-shaped deformed elements. The drawbacks of this method are the computational expense of re-meshing the domain, and the induced inaccuracies by transferring solutions from the degenerated mesh to the new one [75,76].

2.4. Growth and Remodeling Models by the Constrained Mixture Theory

One of the most relevant characteristics of living tissue is its capability to adapt in response to chemical and mechanical stimuli. This adaptation comes with microstructural reconfigurations, which alter the mass composition and the resulting contributions and properties of the tissue constituents. The understanding of the effect of mechanical stimulation on normal and pathological growth and remodeling of soft tissues is an active field of study that bridges biomechanics and mechanobiology.

In 1994, Rodriguez et al. proposed a general continuum formulation for the finite volumetric growth modulated by mechanical stress [77]. The theory of adaptation of living tissues was further developed by Humphrey and Rajagopal who proposed the constrained mixture theory, a mathematical framework to predict not only the growth but also the remodeling of biological tissues under transient mechanical and chemical stimulation [78,79]. The constrained mixture theory is based on the continuum theory of mixtures; that is, each component complies with a modified version of the governing principles of motion in the Eulerian formulation. The modification involves the addition of mass source/sink terms that account for the rate of synthesis or degradation of the constituent in its respective mass balance equation and the component-to-component interaction forces in the momentum balance equation. These source/sink terms respond to a series of constraints of physical and chemical nature and are dependent on the local distribution of strains, stress, and current composition. Constitutive equations must be defined for each constituent, and the overall properties of the construct can be calculated as a combination of its constituents, where simplified linearized forms weighted by their volume fraction are typically chosen [32,80].

2.5. Summary

In Table 1 we summarize the highlights of the governing principles of cardiovascular biomechanics through a continuum mechanics approach.

Table 1. Summary of governing principles of biomechanics.

3. Numerical Methods

The above-described system of governing and constitutive equations can only be solved analytically for a reduced group of oversimplified cases. Thus, mechanical analyses of complex biological systems require the application of numerical methods to obtain approximate solutions. It has been claimed that the development of numerical methods was key to the foundations of modern biomechanics [81,82]. Many of the early simulation analyses of the cardiovascular system and components were developed with in-house codes, but the popularization of commercial software boosted the production of computational research in biomechanics [82]. More recently, open-source specialized software for numerical biomechanics, such as SimVascular and FEBio [83,84], have risen from the collaborative effort of academic groups aiming to incorporate relevant bio-chemo-mechanical models of biological systems into simulation pipelines. Many different options exist for the numerical solution of time-dependent 3D problems. Mesh-based methods are the most popular approaches, particularly the finite volume and finite element methods which are reviewed in the following pages.

3.1. Finite Volume Method

The finite volume method (FVM) is conceptually straightforward. The domain of study is discretized in a series of non-overlapping finite volumes, and the governing equations, usually expressed in Eulerian formulation, are converted into algebraic expressions by integrating them over each discrete volume. The balance equations are applied on a node located in the center of the finite volume, while the flux terms are calculated at its faces (Figure 4a). This allows first and second-order approximations of derivatives. The surface flow for a given shared face is set identical and in opposite direction for the adjacent discrete volumes, and equal to a boundary condition at the edge of the domain. By doing so, the balance equations are held at the whole domain and within each finite volume, which is one of the most attractive features of the FVM. Additionally, since the calculation of properties happens in the center of each volume, it is relatively easy to implement boundary conditions of a higher order [85]. Numerical implementation of this method is also straightforward in the case of structured meshes, becoming more complex for unstructured meshes due to the bookkeeping necessary for the calculations of interface flux balances.

Figure 4. Diagram of finite volume (FVM) and finite element methods (FEM) approximation principles. (a) In FVM, the domain is discretized in finite volumes, and balance equations are solved at the center of each volume. (b) In FEM, the domain is discretized in finite elements, and the variables distribution is assumed to follow a prescribed shape function within each finite element.

The use of FVM for the solution of convection-diffusion problems was first introduced in the early 1960s by Tikhonov and Samarskii [86]. Since then, FVM has been particularly successful in its application to computational fluid dynamics, as many of the current commercial computational fluid dynamics (CFD) software suites are based on this method. Biomechanical applications of this method mostly focus on hemodynamic and tracheobronchial airway simulations. However, this method can be applied to other boundary value problems such as electromagnetics and structural mechanics [87,88].

3.2. Finite Element Method (FEM)

The finite element method (FEM) consists of the discretization of the domain of study on simple geometrical elements (or finite elements), where the unknown fields are discretized as linear combinations of shape functions of any order, linear and quadratic being the most common. The shape functions are typically defined at each element depending on local and normalized coordinates (Figure 4b). The local governing equations for each element are then assembled and organized in a matricial system of algebraic equations. Finally, the solution is approximated by minimizing the weighted error associated with each element. Several weighting rules have been proposed, the Galerkin method and its variations being the most widely used [89]. By converging to the solution through the minimization of an error function and not through the exact solution of balance equations, FEM is said to be formulated in a “weak” form. However, the weak form is equivalent to the exact solution in the limit of refining the domain discretization. In fact, it has been widely shown that mesh-independent FEM solutions do not show any practical difference from the output of more conservative numerical methods such as FVM [90,91].

FEM was developed in the early 1950s to perform structural analysis for the aerospace industry and was soon applied to study the biomechanics of musculoskeletal and cardiovascular tissue [81,82]. As early as 1968, FEM was used to study the non-linear viscoelastic behavior of arteriole tissue [92]. This technique has been traditionally used for the solution of solid mechanics problems; however, it has also been used to solve the governing equations of other physical phenomena, including fluid mechanics [81]. Regarding the convenience of relying on a single solver engine, many multiphysics simulation software suites have introduced FEM formulations for fluid mechanics, [93] which also facilitates the implementation of FSI simulations [94].

3.3. Summary

In Table 2 we summarize the highlights of the principles and formulations of the numerical methods typically applied on cardiovascular biomechanics research.

Table 2. Summary of numerical methods.

4. Inverse Problems

Modeling physical phenomena can be thought of as a mapping operation, where a set of inputs (B) is transformed into a set of outputs (E) by applying the model operator (M) such that M(B)=E. In the realm of physics, there must be a cause–effect relation between the inputs and outputs, and forward modeling consists of designing and applying a mapping function capable of producing outputs that closely follow experimental measurements (e), meaning that the difference E−e should be close to zero (Figure 5a) [95].