Gaussian Process Regressions for Inverse Problems and Parameter Searches in Models of Ventricular Mechanics

Patient specific models of ventricular mechanics require the optimization of their many parameters under the uncertainties associated with imaging of cardiac function. We present a strategy to reduce the complexity of parametric searches for 3-D FE models of left ventricular contraction. The study employs automatic image segmentation and analysis of an image database to gain geometric features for several classes of patients. Statistical distributions of geometric parameters are then used to design parametric studies investigating the effects of: (1) passive material properties during ventricular filling, and (2) infarct geometry on ventricular contraction in patients after a heart attack. Gaussian Process regression is used in both cases to build statistical models trained on the results of biophysical FEM simulations. The first statistical model estimates unloaded configurations based on either the intraventricular pressure or the end-diastolic fiber strain. The technique provides an alternative to the standard fixed-point iteration algorithm, which is more computationally expensive when used to unload more than 10 ventricles. The second statistical model captures the effects of varying infarct geometries on cardiac output. For training, we designed high resolution models of non-transmural infarcts including refinements of the border zone around the lesion. This study is a first effort in developing a platform combining HPC models and machine learning to investigate cardiac function in heart failure patients with the goal of assisting clinical diagnostics.

1. Introduction

Multi-scale models of cardiac mechanics, although are promising (e.g., Kerckhoffs et al., 2007; Nordsletten et al., 2011; Gurev et al., 2015; Land et al., 2017), have found limited applications for diagnosis and treatment. To reach the levels of accuracy needed to assist clinical decisions, models need to overcome major complications related to accessing clinical data, constraining unknown parameters, and coping with computational complexity. Some of the uncertainties associated to patient-specific cardiac models can be partially addressed with increased public access to large clinical datasets (Fonseca et al., 2011) and to high performance computing resources (Towns et al., 2014). Sophisticated finite element (FE) biomechanical simulations can be combined with machine learning techniques to translate parametric studies into efficient statistical models of virtual patient populations. Once an upfront computational cost is paid for training, the coupled effects of varying model parameters can be explored almost in real time, facilitating the solution of the optimization and inverse estimation problems that are required to personalize models for specific patients.

This paper discusses statistical models based on a machine learning technique called Gaussian Process (GP) regression, also known as kriging (Rasmussen and Williams, 2006). After training a “surrogate” of the more expensive FE models, GP regression can be used to assist optimization algorithms, even in complex cases where objective functionals cannot be easily differentiated (Booker et al., 1999; Abramson et al., 2009). More recently, GP regression has also been used in cardiovascular modeling, where it has found application in both fluid and solid mechanics (Marsden et al., 2008; Sankaran and Marsden, 2011; Pérez et al., 2016).

Recent developments in medical imaging techniques have opened new opportunities for cardiac modeling to augment image-based biomarkers from CT, MRI, and ultrasound scans (Lamata et al., 2014). As accuracy and availability of imaging modalities continues to improve, there is a growing need for novel strategies that exploit the capabilities of multi-scale models to enhance diagnostic tools. We present a systematic analysis of the Sunnybrook Cardiac MRI database, a public collection of cine-MRIs (Radau et al., 2009). Statistics gathered from the database were used to design two parametric studies investigating the passive behavior of the myocardium upon inflation and the effects of infarct on cardiac performance.

In the first parametric study, we developed a novel strategy to estimate the unloaded configuration (needed to initialize both passive and active FEM simulations) given either the end-diastolic intraventricular pressure, or the end-diastolic fiber strain. The new method relies on solving multiple forward problems to train a regression model from which unloaded configurations can be inferred for ventricles with arbitrary shapes. Despite such a problem could be alternatively solved with the fixed point iteration method (Sellier, 2011; Genet et al., 2015), our approach has some advantages. Specifically, our method can be easily applied in situations where the intraventricular pressure is not directly known (but could be inferred, for example, from the fiber strain), or where the unloaded geometry is one of the unknown parameters of an optimization problem.

The second example integrates machine learning and multi-scale modeling in a systematic parametric study investigating the effects of infarct on simulated cardiac performance. Location, size, and transmural depth of the infarct were chosen as input variables of a GP regression model predicting changes in simulated stroke volume due to the scar. This work exploited the capabilities of our in-house solver and an automatized workflow to run 40 simulations of infarct with varying shapes and locations. After training on results of FE simulations, the GP regression model provides a useful representation for the analysis of complex effects. Non-transmural infarcts were simulated with a high numerical accuracy.

2. Methods

2.1. Cine-MRI Segmentations and Parameterization via Idealized Models

Publicly available imaging datasets from the Sunnybrook Cardiac MRI database (Radau et al., 2009) were systematically processed to establish boundaries and proper feature distribution for parametric exploration. The Sunnybrook database gathered 45 cine-MRI scans collected from healthy subjects (N, n = 9), patients with ventricular hypertrophy (HYP, n = 12), and patients affected by heart failure both in presence and absence of myocardial infarction (HF-I, n = 12 and HF-NI, n = 12, respectively). For each scan, we considered only the short axis stack series, which provided ~10–15 axial slices per left ventricle (LV) and 20 frames per cardiac cycle. Average voxel sizes were (1.36 ± 0.057 mm) × (1.36 ± 0.057 mm) × (8.8 ± 1.0 mm) in the left-right, anterior-posterior, and apical-basal directions, respectively.

An in-house multi-atlas image processing technique (Xie et al., 2015) was used to co-register the axial slices of each dataset and then segment the LV boundaries. The first 2 columns of Figure 1A show the procedure applied to a representative 3-D image from the database. Outputs were labeled voxels marking the LV blood pool (shown in white semi-transparent overlay) and the ventricular wall (shown in red). The low resolution in the apical-basis direction typical of cine-MRI short axis views introduced segmentation artifacts that prevented direct use in FEM models. We therefore performed a further parameterization step (see third column) to approximate LV geometries as truncated prolate spheroids, as initially proposed by Streeter and Hanna (1973) and more recently revisited by Pravdin et al. (2014). According to such a scheme, the endocardial and epicardial profiles of an idealized axisymmetric LV were described by the following relations

$\begin{array}{l}{\rho }_{epi}=R_b\left[ecos\psi +(1-e)(1-sin\psi )\right]\\ {\zeta }_{epi}=Z(1-sin\psi )\\ {\rho }_{end}=(R_b-L)\left[ecos\psi +(1-e)(1-sin\psi )\right]\\ {\zeta }_{end}=(Z-H)(1-sin\psi )+H& \left(1\right)\end{array}$

linking the radial (ρ) and axial (ζ) coordinates of the epicardial and endocardial boundaries to the angle variable ψ ∈ [ψ₀, π/2]. In the equations above, the idealized geometry is defined by 6 parameters: the outer radius at base, R_b; the length of the longitudinal semi-axis of the outer spheroid, Z; the ventricular wall thicknesses at base and apex, L and H, respectively; the sphericity/conicity of the spheroid, e ∈ [0, 1]; and, finally, the truncation angle, ψ₀. Figure 1B shows a schematic of an idealized LV annotated with geometric descriptions of the parameters.

FIGURE 1

Figure 1. Automatic processing of cine-MRI images from the Sunnybrook Cardiac MRI database and fitting of idealized geometric model. (A) Complete processing of a representative short axis view frame from patient I-01 in the database. This Cine-MRI modality showed sufficient in-plane resolution, but significantly lower detail in the long axis view (e.g., compare first and second row of the first column). An atlas-based image processing algorithm was employed to extract LV boundaries for each patient. Segmented pixels are shown marked in red in the second and third column. Finally, an idealized 6-parameter model of LV geometry was fitted to the segmentation results, partially correcting for the artifacts introduced by the low resolution in the long axis (see model cross-section rendered in white in the third column). (B) Geometric meaning of idealized LV model parameters. Radial and axial coordinates are indicatd by ρ and ζ, respectively. R_b = outer radius at base, L = wall thickness at base, Z = distance from center of the ventricle to apex of outer wall, H = wall thickness at apex, Ψ₀ = truncation angle. Not shown is the e parameter, which governs the curvature of the LV external and internal walls. More details on analytical expressions of the LV geometric profile are provided in the text. (C) Top and lateral 3-D views of overlapped segmentation (rendered as a red surface) and best-fit idealized model (rendered as a gray transparent overlay) for a representative case (I-01 at beginning of diastole).

