Non-invasive Model-Based Assessment of Passive Left-Ventricular Myocardial Stiffness in Healthy Subjects and in Patients with Non-ischemic Dilated Cardiomyopathy

Patient-specific modelling has emerged as a tool for studying heart function, demonstrating the potential to provide non-invasive estimates of tissue passive stiffness. However, reliable use of model-derived stiffness requires sufficient model accuracy and unique estimation of model parameters. In this paper we present personalised models of cardiac mechanics, focusing on improving model accuracy, while ensuring unique parametrisation. The influence of principal model uncertainties on accuracy and parameter identifiability was systematically assessed in a group of patients with dilated cardiomyopathy (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n=3$$\end{document}n=3) and healthy volunteers (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n=5$$\end{document}n=5). For all cases, we examined three circumferentially symmetric fibre distributions and two epicardial boundary conditions. Our results demonstrated the ability of data-derived boundary conditions to improve model accuracy and highlighted the influence of the assumed fibre distribution on both model fidelity and stiffness estimates. The model personalisation pipeline—based strictly on non-invasive data—produced unique parameter estimates and satisfactory model errors for all cases, supporting the selected model assumptions. The thorough analysis performed enabled the comparison of passive parameters between volunteers and dilated cardiomyopathy patients, illustrating elevated stiffness in diseased hearts. Electronic supplementary material The online version of this article (doi:10.1007/s10439-016-1721-4) contains supplementary material, which is available to authorized users.

volume is expected to vary by no more than 10%, in the preliminary tests it was reduced by as much as 30% at systole compared to end diastole, similar to the reduction observed in Shi et al. 1 . This significant variation in wall volume was considered non-physiological and was likely an artefact of the tracking algorithm, which did not ensure displacements were volume conserving. Additionally, the wall thickening over systole was significantly below physiological values for all cases (Table S1).
Accordingly, the motion tracked deformation field was modified to ensure that wall volumes were kept constant throughout the cycle and equal to the wall volume at end diastole.
In order to ensure a close match to the original data, at each time-step the processed displacement field (u d ) was obtained through a volume preserving H 1 projection of the original tracking data (u o ). In this process, the endocardial boundary was kept weakly the same through a Lagrange multiplier λ , and a relaxed boundary condition was introduced on the base boundary 3 . The volume-preserving displacement field u d was obtained by the solution of the following saddle-point problem: where the functional Π at any point in time is given by: The first integral term deals with the data projection, and α = 100 was selected to provide equal weighting between the L 2 norm and H 1 semi-norm terms, over simulations in millimetres. The second and third terms are enforcing the original data motion on endocardial and basal boundaries, respectively. The final term is enforcing the incompressible deformation using Lagrange multiplier λ, where J d is the determinant of the deformation gradient related to the projected displacement u d . Here, Ω ED is the myocardial wall domain at end diastole, and Γ ED and Γ b ED denote the endocardial and basal boundary domains at end diastole, respectively. The employed boundary conditions were selected to ensure consistency between processed and unprocessed quantities of interest. In particular, the endocardial boundary displacements were enforced from the original extracted motion, ensuring that metrics such as stroke volume and ejection fraction were maintained.
The resulting processed meshes were compared against images, demonstrating good agreement over the cardiac cycle as illustrated in a representative example in Figure S1.
Interestingly, the post-processing step enabled physiological wall thickening during systole (Table S1). Additionally, parameter identifiability and reasonable model errors were observed using both processed and unprocessed data. While modest differences were observed in the obtained parameters (relative change up to 20% for a and up to 40% for a f ) the variation between volunteers and patients is consistent (Figures 7b and S2), suggesting that the volume preserving processing step did not fundamentally alter model-derived stiffness.

Figure S1
Unprocessed (blue mesh) and processed (red) extracted motion at end-systole (left) and end diastole (right), compared against cine images in long-axis and short-axis views, for V5.  Figure S2 Comparison of (a) isotropic parameter a and (b) fibre parameter a f between the volunteer and patient groups, when the unprocessed extracted motion was used in model personalisation.

