In vivo estimation of elastic heterogeneity in an infarcted human heart

In myocardial infarction, muscle tissue of the heart is damaged as a result of ceased or severely impaired blood flow. Survivors have an increased risk of further complications, possibly leading to heart failure. Material properties play an important role in determining post-infarction outcome. Due to spatial variation in scarring, material properties can be expected to vary throughout the tissue of a heart after an infarction. In this study we propose a data assimilation technique that can efficiently estimate heterogeneous elastic material properties in a personalized model of cardiac mechanics. The proposed data assimilation is tested on a clinical dataset consisting of regional left ventricular strains and in vivo pressures during atrial systole from a human with a myocardial infarction. Good matches to regional strains are obtained, and simulated equi-biaxial tests are carried out to demonstrate regional heterogeneities in stress–strain relationships. A synthetic data test shows a good match of estimated versus ground truth material parameter fields in the presence of no to low levels of noise. This study is the first to apply adjoint-based data assimilation to the important problem of estimating cardiac elastic heterogeneities in 3-D from medical images.


Introduction
Myocardial infarction (MI) is a condition in which muscle tissue in the heart is damaged due to a loss of blood supply. After an infarction, there is an increased risk of further complications, such as rupture, infarct expansion, ventricular remodelling, hypertrophy, and heart failure (Holmes et al. 2005). Post-MI, the elastic properties of the myocardium Dedicated to the memory of Hans-Petter Langtangen 1962-2016 B Gabriel Balaban gabriel.balaban@kcl.ac.uk 1 have been shown to play a large role in determining the outcome (Morita et al. 2011;Fomovsky et al. 2012).
A promising way to study the elastic properties of in vivo myocardium is by mathematical modelling and computer simulation. With simulation it is possible to create an in silico representation of a patient's heart after an infarction. This opens up new possibilities for quantification of elasticity, beyond what is available in medical imaging today. Additionally, an in silico model that is personalized to a patient can potentially simulate the effects of treatments or therapies on the patient, thereby improving the outcome and reducing risks after MI.
Previous studies have estimated elasticity as a global value in models of infarcted hearts (Chabiniok et al. 2009;Gao et al. 2014;Fan et al. 2016). This resulted in simulated pressurevolume relations that matched well to those observed in vivo. Additionally, estimated elasticity values were shown to be significantly higher in patients with infarction compared to healthy controls (Fan et al. 2016). While these results are intriguing, the use of global parameters neglects the fact that infarction is a local phenomenon.
A more detailed approach has been to identify infarcted and healthy regions a priori and then define separate param-eters for these regions in a model (Walker et al. 2005;Mojsejenko et al. 2015;. The resulting regional parameters were shown to be higher in the infarcted region as compared to the healthy remote myocardium. This demonstrates the potential of modelling to quantify differences in tissue stiffness within the same heart. However, the infarctions that caused the stiffness differences were induced in otherwise healthy animals, leading to clearly demarcated regions of myocardial infarction. In the general clinical setting, however, patients may suffer from multiple infarctions, possibly occurring at different times and locations, and/or may be suffering from other cardiac pathologies. Such conditions may lead to substantial heterogeneity in elastic properties, not known a priori. To address the issue of spatial heterogeneity in cardiac elasticity, we here present a novel 3-D data assimilation procedure. This procedure employs an adjoint gradient-based optimization method which can efficiently handle highdimensional parameter sets. In turn, this allows for the spatial resolution of heterogeneous elastic parameters throughout the myocardium. Previous studies on the topic of soft tissue elastography have proposed the adjoint gradient approach with 2-D models and synthetic data (Oberai et al. 2003) for both compressible (Gokhale et al. 2008) and incompressible mechanics (Goenezen et al. 2011). Furthermore, we applied adjoint gradient optimization to the problem of estimating local cardiac contraction Finsberg et al. 2017), but did not consider spatially resolved elastic parameters, which we now address.
We demonstrate the utility of our method by personalizing an in silico model of cardiac elasticity to data collected from a patient in heart failure with a previous myocardial infarction and a heterogeneous distribution of fibrotic tissue. Input data consist of regional strains, which are computed by speckle tracking echocardiography, and a pressure transient obtained from a catheter. Additionally, we quantify the patient's cardiac scar burden from late gadolinium enhanced magnetic resonance images (LG-MRI) to provide a context for the modelling results.