In order to describe the segmentation results in terms of the idealized models described above, we implemented an ad hoc optimization procedure to find sets of parameters ξ = {R_b, Z, L, H, e, Ψ₀} that would best match the MRI segmentations (I_MR). Each iteration involved first generating a binary 3-D image I_ξ marking the LV volume defined by ξ, and then evaluating an objective function J defined as

$\begin{array}{ll}J(I_{\mathrm{\xi }},I_{\text{MR}})=1-\frac{1}{2}\left(\frac{C_{\mathrm{\xi }}\cap C_{\text{MR}}}{C_{\mathrm{\xi }}\cup C_{\text{MR}}}+\frac{W_{\mathrm{\xi }}\cap W_{\text{MR}}}{W_{\mathrm{\xi }}\cup W_{\text{MR}}}\right),& \left(2\right)\end{array}$

where C_ξ and C_MR indicate the ventricular cavity regions in the idealized and MR segmentation images, respectively; and W_ξ and W_MR similarly indicate corresponding ventricular wall volumes. In other words, J ∈ [0, 1] provides a measure of similarity between a “synthetic” segmentation I_ξ generated for any given ξ and the actual MRI processing results I_MR. The “Nelder-Mead” algorithm available in SciPy was used to carry out the optimization up to convergence for every image dataset included in the database.

The relations in (1) do not include any parameters accounting for the rigid translation and rotations that LVs normally experience during a cardiac cycle. To overcome such limitation and to improve fitting results, each objective function evaluation was preceded by a rigid transformation step aimed at aligning the idealized model to the target segmented geometry. Specifically, we first estimated the main longitudinal axis of the segmented ventricle as the best-fit direction aligning the centers of gravity of the LV segmented axial slices. We then rigidly transformed the idealized models to let the longitudinal axes and the centers of gravity of the two geometries coincide. Figure 1C shows overlapped optimization results and corresponding MRI segmentation for a representative cine-MRI frame after rigid motion correction.

2.2. Passive Material Properties

To assess whether the inverse esimation method presented in this work would generalize to describe other constitutive behaviors (e.g., from future experiments on animal and human tissues, or from novel modeling developments), we considered 3 sets of material parameters (and related functional formulations) from the literature that describe experimental findings on canine, swine, and human ventricle biomechanics. Usyk et al. (2000) fitted a Fung-type orthotropic strain energy function to experiments on canine models

$\begin{array}{l}{W}_U=\frac{C}{2}(\text{exp}(Q)-1),Q=b_{\text{ff}}{E}_{\text{ff}}^2+b_{\text{ss}}{E}_{\text{ss}}^2+b_{\text{nn}}{E}_{\text{nn}}^2\\ +b_{\text{fs}}(E_{\text{fs}}^2+E_{\text{sf}}^2)+b_{\text{fn}}(E_{\text{fn}}^2+E_{\text{nf}}^2)\\ +b_{\text{ns}}(E_{\text{ns}}^2+E_{\text{sn}}^2),& \left(3\right)\end{array}$

where E_ij (i, j = f, s, n) are components of the Green-Lagrange strain tensor expressed in a reference frame locally aligned along the fiber direction (f), the orthogonal direction spanning the myocardial sheet (s), and the cross-fiber direction (n). Values for the C and b_ij (i, j = f, s, n) coefficients are reported in Table 1.

TABLE 1

Table 1. Sets of material properties considered in the study.

The remaining 2 constitutive behaviors here considered followed the constitutive law based on the invariants of the right Cauchy-Green strain tensor C proposed by Holzapfel and Ogden (2009),

$\begin{array}{l}{W}_{HO}=\frac{a}{2b}\left\{\exp \left[b(I_1-3)\right]\right\}+{\displaystyle \sum _{i=\text{\hspace{0.17em}ff},\text{ss}}\frac{a_i}{2b_i}}\left\{\exp \left[b_i{(I_{4i}-1)}^2\right]-1\right\}\\ +\frac{a_{\text{fs}}}{2b_{\text{fs}}}\left\{\exp \left[b_{\text{fs}}{I}_{\text{8fs}}^2\right]-1\right\},& \left(4\right)\end{array}$

where I₁=tr C is the first invariant of C, here applied as the argument of an exponential term; I_4i = v_i · (C · v_i), i = ff, ss is the fourth invariant of C, which corresponds to the squared stretch of a line element oriented along the fiber (v_ff) or sheet (v_ss) directions; finally, I_8fs = f₀ · (C · s₀) is the eighth invariant of C, which captures the effects of strain coupling. Equation (4) has been shown to describe well experiments on pig ventricles (Dokos et al., 2002), and more recently the biaxial and triaxial tests conducted on human myocardial tissue by Sommer et al. (2015). Among best-fit values reported in literature, we selected materials parameters for (4) from Wang et al. (2013) ($W_{HO}^W$, fitted to experiments on swine models) and Gültekin et al. (2016) ($W_{HO}^G$, fitted to experiments on human tissue). The coefficients for all considered material properties are reported in Table 1.

2.3. FEM Models of LV Passive Biomechanics

High-resolution FEM simulations of LV biomechanics are at the core of the parameter exploration and inverse estimation strategies presented in this work. To cope with the complexities of the mechanical behavior of the myocardium, we employed a recently validated numerical solver suitable for dealing with incompressible hyperelastic material laws such as those in (3) and (4) (Gurev et al., 2015), and extended to use stabilized P1/P1 finite elements. The capabilities are necessary for infarct simulations, where capturing sufficient detail at the border zone region around the lesion is pivotal (see section 2.6). The solution algorithm also allows multi-scale effects, and we used the TriSeg ODE-based model with parameters for human to drive myofilament active contraction (Lumens et al., 2009; Gurev et al., 2015). Coupling between cellular and tissue mechanics occurred at the Gauss point level.

To handle the relatively large number of simulations needed to train statistical models, we developed an automatic workflow to construct high-resolution computational domains from any given sets of geometric parameters ξ describing LV anatomy. In this pipeline, analytical models built according to (1) were first converted to 3-dimensional triangulated surfaces, and then to solid meshes of several hundred thousands of tetrahedral elements. Nodes at the base of the ventricle were prevented to move axially, while epicardial nodes in the vicinity of the base (i.e., closer than 3 mm) were fully locked to prevent rigid motions. Boundary traction effects from the pericardial membrane and the right ventricle were neglected, and intraventricular pressure was uniformly applied at the endocardial surface in quasi-static steps. The vector v_ff of alignment of myocardial fibers varies heterogeneously along the radial direction of the myocardium (McCulloch, 1999; Humphrey, 2002). Without specific measurements for the patients in the database, we relied on a rule-based approach to assign fiber directions linearly varying their angle with respect to the circumferential direction from 90° at the endocardial surface (i.e., longitudinally aligned) to -60° at the epicardium.

The mechanical equilibrium equations were solved in parallel on the Cognitive Computing Cluster (CCC), a hybrid high performance shared resource developed at IBM Research deploying both Intel and Power8 nodes. Active infarct simulations required ~10 times more resources than passive models, and were run on the Uran Supercomputer hosted by UB RAS and Ural Federal University. Outputs of the simulations were nodal displacement vector fields, and components of stress and strain tensors defined at the element Gauss points. To relate predictions also to strain dependent length activation of the sarcomere, we also evaluated stretch in the fiber direction, defined as

${\lambda }_{\text{ff}}=\sqrt{v_{\text{ff}}·C·v_{\text{ff}}},\left(5\right)$

where v_ff is the vector aligned along the myofiber direction (as described above), and C is the right Cauchy-Green strain tensor. As a representative scalar of each loading state, we also averaged λ_ff at midwall, which we defined as a tissue slab located between 40 and 60% of the LV wall thickness and between 45 and 55% of the apex-base distance.

2.4. Parameterization of FEM Results

A key aspect of the inverse unloading method presented in this work is the re-parameterization of FEM simulation results in terms of the same geometric parameters employed to process the Sunnybrook database. A 2-step optimization procedure was implemented to fit idealized models of LV anatomy to the deformed configurations predicted by the FEM analyses upon varying loading conditions. First, optimal values for R_b, Z, e, and Ψ₀ were found to minimize average nodal distance between the profile of an idealized epicardium and the corresponding boundary obtained from a FE mesh warped according to the simulations results. Second, a similarly defined nodal distance measure was used to quantify discrepancies between endocardial profiles in order to adjust the remaining L and H parameters. The 2 steps were re-iterated until reaching convergence. An alternative monolithic approach where the 6 parameters were optimized at the same time was also evaluated, but proved to be less computationally efficient.