S2 Evaluation of extracted / processed motion accuracy
The accuracy of both the extracted and processed motion was evaluated against manually tracked points, following the methodology described in Shi et al. 1 . In particular, for each case 16 landmark points were selected, lying mid wall in each AHA region 2 . The position of the landmark points was manually tracked through the cycle, p AHA manual (t), along with their displacement, u AHA manual (t). The relative tracking error ρ AHA (t) for each AHA landmark point through the cardiac cycle was assessed as proposed in Shi et al. 1 : Here p AHA tracked (t) refers to the position of the landmark point, based on the IRTK algorithm or the processed extracted motion. Figure S3a presents the average tracking error over all cases, for both the IRTK extracted motion and the processed motion. These results show that the relative tracking errors for the IRTK extracted motion are consistent with Shi et al. 1 . Additionally, for a volunteer case manual tracking was performed by two individuals, showing that the error between the two operators (∼ 0.08±0.02) was close to the error between manual tracking and IRTK tracking.
Across the tests considered, it was observed that the error metric in Eq. S3 tends to be biased by relative errors in regions with reduced motion. This is due to the fact that the metric expects precision proportional to the motion observed in each AHA segment. Regions with little motion introduce much larger relative errors ( Figure S3b), leading to an increase in the mean errors. However, for the purposes of this work we are more concerned with the overall motion than relative motion when comparing to the personalised models. Therefore, we have additionally estimated tracking errors using a maximal displacement metric: In this case, the error is divided by the maximum displacement over all AHA regions u max manual , 0.07 ± 0.02 0.14 ± 0.03 0.15 ± 0.03 0.11 ± 0.03 0.13 ± 0.04 0.14 ± 0.04 0.09 ± 0.02 0.16 ± 0.06 ρ± std 0.03 ± 0.009 0.07 ± 0.02 0.07 ± 0.02 0.05± 0.01 0.05 ± 0.02 0.07 ± 0.02 0.04 ± 0.006 0.06 ± 0.02 Processed data ρ± std 0.12 ± 0.04 0.20 ± 0.08 0.22 ± 0.09 0.14 ± 0.04 0.17 ± 0.05 0.16 ± 0.06 0.14 ± 0.05 0.18 ± 0.06 ρ± std 0.05 ± 0.02 0.08 ± 0.02 0.09 ± 0.03 0.06± 0.02 0.06 ± 0.02 0.08 ± 0.03 0.05 ± 0.01 0.06 ± 0.02 (d) Figure S3 (a) Average tracking error over all cases, assessed using the error metric in Eq. S3. The error compares IRTK extracted motion (black) and processed motion (red) against manually tracked motion. The shaded regions present the standard deviation of the error between cases. (b) Error of every landmark at every timestep, against the maximum displacement of the landmark through the cycle (manual tracking). (c) Average tracking error over all cases, assessed using the error metric in Eq. S4. (d) Tracking error for each case, averaged over the cardiac cycle. Both processed and unprocessed data are considered, and the error is assessed using ρ andρ.
avoiding the requirement that precision is higher in regions with little to no motion. Tracking errors usingρ for both the unprocessed and processed extracted motion are presented in Figure S3c. Usingρ also reduced the discrepancy between operators to 0.04 ± 0.007. Additionally, average tracking errors over the cardiac cycle for all cases are presented in Figure   S3d, using both error metrics.
While the processed extracted motion presented higher tracking errors than the unprocessed IRTK extracted motion, the errors remain comparable to the tracking errors reported in Shi et al. 1 . Moreover, examining errors using the revised error metricρ ( Figures S3c,   S3d), the tracking errors of the unprocessed and unprocessed data are more consistent.
While the processed motion results in slight increases in error, it also enables physiological radial thickening and wall volumes through the cardiac cycle (see Table S1).