Clinical data
Clinical data were obtained with the permission of the Oslo University Hospital in the context of the Impact study (Hospital 2016). Specifically, we consider the case of a 64-year-old man in systolic heart failure, with left bundle branch block, coronary artery disease, and chronic infarction predominantly in the inferior and inferolateral sections of the left ventricular wall.
Prior to treatment, the patient had echocardiography, LG-MRI, and left ventricular (LV) pressure measurements taken, which are the basis for the clinical data used in this study. Pressure recordings were carried out with an intra-vascular pressure sensor catheter (Millar micro catheter: precision 1 mmHg, accuracy 1.5 mmHg Millar 2017); that is, a pressure catheter that was positioned in the LV via the right femoral artery. Pressure data were obtained automatically and digitized (Powerlab system, AD Instruments) before offline analyses were performed with a low-pass filter of 10Hz.
A 4-D echocardiography examination of the patient's LV was performed using a GE Vingmed E9 machine. Speckle tracking motion analysis was carried out with GE's software package EchoPac. Data from 6 heartbeats were combined in order to obtain a single sequence of images for a single heartbeat. Example short-and long-axis slices taken from the image sequence are shown in Fig. 1. Seven separate measurement points of left ventricular strain during atrial systole were obtained from the echo images. The strains were given as regional averages defined for a standard 17 segment AHA representation (Cerqueira et al. 2002) and measured in the local longitudinal, radial, and circumferential directions, without any off-diagonal shear components.
Strain and pressure data were synchronized using begin of atrial systole (BAS) as the first point of registration. In the pressure data, BAS was located by a deflection in a simultaneously acquired left atrial electrogram. In the strain data, BAS was identified by the onset of longitudinal stretching following diastasis. Pressures corresponding to strains were registered by the use of image acquisition times until just before ventricular systole, which was identified in the strain data by the onset of longitudinal contraction.
Pressure increases in late diastole are generally very small in magnitude, and for our patient strain points 2 and 3 shared the same pressure measurement. In order to give each strain point a unique pressure, an additional cubic polynomial smoothing was carried out. Both smoothed and original pressure data are illustrated in Fig. 2.
Cardiac magnetic resonance imaging was performed with a 3.0 Tesla scanner (Skyra, Siemens, Erlangen, Germany). We quantified the amount of myocardial fibrosis on a per region basis from short-axis late gadolinium enhancement images acquired 10-20 min after intravenous injection of 0.2 mmol/kg of gadoterate meglumine (Guerbet, Villepinte, France). This resulted in an estimated volume ratio of fibrotic to healthy tissue for each myocardial segment (scar burden). In this analysis the apex region was merged into the neighbouring apical regions, giving a 16 segment division. Example LG-MRI images and the scar burden data are displayed in Fig. 3.

Mesh and fibre generation
We created a computational mesh based on a 3-D ultrasound image to capture the details of the patient ventricular geometry in an in silico model. The image was taken at the start of atrial systole, when the pressure was at a minimum. Using GE's EchoPac software, we extracted triangulated data points for the left ventricular endocardial and epicardial surfaces. These surfaces were cut by a plane fitted to the basal points of the surfaces, and adjusted so that the ventricular volume of the computational mesh was within 1 mL of the volume measured in the image. Using the epicardial, endocardial, and basal surfaces as boundaries, we created a volumetric mesh using Gmsh (Geuzaine and Remacle 2009). This mesh contained 741 vertices and 2214 tetrahedra. AHA zones were delineated on this volumetric mesh based on data provided by EchoPac, so that our AHA zones were consistent with those used to calculate image-based strains.
Local myocardial fibre orientations were assigned with a helix angle of 40 degrees on the endocardium rotated clockwise throughout the ventricular wall to −50 degrees on the epicardium using a rule-based method (Bayer et al. 2012). Snapshots of the image-based geometry, along with AHA segments and fibres, are shown in the bottom row of Fig. 1.