2.5. Statistical Learning of LV Unloading

Bulk processing the Sunnybrook cine-MRI image datasets provided information on expected anatomical variability among patients. As part of our inverse unloading estimation strategy, we leveraged database statistics to define a 6-D parameteric space that enclosed all likely LV unloaded configurations. More specifically, we reasoned that the parametric study should conservatively admit and explore large variations in ventricle geometries, since the unloaded state might differ significantly from any of the imaged configurations. Limits of the parametric space were therefore defined to encompass variations of more than 3 standard deviations from the average beginning of diastole (BoD) state, which we chose as most reasonable guess lacking the measurements needed for better estimates (e.g., Xi et al., 2013). More details on the subdivision of the cardiac cycle into its phases are reported in the Supplemental Material. Figure 2A shows pairs of limit parameter values and corresponding LV cross-sections representing maximum allowed variations of each of the 6 geometric features. In drawing the profiles, only one of the 6 parameters was changed while keeping the remaining 5 at corresponding mid-range values. Unloaded configurations admitted to our study were, therefore, intermediate states of the low- and high-parameter geometries shown in Figure 2A in gray and black tones, respectively. The statistical distribution of BoD states was also used to design an efficient probing scheme for the parametric space defined above. As is common now (Marsden, 2014), we used latin hypercube sampling to select 600 points (i.e., 100 times the number of parameters) from a normal distribution centered on the average BoD state and with doubled standard deviation compared to that of the Sunnybrook database. Samples falling beyond the limits defined in Figure 2A were projected onto the closest admissible point. A cloud of chosen probing locations is shown within a unitary 3-D projection of the parameter space in Figure 2B. For this representation, the H, e, and Ψ₀ dimensions were neglected. Figure 2C further shows a mid-slice of the parametric cube exploring geometries corresponding to coupled variations of the R_b and Z parameters.

FIGURE 2

Figure 2. Design of training sets for the 2 statistical models: LV unloading (A–C) and infarct shape effects (D–F). (A) Pairs of LV cross-sections representing extreme geometries limiting parameter space dimensions. Gray (black) cross-sections correspond to extreme negative (positive) variations of one of the geometric parameters, with the remaining 5 parameters kept at mid-range values. (B) Projection of the 6-D parametric space onto a 3-D cube obtained by neglecting the last 3 dimensions (H, e, and Ψ₀). Spherical glyphs indicate locations of 600 sampling points chosen via latin hypercube sampling from a normal distribution centered on the average LV geometry and with a doubled standard deviation compared to that of the complete Sunnybrook database. (C) Cross-section of the parameter space for LV unloading showing combined variations of R_b and Z parameters. (D) Similar to (A), but showing pairs of FE meshes including infarct regions with extreme shapes. The lightest tone of gray indicates the healthy region, the darkest tone indicates the infarct, and the intermediate one marks the refined border zone. (E) 3-D projection of the 4-D parameter space defining infarct shape obtained neglecting the ΔLong. dimension. Similarly to (B), spheres indicate locations of 40 sampling points chosen uniformly in the allowed range parameters. (F) Mid-range slice of the 3-D projection showing representative FE meshes accounting for combined variations of longitudinal location and transmural depth of the infarct.

For each of the 600 sampled ventricle geometries we ran passive inflation simulations for inner LV pressures ranging between 0 and 5 kPa. Results were processed as described in section 2.4 to find optimal geometric parameters for 100 intermediate loading configurations (i.e., differing by 0.05 kPa). These best-fit parameters constituted the training set for GP regression models mapping loaded configurations to their corresponding unloaded state. Overall, we optimized 100 statistical models (one for each considered inner pressure), and fitted 2 additional GP regressions for unloading the inflated configurations for which the midwall fiber stretch reached the values of 10 and 15%.

2.6. Statistical Learning of Infarct Shape and Size on LV Performance

With our solver capable of handling high-resolution tetrahedral meshes, we explored the effects of different infarct shapes and locations on simulated LV cardiac cycles. The lesions were parameterized according to 4 features: longitudinal position (Long. ∈ [0, 1]), indicating whether an infarct was closer to the base (Long.=0) or apex (Long.=1); circumferential extension (ΔCirc. ∈ [0, π]), indicating the portion of circumference (measured in radians) occupied by the infarct; longitudinal extension (ΔLong. ∈ [0, 1]) indicating the fraction of longitudinal cross section harboring a lesion; and wall depth (Depth ∈ [0, 1]), indicating the transmural extension of the infarct, with the maximum value of 1 indicating a fully transmural lesion. Figure 2D shows computational domains reconstructed from limit cases of the infarct parameterization. Similar to that presented in section 2.5, latin hypercube sampling was used to efficiently probe the parameter space. Our sample size was of 40 points, (i.e., 10 times the number of parameters), and we assumed a uniform distribution of parameters across the admissible range. Also, to restrict our attention to the effects of infarct without the added complications introduced by changing LV geometry, all lesions were inserted into the same baseline LV from patient I-02. Infarct effects were simulated by deactivating active contraction in the lesion regions, while maintaining the same passive material properties. Similar to Figures 2B,C,E,F show projections of selected samples onto the considered parameter space of infarct lesions. More details on the general procedure followed to mesh infarcted regions of arbitrary shapes are available in the Supplemental Material.

3. Results

Once enhanced with rigid motion correction, the 6-parameter description of LV geometry was able to capture anatomical and kinematic features from the Sunnybrook MRI scans. Median values of the similarity functional J(I_ξ, I_MR) averaged for each category of patients were 0.29 for N, 0.30 for HYP, 0.23 for HF-NI and 0.19 for HF-I, respectively. Figure 3 shows average group traces of best-fit geometric parameters (see Equation 1) over the course of a normalized cardiac cycle. Certain trends agreed well with known morphologic features of cardiac disease. Patients affected by heart failure (i.e., from the HF-I and HF-NI categories) presented on average the most dilated ventricles, as indicated by the largest R_b values, and the smallest cyclical variations in both e and Ψ₀, probably due to myocardial dysfunction. Hypertrophic patients, on the other hand, maintained highest L values throughout the cycle (L = 12 mm on average) and showed a large systolic thickening (L = 15 mm at peak systole for HYP patients). Only N subjects contracted more visibly, with an average 54% increase in L from diastole to systole. N and HYP subjects overall exhibited the largest changes in truncation angles. Other parameterization findings were less intuitive. For all LVs, contraction in the longitudinal direction was captured mainly by varying Ψ₀ rather than Z, which instead remained close to constant throughout the cycle. Also, the dynamic pattern of e observed in HF patients was peculiar. For example, 10 out of 12 HF-I subjects exhibited increased e at systole compared to diastole, while the opposite was typically observed in the N and HYP categories of patients. Combined behavior of e and Ψ₀ differed also among the 2 classes of HF patients: in presence of an infarct, both e and Ψ₀ were smaller in magnitude, indicating that HF-NI ventricles tended to be more spherical than the HF-I ones. Table 2 reports best fit sets of geometric parameters for all 45 patients at both beginning and end of diastole (BoD and EoD, respectively).

FIGURE 3

Figure 3. Kinematics of LV motion during a normalized cardiac cycle described by changes over time to the 6 parameters of an idealized model of LV geometry. Traces show averages (marked solid line) over the 4 patient groups (N, normal; HYP, hypertrophic; HF-NI, heart failure without infarct; HF-I, heart failure with infarct) plus or minus standard deviations (shown as semi-transparent overlays). Each subplots shows results for one of the geometric parameters. Results were obtained by custom procedure to fit idealized model to segmentation results. See text for more detail.

TABLE 2

Table 2. Best-fit geometric parameters for all patients in correspondence of beginning and end of diastole configuration (BoD, and EoD, respectively).

