Modeling left ventricular dynamics with characteristic deformation modes

A computationally efficient method is described for simulating the dynamics of the left ventricle (LV) in three dimensions. LV motion is represented as a combination of a limited number of deformation modes, chosen to represent observed cardiac motions while conserving volume in the LV wall. The contribution of each mode to wall motion is determined by a corresponding time-dependent deformation variable. The principle of virtual work is applied to these deformation variables, yielding a system of ordinary differential equations for LV dynamics, including effects of muscle fiber orientations, active and passive stresses, and surface tractions. Passive stress is governed by a transversely isotropic elastic model. Active stress acts in the fiber direction and incorporates length–tension and force–velocity properties of cardiac muscle. Preload and afterload are represented by lumped vascular models. The variational equations and their numerical solutions are verified by comparison to analytic solutions of the strong form equations. Deformation modes are constructed using Fourier series with an arbitrary number of terms. Greater numbers of deformation modes increase deformable model resolution but at increased computational cost. Simulations of normal LV motion throughout the cardiac cycle are presented using models with 8, 23, or 46 deformation modes. Aggregate quantities that describe LV function vary little as the number of deformation modes is increased. Spatial distributions of stress and strain change as more deformation modes are included, but overall patterns are conserved. This approach yields three-dimensional simulations of the cardiac cycle on a clinically relevant time-scale. Electronic supplementary material The online version of this article (10.1007/s10237-019-01168-8) contains supplementary material, which is available to authorized users.


S2 Prolate spheroidal bicubic splines
In this section we define bicubic splines in the prolate spheroidal coordinate system, providing a general framework for constructing functions in these coordinates. We define the reference endocardial and epicardial surfaces using these splines, illustrated in Figure S1. Although not shown here, these functions may also be used to construct arbitrary deformation modes within the general kinematic framework. While we have termed our construction "bicubic splines," we do not  We define bicubic functions on the ν × φ plane on the domain φ ∈ [0, 2π) and ν ∈ [ν min , π].
(S1) Bicubic splines are constructed by partitioning the domain into N ν sections in the ν direction and N φ sections in the φ direction. We denote the full bicubic spline function f (ν, φ ) which is composed of N ν ·N φ local bicubic functions f i, j (ν, φ ) indexed by i = 1, . . . , N ν and j = 1, . . . , N φ . There are a total of N ν + 1 nodes in the ν direction (although the last node is degenerate) and N φ nodes in the φ direction (due to the periodicity). The local function domains are illustrated in Figure  S2. We construct the local bicubic functions on rescaled coordinates u ν , u φ ∈ [0, 1] defined by where ∆ ν i and ∆ φ are the widths of each local spline domain. The φ spacing is fixed to The ν spacing is variable to allow for more consistent spacing in Cartesian coordinates, as illustrated in Figure S2. The bicubic functions are defined locally as The full spline function f (ν, φ ) is evaluated by identifying the correct region of the domain and then computing f from (S4). The nodes at the apex ν = π are degenerate. We can directly specify the value of f at these nodes, but not the derivatives. It is possible, however, to define nonaxisymmetric derivative variations in the spline function at ν = π while preserving continuity of the derivatives through the apex. We define a parameter vector wherex andỹ arẽ Writing these conditions in terms of prolate-spheroidal coordinates and evaluating at the apex implies To maintain continuity of the derivatives through the apex the local bicubics must also satisfy The free values that determine the shape of the prolate spline function f are the function values, derivatives, and first cross derivative at the nodes between spline regions. In other words, at each node we assign We define the full set of free parameters This implies a total of N s = 3 + 4N ν N φ free parameters that define f . The values of the function and derivatives at the degenerate nodes are fully determined by the first three values of s s s together with the conditions (S7) and (S8).
The local bicubic functions f i, j defined in (S4) have 16 coefficients p m,n . Each local function is bordered by 4 nodes where the values of are assigned by the input parameter vector s. The 16 coefficients p m,n are exactly determined by the 16 conditions at each node and may be computed by inverting the matrix that prescribes these conditions.