Elastic wall motion model
We adopt a quasi-static continuum mechanics framework to simulate the motion of the left ventricle throughout atrial systole. As primary variables, we consider a vector field u giving the displacement map between a reference configuration Ω and a deformed configuration undergoing a pressure load. Furthermore, we define the deformation gradient F = Grad u + I.

Fig. 3
Top row: two example short-axis late enhancement gadolinium MRI images used for regional scar quantification. Fibrotic sections of the myocardium appear in white. Bottom row: regional quantification of myocardial scar burden based on LG-MRI. The inner, middle, and outer rings represent apical, midwall, and basal sections, respectively. The RV insertion points are marked by two horizontal lines extending to the left of the bullseye In our wall motion model the myocardium is considered to be a hyperelastic material with strain energy given by a transversely isotropic simplification of the Holzapfel-Ogden law (Holzapfel and Ogden 2009), The energy density ψ in (1) defines the amount of elastic energy stored per unit volume myocardium, given the values of the right Cauchy-Green tensor C = F T F. The notation (·) + refers to max{·, 0}, and the mechanical invariants I 1 and I 4 f are defined as with e f indicating the local myocardial fibre direction field. The material parameters a, a f , b, b f are scalar-valued quantities which influence the stiffness of the material. We allow these material parameters to vary spatially with a piecewise linear representation, so that each material parameter has a separate value at each vertex of the mesh. For the sake of improved numerical stability (Land et al. 2015), we employ a modified strain energy densityψ in place of ψ with where J = det F is the deformation gradient. The elastic energy (3) is embedded into a standard pressuredisplacement variational formulation of incompressible hyperelasticity [Chapter 8.5 of Holzapfel (2000)]. Displacements are set to 0 in the longitudinal direction at the base of the ventricular geometry by a Dirichlet boundary condition. Movement in the other directions at the base is restricted by a linear spring with constant k = 1.0 kPa as in our previous study ). The total variational equation, including the effects of blood pressure, p blood , and the basal spring, is given by: find the displacement u and the hydrostatic pressure p such that for all admissible variations δu, δ p in the displacement and pressure respectively. In (4), P is the first Piola-Kirchhoff tensor: P = ∂ψ ∂F , ∂Ω endo represents the endocardium and ∂Ω base the ventricular base, and N is the unit outward facing boundary normal. We discretize (4) by a mixed finite element method with Taylor-Hood interpolation (Hood and Taylor 1974); that is, a piecewise quadratic representation of the displacement field and a piecewise linear representation of the pressure.
The software implementation of the finite element vector and matrix assembly code is based on the software package FEniCS (Logg et al. 2011). Nonlinear systems are solved using the PETSc SNES implementation of a Newton line search algorithm (Balay et al. 2015), while the inner linear solves are handled by a distributed memory parallel LU solver (Li and Demmel 2003).