The distribution of LV shapes at BoD (see Table 2) was pivotal to design our admissible parameter space, both for establishing range limits and for choosing the frequency of allowed variations. Figure 4A shows 3 representative unloaded configurations out of the 600 selected to probe the space. Each geometry was first discretized into a computational domain (see meshes below the idealized profiles) and then inflated with inner pressures up to 5 kPa. Shown also are color coded distributions of the first invariant of the Green-Lagrange strain tensor (I_1E). Strain fields were visibly larger in the LVs endowed with W_U material properties (i.e., those on the first row of each subgroup) than in those endowed with $W_{HO}^G$ (i.e., those on the second row of each subgroups). While the parameteric study extensively explored combined effects of LV geometric features on deformation, the subsequent processing step converted results back to the 6-parameter description (see profiles in gray above and below strain results). Out of the chosen 600 probing profiles, 67 exhibited incompatible features that prevented completion of the FEM simulations (e.g., a disproportionately large L and negative Ψ₀ in a ventricle with minimum R_b), and were therefore excluded from the analysis. Figure 4B shows violin plots of geometric parameter distributions for ventricles at the BoD configuration from the database, at the assumed unloaded configuration, and at 10 deformed configurations for pressures ranging from 0.5 to 5.0 kPa. The BoD distributions (see plots in black tone, leftmost sector) clearly reflected the categories of the database. For example, the violin plot of the R_b parameter (first row) indicated a bimodal distribution, as expected given the sharp differences in ventricle radius between HF patients and the others. By design, the sampled unloaded configurations followed a normal distribution allowing large variations (see plots in lightest gray tone, second sector from the left). Some hard limits on admissible parameter values were enforced to reduce the number of incompatible geometries selected (see section 2.5). The effects of these hard limits were noticeable particularly within the L, e, and Ψ₀ distributions (see last 3 rows), the tails of which were thickened by assimilating parameters beyond allowed range. Finally, the distributions of loaded configurations (see plots in intermediate gray tones, three rightmost sectors) showed the evolution of geometric parameters upon pressurization, which followed the expected behavior for incompressible hyperelastic materials. For example, the R_b parameter increased relatively fast at low pressures, but then dilation progressively stopped accounting for the exponential increase in stiffness. The thickness parameters L and H decreased upon pressurization (also ensuring incompressibility), while the Ψ₀ parameter distributions were the most sensitive to pressure. Finally, the material properties could be ranked in order of increasing stiffness as W_U, $W_{HO}^G$, and $W_{HO}^W$, as shown by changes in mean values from the distributions (see white lines inside the violin plots).

FIGURE 4

Figure 4. Generation of training dataset for the unloading problem. (A) Passive inflation and subsequent parameterization results for selected sample of 3 out of 600 left ventricle geometries considered to build the training datasets for the unloading problem. Idealized geometries, chosen via latin hypercube sampling to probe the parameter space, were discretized and subjected to passive inflation using 3 different sets of material properties. Shown in the panel are results for 2 sets of material properties (W_U and $W_{HO}^G$) and 2 loading pressures (1 and 5 kPa). Shown also are color-coded distributions of strain expressed as the first invariant of the Green-Lagrange tensor. (B) Violin plots depicting changes in geometric parameter distributions upon inflation for all the 533 LV geometries included in the training datasets, and for the 3 sets of material properties. Black tone plots indicate distributions of geometric parameters at the BoD configuration. Lightest gray tone plots correspond to distributions synthesized via latin hypercube sampling from a normal distribution constructed based on the BoD configuration, but allowing 2-times larger variations. White segments close to the center of the distribution indicate mean values. See text for more details.

The computational cost of optimizing a GP regression to a few hundred training points (~1 CPU min) is negligible compared to that of running even only a single passive high resolution simulation. To optimize the use of computational resources, we sought, therefore, the minimum training set size that ensured satisfactory accuracy in estimating the unloaded configurations for all patients in the dataset. Figures 5A,B show cases where predictions by GP regression compared best (see Figure 5A) and least well (see Figure 5B) to the configurations predicted via fixed point iteration for a relatively small training size (n_train=75). As starting (loaded) configurations, we chose geometries from the database at EoD (see first column in both panels), and from these we inferred corresponding unloaded configurations assuming inner LV pressures of either 1 or 2 kPa. Comparison between results from the 2 methods were evaluated in terms of Dice score between unloaded profiles (see Supplemental Material for details on Dice score computations). According to our analysis, n_train=75 was the minimum training set size ensuring Dice scores larger or equal than 0.90 for all cases considered (i.e., including all the LV geometries, both EoD inner pressures, and the 3 sets of material properties). From the last column of Figure 5B one can appreciate how even a Dice score of 0.90 corresponds to a visibly good match between the GP regression prediction (see LV in gray tone) and corresponding geometry obtained via fixed point iteration (see overlapped dashed line). Small mismatches could be observed even in cases with high Dice score in regions close to the base of the LV (e.g., see last column of Figure 5A). These artifacts could be attributed to the zero-displacement boundary condition applied to epicardial elements within 3 mm from the base in the fixed point optimizations. Note that the fixed point iteration method required discretization of the EoD domain and repeated mesh deformation steps (see middle row in both panels). In contrast, after GP regression training the unloaded configurations could be inferred almost in real time, and as another advantage, the GP regression method eliminates potential issues introduced by iteratively warping the mesh (e.g., element degeneration) in the fixed point iteration method. The top row in both panels shows unloaded profiles inferred from GP regressions trained on the full set of simulation results. In the best match case shown in panel A results were essentially the same for n_train=75, while the Dice score increased by 0.06 when we expanded the training dataset from n_train=75 to 533 in the worst match case (see Figure 5B).

FIGURE 5

Figure 5. Unloading via kriging and comparison to the fixed point iteration method. (A) Unloading procedure is shown applied to a representative case (NI-14, unloading pressure P = 1 kPa, and $W_{HO}^G$ material properties) for which a statistical model trained on 75 arbitrary ventricles matched best unloading results via fixed point iteration method. While the fixed point iteration method required meshing of the ventricles in the loaded configuration and iterative updates (middle row), the statistical method allowed to infer the unloaded geometry directly from the 6 parameters describing the end-diastolic (loaded) configuration (bottom row). Top row is similar to bottom row, but shows result obtained after training a statistical model on results from the full parametric study of 500+ LVs. The rightmost column shows overlapped cross-sections of unloaded LVs obtained via the fixed point iteration method (dashed boundary) and 2 statistical models (solid gray tones). (B) Similar to (A), but applied to another representative case (I-07, unloading pressure P = 2 kPa, and W_U material properties) for which the statistical learning method (with n_train = 75) yielded the worst overlap to fixed point iteration results (Dice score of 0.90). In this case, increasing the training set size led to improved results (Dice score of 0.96).

Figure 6 plots average Dice scores comparing GP regressions to fixed point iteration. The 3 rows in Figure 6A show results for different sets of material properties at an unloading pressure of 1 kPa (first column) or 2 kPa (second column). As expected, increasing training sizes generally yielded better Dice scores, although little improvement was observed beyond n_train=75. Also reported are average Dice scores quantifying the overlap between fixed point iteration and the BoD or EoD configurations, as well as the average overlap with the OptD configuration, which was chosen as the imaged diastolic configuration that matched best the unloaded geometry. White dashed lines overlapped to the bars indicate the lower 10th-percentile Dice score observed for predictions from GP regressions.

FIGURE 6

Figure 6. Accuracy of unloading via GP regression compared to unloading via fixed point iteration method. (A) Barplots of average Dice scores comparing beginning of diastole (BoD), end of diastole (EoD), the best-matching diastolic configuration (OptD), and unloaded configurations obtained via kriging with different n_train to unloaded configurations predicted via the fixed point iteration method. Subplots show results for combinations of considered material properties (W_U , $W_{HO}^G$, or $W_{HO}^W$) and unloading pressures (P = 1 kPa or 2 kPa). (B) Similar to (A), but unloaded configurations are estimated prescribing average midwall strain at end diastole (${\lambda }_{\text{f}{\text{f}}^{10\text{\%}}}$ on the left, ${\lambda }_{\text{f}{\text{f}}^{15\text{\%}}}$ in the right column). In both panels the dashed white lines drawn on kriging-related bars indicate lower 10th-percentile Dice score for each subcategory.

Additional GP regression models were trained to handle situations where intraventricular pressure is unknown, but can be estimated by indirect measurements such as the fiber strain at midwall (see section 2.3). Table 3 reports unloaded geometries for all patients in the Sunnybrook database under the assumptions of W_U material properties and end-diastolic fiber strains of either 1.10 (${\lambda }_{\mathrm{\text{ff}}}^{10\%}$) or 1.15 (${\lambda }_{\mathrm{\text{ff}}}^{15\%}$). Outputs of the procedure included end-diastolic LV pressure values corresponding to the target fiber strains in the loaded configurations. Figure 6B reports accuracy of GP regression predictions measured in terms of Dice score with predictions via fixed-point iteration method.

TABLE 3