S3 Incompressibility condition solution
Equation (9) is completed with f c , the values of which are determined by displacements of the endocardial wall µ in through The mapping equation (9) is an implicit cubic equation of cosh µ. We write the solution in terms of convenience variables that provide a computationally efficient code implementation: The mapping equation (9) is equivalent to which has the solution The solution to the mapping equation is therefore

S4 Apex cut-off function
We define a simple C 1 function for sending functions to zero near the apex: where ν T is a fixed transition point near the apex.

S5 Deformation function derivatives
Computation of the deformation gradient tensor (5) requires derivatives of the displacement functions. Here, we give the analytic form of these derivatives using convenience variables to simplify the expressions. We define convenience variables The partial derivatives of these functions are The partial derivatives of the mapping equation (16) yield the derivatives of µ:

S6 Tagged MRI data
We used tagged cardiac MRI data (illustrated in Figure S3) to evaluate the kinematic model. These data were published previously by Kar et al. (2014) and are used here with permission. As noted in that work, subjects signed informed consents in accordance with the university's Institutional Review Board (IRB) guidelines.

Fig. S3
A short-axis tagged cardiac MR imaging plane. The left panel shows the initial tagging near end-diastole, while the right panel illustrates myocardial deformations recorded at end-systole.
The myocardial displacements were registered from the tagged MRI data using a deformable image registration algorithm. The myocardial walls were outlined manually in each imaging plane before the image registration was performed. We described this image registration method previously in Hong (2018). This method registers LV displacements in each plane using a deformable mesh. Although tag fading exists, the registration of consecutive frames is largely unaffected as the intensity between consecutive images is reasonably consistent, despite substantial fading over the full time-course. These 2D registrations are combined to generate a 3D representation of the myocardial deformation using a cubic mesh. These registered displacements are illustrated in Figure S4.
Displacement vectors between end-systole and end-diastole at the endocardial wall registered from the tagged MRI. Vertical motion is registered using long-axis views while the torsion is registered with short-axis views. Greater expansion is seen at the lateral wall compared to the septal wall.

S7 Fiber direction basis vectors
In this section, we develop a method for defining fiber angles that vary linearly from the endocardium to the epicardium. The angles at the surfaces are typically chosen to match average values from the literature. We define the muscle fiber angles in the reference frame. The fiber angles are defined relative to continuous surfaces within the myocardial wall. These surfaces are defined in terms of an auxiliary variable which has been constructed so that it is zero at the inner wall and one at the outer wall. The µ 0 coordinate is computed from u as Internal surfaces are defined by constant values of u. The Cartesian coordinates of the internal surfaces are Surface tangent vectors are found by computing the cross product of the ν 0 and φ 0 derivatives: Because u is a constant on these surfaces, the derivatives of µ 0 are where the partial derivatives of µ in0 and µ out0 are the exact derivatives of the bicubic functions. The local basis vectors are computed by dividing by the Euclidean norm: Note that at the apex (ν = π) |w 3 | = 0 and the expression above is singular. A suitable alternative expression at the apex is We construct a vector e 4 that lies in the x-y and e 2 -e 3 planes. Rodrigues' formula gives a vector e 4 that has been rotated an angle θ about the vector e 1 : This condition has the form The solution is The resulting vector e 4 lies in the x-y plane and is normal to the internal constant u surface. We define fiber coordinates and fiber angles relative to this direction. We define fiber angles that vary linearly through the wall where w increases from zero at the endocardium to one at the epicardium. A simple choice for w is Alternatively, w may be scaled according to the distance along lines of constant ν 0 and φ 0 . The scale factor in the µ 0 direction is Integrating the scale factor along these lines gives the distance: We approximate this integral by the trapezoid rule because it lacks an analytical solution. Thus, We use the cutoff function (S15) defined in the appendix to send the fiber angles to 90 • at the apex. Thus the final form of ψ is The local fiber directions are defined in terms of the local fiber coordinates(s, n, f ). The fiber direction is defined by rotating the reference vector e 4 by an angle ψ in the surface tangent plane. The fiber direction is e f = e 4 cos ψ + (e 1 × e 4 ) sin ψ + e 1 (e 1 · e 4 )(1 − cos ψ).