Elastic parameter estimation via constrained minimization and adjoint gradient calculations
We consider a least squares minimization of the mismatch between model derived and measured strains, to personalize the elastic material properties of our computational mechanics model. We compute both model and measured strains in terms of the deformation gradient tensor F and multiply the model strains by F −1 0 , which is the strain at the smallest measured in vivo pressure. This allows for the simulated strains to be calculated from a reference that is at the same pressure as that used for the image-based strains (Nikou et al. 2016). For a given echo image number i, AHA region Ω j , and strain direction k, we compute the model strain as where |Ω j | is the AHA segment volume, and e k the unit vector pointing in the direction k.
The image-based strain measurements are given as regional engineering strains, which we relate to a diagonal component of the deformation gradient by where ε is the engineering strain and F measured the corresponding measured deformation gradient diagonal component. We note that this implies the linear approximation ε k ≈ ∇u · e k . We quantify the mismatch between model and measured strains with the following functional Here N m = 7 is the number of strain measurements available in atrial systole and N r = 16 the number of AHA regions, with the apex segment excluded for compatibility with the LG-MRI data. Finally, the direction set of index k refers to the circumferential (c), radial (r), or longitudinal (l) directions.
We allow each of the four elastic material parameters a, b, a f , b f in (1) to vary in space, and more precisely, to vary as a continuous piecewise linear function defined relative to the computational mesh. This allows us to resolve spatially heterogeneous material parameters, at the cost of greatly increasing their dimensionality. To constrain the minimization problem at hand, we introduce a first-order Tikhonov regularization functional favouring more smooth material parameter sets. This regularization functional is defined as: where |Ω| is the volume of the simulated myocardium.
In total, we consider the optimization problem of minimizing a combined data and smoothness functional over the admissible material parameter fields a, b, a f , b f : with regularization parameter λ. The total functional (9) is minimized by simultaneously optimizing all of the degrees of freedom of the 4 elastic parameters. This optimization is carried out by a sequential quadratic programming (SQP) algorithm (Kraft 1988). Each iteration of the SQP algorithm requires one or more evaluations of the functional (9), and the gradient of the functional with respect to all of the material parameter variables. This gradient is calculated efficiently by the adjoint gradient method [Eq. 13 of Balaban et al. (2016)] symbolically derived by the software package dolfin-adjoint (Farrell et al. 2013). In particular, the computational cost of the adjoint gradient does not significantly depend on the number of optimization parameters, of which there are 2964 in our study. This compares favourably with a one sided finite difference approach to functional gradient calculation, which would require 2964 model realizations, one for each optimization parameter.
We employ a continuity scheme (Gokhale et al. 2008) to reduce the number of nonlinear solves needed to evaluate the functional (9). In this scheme the first time the functional is evaluated the cavity pressure is applied in small increments, and the displacement-pressure solutions are saved at the seven recorded in vivo pressures. For further functional evaluations, the hyperelastic equation is solved directly at the seven pressure levels, with the previously stored displacement-pressure solution as the initialization point. If convergence in the Newton solver is not achieved, then the difference between the previous and current material parameter vector is divided into smaller increments, which are then applied. In our implementation the number of divisions is doubled every time that convergence is not achieved. Using such divisions we obtained convergence in all cases in our study.

Numerical results
The main results of this study are heterogeneous elastic material parameters optimized to match clinical data, presented in Section 3.2. We also present simulated equi-biaxial extensions tests based on regional averages of estimated elastic parameters. Prior to the presentation of the main results, we present results for a synthetic data test for the purpose of verification and inspection of algorithm performance in Section 3.1.
Optimizations were carried out until the norm of the projected gradient was less than 1.0 × 10 −4 , or 500 iteration of the SQP algorithm had been reached. A lower bound of 0.4 was applied to all material parameter fields pointwise during optimization.

Parameter estimation and evaluation using synthetic data
For the purpose of verification of the model and the optimization procedure, we consider initial trials using synthetically generated data over the ventricular mesh. In these trials, the ground truth elastic parameters were defined as: where y max and z max are the maximum absolute coordinate values in the y and z directions of the computational mesh (and where the yz-plane was defined by the basal plane). Using these ground truth parameters, average regional strains were generated by solving (4) for 6 LV blood pressures: p blood ∈ {0.1, 0.2, 0.3, 0.4, 0.5, 0.6} (kPa). Four sets of strains were generated: one noise-free case and three noisy cases. For the noisy cases, realizations of Gaussian noise with standard deviations of 0.1, 0.2, and 0.3 mm were applied to the displacements from which strains were calculated. We quantified the effect of this noise on the value of the synthetic strains in the second column of Table 1. We note that though the average effect of the noise is small, individual strains have relative errors as high as 24, 25, and 12 per cent for the 0.1, 0.2, and 0.3 mm noise levels, respectively. Optimizations were carried out using the synthetic strains as target data in the total functional (9). All material parameters were initialized to a spatially constant value of 1.5. For each level of noise, the values λ ∈ {1, 10, 100, 1000, 10000} of the regularization parameter were tested, and the case with the lowest relative L 2 -error, averaged across the 4 parameters, was selected. The λ values that were selected are listed in Table 1. As expected, the regularization value increases with the noise level. In the target functional (7), F −1 0 was calculated from the model strains at 0.1 kPa.
We remark that, in order to represent non-trivial material parameters, the synthetic material parameter fields were chosen with a nonzero spatial gradient. In turn, this gave a nonzero contribution from the regularization functional, cf. (8). Thus, even in the case of an exact optimization, we did not expect to obtain an optimal functional value of 0 and did not expect to recover the exact material parameter fields in this test case.
The ground truth and estimated material parameter fields for this noise-free case are presented in Fig. 4. Moreover, the differences between the ground truth and estimated parameters are given in terms of the relative L 2 -errors in Table 1, along with the optimal data and smoothness functional values.
We note that for the noise-free and 0.1 mm noise case, the parameters are accurately reproduced, with all relative errors being less that 6%. For the 0.2 and 0.3 mm noise cases, the errors in the first three material parameters a, b, a f are also less than 6%, but the error in the parameter b f is 19% and 17% for 0.2 mm and 0.3 mm noise respectively. The accuracy of the reproductions is also visible in Fig. 4, where we can see that the linear gradients are reproduced for a, b, a f in all cases and for b f in the noise-free and 0.1 mm cases.