S3 Computation of clinical metrics
Motion extracted from TMRI and processed as described in section S1 was utilised for the computation of typical clinical metrics (Table 1). End-diastolic and end-systolic volumes were computed from the respective deformed meshes of the processed data. Cavity volume was computed following the approach presented in Asner et al.  This section examines parameter identifiability and model fidelity for the personalised models created for four volunteers (V1-V4) and three patients (P1-P3). Patient-specific models were developed following the pipeline described in sections 2.1 and 2.2. The effect of these model assumptions on the obtained parameter ratio is presented in Table S2. Additionally, the influence of the assumed fibre distribution and epicardial boundary condition on model accuracy and parameter identifiability is summarised in Figures S4 and S5.
Based on Figures S4 and S5 as well as the analysis in section 3.1 the combination of the RV BC and a fibre distribution with θ = 50 • consistently resulted in lower errors. Additional assessment of model fidelity for these simulations was provided through the nodal distance between model and data (Table S3). The average nodal model error was ∼ 2mm and the relative errors were between 17% and 32%. These errors suggest satisfactory model accuracy.
Case e m ± std e b ± std e rv ± std e e ± std e rem ± std u d ± std V1  Table S3 Average and standard deviation of the magnitude of nodal displacement error between model and data at end-diastole in mm. Model error in this case is assessed over the entire domain (e m ), the basal surface (e b ), the RV attachment region (e rv ), the endocardial surface (e e ), and the remainder of the domain (e rem ). Also presented is the average and standard deviation of the data displacement (u d ) for each case, in mm.
Nodal model error was also computed over different regions of the LV mesh domain (Table S3). The relaxed basal boundary condition employed enabled excellent agreement (< 0.3mm) between model and data at the base. Low errors were also observed on the RV attachment region that was constrained to deform following the data. While low, the errors on the RV attachment region were non-negligible due to the moderately high value used for the relaxation parameter in Eq. 6. The value was selected to enable a good agreement with the data without causing pressure peaks. The displacement errors were higher on the endocardial surface (< 3mm) where only the cavity volume was prescribed.
Reasonable errors (< 3mm) were observed in the mesh region which was not driven by the data, indicating the ability of the model to reproduce cardiac behaviour.
Additionally, the effect of the data frame assumed as reference domain was assessed on all volunteer and patient cases. Following the description of section 3.1, in Figure S6 parameter identifiability and parameter estimates were examined using three different data-frames (endsystolic, second and fourth after end-systolic frames), for V2-V5 and P2-P3.   Figure S6 Objective function J over the parameter ratio γ, with three different data frames assumed as the reference geometry.

S5 Regional strain analysis
This section provides an additional means of assessing model accuracy through the evaluation of regional strain distributions in both data and models. Analysis of regional cardiac behaviour derived from strains might elucidate differences between DCM and volunteer cases which could not be derived from global metrics. To this end, each personalised geometry was divided in 16 regions according to the American Heart Association (AHA) recommendations 2 . The apex was not included as a 17th region due to difficulties in accurately identifying the apical location.
Based on the Lagrange-Green definition, E = 1 2 (F T F − I), strain E is dependent on spatial derivatives of data-derived displacement and is thus very sensitive to noise. In order to circumvent this data sensitivity, an alternative way of computing strain was employed.
Considering the fundamental definition of the deformation gradient (dx = F dX), F maps infinitesimal vectors dX in the undeformed domain to their corresponding vectors dx in the deformed domain. Instead of computing F at every node in the geometry, an approximate valueF per AHA region was estimated. Specifically, for each mesh node of the AHA region, the "infinitesimal" vectors dx and dX were computed with respect to the centroid of the region in the deformed and undeformed domain, respectively.F for every AHA region was then computed as the tensor that satisfied dx =F dX for all nodes k in the region, in the least-squares sense: Strain was subsequently computed using the approximated deformation gradient per region The use of definition S6 for strain avoids bias introduced by noise-sensitive derivative estimates. Combined with the use of mean quantities per AHA region, this strain definition should reduce noise sensitivity, providing physiological strains.
Using this definition, strain was computed for all AHA regions and reported values correspond to strains at end diastole with respect to end systole. Strains in the radial (E rr =Ēr · r), longitudinal (E ll =Ēl · l) and circumferential (E θθ =Ēθ · θ) directions were also computed for each region, with r, l and θ the unit vectors in the radial, long-axis and circumferential directions, respectively. The spectral norm of the strain tensor was computed as well, as the maximum strain eigenvalue (||Ē|| 2 = λ max (Ē)), providing a global metric.
Initially, average values of strain metrics over all regions were considered, to provide  Table S4 Average data strain at end diastole with respect to end systole for all volunteer and patient cases. Strain in the radial (E rr ), longitudinal (E ll ) and circumferential (E θθ ) directions along with strain norm (||Ē|| 2 ) are presented.
overall deformation characteristics (Table S4) and enable direct comparisons between the volunteer and patient groups. As indicated by the decrease in ||Ē|| 2 in the DCM group (p = 0.02), global strain was decreased in DCM hearts. In fact, strain in the radial and circumferential directions was notably reduced in the patient group (p < 0.1 for both E rr and E θθ ). Decrease was also observed in longitudinal strain, even though not statistically significant.
Strain per AHA region was also computed for the personalised models, to enable regional comparisons between data and model. Data and model strains per AHA region are presented for a representative example (V3) in Figure S7. Notable similarities can be observed between spatial distributions of data and model strains. Despite differences in magnitude, the regions of maximum and minimum strains are -in most cases -consistent between model and data.
The similarity in strain distributions suggests that the personalised models reasonably reflect the true data-derived strains.