Table 3. Unloaded geometries inferred via GP regression assuming EoD fiber stretches at midwall of either 1.10 (${\lambda }_{\mathrm{\text{ff}}}^{10\%}$) or 1.15 (${\lambda }_{\mathrm{\text{ff}}}^{10\%}$) and W_U set of material properties.

Other than being used for inverse problems, GP regressions are ideal as tools for rapidly exploring multi-dimensional parameter spaces. As a proof of concept for the usage, we show preliminary results for a parametric study of the effect of infarct location and shape on cardiac performance as assessed by stroke volume (SV). Figure 7A shows color maps of simulated SV over 2-D slices of the 4-D parameteric space. Also shown, are projections onto each slice of the probing locations composing the full training set (see black dots). Each plot isolates the combined effects of 2 out of the 4 parameters used to define infarct shape and location. As expected, increases in lesion sizes yielded significant drops in SV. Maximum combined effect was reached by increasing both circumferential and transmural extension. Starting from a baseline failing LV with SV = 49 ml, GP regression predicted a drop down to SV = 21 mL at maximum depth and circumferential extension. Figure 7B shows 5-fold cross-validation for evaluating progressive convergence of GP regression for increasing training sizes. Average relative discrepancies between SV values from simulations and corresponding predictions from GP regression progressively decreased to 6% for a maximum training size of 40 simulations.

FIGURE 7

Figure 7. Statistical model of infarct shape and location effects on simulated SV. (A) Color-coded distribution of SV as predicted by kriging on 6 slices of the 4-D parametric hyperspace. Each plot shows combined effects of variations of 2 parameters on simulated SV, as shown by scale bar (values outside the range are truncated). Darker (lighter) color tones indicate stronger (weaker) impairment due to infarct. Dots represent projections of the probing points onto the slice plane. (B) 5-fold cross-validation to assess performance of the statistical model for varying training sizes n_train. Relative error on simulated SV predictions approached 6% for the maximum training set size (n_train = 40).

Figure 8A compares in detail 2 simulations from the training set characterized by different infarct morphologies. While INF₁₆ (on the left) harbored a non-transmural basal infarct, the lesion in INF₃₀ was larger, more apical, and fully transmural. The high level of mesh refinement within and surrounding the infarct (see regions in darkest and intermediate gray tones, respectively) required the capability of our solver of handling high resolution tetrahedral meshes. Figure 8B compares simulated PV loops for the 2 models described above. As expected, INF₃₀ (see dashed line), which harbored a larger lesion, exhibited a stronger impairment in simulated cardiac performance. The PV loops show the weaker contraction generated by INF₃₀ despite an increase in end-diastolic volume (i.e., SV = 40 ml and SV = 32 ml for INF₁₆ and INF₃₀, respectively).

FIGURE 8

Figure 8. Comparison between 2 select simulations (out of the 40 considered). (A) On the left, the INF₁₆ model has a smaller basal infarct (volume of 8.8 ml, Long. = 0.43, ΔCirc. = 1.22, ΔLong. = 0.40, Depth = 0.40). On the right, INF₃₀ presents a larger transmural lesion (volume of 16 ml, Long. = 0.79, ΔCirc. = 1.57, ΔLong. = 0.46, Depth = 1.0). (B) Simulated PV loops showing smaller SV for the largest lesion INF₃₀ (dashed line), as expected.

4. Discussion

Numerous computational models of LV mechanics have been developed over the years to understand better LV function in normal and diseased hearts with the ultimate goal of assisting personalized diagnostics and treatment. Available models differ both in terms of enclosed biophysical detail and of anatomical representation. In the simplest form, left ventricular function can be captured by a time-varying elastance model, where a single time-varying ODE couples the evolution of intraventricular pressure and volume over the course of a cycle (Suga and Sagawa, 1972; Stergiopulos et al., 1996). At the other end of the complexity scale, models of LV mechanics incorporate phenomenological or biophysical descriptions of muscle contraction at the microscopic level, while at the same time capturing in detail the cardiac anatomy on high-resolution computational domains (e.g., Guccione et al., 1995; Kerckhoffs et al., 2007; Göktepe and Kuhl, 2010; Baillargeon et al., 2014; Sundnes et al., 2014; Gurev et al., 2015; Augustin et al., 2016). Although these highly refined 3D models provide valuable information, they entail high computational costs. To improve computational efficiency, models with intermediate levels of complexity have been based on simplifying assumptions on ventricular geometry and structure (Arts et al., 1979; Beyar and Sideman, 1984; Lumens et al., 2009). For prolate spheroid geometries and passive mechanics simulations, distributions of stress in other low order models can match well FEM results despite running faster than in real-time (Moulton and Secomb, 2013, 2014; Moulton et al., 2017).

Significant reductions in computational costs can be similarly achieved by training machine learning models on the results of opportunely sampled biophysical simulations. As a proof of concept, in this paper we applied GP regression, a popular supervised learning technique, to 2 problems of interest in cardiac mechanics modeling. First, 600 LV geometries described by a 6-parameter (R_b, L, Z, H, Ψ₀, e) prolate spheroid were extracted randomly from a conservatively defined parameter space. For each geometry, a forward simulation was run to trace ventricle geometries upon inflation at progressively larger intraventricular pressures. GP regression models then allowed to infer unloaded configurations given sets of 6 parameters defining the loaded geometries and either their corresponding intraventricular pressure or their fiber strain at midwall. For the second statistical model, we built a GP regression between parameters characterizing the location and shape of an infarct and corresponding stroke volumes predicted by high-resolution simulations accounting for the presence of the lesion.

4.1. Ventricular Shape Analysis

The Sunnybrook Cardiac MRI database was the primary source of imaging data for this study. Conventional analyses of the segmentations from such a database have employed methods to either extract features directly from images (e.g., Chumarnaya et al., 2016), or have used finite element models to analyze ventricular shapes and build statistical classifiers of patient disease (e.g., Piras et al., 2017). A geometric description with fewer parameters is better suited for parameterizing the geometry of ventricles in regressions trained on biophysical simulation results. Therefore, instead of finite element models, we adopted a 6-parameter description (Streeter and Hanna, 1973; Pravdin et al., 2014) to approximate ventricular geometry. In spite of its simplicity, this approach was able to capture some of the shape features and biomarkers that have been previously extracted using the conventional finite element models (e.g., Zhang et al., 2014). In particular, ventricular sphericity (e) separated ventricles with and without myocardial infarction in patients with heart failure (see HF-I and HF-NI traces in Figure 3). The 6-parameter model analysis also captured higher average wall thickness in hypertrophic hearts and highest relative dynamic thickening in normal patients. To partially compensate for the limits of considering a fully axisymmetric parameterization, we accounted for eventual rigid rotations and translations to better align parameterized and segmented ventricles throughout the cardiac cycle. This ensured us overall good fitting results, especially for the failing hearts, which proved to be more symmetric. Nonetheless, the methods here presented could be promptly extended also to non-axisymmetric parameterizations such as those based on non-uniform rational B-splines at the expense of extending the parameter space to additional dimensions.

Out of the several field views provided in the Sunnybrook database, we restricted our analyses to short-axis stack series, which have the disadvantage of providing relatively low resolution in the coronal planes. As a result, some artifacts were particularly evident close to the apex of the ventricle, where the segmentation and subsequent parameterization were sometimes not able to resolve correctly the apical thickness, especially in the thinner failing LVs. Not surprisingly then, the H parameter showed the largest relative standard deviations within the same cardiac cycle for all patients, indicating that apex parameterization accuracy could be likely corrected by registering and merging multiple MRI views.

4.2. Ventricular Unloading

Standard FE simulations need to be initialized from an unloaded state, which cannot be directly extracted from images because ventricles are pressurized in all of the configurations imaged by cine-MRI or CT scans. Given material properties and inner LV pressure, iterative approaches such as the fixed point iteration method allow to estimate the unloaded configuration by progressively correcting a loaded state (Sellier, 2011; Genet et al., 2015). Nonetheless, due to their large computational cost and added complexity, these techniques are not typically incorporated into sophisticated optimization schemes proposed to estimate model parameters from images (Asner et al., 2016, 2017; Nasopoulou et al., 2017). To ensure feasibility, many modeling studies tend instead to use representative loaded configurations (i.e., at beginning or end of diastole) as approximations for the unknown unloaded state. As shown by our analyses, this could significantly bias results, since BoD and EoD configurations tend to match poorly to the profiles of unloaded geometries (see Figure 6). GP regression models of unloading can help circumvent some of the limitations associated with iterative methods and enable larger parameter search studies. Somewhat surprisingly, even a training set of n_train = 75 forward simulations was sufficient to ensure good inverse estimation results. LV profiles inferred from the statistical model matched those obtained via fixed point iteration with Dice scores always larger than 0.90 under 2 loading pressures and for 3 different sets of material properties. Considering that in our experience 7–10 iterations are needed to reach convergence via fixed point iteration, the preparation of an accurate statistical model might then require a computational cost comparable to unloading 7–10 ventricles with the standard method. Unlike fixed point iteration our strategy requires also an additional step of re-parameterizing simulation results in a format that can be handled by the machine learning model. The computational cost of reparameterizing is often negligible (on the order of few CPU mins), and after training the statistical model can be further interrogated to unload additional geometries at essentially no computational cost.