Parameter estimation using patient-specific strain data
As a first step towards creating a patient-specific model of the infarcted left ventricle in atrial systole, we identified suitable values for the regularization parameter λ.
We tested a series of trial material parameter optimizations with the patient strains as target data using λ ∈ Synthetic data are corrupted by varying levels of Gaussian noise applied to the ground truth displacement field. The standard deviations of these noise realizations are given in the first column. The relative strain error due to noise is given in the second column. The numbers presented are the mean, and max of the quantities  Table 5, case P2). Optimal data and regularization functional values I data and I smooth were obtained for each of the λ values tested. These are shown in Fig. 5. For the subsequent experiments, we selected λ = 5 as the corresponding optimal functional values lay in a corner of the trade-off curve, and therefore represented a good compromise between smoothness and data fit. This choice of λ is inspired by the so called L-criterion (Hansen and O'Leary 1993).
With the value of the regularization parameter λ fixed, we carried out a series of optimizations using various initializations for the elastic parameters. These initializations  Table 2 along with the spatial mean and standard deviation for each elastic parameter. We note that there is great variability in the optimal parameter sets calculated, and that there is a clear best fitting parameter set. Furthermore, the values of the smoothness functional are similar among all parameter sets and small in comparison to the total functional values. This indicates that the parameter sets differ in their ability to fit the model to the data, but are similar in their smoothness.
The best fitting parameter fields (corresponding to the first row of Table 2) are visualized in the top row of Fig. 6. We note that these fields are fairly smooth, yet show large variation across the ventricle. We also compare strain curves generated by the optimized model to the patient strains in Fig. 7. Optimized material parameter fields are given as a spatial mean followed by the standard deviation in brackets. Standard deviations are calculated as the L 2 (Ω) distance of a parameter field to its mean, normalized by the mesh volume. The first three columns show the value of the optimized total, data, and smoothness functionals. The last column shows the number of least squares function evaluations needed for convergence. A limit of 500 SQP iterations was used, with each SQP iteration consisting of 1 or more function evaluations

Stability of optimized material parameters under mesh refinement
In order to test the effect of mesh refinement on the estimated material parameters, we have carried out a parameter estimation with a slightly finer mesh (1117 vertices, 3373 elements). This estimation was initialized with the same constant values that were used previously in the best fitting optimization to clinical data. The target data were the clinical strains. The resulting optimal material parameter fields are shown in the bottom two rows of Fig. 6. We note that the corresponding original and higher resolution parameters appear to be very similar. We note that the higher resolution parameters came with an increased computational cost, as the time required for an average evaluation of the total functional increased from 30 s to 46 s as compared to the original resolution.

Regional stress-strain relationships
The personalization of the mechanics model to the patient data resulted in four material parameters fields that were resolved in space over the ventricular geometry. We combined these four parameters into a more intuitive visualization of stiffness by considering regional stress-strain relationships. This allows for regional comparisons to be made for a given level of strain, as stiffer materials give higher stresses given the same strain. Regional stress-strain curves were calculated with in silico equi-biaxial extension tests, using analytical values for the stresses based on [Eqs. 17, 18 of Holzapfel and Ogden (2008)]. A test was conducted per AHA region using the average of the material parameter fields over the corresponding region. The resulting stress-strain relations along the fibre and cross-fibre directions are presented in Fig. 8.

Discussion
By applying an adjoint gradient-based data assimilation method, we were able to estimate spatially heterogeneous material properties in an infarcted left ventricle with a good match of simulated to measured strains. This has important implications for the use of computational mechanics models in planning and optimizing therapies in silico. Conditions such as myocardial infarction are local and lead to elastic heterogeneities which should be accounted for in a personalized model. This study presents a general and flexible method to account for these elastic heterogeneities. Fig. 6 View of optimal material parameters at two different mesh resolutions estimated from patient strain data. The first and third rows show the inferior view and the second and fourth rows the anterior view Our experiments with synthetic data indicate that fairly accurate reproductions of spatially varying parameters are possible in the absence of noise, as the relative L 2 errors were less than 5% for all parameters in this case. In the presence of Gaussian noise the relative errors in the a-type parameters increased slightly, but were still below 6%. The effect of the noise on the reproduced b-type parameters was more pronounced, and in particular, errors in the b f parameter for the 0.02 and 0.03 mm noise cases were large enough that the spatial gradient present in the ground truth b f parameter could no longer be reproduced. These results suggest that spatial heterogeneities can be more robustly estimated with a-type parameters rather than b-type exponential parameters. Indeed, several recent data assimilation studies have limited parameters estimations to the a-type parameters only Asner et al. 2017).
We note that the optimal I data values were two orders of magnitude higher in the clinical case than in the synthetic case. This could be due to higher noise in the clinical data and or modelling error in the representation of in vivo cardiac motion (1). Similarly, the optimal I smooth values were several orders of magnitude higher in the clinical case than in the synthetic case. This is due to relatively higher gradients in the optimal material parameters fitted to the clinical data.
In both simulated and measured patient data, we noticed that the heavily infarcted region encompassing inferior to inferolateral segments at the base and the mid-posterior segment differed in several ways from healthy segments. In these infarcted segments strains were smaller, and the simulated equi-biaxial stress-strain relationships showed greater fibre stresses. Additionally, the optimal a and a f material parameters are larger in the infarcted anterolateral segment, and there is a band of high a parameter value running through the infarcted inferior segments. These observations indicate increased myocardial stiffness. This is consistent with the increased stiffness observed in healing infarcts during an ex vivo tissue experiment (Gupta et al. 1994) and in previous computational modelling with in vivo data .
We also observed that the mid-anterolateral segment was identified as free from scar in the late enhancement MRI analysis, yet showed signs of stiffness similar to the heavily infarcted segments described above, that is both low strain and high simulated stress. Such apparent stiffness in a healthy segment is consistent with an infarction impairing the mechanics of neighbouring healthy tissue (Holmes et al. 2005) or could be an effect of myocardial border zone tissue.
Ideally, model material parameters should be uniquely identifiable from in vivo data in order to produce potentially useful biomarkers for clinical practice. Recently, it has been shown that the linear parameters a and a f of a reduced Holzapfel-Ogden law (1), are structurally identifiable . Structural identifiability means that there exist sets of model loaded states such that only one set of parameters produces them, making it theoretically possible to uniquely identify the parameters. Our in vivo data are corrupted by noise, which makes the question of the unique identifiability of parameters more complex. Additionally, we have optimized the exponential b and b f parameters in our in vivo experiment, for which possible structural identifiability is still unknown. Last but far from least, we have spatially resolved all of the parameters, thereby greatly increasing their dimensionality. Under such circumstances the theoretical identifiability of material parameters is an open question.
To improve the identifiability of material parameters in our estimations, we have added regularization to the optimized functional. Indeed, Fig. 5 confirms the existence of several material parameter sets that fit the model to the data very similarly, but differ in their smoothness. By choosing a corner point in the space of optimized data and smoothness functionals our aim was to pick the smoothest set of elastic parameters that still fit the data well. However, even with the regularization, our parameter estimation still showed a dependency on the choice of initial parameters, and a variety of results were obtained (Table 2). Nevertheless only one parameter set fit the best, allowing us to choose it from among the others.  . 8 Regional fibre and cross-fibre stress-strain curves generated from simulated equi-biaxial extension tests. In each AHA region the spatial average of the optimal material parameters is used in the simulated extension experiment. The colour of each line indicates the corresponding regional scar burden value

