ORIGINAL RESEARCH article
An Implementation of Patient-Specific Biventricular Mechanics Simulations With a Deep Learning and Computational Pipeline
- 1School of Biomedical Engineering and Imaging Sciences, King's College London, London, United Kingdom
- 2Auckland MR Research Group, University of Auckland, Auckland, New Zealand
- 3Department of Biomedical Engineering and Cardiac Surgery, University of Michigan, Ann Arbor, MI, United States
Parameterised patient-specific models of the heart enable quantitative analysis of cardiac function as well as estimation of regional stress and intrinsic tissue stiffness. However, the development of personalised models and subsequent simulations have often required lengthy manual setup, from image labelling through to generating the finite element model and assigning boundary conditions. Recently, rapid patient-specific finite element modelling has been made possible through the use of machine learning techniques. In this paper, utilising multiple neural networks for image labelling and detection of valve landmarks, together with streamlined data integration, a pipeline for generating patient-specific biventricular models is applied to clinically-acquired data from a diverse cohort of individuals, including hypertrophic and dilated cardiomyopathy patients and healthy volunteers. Valve motion from tracked landmarks as well as cavity volumes measured from labelled images are used to drive realistic motion and estimate passive tissue stiffness values. The neural networks are shown to accurately label cardiac regions and features for these diverse morphologies. Furthermore, differences in global intrinsic parameters, such as tissue anisotropy and normalised active tension, between groups illustrate respective underlying changes in tissue composition and/or structure as a result of pathology. This study shows the successful application of a generic pipeline for biventricular modelling, incorporating artificial intelligence solutions, within a diverse cohort.
Cardiovascular disease causes changes in cardiac anatomy, structure, and function—all resulting in changes to the active and passive biomechanics of the myocardium. However, it is difficult to assess intrinsic properties from imaging data alone. Patient-specific computational models can be used to simulate cardiac mechanics and measure quantities such as stress and strain and have the potential to augment current steps in therapy planning, allowing clinicians to test devices, such as left ventricular assist devices (Sack et al., 2018b), and therapies, such as septal myectomy (Huang et al., 2021). Personalised models can also be used to create “virtual cohorts,” running large-scale trials on large numbers of realistic heart models in concert with animal and human studies (Peirlinck et al., 2021).
Personalised models have been used to estimate both passive (e.g., Augenstein et al., 2006) and active (e.g., Marchesseau et al., 2013) parameters. Differences between global stiffness parameters have been identified between healthy and diseased cohorts through the use of patient-specific modelling (Hadjicharalambous et al., 2017; Wang et al., 2018). These passive parameters could be used as an additional diagnostic tool or to track disease progression. Additionally, estimation of heterogeneous stiffness parameters demonstrate the feasibility to identify local differences in tissue properties (Balaban et al., 2018), which can give an indication of regional changes. The optimisation of material parameters has been formulated as a nonlinear optimisation problem which aims to minimise an objective function based on the observation error, typically using displacement (Wang et al., 2018), strain (Augenstein et al., 2005; Wang et al., 2009), or geometric metrics (Nasopoulou et al., 2017). Filtering approaches, such as the use of Kalman filters, have also demonstrated robust and accurate estimation of passive parameters in the presence of noise (Xi et al., 2011). Regional contractility parameters estimated in personalised models have been shown to decrease in infarcted regions (Chabiniok et al., 2012). Although personalised modelling has been demonstrated to offer insights into intrinsic properties of the heart in health and disease, key challenges remain including automation of many cumbersome steps in model development as well as integration of key biomechanical information. For example, most studies have utilised manual segmentation for model development. Additionally, many studies developing personalised models have relied on data which is not typically acquired in a clinical scan (e.g., tagged MRI or intraventricular pressure measurements), thus limiting the size of their cohorts.
The process of generating patient-specific models was once a time-consuming task, requiring manual annotation and segmentation of images to construct an accurate geometric model (Heijman et al., 2008). The advent of and advances in machine learning have enabled automation of many of these tasks, with results ranging in accuracy and reliability (Henglin et al., 2017; Leiner et al., 2019). Deep learning is a subset of the machine learning field of techniques focusing on artificial neural networks which are constructed as deeply interconnected neural structures (Zhang et al., 2018). Neural networks have greater capacity to learn more complex problems than other machine learning methods with a greater ability to generalise to unseen data. These can be applied directly to the labelling of anatomical structures by assigning each pixel/voxel of an image a category probability which associates them with one or more structures. This allows the automation of cardiac segmentation such that an entire short-axis cine dataset can be labelled in seconds without manual initialisation or intervention, rather than hours (e.g., Bai et al., 2018; Chen et al., 2020). Cardiac segmentations can then be used to automatically calculate clinical metrics such as ejection fraction and long-axis strain (Ruijsink et al., 2019), and can be combined with other networks to perform disease classification (Martin-Isla et al., 2020) and feature detection (Bizopoulos and Koutsouris, 2018). In addition to expediting segmentation, trained neural networks can also lead to consistent and standardised results improving reliability and reproducibility. Neural network segmentations can then be used to automate the process of generating patient-specific geometric models.
In order to extend personalised modelling into the clinical domain, there is a need to develop a robust pipeline to not only generate models for diverse cardiac morphologies, but also to run biomechanical simulations using data acquired within a clinical scan. This study presents an AI-driven pipeline for the development of personalised biventricular mechanical models which were used to simulate passive and active mechanics. Novel boundary conditions, driven by neural network derived landmarks, were used to constrain valve motion and cavity volumes. The pipeline was tested in a diverse cohort which included healthy volunteers, patients with DCM and hypertrophic cardiomyopathy (HCM), using only a short-axis cine stack and three long-axis image planes, equivalent to images that would be collected in a standard clinical MR scan.
2. Materials and Methods
For all cases in this study, a balanced steady-state free precession sequence was used to collect cine images at short-axis slice locations and three long-axis imaging planes including two (2CH), three (3CH), and four-chamber (4CH) views. Between 20 and 40 images were acquired per cardiac cycle depending on the individual's heart rate. All images were acquired on a Philips Achieva 1.5 T scanner at St. Thomas' Hospital in London. Written informed consent was obtained from all participants prior to scanning. The study protocols for the DCM patients and healthy volunteers (study number 12/LO/1456) and HCM patients (study number 15/NS/0030) were approved by the London Bridge National Research Ethics Service. This initial study includes a cohort of patients with HCM (n = 4), DCM (n = 4), and healthy volunteers (n = 4). Patients with HCM demonstrated heterogeneous patterns of wall thickening, in keeping with the underlying diagnosis. The entire processing pipeline for each case is shown in Figure 3. Each block within the figure will be discussed in greater detail in the following sections.
2.1. Neural Network Image Labelling
Cine images were passed to two neural networks. The first of which labelled the left ventricular (LV) blood pool, LV myocardium and right ventricular (RV) blood pool in all short and long-axis images, whereas the second returned labels of 10 valve landmarks identifying leaflet insertion points in all long-axis images.
2.1.1. Cine Image Labelling
Full-cycle three-label segmentation was accomplished using a UNet-derived (Ronneberger et al., 2015; Kerfoot et al., 2018, 2019) neural network. Left-ventricular blood pool, left-ventricular myocardium, and right-ventricular blood pool were identified by this network by analysing each two-dimensional slice from a full short-axis stack individually.
The network architecture is composed of a stack of blocks incorporating the encode and decode paths of the UNet structure (Figure 1). Data flows through the encode side on the left where it passes through a residual unit (He et al., 2016) of convolution/normalisation/activation layers. The output from this unit passes to the next layer in the stack, which is either a further layer of such encode/decode pathways or a final residual unit. The data from the encode path is concatenated with the output from the layer below before being passed through another residual unit in the decode side.
Figure 1. The segmentation network is implemented as a stack of blocks illustrated here. The encode and decode paths along with the skip connection are defined in the same block. The “Next Layer” is either another such block or the bottom encoding block comprised of convolution/normalisation/activation sequences. The overall structure of the network is shown on the right, with the dimensions of tensors passing between layers given relative to an input of shape (1, N, W).
The dataset used for training consisted of 9,095 segmented MR short-axis images (Kerfoot et al., 2019). These were derived from the ACDC challenge dataset (Bernard et al., 2018) of 100 cases and 175 UK Biobank healthy cases. Of the latter, an expert clinician at St. Thomas' Hospital in London segmented 100 healthy cases, 50 cardiomyopathy cases, and 25 randomly selected cases that exhibited sufficient image quality for use as input. Additionally, 215 cases were acquired on a 1.5 T Philips Ingenia scanner at St. Thomas' Hospital in London, and 116 cases from a Siemens Trio 3T scanner (Siemens Healthineers, Erlangen, Germany), and were also segmented by an expert clinician at St. Thomas' Hospital. These cases consisted of healthy volunteers, HCM patients, and patients with cardiac resynchronisation therapy (CRT).
The network was trained for 10,000 iterations. For each iteration, a mini-batch was created by selecting 250 randomly selected images from the dataset. A random selection of flip, transpose, 90° rotation, shift and non-rigid deformation operations were applied to the image and segmentation pairs. The loss function used was a simple Dice loss (Dice, 1945).
2.1.2. Valve Landmark Identification
Landmark coordinates in the three long-axis views were used to identify the locations of the leaflet insertions into the myocardium. Ten landmark locations in total were estimated: six mitral valve locations (two from each view), two aortic locations in the three-chamber view, and two tricuspid locations in the four-chamber view. These landmarks were estimated using a convolutional neural network implemented as a regression from two-dimensional images to a landmark coordinate array (Kerfoot et al., 2021). See Figure 3: Valve Landmark Identification. Briefly, the network was trained on 8,574 long-axis images collected from HCM (n = 3,069) and myocardial infarction (MI, n = 5,505) patients. Further details of the dataset used for training can be found in Kerfoot et al. (2021).
Figure 2C illustrates the general architecture of the network composed of a sequence of densely-connected blocks of convolutions. The output data from these blocks is then passed to a series of small neural networks trained to recognise the 10 different landmark coordinates. Having condensed the information from the input image to a deep representation, each sub-network learns to recognise which view is represented and determine a location from this representation. From each long-axis image, all ten landmarks are identified. However, the landmarks which do not occur in the input image are inferred to be in the top-left corner at coordinate [0,0].
Figure 2. The valve estimation network is composed primarily of a series of densely-connected convolutional layers. Each dense block is composed of three residual units containing 2D convolutions using progressively larger dilation rates. A final convolution reduces the spatial dimension of the volume by 2. The regression network is implemented as a sequence of densely-connected blocks followed by a series of small fully-connected networks relating the final output volume to each landmark coordinate. (A) Residual unit, (B) dilated dense block, (C) network definition.
Figures 2A,B illustrates the architecture of the densely-connected blocks (Huang et al., 2017). Within each block is a residual unit composed of two sets of convolution/normalisation/regularisation layers. The dilation value for the convolutions increments for each succeeding unit, which allows convolutions to recognise features of different scales in the input volume. The output from each unit is concatenated with outputs from previous units. This combined volume is used as the input to the next unit. All such outputs, plus the original input, are concatenated into the final output volume. A final convolution/normalisation/regularisation reduces the output volume in the spatial dimensions by a factor of two.
During training, data augmentation was applied to the images from the manually-annotated dataset. A randomised combination of flip, transpose, zoom, rotate, shift, and non-rigid deformation operations were applied to the image and ground-truth landmark pairs to be fed into the network during training. The images were further augmented with added noise, smooth image intensity variation and k-space dropout to simulate a poor-quality acquisition. The objective of these augmentations was to vary the data the network is trained with to reduce overfitting and improve its generalisation to unseen image types.
2.1.3. Label Quality Control
For the short and long-axis segmentations, labels were cleaned (a) by removing labelled regions with fewer than 50 pixels, disregarding improperly labelled “islands” far from the heart as well as (b) filling holes in the labelled regions. Since valve landmarks were identified for each 2D long-axis image independently (not incorporating temporal continuity), an additional step was implemented to automatically identify landmarks which were incorrectly labelled in order to omit these points. Then, a linear interpolation step was used to interpolate missing points before applying a low-pass filter to temporally smooth landmark displacements. Valve landmarks were used both as input to the model fitting step as well as boundary conditions to constrain valve annuli motion throughout simulations of the cardiac cycle.
2.2. Segmentations to Models
Short-axis alignment was performed using the IRTK toolbox (Schnabel et al., 2001) which applies a rigid transformation to individual short-axis planes in order to optimise the overlap between short-axis masks and a model template. Starting with the rigid registration tool's default initialisation, the short-axis images are rigidly moved in the in-plane dimension to account for misalignment during acquisition.
Subsequently, long-axis images were rigidly registered to the short-axis aligned images using the rigid registration algorithm in the IRTK package, also using the tool's default initialisation. Dice scores along the line of intersection between each short and long-axis mask were used to (a) automatically determine which short-axis slices would be used for the model fitting by omitting slices with a dice score <0.5 and (b) to assign weights to each contour point based on their overlap with other data. In this way, long-axis slices which were poorly registered to the short-axis data, even after running the registration step, did not skew or greatly impact the final fitted model.
A biventricular template was then fitted to the segmentations using the two-step iterative method developed in Mauger et al. (2018). In order to do this, contours were automatically generated from the short and long-axis labels (i.e., LV endocardium, RV septum, RV free wall, epicardium, RV insertion points, apex, etc.) in order to fit model surfaces to the contour points. Locations of the mitral, tricuspid, and aortic valve annuli were obtained from the annotated valve landmarks. Due to a lack of segmentation of the RV myocardium, RV epicardial contours were automatically generated by projecting the RV free wall contours in the normal direction at a distance of 3 mm. Briefly, a series of stiff linear least squares fits with a high D-affine regularisation weight was performed to provide an adequate first solution. For each iteration, the Jacobians on 4 × 4 × 4 Gaussian quadrature points were calculated. If all were positive, the model was updated, the regularisation weight was decreased and another iteration was performed. If not, the model was not updated and another optimisation step was performed using diffeomorphic constraints based on the magnitude of the displacement. Models were fit to segmentations at all frames of the cardiac cycle. Surface meshes were used to construct cavity volume curves throughout the cardiac cycle as well as quantify metrics such as wall thickness and ejection fraction.
From the fitted surface meshes, tetrahedral meshes were generated for the end-systolic time point using SimModeler (Simmetrix1). Mesh metrics, including number of nodes and element quality, can be found in Supplementary Table 1. Biventricular fibre fields were created using a rule-based method adapted from Doste et al. (2019) and Bayer et al. (2012). Fibre angles varied from −60 to 60° and −25 to 90° from the epicardium to endocardium in the LV and RV, respectively. Fibre angles at the valve annuli were determined based on high-resolution DTI measurements from ex-vivo porcine hearts. Specific angles at each boundary can be found in Supplementary Table 1. An example fibre field can be seen in Figure 3, Rule-based Fibres.
Figure 3. Short and long-axis cine MR images are simultaneously fed into two neural networks, one for labelling the LV blood pool (red), LV myocardium (green), and RV blood pool (blue) and the second labelling ten valve landmarks throughout the cardiac cycle. The segmentations are converted to labelled contours and a biventricular template surface mesh is fitted to the labelled contours. Volumes, derived from the network generated cavity labels, as well as valve annuli motion are used as boundary conditions in the biomechanical simulations. Passive parameters are optimised by minimising the difference between the model and imaged geometries at end-diastole.
2.3. Biventricular Modelling
The personalised mechanical models were solved using energy potential minimisation, following (Asner et al., 2017; Hadjicharalambous et al., 2017). In brief, the myocardium is defined by the reference domain with initial coordinates X ∈ Ω0. The biventricular domain, Ω0, consists of boundaries on the endocardial sides of the LV and RV (denoted and ), the wall marking the rings for the mitral (), aortic (), tricuspid (), and pulmonary valves (), as well as the epicardium (). The orientation of local tissue microstructure across the myocardial wall is given by the fibre, sheet and sheet normal vector fields, (f0, s0, n0). Similarly, at each valve boundary, a circumferential vector field is defined (denoted for , ) which describes the local orientation of connective tissue that comprises each valve orifice. Finally, to enable variations between the LV/septum and the RV, we define a labelling field, ϕ, where ϕ = 1 in the LV/LV septum and ϕ = 0 in the RV/RV septum.
For simulating myocardial function, imaging data is extracted to describe functional changes through time. The change in LV and RV luminal volumes is extracted from images and interpolated to provide describing the mean volume trace as computed over a truncated region of each endocardial lumen, . The truncation planes are similarly defined by normal vectors across both LV and RV lumens. The pressure is given over the cardiac cycle by and can be defined either via invasive measures, coupled via a full-circulation model (Arts et al., 2005), or estimated from noninvasive data (Asner et al., 2015). Finally, the motion of each valve plane is encapsulated by the estimated motion of the centre of mass, , interpolated over time for each valve, .
As the biventricular model deforms, the physical domain at time t, Ω(t), is described using coordinates of its current position x = X + u(t), where u denotes displacement. Typically the displacement is used to describe the deformation gradient tensor F = ∇0u + I, its determinant J = detF > 0, as well as the material stretch described by the right Cauchy-Green strain tensor C = FTF. The displacement of the heart is solved by considering either the quasi-static (Asner et al., 2017) or dynamic (Chabiniok et al., 2012; Sermesant et al., 2012) principle of virtual work, with the additional state variables of pressure (p), activation state in LV / RV (αlv, αrv), and the forces present at each valve plane (). In this study, at each time point, t ∈ [0, T], we seek to find the state variables , , αlv(t), αrv(t) ∈ ℝ, and λmv(t), …λpv(t) ∈ ℝ3 satisfying the quasi-static virtual work equation,
The virtual work equation describes the internal myocardial stresses and balance of volumetric/pressure change (), the stresses induced by the collagenous valve tissue (), the internal pressure exerted by the blood (), the constraint on chamber volumes (), and the added forces required to ensure motion of the valve orifices (). The specifics of these terms are detailed below.
The internal stresses () are given by the first Piola-Kirchhoff tensor, Pmyo, which is described by the hyperelastic-strain energy, Ψ, that can be broken into passive, volumetric, and active strain energy components,
and f = Ff0 describes the deformed fibre direction. Here the passive component (Equation 3a) follows the reduced form of the Holzapfel-Ogden model (Holzapfel and Ogden, 2009; Hadjicharalambous et al., 2014a, 2017; Asner et al., 2015) adapted for appropriate use within a nearly-incompressible framework (Nolan et al., 2014) (though numerous alternative models exist, see Chabiniok et al., 2016). The parameters a0 and af linearly scale the stiffness of the ground substrate and fibre direction, respectively, and have units of stress whereas b0 and bf scale the exponential behaviour of the isotropic and fibre components and are unitless. is the isochoric form of the first invariant of the right Cauchy-Green strain tensor (I1 = C : I) and If is the fibre pseudo-invariant If = C:(f0 ⊗ f0). To ensure unique parameters, the number of personalised parameters was reduced to two: a and af. The values of b0 and bf were each set to 5.0 to ensure physiological pressure-volume response (Hadjicharalambous et al., 2014a).
The resulting stress from volumetric effects (Equation 3b) results from the nearly-incompressible strain energy,
which also governs the relation between volume change and hydrostatic pressure (second part of the ). The parameter, K, denotes the bulk modulus of the tissue (in this study K = 1, 000 kPa).
The active stress, given in Equation (3c), defines the amount of contraction as well as the length dependent mechanisms (Kerckhoffs et al., 2003). Here, stresses were also applied both along fibres as well as across fibres based on known myofibre dispersion (Tangney et al., 2013; Krishnamurthy et al., 2016). Note that the active scalings, αlv, αrv, are dotted with the region identifier, ϕ, in order to selectively activate LV and RV chambers. While regional activation can be defined based on eikonal activation times (Tomlinson et al., 2002) or monodomain/bidomain simulations (Potse et al., 2006); here, the contraction of the chambers was approximated by uniform contraction parameters, αlv and αrv.
Additional stresses were added along the surface of each valve orifice (), reflecting the fact that each valve annulus is comprised of thin cartilaginous tissue (Hamdan et al., 2012; Gunning and Murphy, 2014). This tissue, comprised of circumferential collagen fibres, is extremely flexible but exhibits strong resistance to annular dilation (e.g., stretch of the collagen fibres). This was incorporated into the model by adding stresses, , applied over each annular plane, , where
Here, describes the deformed circumferential direction of collagen fibres in the kth−annulus and is the pseudo-invariant along collagen fibres. As the stresses induced are exerted along an extremely thin area, and were principally oriented along fibres, the added stresses were incorporated over the annular surfaces. The parameter c1 accounts for the collagen stiffness scaled by the thickness while the parameter c2 allows for exponential growth in the fibre stresses. Here, a value of c1 = 0.1 kPa and c2 = 0.5 were selected for all valves based on achieving a consistent qualitative annular dilation as typically found in vivo and were unchanged across patients (assuming the collagenous structures around the valves was consistent).
Instead of using parameter estimation techniques (Chabiniok et al., 2012; Marchesseau et al., 2013; Asner et al., 2015) to determine the activation of the myocardium, the LV/RV activation was solved for as part of the forward model problem. In this context, the active parameters αlv, αrv act as Lagrange multipliers with the constraint held being that both chambers follow the volume trends observed in the data (). Here, the first term provides the model predicted volumes which must be equal to the volumes prescribed by Vk(t) (where and β(X) is a binary variable taking the value of 1 below and 0 above the truncation plane) (Asner et al., 2017). As the pressure at each time point is given (), the active tension scalings are found which enable matching between the model/data. Note, for consistency and stability, the applied pressure Pk(t) must be greater or equal to the passive pressure at the specified volume Vk(t).
Valve plane motion was prescribed () using the valve landmark displacements, , extracted from the points predicted by the neural network (discussed in section 2.1). For each valve, the average position over the cardiac cycle was enforced using Lagrange multipliers, λk. The pulmonary annulus is not visible in any of the long-axis images acquired in this study. It can be viewed in a right ventricular outflow tract (RVOT) view, which is not always acquired in clinical scans. Therefore, in this study, an average displacement, computed from the displacements of the other three valves, was applied to the pulmonary valve. The introduce the multipliers (that can be thought of as reference tractions) which constrain the motion of the centre of mass to move as observed in the data.
Since the unloaded state of myocardium is unknown, the end-systolic geometry was used as the reference geometry. Some studies have used inverse methods to estimate the reference geometry (Krishnamurthy et al., 2013; Wang Y. et al., 2020). These methods are dependent on the choice of material law, constitutive parameters and boundary conditions. An analysis illustrating the impact of these choices and ramifications of boundary conditions is presented in Hadjicharalambous et al. (2021).
Personalised models were solved in a finite element framework with displacements and pressure defined using linear ℙ1 elements, see Supplementary Material for further details. Endocardial and valve Lagrange multipliers were scalars. All problems were solved in Heart, a multi-physics finite element solver (Lee et al., 2016).
2.4. Diastolic Inflation and Passive Parameter Estimation
Estimation of both a0 and af is not feasible using displacements alone when driving simulations with cavity volumes. However, due to the linear parameter dependence of a0 and af in the reduced Holzapfel-Ogden law, both passive parameters scale linearly with pressure. Therefore, each model was personalised by estimating, γ = a0/af, describing the anisotropy of the tissue from displacements and then using an end-diastolic pressure value to obtain a0 and af. To do this, simulations of diastolic inflation were first run starting from the end-systolic geometry, prescribing cavity volumes, and valve motion. This was done by solving Equation (1), assuming αlv, αrv were zero, and considering Plv, Prv as state variables (Asner et al., 2017). In this first phase, the LV/RV volumes were inflated to their end diastolic state, after which the volume was kept constant and values of γ were swept between 1.0 and 0.1. Practically, this was done by setting and to 1.0 kPa during inflation. Then, during the sweep, was kept constant and asim was varied. Absolute values of a0 and af were then retrieved by scaling them by the ratio between an end-diastolic pressure value appropriate for each patient group () and the simulated end-diastolic pressure (). The for each individual case was then found by using the values of a0 since the stiffness scales with both EDP in the left and right ventricles.
For each value of γ, the objective function, was calculated as the root mean square of the distance between contour points obtained from the neural network labels and the deformed model surface, Γlv (see Equation 6). The objective function utilised only contour points from the LV epicardium, LV endocardium and RV septal wall, omitting RV free wall points.
Since all data used were acquired from a standard clinical scan, no catheter pressure measurements were acquired. Filling pressures have previously been estimated from the E/A ratio measured from echocardiography (Nagueh et al., 1997) which requires blood flow measurements through the mitral valve. In the absence of echo and 4D flow MRI data, end-diastolic and end-systolic pressure values were found from literature in studies which obtained invasive catheter pressure measurements from within the LV in each patient group (see Table 1). Additionally, normal end-diastolic and end-systolic pressures were taken from Klingensmith et al. (2012). Taking the mean (weighted by sample size) of pressure values from literature, LV end-diastolic pressures () were assigned to be 8 mmHg (1.1 kPa), 20.2 mmHg (2.7 kPa), and 24.2 mmHg (3.2 kPa) for the healthy volunteers, DCM patients and HCM patients, respectively. Similarly, values were set to 120.0 mmHg (16.0 kPa), 120.0 mmHg (16.0 kPa), and 183.1 mmHg (24.4 kPa) for each group, respectively. A representative pressure trace (Russell et al., 2012) was scaled to group pressure values at end-diastole and end-systole for both the LV and RV. Then, each segment of the pressure trace (i.e., ED to eIVC, eIVC to ES, ES to eIVR, eIVR to diastasis, and diastasis to ED) was temporally scaled for each individual based on valve opening and closing times in the cine images.
Table 1. Left ventricular pressure measurements from literature denoting mean pressure ± one standard deviation as well as sample sizes in each study.
3.1. Neural Network Segmentation and Landmark Labelling
All short and long-axis images were manually segmented at the end-diastolic state in order to measure accuracy of the network. Boxplots in Figure 4 plot dice scores measuring similarity between manual and neural network segmentations for each label and group. Results show that the largest errors occur in the segmentation of the myocardium, with HCM cases having the lowest dice scores for this label.
Figure 4. Boxplots illustrate dice scores for each label: LV blood pool (LV), LV myocardium (Myo), and RV blood pool are shown for the healthy volunteers (V), DCM patients (D), and HCM patients (H). The centre line of each boxplot represents the median and the whiskers denote the 25th and 75th percentiles.
Errors (in mm) between predicted and manually annotated valve landmarks are shown in Figure 5 and were generally within 3 mm of the manually annotated position (~2 pixels). Six landmarks are labelled for the mitral valve (in the 2CH, 3CH, and 4CH images) whereas only two landmarks are labelled for the aortic and tricuspid valves each. Generally, landmarks were more accurately identified in the HCM group.
Figure 5. Boxplots of valve annotation errors for all 10 valve landmarks are shown for selected healthy volunteers (V), DCM patients (D), and HCM patients (H) in which valve landmarks were manually identified over the entire cardiac cycle. Each boxplot represents errors throughout the cardiac cycle. The centre line of each boxplot represents the median and the whiskers denote the 25th and 75th percentiles. The 10 valve landmarks correspond to those shown in Figure 3, Valve Landmark Identification.
3.2. Model Fitting and Geometric Measurements
The model fitting algorithm was able to accurately represent the various morphologies within the diverse cohort. Figure 6 illustrates models fit to a healthy volunteer, a DCM patient and an HCM patient with septal hypertrophy at the end-systolic time point. Clinical metrics, such as end-diastolic volumes, end-systolic volumes and wall thickness are reported for each case in Table 2. Compared to the healthy volunteers, the DCM cases have a higher mean LV end-diastolic (EDV) and end-systolic volume (ESV). Conversely, HCM patients have both reduced EDV and ESV in the left ventricle. Healthy volunteers and HCM patients demonstrate LV ejection fractions (EF) in a normal range (50% < EF < 70%) whereas DCM patients exhibit a depressed EF by definition. In all groups, RVEF values fell between (37.9% < EF < 55.9%) without discernible differences between group means. Wall thickness was greatest in the HCM cohort and showed minimal changes in the DCM group between ED and ES.
Figure 6. Three representative cases for the healthy volunteer, DCM and HCM groups are shown below with models fit to the neural network segmentations. Model surfaces at end-systole (purple) are overlayed on a single long-axis and short-axis image.
Table 2. Functional and geometric indices: end-diastolic volume (EDV), end-systolic volume (ESV), ejection fraction (EF), wall thickness (WT).
3.3. Passive and Active Parameterisation
The fibre stiffness ratio, γ, was estimated for all cases using parameter sweeps and the optimal values are listed in Table 3, along with the end-diastolic pressure values used to scale a0 and af to meaningful stiffness estimates. Values of γ close to 1 indicate that the material is more isotropic whereas a value of 0.1 would indicate a highly anisotropic material. The value of γ also influences the final inflated geometry where larger values result in a more spherical shape. The mean value of γ for the volunteers is less than that for the DCM (p = 0.2) and HCM (p = 0.05) patient groups, indicating that healthy myocardium may be slightly more anisotropic in this small cohort. Additionally, larger γ values are in line with more spherical shapes observed in DCM hearts.
Two active tension scaling parameters, for the LV and RV, were estimated throughout the cardiac cycle for each case. Normalised time-to-peak ( and ) as well as peak scaling parameters ( and ) are listed in Table 3. It should be noted that the traces of αlv and αrv are dependent on volume changes as well as pressures. Since higher ESP values were assigned in the HCM cases, it can be seen that the peak values of and are greater in this group. There were no significant differences between time to peak activation. Figure 7 shows mean active fibre stress over the cardiac cycle in the LV and RV as well as panels showing fibre stress patterns throughout the model for a single case (v1) at three time points during systole. The highest stresses are seen near the base of the left and right ventricles. Fibre stretch with respect to the end-diastolic state is plotted for a representative case illustrating model deformation and regional stretch patterns over the cardiac cycle (Figure 8). A bullseye plot of fibre stretch at end-systole illustrates that the largest values (<0.55) are seen in the LV free wall whereas fibre stretch is restricted in the basal septal region near the valves.
Figure 7. Mean active fibre stress over the cardiac cycle in both the LV and RV for a single case (v1) illustrating the active stress distribution over the entire heart at three points during active contraction: early systole (t = 50 ms), peak active contraction (t = 230 ms), and end-systole (t = 360 ms). The bullseye plot shows the regional distribution of mean active fibre stress over the 17 AHA regions at end-systole.
Figure 8. Mean fibre stretch over the cardiac cycle is shown for both the LV and RV for a single case (v1). Fibre stretch over the cardiac cycle is also plotted for nine time points with the reference state model (ES) shown as a wireframe mesh. The bullseye plot illustrates differences in regional stretch at the end-systolic state in the 17 AHA regions.
Both mean active fibre stress and mean fibre stretch in the LV are plotted for 16 AHA segments for all cases in Figure 9, illustrating group differences. Peak fibre stretch is smaller in DCM cases when compared to healthy volunteers in 12 out of 16 AHA regions (p < 0.05). Regional fibre stretch demonstrates that, in some DCM cases, some segments exhibit further stretching of fibres (values >1.0) in early phases of systolic contraction. Circumferential and longitudinal stretch, common clinical metrics, are also plotted for 16 AHA segments for all cases in Supplementary Figure 8. The mean for each segment and group are listed in Supplementary Table 4.
Figure 9. (A) Mean active fibre stress and (B) fibre stretch over the cardiac cycle in 16 AHA regions of the LV for healthy (black), DCM (blue), and HCM (red) groups.
The primary goal of this study was to implement a pipeline for running full-cycle simulations using personalised biventricular models generated entirely from neural-network labels. In this process, no manual segmentation was done for the cases presented, other than for analysis of the accuracy of each network. The time and computational resources used for the pipeline are given in Supplementary Table 5 and demonstrate a clear advantage over manual methods. Additionally, all data used in this study was obtained using sequences common to any standard clinical MR scan. Due to its ability to be applied to diverse datasets, this pipeline could be used to develop an in silico cohort based on true patient data. This virtual cohort would be invaluable for testing novel therapies and devices alongside human and animal studies. Additionally, personalised metrics obtained from the models (e.g., anisotropy, material stiffness) could be further used to either classify patients or mark disease progression. However, larger sample sizes are needed in order to better understand differences between patient classes. Additional data, where available, could be used to augment the robustness of the personalised models, such as the use of tagged MR data for passive parameterisation (Asner et al., 2015). The use of a biventricular template along with fitting weights assigned to contours based on data fidelity enabled the generation of high-quality meshes suitable for biomechanical simulations with minimal user intervention. Neural-network identified leaflet landmarks were used to prescribe average valve motion on each valve in the model, allowing for physiological basal motion of both ventricles.
4.1. Neural Networks
The neural network was able to accurately label the left and right ventricles from standard clinical images in a diverse cohort. The network captured the varied morphology of heart shapes in both DCM and HCM patients. Dice scores for labelling the LV blood pool were comparable to those from other segmentation networks (Wang et al., 2021) and the RV dice scores demonstrated greater accuracy than previous studies (Luo et al., 2016; Tran, 2016). However, the largest errors arose in labels of the myocardium. Although a comparison to inter-observer error was not done as part of this study, previous groups have compared annotations from multiple observers using the UK Biobank (Attar et al., 2019) as well as ACDC (Bernard et al., 2018) data sets. Similar to results shown in Figure 4, inter-observer errors for myocardium are greater than those for both the LV and RV blood pools in both data sets. Dice scores observed in this study are higher than the inter-observer dice scores reported for the LV and RV blood pools in Attar et al. (2019). One possible reason for the larger errors in the segmentation of the myocardium could be due to its annular shape which has a larger perimeter. Any equal overlap shifts would produce a greater error when compared to any shift in segmentation of the blood pools.
The second neural network labelled 10 different valve landmarks in each long-axis image to within 2–3 pixels of accuracy. These errors are similar to those encountered using common tracking algorithms and the method does not require manual initialisation (Kerfoot et al., 2021). Of the 10 valve landmarks, four demonstrated lower neural network predicted errors than interobserver errors (Kerfoot et al., 2021). The error can vary considerably throughout the cycle and between patients as each image is treated individually—i.e., no temporal consistency is taken into account in the neural network. There was no single patient that performed worse than others. Higher errors seen in the identification of landmark 8, on the septal side of the aortic valve can be attributed to image artefacts during systolic blood flow through the aortic outflow tract. Although improvements can be made in future work to increase the accuracy of both neural networks, the study focused on demonstrating their utility in driving model generation and biomechanical personalisation for a diverse set of patients.
In order to have a pipeline that is robust to the presence of noise and artefacts in the imaging data, the neural network training process introduces noise to the images in various ways (e.g., dropout in k-space) so that it learns to account for the noise it may encounter in the imaging data. Additionally, due to the use of a model template fit to all short- and long-axis contours simultaneously, the pipeline is robust to the presence of a single or even multiple poor-quality images within a dataset. This, however, can result in a smooth surface that does not conform to small features. To demonstrate the robustness of the neural network segmentation and resulting pipeline, selected poor-quality images are shown in Figure 10 which were part of the 12 datasets used in this study. The poor image quality resulted in deteriorated segmentations. However, the final fitted model, which also takes into account all short-axis information, produced an adequate estimate of the long-axis shape. The pipeline could be improved by further augmenting the neural network with images that mimic typical artefacts found in MR images.
Figure 10. Selected long-axis images which were part of the 12 cases (h1, h4, d2, d4) used in this study which demonstrate poor image quality due to imaging artefacts. In the first row, the LV blood pool and myocardial segmentations obtained using the neural network are overlain on top of each image to qualitatively show the impact of the image quality on the network segmentation. In the second row, the final fitted model surface is shown on top of the image.
4.2. Clinical Metrics
Beyond improving the generation of computational models, trained neural networks provide a mechanism for automatically characterising common clinical metrics. In the DCM group, the mean EDV and ESV were greater than those measured in the healthy volunteers. Similarly, the mean LVEF was less than that in both the healthy group, marking the deteriorated contractile and diastolic filling function typically clinically associated with DCM (Rihal et al., 1991). Conversely, the mean EDV and ESV values were slightly smaller in the HCM group when compared to the healthy volunteers. As commonly reported in HCM patients, the LVEF in this group was slightly elevated when compared to the healthy group (Haland et al., 2017). HCM patients exhibited greater wall thickness at both end-diastole and end-systole when compared to both the DCM patients and healthy volunteers. Typically, HCM is characterised by a wall thickness >12 mm during diastole. Although the wall thickness values shown in Table 2 report mean values <12 mm, isolated hypertrophic regions in each patient demonstrate areas of hypertrophy >12 mm. Regional plots of wall thickness averaged over each cohort are shown in Supplementary Figure 4 using the 17-segment AHA model. In each of these three groups, no significant differences were observed in the mean RVEF. However, the HCM patients demonstrated lower values of EDV and ESV in the RV than the other two groups. This pipeline has demonstrated the ability to rapidly generate common clinical metrics such as EF and wall thickness as well as cavity volumes over the entire cardiac cycle without the need for manual processing. Aside from using these values in clinical decision making, they can also be used as input into personalised biomechanical models.
4.3. Valve Motion
This study presents a novel means of constraining valve motion. Displacement was prescribed to valve centroids based on the motion of the identified landmarks from the neural network. In other cardiac modelling studies, basal motion is often constrained by restricting longitudinal motion (e.g., Sack et al., 2018b; Finsberg et al., 2019; Wang Z. J. et al., 2020) or applying an average motion measured from imaging data (e.g., Hadjicharalambous et al., 2017). In another study, basal motion was constrained by tethering the pulmonary outflow tract to a fixed point (Sack et al., 2018a). In truncated models, without the inclusion of anatomical landmarks, tagged magnetic resonance imaging (MRI) data is necessary to measure longitudinal motion, which may not be available in all clinical scans. The use of a biventricular model with all four valve annuli along with the neural network-defined leaflet insertion points allowed for the integration of longitudinal motion measured from imaging data into the computational model. Further, to impose a constraint similar to the stiff valve annulus, an additional stiffness term was used to restrict annular dilation.
4.4. Model Personalisation
Integration of imaging data with personalised biomechanical models enables estimation of intrinsic material stiffness parameters, providing important information about the mechanical state of the myocardium. In this study, we focused on determination of bulk and fibre material parameters, fit by adjusting their ratio, γ. The mean value of γ, which is independent of pressure, was smaller for the healthy volunteers than those estimated for the DCM and HCM patient groups. These weak differences may indicate that the myocardium in healthy individuals is more anisotropic than in pathological hearts. However, these differences were not statistically significant. A power analysis suggests that using eight samples would enable these differences to reach significance. In order to demonstrate statistically significant differences between αmax values in the LV for DCM and HCM groups, 11 samples would be needed. To distinguish differences between the time to peak contraction in diseased patients and healthy volunteers, 46 samples are needed. Therefore, future studies will aim to expand the sample size to demonstrate the pipeline's utility in providing metrics which distinguish between patient groups.
Simulation outcomes such as strain are relatively independent of the estimated value of γ since the simulations are driven by cavity volumes. However, stress would be more affected by a change in passive material properties. Fibrosis, common in both DCM and HCM patients (Aurigemma et al., 2006; Marian and Braunwald, 2017), results from the growth of collagen within cardiac tissue and may impact tissue anisotropy. The estimated value of γ is also strongly influenced by the angles defined in the rule-based fibre field (Asner et al., 2015; Hadjicharalambous et al., 2017; Campos et al., 2020). As rapid in vivo diffusion tensor MRI sequences improve (Stoeck et al., 2018), personalised fibre fields will augment the robustness of the presented pipeline. Methods of using low-resolution in vivo data along with statistical models of population fibre fields may provide a new means of personalisation without significantly adding to the clinical scan time (Stimm et al., 2021).
The objective function used for determining the optimal value of γ utilised LV contour points only. The RV deformation is impacted significantly by epicardial boundary conditions due to its thin wall. Various approaches have been used in previous studies to constrain epicardial dilation, such as a spring force acting in the normal direction (Levrero-Florencio et al., 2020; Strocchi et al., 2020) or parallel spring and dashpot forces (Pfaller et al., 2019). However, there remains a lack of clear consensus on the role of the pericardium in restricting myocardial motion and whether or not the inclusion of epicardial constraints improves model personalisation. Therefore, simulations in this study were run without the addition of boundary conditions on the epicardium. However, without these constraints, the right ventricular deformation did not sufficiently match the imaging data (Supplementary Figure 6). Objective function curves with and without the inclusion of RV free wall points are shown in Supplementary Figure 7. Including the RV in the objective function resulted in larger errors and, in some cases, resulted in curves with no unique minimum. The inclusion of the RV in the mechanics problem, however, plays a vital role in restricting motion of the septum (Hadjicharalambous et al., 2017). In future studies, RV epicardial boundary conditions should be tested which result in accurate RV deformation.
This method also presents an elegant solution for estimating dynamically varying active scaling parameters in both the LV and RV in the forward model problem. It avoids data assimilation methods which often require repeated simulations and thus, the presented method reduces computation time. In some diseased states, contractility can vary regionally over the whole heart, such as the case in patients with a myocardial infarction (Chabiniok et al., 2012). In this case, utilising additional constraints on regional displacements, the current method could be adapted to have a spatially varying activation parameter. The method could be also be augmented by adding a time-varying activation model, such as an Eikonal model (Keener, 1991), to specify the spatially varying sequence of activation.
4.5. Regional Stress and Strain
From full-cycle simulations, active fibre stress and fibre stretch were plotted over 16 AHA regions in Figure 9 illustrating regional differences. In some DCM cases, fibre stretch in some regions was greater than one in early systole, signifying dilation. This may be due to regional systolic dysfunction in these cases. In general, DCM cases showed smaller magnitudes of fibre stretch than healthy and HCM groups, in line with typical systolic dysfunction marked in DCM (Hayashida et al., 1990). Fibre stretch, as opposed to circumferential and longitudinal stretch, could provide more direct measurements of how muscle fibres change with disease. Circumferential stretch measured from the models, plotted in Supplementary Figure 8, were comparable to circumferential strain measured from ultrasound in healthy individuals (Hurlburt et al., 2007; Leitman et al., 2010; Duan et al., 2012). However, model-derived longitudinal stretch was underestimated compared with longitudinal strain from ultrasound. Longitudinal strain is largely dependent on the defined fibre orientation and the model used to describe active contraction. In future, patient-specific fibres as well as constitutive models should be adapted to achieve physiological longitudinal strains. Stress and strain from personalised simulations such as these can provide valuable insights into cardiac function on an individual basis.
In previous cardiac modelling studies, inverse methods have been used to estimate the unloaded geometry of the heart (Krishnamurthy et al., 2013; Wang Y. et al., 2020) which are dependent on the choice of material law, stiffness parameters and boundary conditions. In other studies, various points in the cardiac cycle have been used as the reference geometry including end-systole (Wang et al., 2009), early-diastole (Xi et al., 2013), and diastasis (Wang et al., 2018). However, physiologically, the heart is never in an entirely unloaded state. In early diastole, residual active stress may be present and in all phases of diastolic filling, the cavity pressure is never zero. Although passive parameter estimates have been shown to be minimally affected by changing the reference state from end-systolic to early-diastolic geometries (Hadjicharalambous et al., 2014b), the impact of the choice of reference state is examined further in Hadjicharalambous et al. (2021) and should be assessed in biventricular patient-specific modelling.
The personalised parameters in this study, e.g., a0, af, αlv, and αrv, are all dependent on pressure estimates. If available, catheter measurements from within the LV cavity would enable accurate scaling of these parameters for each individual, and better certainty on model data. LV filling pressures can also be approximated with knowledge of the peak blood flow through the mitral valve as well as the mitral valve peak annular velocity (Nagueh et al., 1997). Given that the mitral valve annular velocity can be obtained in the current pipeline using the landmark predicted valve points, the additional acquisition of 4D flow MR imaging could provide measurements of peak blood flow, enabling appropriate personalisation of all parameters through non-invasive imaging. New methods, such as the use of microbubbles within the LV (Forsberg et al., 2005; Dave et al., 2012), may soon enable more accurate non-invasive cavity pressure measurements. Here, we demonstrate the feasibility of the personalised modelling method using standardised pressure data. If available, pressure data can easily be incorporated into the current pipeline.
This work presents a pipeline using neural networks for generating high quality biventricular models from standard MR cine data. A cohort of 12 individuals were used to demonstrate the pipeline in three different groups: healthy volunteers, DCM patients and HCM patients. Despite the varied morphology and motion of each case, the automated pipeline robustly allowed for determination of a unique passive material parameter describing the tissue anisotropy (γ) as well as two active scaling parameters controlling systolic contraction in the LV and RV (αlv and αrv). The entire pipeline was run using only images from a typical clinical scan, demonstrating its potential to be applied to a large cohort of retrospective data. The use of neural networks along with the model fitting step significantly sped up the process for creating high-quality finite element models. cardiac cycle. This study demonstrates a pipeline that is suitable to model cardiac mechanics and estimate personalised parameters in a diverse cohort of individuals including healthy volunteers, DCM and HCM patients with varying morphologies.
Data Availability Statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.
The studies involving human participants were reviewed and approved by London Bridge National Research Ethics Service. The patients/participants provided their written informed consent to participate in this study.
RM implemented the pipeline, created the tool for fibre field generation, and ran all simulations to estimate patient-specific material properties. EK developed both networks for MR labelling and landmark identification. CM developed the biventricular model fitting tool. TI collected the MR imaging data from HCM patients. AY gave input on the landmark network and revised the manuscript. DN was involved in supervision and collection of DCM and healthy volunteer imaging data. DN, RM, and EK were involved in writing the manuscript. All authors critically reviewed the manuscript.
DN would like to acknowledge funding from Engineering and Physical Sciences Research Council (EP/N011554/1 and EP/R003866/1). This work was also supported by the Wellcome ESPRC Centre for Medical Engineering at King's College London (WT 203148/Z/16/Z) and the British Heart Foundation (TG/17/3/33406). AY would like to acknowledge funding from the National Heart, Lung and Blood Institute (NIH R01HL121754).
Conflict of Interest
The 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.
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2021.716597/full#supplementary-material
Arts, T., Delhaas, T., Bovendeerd, P., Verbeek, X., and Prinzen, F. W. (2005). Adaptation to mechanical load determines shape and properties of heart and circulation: the CircAdapt model. Am. J. Physiol. 288, 1943–1954. doi: 10.1152/ajpheart.00444.2004
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
Asner, L., Hadjicharalambous, M., Chabiniok, R., Peresutti, D., Sammut, E., Wong, J., et al. (2015). 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
Attar, R., Perea nez, M., Gooya, A., Albá, X., Zhang, L., de Vila, M. H., et al. (2019). Quantitative CMR population imaging on 20,000 subjects of the UK Biobank imaging study: LV/RV quantification pipeline and its evaluation. Med. Image Anal. 56, 26–42. doi: 10.1016/j.media.2019.05.006
Augenstein, K. F., Cowan, B. R., LeGrice, I. J., Nielsen, P. M. F., and Young, A. A. (2005). Method and apparatus for soft tissue material parameter estimation using tissue tagged magnetic resonance imaging. J. Biomech. Eng. 127, 148–157. doi: 10.1115/1.1835360
Augenstein, K. F., Cowan, B. R., LeGrice, I. J., and Young, A. A. (2006). “Estimation of cardiac hyperelastic material properties from MRI tissue tagging and diffusion tensor imaging,” in Proceedings of the International Conference on Medical Image Computing and Computer-Assisted Intervention, Vol. 9 (Berlin; Heidelberg), 628–635. doi: 10.1007/11866565_77
Aurigemma, G. P., Zile, M. R., and Gaasch, W. H. (2006). Contractile behavior of the left ventricle in diastolic heart failure: with emphasis on regional systolic function. Circulation 113, 296–304. doi: 10.1161/CIRCULATIONAHA.104.481465
Bai, W., Sinclair, M., Tarroni, G., Oktay, O., Rajchl, M., Vaillant, G., et al. (2018). Automated cardiovascular magnetic resonance image analysis with fully convolutional networks. J. Cardiovasc. Magn. Reson. 20, 1–12. doi: 10.1186/s12968-018-0471-x
Balaban, G., Finsberg, H., Funke, S., Håland, T. F., Hopp, E., Sundnes, J., et al. (2018). In vivo estimation of elastic heterogeneity in an infarcted human heart. Biomech. Model. Mechanobiol. 17, 1317–1329. doi: 10.1007/s10237-018-1028-5
Bayer, J. D., Blake, R. C., Plank, G., and Trayanova, N. A. (2012). A novel rule-based algorithm for assigning myocardial fiber orientation to computational heart models. Ann. Biomed. Eng. 40, 2243–2254. doi: 10.1007/s10439-012-0593-5
Bernard, O., Lalande, A., Zotti, C., Cervenansky, F., Yang, X., Heng, P. A., et al. (2018). Deep learning techniques for automatic MRI cardiac multi-structures segmentation and diagnosis: is the problem solved? IEEE Trans. Med. Imaging 37, 2514–2525. doi: 10.1109/TMI.2018.2837502
Campos, J. O., Sundnes, J., Dos Santos, R. W., and Rocha, B. M. (2020). Uncertainty quantification and sensitivity analysis of left ventricular function during the full cardiac cycle: UQ and SA in cardiac mechanics. Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 378:2173. doi: 10.1098/rsta.2019.0381
Chabiniok, R., Moireau, P., Lesault, P. F., Rahmouni, A., Deux, J. F., and Chapelle, D. (2012). Estimation of tissue contractility from cardiac cine-MRI using a biomechanical heart model. Biomech. Model. Mechanobiol. 11, 609–630. doi: 10.1007/s10237-011-0337-8
Chabiniok, R., Wang, V. Y., Hadjicharalambous, M., Asner, L., Lee, J., Sermesant, M., et al. (2016). Multiphysics and multiscale modelling, data-model fusion and integration of organ physiology in the clinic : ventricular cardiac mechanics. Interface Focus 6, 1–24. doi: 10.1098/rsfs.2015.0083
Dave, J. K., Halldorsdottir, V. G., Eisenbrey, J. R., Raichlen, J. S., Liu, J. B., McDonald, M. E., et al. (2012). Noninvasive LV pressure estimation using subharmonic emissions from microbubbles. JACC Cardiovasc. Imaging 5, 87–92. doi: 10.1016/j.jcmg.2011.08.017
Doste, R., Soto-Iglesias, D., Bernardino, G., Alcaine, A., Sebastian, R., Giffard-Roisin, S., et al. (2019). A rule-based method to model myocardial fiber orientation in cardiac biventricular geometries with outflow tracts. Int. J. Numer. Methods Biomed. Eng. 35, 1–17. doi: 10.1002/cnm.3185
Duan, F., Xie, M., Wang, X., Li, Y., He, L., Jiang, L., and Fu, Q. (2012). Preliminary clinical study of left ventricular myocardial strain in patients with non-ischemic dilated cardiomyopathy by three-dimensional speckle tracking imaging. Cardiovasc. Ultrasound 10:8. doi: 10.1186/1476-7120-10-8
Finsberg, H., Xi, C., Zhao, X., Tan, J. L., Genet, M., Sundnes, J., et al. (2019). Computational quantification of patient-specific changes in ventricular dynamics associated with pulmonary hypertension. Am. J. Physiol. Heart Circul. Physiol. 317, H1363–H1375. doi: 10.1152/ajpheart.00094.2019
Hadjicharalambous, M., Asner, L., Chabiniok, R., Sammut, E., Wong, J., Peressutti, D., et al. (2017). Non-invasive model-based assessment of passive left-ventricular myocardial stiffness in healthy subjects and in patients with non-ischemic dilated cardiomyopathy. Ann. Biomed. Eng. 45, 605–618. doi: 10.1007/s10439-016-1721-4
Hadjicharalambous, M., Chabiniok, R., Asner, L., Sammut, E., Wong, J., Carr-White, G., et al. (2014a). Analysis of passive cardiac constitutive laws for parameter estimation using 3D tagged MRI. Biomech. Model. Mechanobiol. 14, 807–828. doi: 10.1007/s10237-014-0638-9
Hadjicharalambous, M., Lee, J., Smith, N. P., and Nordsletten, D. A. (2014b). A displacement-based finite element formulation for incompressible and nearly-incompressible cardiac mechanics. Comput. Methods Appl. Mech. Eng. 274, 213–236. doi: 10.1016/j.cma.2014.02.009
Hadjicharalambous, M., Stoeck, C., Weisskopf, M., Cesarovic, N., Loannou, E., Vavourakis, V., et al. (2021). Investigating the reference domain influence in personalised models of cardiac mechanics: effect of unloaded geometry on cardiac biomechanics. Biomech. Model. Mechanobiol. 20, 1579–1597. doi: 10.1007/s10237-021-01464-2
Haland, T. F., Hasselberg, N. E., Almaas, V. M., Dejgaard, L. A., Saberniak, J., Leren, I. S., et al. (2017). The systolic paradox in hypertrophic cardiomyopathy. Br. Med. J. 4:e000571. doi: 10.1136/openhrt-2016-000571
Hamdan, A., Guetta, V., Konen, E., Goitein, O., Segev, A., Raanani, E., et al. (2012). Deformation dynamics and mechanical properties of the aortic annulus by 4-dimensional computed tomography: insights into the functional anatomy of the aortic valve complex and implications for transcatheter aortic valve therapy. J. Am. Coll. Cardiol. 59, 119–127. doi: 10.1016/j.jacc.2011.09.045
Hayashida, W., Kumada, T., Nohara, R., Tanio, H., Kambayashi, M., Ishikawa, N., et al. (1990). Left ventricular regional wall stress in dilated cardiomyopathy. Circulation 82, 2075–2083. doi: 10.1161/01.CIR.82.6.2075
He, K., Zhang, X., Ren, S., and Sun, J. (2016). “Deep residual learning for image recognition,” in 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 770–778. doi: 10.1109/CVPR.2016.90
Heijman, E., Aben, J.-P., Penners, C., Niessen, P., Guillaume, R., van Eys, G., et al. (2008). Evaluation of manual and automatic segmentation of the mouse heart from cine MR images. J. Magn. Reson. Imaging 27, 86–93. doi: 10.1002/jmri.21236
Holzapfel, G. A., and Ogden, R. W. (2009). Constitutive modelling of passive myocardium: a structurally based framework for material characterization. Philos. Trans. R. Soc. Ser. A Math. Phys. Eng. Sci. 367, 3445–3475. doi: 10.1098/rsta.2009.0091
Huang, G., Liu, Z., and Weinberger, K. Q. (2017). “Densely connected convolutional networks,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. doi: 10.1109/CVPR.2017.243
Huang, X., Deng, L., Zuo, H., Yang, C., Song, Y., Lesperance, M., et al. (2021). Comparisons of simulation results between passive and active fluid structure interaction models for left ventricle in hypertrophic obstructive cardiomyopathy. BioMed. Eng. Online 20, 1–15. doi: 10.1186/s12938-020-00838-4
Hurlburt, H. M., Aurigemma, G. P., Hill, J. C., Narayanan, A., Gaasch, W. H., Vinch, C. S., et al. (2007). Direct ultrasound measurement of longitudinal, circumferential, and radial strain using 2-dimensional strain imaging in normal adults. Echocardiography 24, 723–731. doi: 10.1111/j.1540-8175.2007.00460.x
Kass, D. A., Chen, C. H., Curry, C., Talbot, M., Berger, R., Fetics, B., et al. (1999). Improved left ventricular mechanics from acute VDD pacing in patients with dilated cardiomyopathy and ventricular conduction delay. Circulation 99, 1567–1573. doi: 10.1161/01.CIR.99.12.1567
Kerckhoffs, R. C., Bovendeerd, P. H., Kotte, J. C., Prinzen, F. W., Smits, K., and Arts, T. (2003). Homogeneity of cardiac contraction despite physiological asynchrony of depolarization: a model study. Ann. Biomed. Eng. 31, 536–547. doi: 10.1114/1.1566447
Kerfoot, E., Anton, E. P., Ruijsink, B., Clough, J., King, A. P., and Schnabel, J. A. (2018). “Automated CNN-based reconstruction of short-axis cardiac MR sequence from real-time image data,” in Image Analysis for Moving Organ, Breast, and Thoracic Images (Cham: Springer), 32–41. doi: 10.1007/978-3-030-00946-5_4
Kerfoot, E., King, C. E., Ismail, T., Nordsletten, D., and Miller, R. (2021). “Estimation of cardiac valve annuli motion with deep learning,” in Statistical Atlases and Computational Models of the Heart, eds E. Puyol-Anton, M. Pop, M. Sermesant, V. Campello, A. Lalande, K. Lekadir, et al. (Cham: Springer), 146–155. doi: 10.1007/978-3-030-68107-4_15
Kerfoot, E., Puyol-Antón, E., Ruijsink, B., Ariga, R., Zacur, E., Lamata, P., et al. (2019). Synthesising Images and Labels Between MR Sequence Types With cycleGAN. Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), Springer, 45–53. doi: 10.1007/978-3-030-33391-1_6
Krishnamurthy, A., Coppola, B., Tangney, J., Kerckhoffs, R. C. P., Omens, J. H., and McCulloch, A. D. (2016). “A microstructurally based multi-scale constitutive model of active myocardial mechanics,” in Structure-Based Mechanics of Tissues and Organs (Boston, MA: Springer), 439–460. doi: 10.1007/978-1-4899-7630-7_22
Krishnamurthy, A., Villongco, C. T., Chuang, J., Frank, L. R., Nigam, V., Belezzuoli, E., et al. (2013). Patient-specific models of cardiac biomechanics. J. Comput. Phys. 244, 4–21. doi: 10.1016/j.jcp.2012.09.015
Leiner, T., Rueckert, D., Suinesiaputra, A., Baeßler, B., Nezafat, R., Išgum, I., et al. (2019). Machine learning in cardiovascular magnetic resonance: basic concepts and applications. J. Cardiovasc. Magn. Reson. 21:61. doi: 10.1186/s12968-019-0575-y
Leitman, M., Lysiansky, M., Lysyansky, P., Friedman, Z., Tyomkin, V., Fuchs, T., et al. (2010). Circumferential and longitudinal strain in 3 myocardial layers in normal subjects and in patients with regional left ventricular dysfunction. Heart Failure 23, 64–70. doi: 10.1016/j.echo.2009.10.004
Levrero-Florencio, F., Margara, F., Zacur, E., Bueno-Orovio, A., Wang, Z. J., Santiago, A., et al. (2020). Sensitivity analysis of a strongly-coupled human-based electromechanical cardiac model: effect of mechanical parameters on physiologically relevant biomarkers. Comput. Methods Appl. Mech. Eng. 361:112762. doi: 10.1016/j.cma.2019.112762
Marchesseau, S., Delingette, H., Sermesant, M., Cabrera-Lozoya, R., Tobon-Gomez, C., Moireau, P., et al. (2013). Personalization of a cardiac electromechanical model using reduced order unscented Kalman filtering from regional volumes. Med. Image Anal. 17, 816–829. doi: 10.1016/j.media.2013.04.012
Marian, A. J., and Braunwald, E. (2017). Hypertrophic cardiomyopathy: genetics, pathogenesis, clinical manifestations, diagnosis, and therapy. Circul. Res. 121, 749–770. doi: 10.1161/CIRCRESAHA.117.311059
Martin-Isla, C., Campello, V. M., Izquierdo, C., Raisi-Estabragh, Z., Baeßler, B., Petersen, S. E., et al. (2020). Image-based cardiac diagnosis with machine learning: a review. Front. Cardiovasc. Med. 7:1. doi: 10.3389/fcvm.2020.00001
Mauger, C., Gilbert, K., Suinesiaputra, A., Pontre, B., Omens, J., McCulloch, A., et al. (2018). “An iterative diffeomorphic algorithm for registration of subdivision surfaces: application to congenital heart disease,” in Proceedings of the Annual International Conference of the IEEE Engineering in Medicine and Biology Society, EMBS, 596–599. doi: 10.1109/EMBC.2018.8512394
Nagueh, S. F., Lakkis, N. M., Middleton, K. J., Spencer, W. H., Zoghbi, W. A., and Zuinones, M. A. (2005). Doppler estimation of left ventricular filling pressures in patients with hypertrophic cardiomyopathy. Circulation 111, 3281–3289. doi: 10.1161/CIRCULATIONAHA.104.508812
Nagueh, S. F., Middleton, K. J., Kopelen, H. A., Zoghbi, W. A., and Qui nones, M. A. (1997). Doppler tissue imaging: a noninvasive technique for evaluation of left ventricular relaxation and estimation of filling pressures. J. Am. Coll. Cardiol. 30, 1527–1533. doi: 10.1016/S0735-1097(97)00344-6
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
Nishimura, R. A., Appleton, C. P., Redfield, M. M., Ilstrup, D. M., Holmes, D. R., and Tajik, A. J. (1996a). Noninvasive Doppler echocardiographic evaluation of left ventricular filling pressures in patients with cardiomyopathies: a simultaneous Doppler echocardiographic and cardiac catherization study. J. Am. Coll. Cardiol. 28, 1226–1233. doi: 10.1016/S0735-1097(96)00315-4
Nishimura, R. A., Hayes, D. L., Ilstrup, D. M., Holmes, D. R., and Tajik, A. J. (1996b). Effect of dual-chamber pacing on systolic and diastolic function in patients with hypertrophic cardiomyopathy: acute Doppler echocardiographic and catheterization hemodynamic study. J. Am. Coll. Cardiol. 27, 421–430. doi: 10.1016/0735-1097(95)00445-9
Nolan, D. R., Gower, A. L., Destrade, M., Ogden, R. W., and McGarry, J. P. (2014). A robust anisotropic hyperelastic formulation for the modelling of soft tissue. J. Mech. Behav. Biomed. Mater. 39, 48–60. doi: 10.1016/j.jmbbm.2014.06.016
Opherk, D., Schwarz, F., Mall, G., Manthey, J., Baller, D., and Kübler, W. (1983). Coronary dilatory capacity in idiopathic dilated cardiomyopathy: analysis of 16 patients. Am. J. Cardiol. 51, 1657–1662. doi: 10.1016/0002-9149(83)90205-9
Peirlinck, M., Costabal, F. S., Yao, J., Guccione, J. M., Tripathy, S., Wang, Y., et al. (2021). Precision medicine in human heart modeling. Biomech. Model. Mechanobiol. 20, 803–831. doi: 10.1007/s10237-021-01421-z
Pfaller, M. R., Hörmann, J. M., Weigl, M., Nagler, A., Chabiniok, R., Bertoglio, C., et al. (2019). The importance of the pericardium for cardiac biomechanics: from physiology to computational modeling. Biomech. Model. Mechanobiol. 18, 503–529. doi: 10.1007/s10237-018-1098-4
Potse, M., Dubé, B., Vinet, A., and Cardinal, R. (2006). “A comparison of monodomain and bidomain propagation models for the human heart,” in Annual International Conference of the IEEE Engineering in Medicine and Biology - Proceedings, Vol. 53, 3895–3898. doi: 10.1109/IEMBS.2006.259484
Rihal, C. S., Nishimura, R. A., Hatle, L. K., Bailey, K. R., and Tajik, A. J. (1991). Systolic and diastolic dysfunction in patients with clinical diagnosis of dilated cardiomyopathy relation to symptoms and prognosis. Circulation 90, 2772–2779. doi: 10.1161/01.CIR.90.6.2772
Romeo, F., Pelliccia, F., Cianfrocca, C., Barilla, F., Reale, A., Gallo, P., et al. (1989). Determinants of end-stage idiopathic dilated cardiomyopathy: a multivariate analysis of 104 patients. Clin. Cardiol. 12, 387–392. doi: 10.1002/clc.4960120708
Ronneberger, O., Fischer, P., and Brox, T. (2015). “U-Net: convolutional networks for biomedical image segmentation,” in MICCAI 2015, Vol. 9351 of LNCS (Cham: Springer), 234–241. doi: 10.1007/978-3-319-24574-4_28
Ruijsink, B., Puyol-Antón, E., Oksuz, I., Sinclair, M., Bai, W., Schnabel, J. A., et al. (2019). Fully automated, quality-controlled cardiac analysis from CMR. JACC Cardiovasc. Imaging 13, 684–695. doi: 10.1016/j.jcmg.2019.05.030
Russell, K., Eriksen, M., Aaberge, L., Wilhelmsen, N., Skulstad, H., Remme, E. W., et al. (2012). A novel clinical method for quantification of regional left ventricular pressure strain loop area: a non-invasive index of myocardial work. Eur. Heart J. 33, 724–733. doi: 10.1093/eurheartj/ehs016
Sack, K. L., Aliotta, E., Ennis, D. B., Choy, J. S., Kassab, G. S., Guccione, J. M., et al. (2018a). Construction and validation of subject-specific biventricular finite-element models of healthy and failing swine hearts from high-resolution DT-MRI. Front. Physiol. 9:539. doi: 10.3389/fphys.2018.00539
Sack, K. L., Dabiri, Y., Franz, T., Solomon, S. D., Burkhoff, D., and Guccione, J. M. (2018b). Investigating the role of interventricular interdependence in development of right heart dysfunction during LVAD support: A patient-specific methods-based approach. Front. Physiol. 9:520. doi: 10.3389/fphys.2018.00520
Schnabel, J. A., Rueckert, D., Quist, M., Blackall, J. M., Castellano-Smith, A. D., Hartkens, T., et al. (2001). “A generic framework for non-rigid registration based on non-uniform multi-level free-form deformations,” in Medical Image Computing and Computer-Assisted Intervention-MICCAI 2001, eds W. J. Niessen and M. A. Viergever (Berlin; Heidelberg: Springer Berlin Heidelberg), 573–581. doi: 10.1007/3-540-45468-3_69
Sermesant, M., Chabiniok, R., Chinchapatnam, P., Mansi, T., Billet, F., Moireau, P., et al. (2012). Patient-specific electromechanical models of the heart for the prediction of pacing acute effects in CRT: a preliminary clinical validation. Med. Image Anal. 16, 201–215. doi: 10.1016/j.media.2011.07.003
Stimm, J., Buoso, S., Berberoǧlu, E., Kozerke, S., Genet, M., and Stoeck, C. T. (2021). A 3D personalized cardiac myocyte aggregate orientation model using MRI data-driven low-rank basis functions. Med. Image Anal. 71:102064. doi: 10.1016/j.media.2021.102064
Stoeck, C. T., Deuster, C. V., Fleischmann, T., Lipiski, M., Cesarovic, N., and Kozerke, S. (2018). Direct comparison of in vivo versus postmortem second-order motion-compensated cardiac diffusion tensor imaging. Magn. Reson. Med. 79, 2265–2276. doi: 10.1002/mrm.26871
Strocchi, M., Gsell, M. A., Augustin, C. M., Razeghi, O., Roney, C. H., Prassl, A. J., et al. (2020). Simulating ventricular systolic motion in a four-chamber heart model with spatially varying robin boundary conditions to model the effect of the pericardium. J. Biomech. 2020:109645. doi: 10.1016/j.jbiomech.2020.109645
Tangney, J. R., Chuang, J. S., Janssen, M. S., Krishnamurthy, A., Liao, P., Hoshijima, M., et al. (2013). Novel role for vinculin in ventricular myocyte mechanics and dysfunction. Biophys. J. 104, 1623–1633. doi: 10.1016/j.bpj.2013.02.021
Tomlinson, K. A., Hunter, P. J., and Pullan, A. J. (2002). A finite element method for an eikonal equation model of myocardial excitation wavefront propagation author. SIAM J. Appl. Math. 63, 324–350. doi: 10.1137/S0036139901389513
Wang, V. Y., Lam, H. I., Ennis, D. B., Cowan, B. R., Young, A. A., and Nash, M. P. (2009). Modelling passive diastolic mechanics with quantitative MRI of cardiac structure and function. Med. Image Anal. 13, 773–784. doi: 10.1016/j.media.2009.07.006
Wang, Y., Guo, H., Wang, S., and Lai, Y. (2020). Effectiveness of a patient-specific 3-dimensional printed model in septal myectomy of hypertrophic cardiomyopathy. Pakistan J. Med. Sci. 36, 1678–1682. doi: 10.12669/pjms.36.7.2620
Wang, Y., Zhang, Y., Wen, Z., Tian, B., Kao, E., Liu, X., et al. (2021). Deep learning based fully automatic segmentation of the left ventricular endocardium and epicardium from cardiac cine MRI. Quant. Imaging Med. Surg. 11, 1600–1612. doi: 10.21037/qims-20-169
Wang, Z. J., Santiago, A., Zhou, X., Wang, L., Margara, F., Levrero-florencio, F., et al. (2020). Human biventricular electromechanical simulations on the progression of electrocardiographic and mechanical abnormalities in post-myocardial infarction. Europace 23, 143–152. doi: 10.1093/europace/euaa405
Wang, Z. J., Wang, V. Y., Bradley, C. P., Nash, M. P., Young, A. A., and Cao, J. J. (2018). Left ventricular diastolic myocardial stiffness and end-diastolic myofibre stress in human heart failure using personalised biomechanical analysis. J. Cardiovasc. Transl. Res. 11, 346–356. doi: 10.1007/s12265-018-9816-y
Xi, J., Lamata, P., Lee, J., Moireau, P., Chapelle, D., and Smith, N. (2011). Myocardial transversely isotropic material parameter estimation from in-silico measurements based on a reduced-order unscented Kalman filter. J. Mech. Behav. Biomed. Mater. 4, 1090–1102. doi: 10.1016/j.jmbbm.2011.03.018
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
Keywords: personalised modelling, biventricular mechanics, parameter identification, automatic segmentation, valve landmark identification
Citation: Miller R, Kerfoot E, Mauger C, Ismail TF, Young AA and Nordsletten DA (2021) An Implementation of Patient-Specific Biventricular Mechanics Simulations With a Deep Learning and Computational Pipeline. Front. Physiol. 12:716597. doi: 10.3389/fphys.2021.716597
Received: 28 May 2021; Accepted: 06 August 2021;
Published: 16 September 2021.
Edited by:Linwei Wang, Rochester Institute of Technology, United States
Reviewed by:Paolo Di Achille, Broad Institute, United States
Aurore Lyon, Maastricht University, Netherlands
Copyright © 2021 Miller, Kerfoot, Mauger, Ismail, Young and Nordsletten. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.