In addition to morphology, estimating the unloaded configuration relies on the knowledge of loading conditions and of the material properties of the myocardium. To fully characterize the material behavior of cardiac tissue, sophisticated experiments are required to reproduce in vitro the principal strain modes experienced by the heart during the cycle. The most extensive dataset on the passive behavior of the human myocardium is provided by Sommer et al. (2015). This work confirms how the micro-architecture of myocardial sheets leads to complex nolinear anisotropic behavior combined to a persisting viscoelastic response. Although viscoelastic effects were neglected in this work, we considered material properties based on the triaxial experiments of Sommer et al. (2015) as well as 2 other sets of constitutive behaviors based on experiments on animal models (Usyk et al., 2000; Wang et al., 2013; Gültekin et al., 2016). Our unloading procedure proved to work well for all of these sets of material properties.

In the form presented herein, our method for ventricular unloading required building a new training dataset and subsequently a new GP regression model for each set of material properties considered. Nonetheless, for future applications, the input parametric space could be extended to additional dimensions to account also for variations in material properties. While more training simulations would likely be needed to reach the desired convergence, the presented approach could still prove to be convenient for material property identification based on strain energy functions with reduced number of parameters (e.g., Nasopoulou et al., 2017), and especially in cases where large high performance machines are available to tackle the required computational cost in a distributed fashion.

The diastolic fiber strain at midwall, the constitutive law, and the shape of the ventricles at end-diastole are sufficient to uniquely unload geometries either via fixed-point iteration or GP regression. In this paper, we proposed to constrain end-diastolic fiber stretch to account for scenarios where diastolic pressure in the ventricles is not known. Animal model experiments suggest that end-diastolic fiber strain varies within a relatively small range in several circumstances (e.g., see Ross et al., 1971). Inspired by studies on inverse stress identification (Miller et al., 2010; Miller and Lu, 2013), we therefore tried to find the unloaded ventricular shape without solving for ventricular pressure. This was also motivated by the fact that unloading by strain would yield the same unloaded configuration independently from a homogeneous scaling of the constitutive law (i.e., predicted end-diastolic pressures would scale accordingly). To illustrate the potential of such approach, we additionally computed Dice scores between unloaded ventricles with 10% diastolic fiber strain using different constitutive laws. Our results (Dice scores of 0.90±0.05 for W_U vs. $W_{HO}^G$, 0.85±0.03 for W_U vs. $W_{HO}^W$ and 0.96±0.03 for $W_{HO}^G$ vs. $W_{HO}^W$, respectively) suggested strong similarity between unloaded ventricles endowed with umlaut Gultekin and Wang material behaviors, which followed the same Holzapfel-Ogden functional formulation.

4.3. Modeling of Infarct Mechanics

Two main factors increase the complexity of ischemia and myocardial infarction models. The first one is the need to account for the progressive changes in passive and active material properties that are triggered by the lesion and driven by tissue damage recovery and remodeling (Holmes et al., 2005). The second one is the more complex numerical framework required to handle the large finite element meshes needed to accurately capture realistic infarct shapes. In the past, only few studies have simulated non-transmural infarcts (Leong et al., 2015; Duchateau et al., 2016; Leong et al., 2017), while most models have either simulated infarct with simplified morphologies, or have allowed infarct/ischemic regions crossing the finite element boundaries (e.g., Mazhari et al., 2000; Jie et al., 2010; Wenk et al., 2011; Mojsejenko et al., 2015). Here, we present a model of non-transmural infarct that has refined elements in the border region of infarct. To handle large finite element meshes that result from such a refinement, we use an iterative solver for the large system of linearized equations with an efficient preconditioner (Gurev et al., 2015). To quickly summarize our results, the 2 main parameters affecting simulated SV were the transmural and circumferential extensions of the lesion, while location of the infarct played a minor role. Our models of infarct and the corresponding statistical model are still at a preliminary stage of development, and were here presented mainly to demonstrate the concept of integration between statistical and physical models.

4.4. Summary

This work shows 2 applications of GP regression in modeling ventricular heart mechanics. First, we present a strategy to estimate the ventricular unloaded configuration given material properties and intraventricular pressure (or alternatively fiber strain at midwall). Once an upfront computational cost (amounting to ~10 applications of a conventional iterative method) is paid for training, GP regression models allow the estimation of unlimited unloaded geometries at no additional cost. The method is therefore suitable to be used in analyses involving large number of patients such as those collected in publicly available databases. Second, we use GP regression as a convenient tool to explore results of a parametric study investigating coupled effects of infarct shape and location. While just a proof of concept study, these preliminary results demonstrate the power of the approach. That is, we were able to characterize a large variation in infarct location and size, including non-transmural infarcts with highly complex meshes that are computationally demanding to solve.

Author Contributions

PD designed the study, analyzed results, and wrote the manuscript, AH provided the image segmentations and supervised the image processing pipeline, SK performed the parametric infarct simulations, OS supervised the infarct simulations, JR wrote the manuscript and supervised the project; VG designed the study, wrote the manuscript, and supervised the project. All authors agree to be accountable for the content of the work.

Funding

Supported by the Program No. 27 of the Presidium of the RAS, the Decree of the Government of the Russian Federation No. 211 of 16/03/2013.

Conflict of Interest Statement

PD, AH, JR, and VG were employed by IBM Research.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2018.01002/full#supplementary-material

References

Abramson, M., Audet, C., Dennis, J., and Digabel, S. (2009). OrthoMADS: a deterministic MADS instance with orthogonal directions. SIAM J. Optim. 20, 948–966. doi: 10.1137/080716980

CrossRef Full Text | Google Scholar

Arts, T., Reneman, R. S., and Veenstra, P. C. (1979). A model of the mechanics of the left ventricle. Ann. Biomed. Eng. 7, 299–318. doi: 10.1007/BF02364118

PubMed Abstract | CrossRef Full Text | Google Scholar

Asner, L., Hadjicharalambous, M., Chabiniok, R., Peressutti, D., Sammut, E., Wong, J., et al. (2017). Patient-specific modeling for left ventricular mechanics using data-driven boundary energies. Comput. Methods Appl. Mech. Eng. 314, 269–295. doi: 10.1016/j.cma.2016.08.002

CrossRef Full Text | Google Scholar

Asner, L., Hadjicharalambous, M., Chabiniok, R., Peresutti, D., Sammut, E., Wong, J., et al. (2016). Estimation of passive and active properties in the human heart using 3D tagged MRI. Biomech. Model. Mechanobiol. 15, 1121–1139. doi: 10.1007/s10237-015-0748-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Augustin, C. M., Neic, A., Liebmann, M., Prassl, A. J., Niederer, S. A., Haase, G., et al. (2016). Anatomically accurate high resolution modeling of human whole heart electromechanics: a strongly scalable algebraic multigrid solver method for nonlinear deformation. J. Comput. Phys. 305, 622–646. doi: 10.1016/j.jcp.2015.10.045

PubMed Abstract | CrossRef Full Text | Google Scholar

Baillargeon, B., Rebelo, N., Fox, D. D., Taylor, R. L., and Kuhl, E. (2014). The living heart project: a robust and integrative simulator for human heart function. Eur. J. Mech. 48, 38–47. doi: 10.1016/j.euromechsol.2014.04.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Beyar, R., and Sideman, S. (1984). A computer study of the left ventricular performance based on fiber structure, sarcomere dynamics, and transmural electrical propagation velocity. Circul. Res. 55, 358–375. doi: 10.1161/01.RES.55.3.358

PubMed Abstract | CrossRef Full Text | Google Scholar