Limitations
The identifiability of material parameters was limited in our study, and all optimizations depended upon their initial guess. This dependence is demonstrated in Table 2 by the variety of minima. In the future it would be of interest to further examine constraints to spatially resolved material parameters, ideally yielding an optimization procedure that yields the same parameters regardless of the initialization. One such possible constraint is the left ventricular chamber volume, which has been previously matched together with strain data Mojsejenko et al. 2015;Sun et al. 2009). Further possible constraints are aggregated geometry measures such as LV twist, and long-and short-axis motion. These have been shown to improve identifiability of elastic parameters in experiments with mouse ventricles (Nordbø et al. 2014).
Further limitations were related to the rule-based fibres, mechanics modelling, computational efficiency, and strain and pressure synchronization

Rule-based fibres
The fibre orientations in our model were generic and not patient specific. As a result, healthy fibre angles were used in infarcted areas. Previous studies have shown that fibre orientations of infarcted areas can be significantly different from healthy tissue Fomovsky et al. 2012). If this effect were incorporated in our parameter estimation, we would expect a change in the optimal material parameters in the infarcted areas, especially in the a f and b f parameters, which control the amount of anisotropy in the model along the fibre direction. In the future, further improvements to diffusion-tensor MRI technology may allow for in vivo identification of local myocardial fibre directions, which would allow for the fibre directions to be directly incorporated into the optimized model without needing to be estimated.