(S45)
The s coordinate direction is the surface normal vector, i.e., To generate a right handed orthogonal coordinate system, the final fiber coordinate direction is

S9 Closing surface Γ
The LV chamber must be closed to define the volume of the cavity and approximate work done at the base of the LV. However, an appropriate lid is difficult to define in prolate spheroidal coordinates. As a result, we define the lid in cylindrical coordinates.

S9.1 Cavity volume
The cavity volume is computed with three terms. This is necessary because of singularities associated with the prolate coordinate system at the axis.
Another volume is computed between the boundary of the first volume at the base and a second surface that has a fixed ν value at the axis. For simplicity we set the value of ν at the axis to ν up0 , the largest value of ν up in the reference configuration. This choice is unimportant because it does not define the final closing surface. This surface is defined in the prolate coordinate system. The µ values are The ν values are The volume V 2 is therefore We define a third surface that determines the actual boundary of the LV cavity at the base. This surface is defined in cylindrical coordinates to avoid singularities associated with prolate spheroidal coordinate system surfaces. The cylindrical coordinates of this boundary can be expressed as The surface is defined on the domain u r ∈ [0, 1]×φ 0 ∈ [0, 2π). The value of the radius is given in terms of the reduced radial coordinate u r by r(φ 0 ) = u r r out (φ 0 ).
The z-coordinate is defined through The value of z mean is the value of z at the axis. This value is averaged from the endocardial z values as This choice provides a nearly flat closing surface for the cavity volume. The third volume term is therefore where z 2 is the z value of the previous surface at the same (r, φ ) locations.

S9.2 Basal work
Conservation of virtual work requires a term for the work done at the base. The basal work is equal to the work done by the fluid through the chamber closing surface Γ . The traction vector is We define surface tangent vectors in terms of the Cartesian coordinate vector The local surface tangent vectors are The surface normal can be computed as the normalized cross product of the two tangent surface vectors: The partial derivatives are where the two auxiliary derivatives are (S65) S10 Lumped parameter system The lumped parameter system includes variable resistance mitral and aortic valves along with 0D lumped parameter models for the left atrium, systemic arteries, and systemic periphery. The circulation is not treated as a closed loop, but is terminated with fixed pressures P pv and P sv at the pulmonary and systemic venous systems, respectively. The variable resistance valves are described by The parameters are defined in Table S2. P la and P pao represent the pressures in the left atrium and proximal aorta, respectively. The lumped parameter models are described by a set of volume-conservation ODEs.
The compliances, resistances, and fixed pressures are defined in Table S2. The variables P sa and P sp represent the pressure in the systemic arteries and systemic periphery, respectively. One additional algebraic equation is included to model the resistance of the proximal aorta:

S11 Numerical solution to the coupled system
Combining the LV ODE system (41), LV blood volume conservation (42), and proximal aorta pressure drop (S69), yields the coupled system: This is a square (n q + 2) × (n q + 2) system. It is linear in all terms except R mv and R aov , which depend exponentially on the pressures in the left ventricle P lv and proximal aorta P pao . We solve this system for the time derivatives of the LV deformation parameters, LV pressure, and proximal aorta pressure using a Newton iteration. We define the residual function H H H to be the right hand side of the coupled system (S70). Defining the Newton iteration variable the Jacobian is α 1,1 α 1,2 · · · α 1,n q −η 1 0 α 2,1 α 2,2 · · · α 2,n q −η 2 0 . . . . . . α n q ,1 α n q ,2 · · · α n q ,n q −η n q 0 The derivatives required by the Jacobian may be computed exactly using These equations have been written in a form that avoids numerical precision errors in the exponentials. The Jacobian entries are consequently A Newton iteration updates the vector X X X by solving for X X X i+1 and iterating to convergence. The full ODE system includes the time derivatives of the kinematic variables q q q which were solved for in the Newton interation and those of the lumped parameter models (S68). This ODE system is solved using a standard explicit time stepping scheme such as 4th order Runge-Kutta (RK4).