Booker, A. J., Dennis, J. E., Frank, P. D., Serafini, D. B., Torczon, V., and Trosset, M. W. (1999). A rigorous framework for optimization of expensive functions by surrogates. Struct. Optim. 17, 1–13. doi: 10.1007/BF01197708

CrossRef Full Text | Google Scholar

Chumarnaya, T. V., Alueva, Y. S., Kochmasheva, V. V., Mihailov, S. P., Revishvili, A. S., Tsyv'ian, P. B., et al. (2016). Features of the left ventricular functional geometry in patients with myocardial diseases with varying degrees of systolic dysfunction. Bull. Exp. Biol. Med. 162, 30–34. doi: 10.1007/s10517-016-3537-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Dokos, S., Smaill, B. H., Young, A. A., and LeGrice, I. J. (2002). Shear properties of passive ventricular myocardium. Am. J. Physiol. Heart Circul. Physiol. 283, H2650–H2659. doi: 10.1152/ajpheart.00111.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

Duchateau, N., Craene, M. D., Allain, P., Saloux, E., and Sermesant, M. (2016). Infarct localization from myocardial deformation: prediction and uncertainty quantification by regression from a low-dimensional space. IEEE Trans. Med. Imaging 35, 2340–2352. doi: 10.1109/TMI.2016.2562181

PubMed Abstract | CrossRef Full Text | Google Scholar

Fonseca, C. G., Backhaus, M., Bluemke, D. A., Britten, R. D., Chung, J. D., Cowan, B. R., et al. (2011). The Cardiac Atlas Project—an imaging database for computational modeling and statistical atlases of the heart. Bioinformatics 27, 2288–2295. doi: 10.1093/bioinformatics/btr360

PubMed Abstract | CrossRef Full Text | Google Scholar

Genet, M., Rausch, M. K., Lee, L. C., Choy, S., Zhao, X., Kassab, G. S., et al. (2015). Heterogeneous growth-induced prestrain in the heart. J. Biomech. 48, 2080–2089. doi: 10.1016/j.jbiomech.2015.03.012

PubMed Abstract | CrossRef Full Text | Google Scholar

GPy (2018). GPy: Gaussian Processes Framework in Python. Available online at: http://github.com/SheffieldML/GPy

Göktepe, S., and Kuhl, E. (2010). Electromechanics of the heart: a unified approach to the strongly coupled excitation–contraction problem. Comput. Mech. 45, 227–243. doi: 10.1007/s00466-009-0434-z

CrossRef Full Text | Google Scholar

Guccione, J. M., Costa, K. D., and McCulloch, A. D. (1995). Finite element stress analysis of left ventricular mechanics in the beating dog heart. J. Biomech. 28, 1167–1177. doi: 10.1016/0021-9290(94)00174-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Gültekin, O., Sommer, G., and Holzapfel, G. A. (2016). An orthotropic viscoelastic model for the passive myocardium: continuum basis and numerical treatment. Comput. Methods Biomech. Biomed. Eng. 19, 1647–1664. doi: 10.1080/10255842.2016.1176155

PubMed Abstract | CrossRef Full Text | Google Scholar

Gurev, V., Pathmanathan, P., Fattebert, J.-L., Wen, H.-F., Magerlein, J., Gray, R. A., et al. (2015). A high-resolution computational model of the deforming human heart. Biomech. Model. Mechanobiol. 14, 829–849. doi: 10.1007/s10237-014-0639-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Holmes, J. W., Borg, T. K., and Covell, J. W. (2005). Structure and mechanics of healing myocardial infarcts. Annu. Rev. Biomed. Eng. 7, 223–253. doi: 10.1146/annurev.bioeng.7.060804.100453

PubMed Abstract | CrossRef Full Text | Google Scholar

Holzapfel, G. A., and Ogden, R. W. (2009). Constitutive modelling of passive myocardium: a structurally based framework for material characterization. Philos. Trans. R. Soc. Lond. A Math. Phys. Eng. Sci. 367, 3445–3475. doi: 10.1098/rsta.2009.0091

PubMed Abstract | CrossRef Full Text | Google Scholar

Humphrey, J. D. (2002). Cardiovascular Solid Mechanics: Cells, Tissues, and Organs. New York, NY: Springer Science & Business Media.

Google Scholar

Jie, X., Gurev, V., and Trayanova, N. (2010). Mechanisms of mechanically induced spontaneous arrhythmias in acute regional ischemia. Circul. Res. 106, 185–192. doi: 10.1161/CIRCRESAHA.109.210864

PubMed Abstract | CrossRef Full Text | Google Scholar

Kerckhoffs, R. C. P., Neal, M. L., Gu, Q., Bassingthwaighte, J. B., Omens, J. H., and McCulloch, A. D. (2007). Coupling of a 3D finite element model of cardiac ventricular mechanics to lumped systems models of the systemic and pulmonic circulation. Ann. Biomed. Eng. 35, 1–18. doi: 10.1007/s10439-006-9212-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Lamata, P., Casero, R., Carapella, V., Niederer, S. A., Bishop, M. J., Schneider, J. E., et al. (2014). Images as drivers of progress in cardiac computational modelling. Prog. Biophys. Mol. Biol. 115, 198–212. doi: 10.1016/j.pbiomolbio.2014.08.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Land, S., Park-Holohan, S.-J., Smith, N. P., dos Remedios, C. G., Kentish, J. C., and Niederer, S. A. (2017). A model of cardiac contraction based on novel measurements of tension development in human cardiomyocytes. J. Mol. Cell. Cardiol. 106, 68–83. doi: 10.1016/j.yjmcc.2017.03.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Leong, C.-N., Al Abed, A., Lim, E., Lovell, N. H., Marasco, S., Hashim, S. A., et al. (2015). “Electromechanics modeling of the effects of myocardial infarction on left ventricular function,” in Engineering in Medicine and Biology Society (EMBC), 2015 37th Annual International Conference of the IEEE (Milan: IEEE), 5684–5687.

Google Scholar

Leong, C.-N., Lim, E., Andriyana, A., Al Abed, A., Lovell, N. H., Hayward, C., et al. (2017). The role of infarct transmural extent in infarct extension: A computational study: the role of infarct transmural extent in infarct extension. Int. J. Numer. Methods Biomed. Eng. 33:e02794. doi: 10.1002/cnm.2794

PubMed Abstract | CrossRef Full Text | Google Scholar

Lumens, J., Delhaas, T., Kirn, B., and Arts, T. (2009). Three-wall segment (TriSeg) model describing mechanics and hemodynamics of ventricular interaction. Ann. Biomed. Eng. 37, 2234–2255. doi: 10.1007/s10439-009-9774-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Marsden, A. L. (2014). Optimization in cardiovascular modeling. Annu. Rev. Fluid Mech. 46, 519–546. doi: 10.1146/annurev-fluid-010313-141341

CrossRef Full Text | Google Scholar

Marsden, A. L., Feinstein, J. A., and Taylor, C. A. (2008). A computational framework for derivative-free optimization of cardiovascular geometries. Comput. Methods Appl. Mech. Eng. 197, 1890–1905. doi: 10.1016/j.cma.2007.12.009

CrossRef Full Text | Google Scholar

Mazhari, R., Omens, J. H., Covell, J. W., and McCulloch, A. D. (2000). Structural basis of regional dysfunction in acutely ischemic myocardium. Cardiovasc. Res. 47, 284–293. doi: 10.1016/S0008-6363(00)00089-4

PubMed Abstract | CrossRef Full Text | Google Scholar

McCulloch, A. D. (1999). “Cardiac biomechanics,” in The Biomedical Engineering Handbook, 2nd Edn. 2 Volume Set (CRC Press).

Google Scholar

Miller, K., and Lu, J. (2013). On the prospect of patient-specific biomechanics without patient-specific properties of tissues. J. Mech. Behav. Biomed. Mater. 27, 154–166. doi: 10.1016/j.jmbbm.2013.01.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Miller, K., Wittek, A., Joldes, G., Horton, A., Dutta-Roy, T., Berger, J., et al. (2010). Modelling brain deformations for computer-integrated neurosurgery. Int. J. Numer. Methods Biomed. Eng. 26, 117–138. doi: 10.1002/cnm.1260

CrossRef Full Text | Google Scholar