Mechanics modelling
The image-based reference geometry contained a pressure load that was not accounted for in the current study as 0 pressure was assumed for the reference geometry. Using recently developed techniques, it is possible to calculate a pressure-free reference geometry simultaneously with material parameter estimation (Nikou et al. 2016). Applying this technique in our study was unfortunately not possible as the unloaded mesh self-intersected partway into the calculation when we attempted it.
Active tension was assumed to vanish in our model. Typically this tension has decayed to 0 in the diastasis phase of a healthy heart, but may extend into atrial systole under pathological conditions. If active tension were present in the diastasis phase of our patient, then it could add additional stiffness to the myocardium. At the same time, the release of active tension could contribute to strain in atrial systole. Missing these effects would lead to potential overestimation of passive tissue stiffness in the first case and an underestimation in the second.
The computational model lacked several relevant physical effects, notably inertia, viscoelasticity, residual stresses in the unloading geometry, mechanical coupling of the LV to the right ventricle and atria, the effect of sheet microstructure, and tissue compressibility due to blood entering and exiting the ventricle via coronary vessels. The spring constant at the base was a rough approximation and could be replaced by displacement data at the basal boundary if it were available. The apex of the computational model was free, while longitudinal motion at the base was fixed. The in vivo situation is the opposite, the base moves longitudinally, and the apex is mostly stationary.

Computational efficiency
The spatial discretizations of the material parameters were not optimized. Instead, the computational mesh used to solve the variational equation of motion (4) was also used for the representation of the spatial parameters due to ease of implementation. It is possible that a coarser representation of the material parameters could have also produced good modeldata fits. Using fewer parameters could potentially improve the identifiability of parameters and reduce the number of SQP iterations needed to find a minimum.
For the sake of computational efficiency, the resolution of the mesh was not increased to the point of obtaining a numerically convergent solution. Errors in the discretization of the hyperelastic variational equation (4) may have affected the optimized elastic parameter values. However, the results of our test optimization with a finer mesh indicate that any errors due to insufficient mesh resolution did not substantially affect the overall pattern of the optimized parameters.

Strain and pressure synchronization
LV pressure and strain measurements were not taken simultaneously and had to be synchronized in our study. Though both strain and pressure measurements were taken when the patient was relaxed and prone, there could have been slight differences in heart rate which would confound the strainpressure synchronization.

Conclusion
Adjoint-based data assimilation has been used to personalize a mechanics model to reflect the heterogeneity in material properties throughout an infarcted human left ventricle. Further trials with more datasets and more methodological development are warranted in order to evaluate the applicability of the technique.