S12 Analytic cylinder solution
The stress distributions evaluated by the prolate spheroidal model are compared to the analytic cylinder stress distributions in Figure S5. The isotropic component of stress σ iso is not computed in the virtual work equations, but is integrated numerically from the computed solution.
Following the analysis of Rivlin (1949), using the convenient formulation written by Nash (1998), we consider an infinitely extended hollow incompressible cylinder of an isotropic Mooney-Rivlin material. The deviatoric part of the elastic strain energy is: where the invariants are The elastic component of the PK2 stress S S S e is computed from the elastic strain energy (S76) using (26). The cylinder is assumed to have undergone a fixed stretch λ in the z direction as well as a fixed torsion ψ. In cylindrical coordinates (r, θ , z), the deformation is Rivlin solved the finite elasticity problem under pressure loading conditions p in on the interior surface and p ext at the exterior surface. The radial deformation r(r 0 ) is computed by solving the strong form of the equilibrium equations (44). The results are outlined here for reference.
Here a in is the undeformed inner radius, a out the undeformed outer radius, µ in = µ(a in ), and µ out = µ(a out ). The analytical Cauchy stress components are with the isotropic component of the stress given by The system may be solved by first numerically solving for µ out using the boundary condition and subsequently solving for µ in in (S80). The remainder of the variables are evaluated explicitly. To compare our model solution to the analytical solution, we set the focal length of the prolate spheroidal coordionate system to a = 100, which produces an extremely elongated ellipsoid. In this configuration, the equatorial region closely approximates a cylinder. We limit the ν coordinate to include only a 1 cm region near the equator, and solve the virtual work equations using the Mooney-Rivlin material law. For this simulation, no active or viscous stresses are assumed, so the total stress (38) is given by S S S = S S S e . At the interior and exterior surfaces, the boundary conditions are as for the analytic solution, while natural boundary conditions are used at the other surfaces. The equilibrium virtual work system (43) is solved for the deformation with the assumed cavity pressure P lv = p in . We use the same parameters as those used by Nash (1998) shown in Table S3. The stress distributions evaluated by the prolate spheroidal model are compared to the analytic cylinder stress distributions in Figure S5. The isotropic component of stress p is not computed in the virtual work equations, but is integrated numerically from the computed solution.  Table S3. σ iso is the isotropic stress component that results from the material incompressibility.

S13 Cauchy stress divergence
We compute the Cauchy stress divergence in Cartesian coordinates. The Cartesian representation of the Cauchy stress may be computed from the prolate coordinate representation of the PK2 stress using where Q Q Q is a rotation matrix that transforms from prolate to Cartesian coordinates. This matrix is defined analogously to (22) as e e e x · e e e µ e e e x · e e e ν e e e x · e e e φ e e e y · e e e µ e e e y · e e e ν e e e y · e e e φ e e e z · e e e µ e e e z · e e e ν e e e z · e e e φ    .
To approximate derivatives of the Cartesian Cauchy stress, we first approximate derivatives with respect to the undeformed prolate spheroidal coordinates Θ Θ Θ = [Θ 1 ,Θ 2 ,Θ 3 ] = [µ 0 , ν 0 , φ 0 ] using finite differences. For example: approximates the derivative with respect to the first prolate spheroidal coordinate. These may be converted to derivatives with respect to the deformed Cauchy stress using To this end, we define the tensor transformation T T T by where θ are the deformed prolate spheroidal coordinates. Note that implying that we may compute the Cartesian divergence using The matrices required to compute T T T are