Mojsejenko, D., McGarvey, J. R., Dorsey, S. M., Gorman, J. H., Burdick, J. A., Pilla, J. J., et al. (2015). Estimating passive mechanical properties in a myocardial infarction using MRI and finite element simulations. Biomech. Model. Mechanobiol. 14, 633–647. doi: 10.1007/s10237-014-0627-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Moulton, M. J., Hong, B. D., and Secomb, T. W. (2017). Simulation of left ventricular dynamics using a low-order mathematical model. Cardiovasc. Eng. Technol. 8, 480–494. doi: 10.1007/s13239-017-0327-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Moulton, M. J., and Secomb, T. W. (2013). A low-order model for left ventricle dynamics throughout the cardiac cycle. Math. Med. Biol. J. IMA 30, 45–63. doi: 10.1093/imammb/dqr024

PubMed Abstract | CrossRef Full Text | Google Scholar

Moulton, M. J., and Secomb, T. W. (2014). A low-order parametric description of left ventricular kinematics. Cardiovasc. Eng. Technol. 5, 348–358. doi: 10.1007/s13239-014-0191-9

CrossRef Full Text | Google Scholar

Nasopoulou, A., Shetty, A., Lee, J., Nordsletten, D., Rinaldi, C. A., Lamata, P., et al. (2017). Improved identifiability of myocardial material parameters by an energy-based cost function. Biomech. Model. Mechanobiol. 16, 971–988. doi: 10.1007/s10237-016-0865-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Nordsletten, D. A., Niederer, S. A., Nash, M. P., Hunter, P. J., and Smith, N. P. (2011). Coupling multi-physics models to cardiac mechanics. Prog. Biophys. Mol. Biol. 104, 77–88. doi: 10.1016/j.pbiomolbio.2009.11.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Pérez, E. A., Rojas-Solórzano, L. R., and Finol, E. (2016). Geometric predictors of abdominal aortic aneurysm maximum wall stress. Chem. Eng. Trans. 49, 73–78. doi: 10.3303/CET1649013

PubMed Abstract | CrossRef Full Text | Google Scholar

Piras, P., Teresi, L., Puddu, P. E., Torromeo, C., Young, A. A., Suinesiaputra, A., et al. (2017). Morphologically normalized left ventricular motion indicators from MRI feature tracking characterize myocardial infarction. Sci. Rep. 7:12259. doi: 10.1038/s41598-017-12539-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Pravdin, S. F., Dierckx, H., Katsnelson, L. B., Solovyova, O., Markhasin, V. S., and Panfilov, A. V. (2014). Electrical wave propagation in an anisotropic model of the left ventricle based on analytical description of cardiac architecture. PLoS ONE 9:e93617. doi: 10.1371/journal.pone.0093617

PubMed Abstract | CrossRef Full Text | Google Scholar

Radau, P., Lu, Y., Connelly, K., Paul, G., Dick, A. J., and Wright, G. A. (2009). Evaluation framework for algorithms segmenting short axis cardiac MRI. MIDAS J. 49, 1–7.

Google Scholar

Rasmussen, C. E., and Williams, C. K. (2006). Gaussian Processes for Machine Learning, Vol. 1. Cambridge: MIT Press.

Google Scholar

Ross, J., Sonnenblick, E. H., Taylor, R. R., Spotnitz, H. M., and Covell, J. W. (1971). Diastolic geometry and sarcomere lengths in the chronically dilated canine left ventricle. Circul. Res. 28, 49–61. doi: 10.1161/01.RES.28.1.49

PubMed Abstract | CrossRef Full Text | Google Scholar

Sankaran, S., and Marsden, A. L. (2011). A stochastic collocation method for uncertainty quantification and propagation in cardiovascular simulations. J. Biomech. Eng. 133:031001. doi: 10.1115/1.4003259

CrossRef Full Text | Google Scholar

Sellier, M. (2011). An iterative method for the inverse elasto-static problem. J. Fluids Struct. 27, 1461–1470. doi: 10.1016/j.jfluidstructs.2011.08.002

CrossRef Full Text | Google Scholar

Si, H. (2015). TetGen, a delaunay-based quality tetrahedral mesh generator. ACM Trans. Math. Softw. 41, 11:1–11:36. doi: 10.1145/2629697

CrossRef Full Text | Google Scholar

Sommer, G., Haspinger, D. C., Andrä, M., Sacherer, M., Viertler, C., Regitnig, P., et al. (2015). Quantification of shear deformations and corresponding stresses in the biaxially tested human myocardium. Ann. Biomed. Eng. 43, 2334–2348. doi: 10.1007/s10439-015-1281-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Stergiopulos, N., Meister, J. J., and Westerhof, N. (1996). Determinants of stroke volume and systolic and diastolic aortic pressure. Am. J. Physiol. Heart Circul. Physiol. 270, H2050–H2059. doi: 10.1152/ajpheart.1996.270.6.H2050

PubMed Abstract | CrossRef Full Text | Google Scholar

Streeter, D. D., and Hanna, W. T. (1973). Engineering mechanics for successive states in canine left ventricular myocardium. Circul. Res. 33, 639–655. doi: 10.1161/01.RES.33.6.639

CrossRef Full Text | Google Scholar

Suga, H., and Sagawa, K. (1972). Mathematical interrelationship between instantaneous ventricular pressure-volume ratio and myocardial force-velocity relation. Ann. Biomed. Eng. 1, 160–181. doi: 10.1007/BF02584205

PubMed Abstract | CrossRef Full Text | Google Scholar

Sundnes, J., Wall, S., Osnes, H., Thorvaldsen, T., and McCulloch, A. (2014). Improved discretisation and linearisation of active tension in strongly coupled cardiac electro-mechanics simulations. Comput. Methods Biomech. Biomed. Eng. 17, 604–615. doi: 10.1080/10255842.2012.704368

PubMed Abstract | CrossRef Full Text | Google Scholar

Towns, J., Cockerill, T., Dahan, M., Foster, I., Gaither, K., Grimshaw, A., et al. (2014). XSEDE: accelerating scientific discovery. Comput. Sci. Eng. 16, 62–74. doi: 10.1109/MCSE.2014.80

CrossRef Full Text | Google Scholar

Usyk, T. P., Mazhari, R., and McCulloch, A. D. (2000). Effect of laminar orthotropic myofiber architecture on regional stress and strain in the canine left ventricle. J. Elast. Phys. Sci. Solids61, 143–164.

Google Scholar

Valette, S., Chassery, J. M., and Prost, R. (2008). Generic remeshing of 3D triangular meshes with metric-dependent discrete Voronoi diagrams. IEEE Trans. Visual. Comput. Graph. 14, 369–381. doi: 10.1109/TVCG.2007.70430

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, H. M., Gao, H., Luo, X. Y., Berry, C., Griffith, B. E., Ogden, R. W., et al. (2013). Structure-based finite strain modelling of the human left ventricle in diastole. Int. J. Numer. Methods Biomed. Eng. 29, 83–103. doi: 10.1002/cnm.2497

PubMed Abstract | CrossRef Full Text | Google Scholar

Wenk, J. F., Eslami, P., Zhang, Z., Xu, C., Kuhl, E., Gorman, J. H., et al. (2011). A novel method for quantifying the in-vivo mechanical effect of material injected into a myocardial infarction. Ann. Thoracic Surg. 92, 935–941. doi: 10.1016/j.athoracsur.2011.04.089

PubMed Abstract | CrossRef Full Text | Google Scholar

Xi, J., Lamata, P., Niederer, S., Land, S., Shi, W., Zhuang, X., et al. (2013). The estimation of patient-specific cardiac diastolic functions from clinical measurements. Med. Image Anal. 17, 133–146. doi: 10.1016/j.media.2012.08.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, L., Sedai, S., Liang, X., Compas, C. B., Wang, H., Yushkevich, P. A., et al. (2015). “Multi-atlas label fusion with augmented atlases for fast and accurate segmentation of cardiac MR images,” in 2015 IEEE 12th International Symposium on Biomedical Imaging (ISBI) (New York, NY), 376–379.

Google Scholar

Young, R. J., and Panfilov, A. V. (2010). Anisotropy of wave propagation in the heart can be modeled by a Riemannian electrophysiological metric. Proc. Natl. Acad. Sci. U.S.A. 107, 15063–15068. doi: 10.1073/pnas.1008837107

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, X., Cowan, B. R., Bluemke, D. A., Finn, J. P., Fonseca, C. G., Kadish, A. H., et al. (2014). Atlas-based quantification of cardiac remodeling due to myocardial infarction. PLoS ONE 9:e110243. doi: 10.1371/journal.pone.0110243

PubMed Abstract | CrossRef Full Text | Google Scholar