Parameter estimation in a Holzapfel–Ogden law for healthy myocardium

A central problem in biomechanical studies of personalized human left ventricular (LV) modelling is to estimate material properties from in vivo clinical measurements. In this work we evaluate the passive myocardial mechanical properties inversely from the in vivo LV chamber pressure–volume and strain data. The LV myocardium is described using a structure-based orthotropic Holzapfel–Ogden constitutive law with eight parameters. In the first part of the paper we demonstrate how to use a multi-step non-linear least-squares optimization procedure to inversely estimate the parameters from the pressure–volume and strain data obtained from a synthetic LV model in diastole. In the second part, we show that to apply this procedure to clinical situations with limited in vivo data, additional constraints are required in the optimization procedure. Our study, based on three different healthy volunteers, demonstrates that the parameters of the Holzapfel–Ogden law could be extracted from pressure–volume and strain data with a suitable multi-step optimization procedure. Although the uniqueness of the solution cannot be addressed using our approaches, the material response is shown to be robustly determined.

is not suitable for modelling the passive myocardial response in simple shear tests. Therefore, to account for the layered micro-structure, orthotropic constitutive laws have been proposed, such as the Fung-type law [33], polezero law [34] and strain-invariant-based law [1]. Generally, there are more parameters in orthotropic constitutive laws compared to transversely isotropic laws, which makes the inverse estimation harder to do. Using a synthetic LV model, Remme et al. [20] inversely estimated parameters for the pole-zero law. However, because of the high intercorrelations among the parameters, only 3 out of 18 parameters were estimated.
Although a number of studies have demonstrated the feasibility of inversely estimating constitutive parameters from in vivo clinical measurements using simpler constitutive relations [10], fewer attempts have been made to estimate parameters from in vivo data using orthotropic laws [20]. The orthotropic constitutive law proposed in [1] (known as the Holzapfel-Ogden law) has strong ellipticity with respect to material stability and fewer parameters than the pole-zero law. Most importantly, the advantage of the Holzapfel-Ogden law is that it can account for a layered myofibre architecture. The Holzapfel-Ogden law is also relatively easy to implement using the FEM. However, to the authors' best knowledge, the feasibility of identifying the parameters of the Holzapfel-Ogden law from non-invasive clinical measurements has not yet been investigated.
In this paper, we carry out such a study for the first time using a previously published LV model with known parameters to provide a set of synthetic strain data and pressure-volume relationships [35]. Once verified, the optimization method is then extended and applied to in vivo models with limited clinical measurements.

Constitutive law for passive myocardium
The Holzapfel-Ogden constitutive law [1] assumes that the strain energy function ψ for the myocardium is where a, b are the parameters for the matrix response; a f , b f are the parameters for the myocardial fibres; a s and b s account for the fibre sheet contribution; and a fs , b fs represent the shear effects in the sheet-plane. I 1 = tr(C), C = F T F is the right Cauchy-Green tensor, and F is the deformation gradient. I 4f = f 0 · (Cf 0 ) and I 4s = s 0 · (Cs 0 ) are the stretch-related invariants along the myocyte and sheet directions, f 0 and s 0 respectively, in the reference configuration. I 8fs = f 0 · (Cs 0 ) reflects the coupling between the fibre and sheet stretches. For a more detailed description of the Holzapfel-Ogden law, please refer to [1] and, for its applications in LV modelling, to [35][36][37][38][39].

Synthetic model
To check the feasibility of the inverse estimation of parameters in the Holzapfel-Ogden law from strain and pressurevolume data, a previously published passive LV model [35] with known material parameters (i.e. the so-called true parameters) is simulated using ABAQUS FEA. The myofibre structure and boundary conditions of the LV model are kept the same as in [35]. This produces a set of synthetic data (i.e. the so-called experimental data) when the endocardial pressure is increased from 0 to 8 mmHg. From these data we extract the LV cavity volume, as well as the first, second and third principal strains from randomly distributed observation points inside the LV wall (excluding the basal plane where the boundary conditions are applied) throughout the LV diastolic loading phase.

Optimization procedure
We define three objective functions for matching the experimental data: where T N exp = 25 is the number of timesteps from the beginning to the end of pressure loading; N point = 20 is the number of observation points randomly distributed within the LV wall; V * i and ε k* j,i are the LV chamber volume and principal strains at timestep i from the forward modelling; V i and ε k j,i are the corresponding data from the inverse procedure; and k indicates the first, second and third principal strains.
A sensitivity analysis of parameters is performed by varying each of the eight parameters in a range of ±10% around their so-called true values, while the seven remaining parameters are chosen to be the true values. A study of the sensitivity of the objective functions to the variation of a parameter is then performed and summarized in Table 1. Notice that because the total number of data points for the volume (25) is much less than the total number of the strains (1,500), the change Δf ε and, hence, Δf obj is much greater than Δf vol . In addition, f obj , f vol and f ε are most sensitive to changes in a f , followed by a, b, a fs , but are less sensitive to changes in a s and b s . This suggests that a s and b s may not be accurately estimated. However, a greater change in f vol due to the variation in b f is observed; thus, using f vol in addition to f ε allows one to extract extra information.
Based on the sensitivity analysis, we propose a multi-step optimization procedure, as shown in Fig. 1. In each step, a forward problem is solved to yield the updated data set V i and ε k j,i . The MATLAB function lsqnonlin, together with a trust-region-reflective algorithm, is used to update the material parameters by minimizing the least-squares difference between the experimental data and the model responses from the forward simulations.

Study 2: Application to in vivo modelling
Having developed the multi-step optimization procedure, we proceed to estimate the parameters of the Holzapfel-Ogden law using in vivo data. An in vivo LV geometry is first reconstructed from a CMR study of a healthy male volunteer, age 28. The CMR study was performed on a Siemens MAGNETOM Avanto (Erlangen, Germany) 1.5 Tesla scanner with a 12-element phased array cardiac surface coil. The imaging protocol included cine sequence with steady-state free precession. Custom MATLAB software was written to segment the endocardial and epicardial boundaries at early diastole when the LV pressure is lowest (Fig. 2a). Following the manual segmentation, prolate spheroidal coordinates and cubic B-splines are used to fit endocardial and epicardial surfaces, smoothness regularization is imposed to minimize the geometric distortion, and constraints are imposed on the fitting parameters to Step 1, all the parameters are estimated using both the volume and strain data, followed by (Step 2) an estimation of a f and b f using the volume data only. Similarly in Step 3, a s , b s , a fs , b fs are estimated using the strain data only. In Step 4, the remaining parameters a, b, a f , b f from Step 3 are optimized using both the volume and strain data start Step 1: optimizing a, b, a f , b f , a s , b s , a fs , b fs with f obj Step 2: optimizing a f , b f with f ol and fixed a, b, a s , b s , a fs , b fs Step 3: optimizing a s , b s , a fs b fs , with f and fixed a, b, a f , b f Step 4: optimizing a, b, a f , b f with f obj and fixed a s , b s , a fs , b fs End maintain C 2 continuity at the apex. Figure 2b shows the unfitted reference LV geometry. Hexahedral elements are generated using a linear interpolation from endocardial to epicardial surfaces, as shown in Fig. 2c. The interested reader is referred to [40] for more details on the LV geometry reconstruction.
An in-house B-spline deformable registration method [41] is used to estimate the regional circumferential strain from early to ED from four positions of short-axis cine images from basal to middle ventricles and six regions for each short-axis position. The LV cavity volumes at ED are also calculated from the cine images. In summary, the data from the in vivo measurements consist of 24 regional circumferential strains and LV cavity volume at ED. Because the ventricular pressure recording is not available, a population-based ED pressure (8 mmHg) is assumed.

Sensitivity analysis
Generally, the in vivo measurements are more noisy and contain significantly fewer data compared to what can be produced by the synthetic model (25 versus 1,525), so a more comprehensive sensitivity analysis is performed.
where N = n + 1, ε i is the strain data from the ith position (i = 1, . . . , n), and V is the LV volume. To measure the model response sensitivity to parameter variations (10% as in Sect. 2.2.2), we introduce the matrix where the relative changes in the strains and volume due to the parameter variations Δκ m are denoted by Δd i = ε i − ε i0 , i ∈ N , 0 indicates the data from the baseline model. If the variations in any of the parameters produce similar responses, then these parameters may be correlated. To quantify the correlations, we further calculate the sensitivity coefficient matrix SCM from the normalized columns of S:  Table 3 Correlation coefficient SCM with normalized LV volume If SCM i, j is close to ±1, then κ i and κ j are closely correlated. The existence of correlated parameters poses nonuniqueness in the inverse problem. A possible remedy is to estimate one parameter only while leaving the correlated ones constant. In addition, the norm reflects the sensitivity of the objective function to the variation of each parameter. Low sensitivity to a parameter makes it difficult to estimate. Table 2 lists the average results of two sensitivity analyses: one is from a 10% increase in each parameter, the other is from a 10% decrease in each parameter; the baseline parameters are from [36]. These results show that a has the highest sensitivity, followed by a f , a fs , b f , b, b fs . a s has a much lower sensitivity, and b s has virtually zero sensitivity. High correlations can also be found among a, b, a f , b f , a s , a fs , b fs . These observations suggest that it is unlikely that the eight parameters can be estimated unambiguously.
Since the absolute LV volume is much greater than the strains, it can overdominate the optimization procedure. Hence, we also perform the same sensitivity study using the normalized LV volume, (V − V 0 )/V 0 , and the results are shown in Table 3. Now the correlation matrix SCM is different from that in Table 2. The four parameters a, b, a fs , b fs are highly correlated, while a f and b f are only correlated to each other, not to a, b, a fs , b fs . Again, a s and b s have very low sensitivity to parameter variations.
To deal with the difficulties caused by the correlation of the parameters, we follow the approach by [2] and divide the parameters into two groups, where a group 0 and b group 0 are the experimentally estimated parameters from the previous study [35] (referred to as the original parameters), and C a , C b are the scaling factors.
In the first optimization step, the optimal C a and C b are found by minimizing in which, ε * i is the measured regional circumferential strain and V 0 is the measured LV cavity volume from cine images. The minimization problem is solved by sweeping the parameter space (C a and C b ) as in [2]. The final parameter set which minimizes f O1 is chosen. If multiple sets of C a and C b exist, the set with the minimized Because f O1 is volume dominated, further steps are needed to match the strains. According to Table 3, a f is only highly related to b f , but not to a, b, a fs and b fs ; therefore, we could set a, b, a fs and b fs at the same values from the first step but optimize a f , b f with the following objective function: However, since the sensitivities of a s and b s are very low, the optimization procedure is still unable to update these without further information.
In a healthy myocardium, we know that the stiffness in the sheet direction should be much smaller than that in the myofibre direction; hence, we introduce the following constraints: in the optimization procedure to allow a s and b s to be updated. The SQP method in MATLAB is then employed to minimize f O2 with the constraints (Eq. 8). A further step is used to ensure that the volume matching is still satisfied by applying the objective function f O1 once more. Since Table 2 indicates that f O1 is highly sensitive to variations in a and a fs , which are also correlated, we may reduce one of these two parameters by introducing C 3 , so that a = C 3 a and a fs = C 3 a fs .
Now we only need to estimate C 3 in the third step. The remaining parameters are kept the same from Step 2. The overall optimization procedure for the in vivo model is illustrated in the flowchart in Fig. 3. Step 1: optimizing C a , C b with f O1 Step 2: optimizing a f , b f with f O2 using the constraints (Eq. (8)), and fixed a, b, a fs , b fs Step 3: optimizing C 3 with f O1 Case 1: estimation from the proposed optimization procedure Case 2: estimation using only the strain data in objective function Case 3: estimation using only the first principal strain and volume data Case 4: estimation terminated after the first step Table 4 shows the inversely estimated eight parameters from the proposed multi-step optimization procedure with both f vol and f ε (Case 1). These are fairly close to the so-called true parameters from [35], though some discrepancy is seen in the values of a s and b s . This is because the objective functions are insensitive to changes in a s and b s .   Table 4 also shows the results when using f ε only for the multi-step optimization procedure (Case 2). Compared to the proposed multi-step optimization (Case 1), the parameters optimized in Case 2 are less accurate, especially for b f . This highlights the importance of using volume measurements.

Synthetic model
Case 3 in Table 4 shows the estimated parameters when we use the proposed multi-step optimization procedure but with only the first principal strains (k = 1st). In this case, the discrepancies in a s , b s and b fs are greater, though the remaining parameters are still close to the so-called true values. Case 4 arises when all the strains and volume data are used but only the first step of the optimization is performed. The values of a f , b f , a s and b s are all less accurate in this case. In summary, the proposed multi-step sequential optimization procedure is the best approach of all cases considered.

In vivo estimation
The parameters for the in vivo human LV model are now estimated using the proposed optimization procedure in Fig. 3 with 8 mmHg ED pressure. The mapping of the objective function f O1 , with C a and C b varying in a range of (0.1, 1), is shown in Fig. 4. From the centre of the valley, we select the point (C a = 0.2, C b = 0.3). The final estimated parameters are summarized in Table 5, with the predicted strains compared with the measurements Fig. 5 Regional circumferential strain at end of diastole after each optimization step (strain is calculated related to early diastole: the reference state) Fig. 6 Distributions of myofibre stress at 8 mmHg endocardial pressure for in vivo LV model using a initial parameters and b optimized parameters following each step. Clearly, the original parameters used by [35] are over-stiff for the human subjects and result in a much smaller ED LV cavity volume compared to the measured ED volume. This over-stiffness in the parameter set is also noticed in [12]. Figure 5 plots the detailed regional circumferential strains at ED after each step. After Step 3, the simulated regional circumferential strain is in good agreement with the measurements. Figure 6 shows the myofibre stress distributions within the LV at ED using both the initial and optimized parameters. Clearly, the initial parameters yield a very different stress pattern compared to the optimized ones.
Since we use a population-based end-diastolic pressure (EDP) in the parameter estimation, we check the change of the parameters when the EDP is varied from 10 mmHg to 16 mmHg; this is summarized in Table 6. Changes in the myofibre stress-strain relationships, when estimated with different EDPs, are plotted in Fig. 7. In general, with increased EDP, the stiffness along the myofibres increases.
In Fig. 7, we further compare the myofibre stress-strain relationship with those from other studies which inversely estimated human myocardial parameters based on different constitutive laws and EDPs. Xi et al. [42] used a healthy subject and 13.6 mmHg EDP; Wang et al. [29] reported the averaged healthy myocardial stiffness from six subjects  [12] estimated the functional myocardial stiffness from one LV dysfunctional patient with an EDP of 15 mmHg. Although there are necessary discrepancies due to the subject variety and different constitutive laws used, the overall trend of the mechanical responses is similar. As for most of the inverse problems, the uniqueness of inversely estimated coefficients cannot be guaranteed. To test the robustness of the procedure, however, we look at seven more cases: (1) using different values of C a and C b ; (2) increasing a ini f by half; (3) increasing b ini f by half; (4) decreasing a ini f by half; (5) decreasing b ini f by half; (6) adding Gaussian noise to the measured circumferential strains with a half standard deviation of the measured values; and (7) as in (6) but using one standard deviation of the measured values. Table 7 summarizes the estimated parameters from all seven cases, and the corresponding myofibre stress-strain relationships for a cubic myocardial  Table 7 under uni-axial tension Table 7 Estimated parameters for uncertainty analysis tissue under uni-axial tension are shown in Fig. 8. We notice that although there are some differences in the estimated parameters for the different cases, there seems to be reasonable agreement in the predicted myofibre stress-strain relationships (Fig. 8). This is confirmed by applying the optimization procedure to two more healthy volunteers (see the appendix for details). The myofibre stress distributions obtained using different sets of parameters in Table 7 are all found to be similar to the stress distribution in Fig. 6b (and hence not shown). Since the myofibre stressstrain relationship largely determines the LV responses, and given the extreme difficulty in the in vivo parameter estimation, the proposed optimization procedure is deemed to be sufficiently robust.

Discussion
Estimating the material parameters of myocardial constitutive laws from limited in vivo data remains a major challenge due to the non-linearity of the model responses and strong intercorrelations between the material parameters [2]. Furthermore, different constitutive laws have different parameter sets and behaviours, and specific treatment in the inverse problem is often required for a given constitutive law. Although the structure-based Holzapfel-Ogden law is gaining popularity in heart modelling, to the best of our knowledge, no studies have been performed to estimate the parameters of the Holzapfel-Ogden law from in vivo measurements. In this study, inspired by the ideas from previous studies [20,27,30], we first developed an optimization procedure for estimating the parameters of the Holzapfel-Ogden law using a synthetic model to produce the so-called true values. The synthetic model allows us to investigate the estimation errors and verify the correctness of the solution.
By performing the sensitivity analysis, we show that even with a large amount of strain and pressure-volume data, it is difficult to accurately estimate the sheet parameters a s and b s , which are usually very small [20]. We also find that by using a multi-step and sequential optimization procedure we can achieve much higher accuracy compared to a single-step optimization (Table 4).
When applying the approach to in vivo LV models, we have encountered difficulties since there are insufficient data from clinical measurements. Various constraints must be introduced to reduce the complexity of the problem. One constraint is to assume that the eight parameters can be grouped into two (one with exponential terms, and one without) and scaled to the corresponding groups of the original parameters [35]. In this way, we only need to estimate two scaling factors, C a and C b , in the first step. In addition, the set which minimizes |C a − C b | is chosen when multiple sets of C a and C b are available from the sweeping procedure. This also avoids over-stiffness responses. A similar approach is used by [12], though for a different constitutive law.
Furthermore, constraints (Eq. 8) must be introduced when estimating a f and b f using normalized volume and strain data. Because the measured circumferential strain from cine images is in the circumferential-longitudinal axis plane, which aligns with the myofibre direction, this procedure can effectively update the myofibre stiffness parameters, a f and b f . However, we are unable to accurately estimate a s and b s , which are related to myofibre sheet strains. Measurements of sheet strains are extremely difficult in vivo. In addition, the synthetic model study shows that even with all the principal strains, a s and b s have very low sensitivities and cannot be estimated properly. Hence, additional constraints are applied to a s and b s in the second optimizing step.
The sensitivity analysis shows that the parameters in the Holzapfel-Ogden law are highly intercorrelated; for example, an increase in a can be compensated by a decrease in b or other correlated parameters. This leads to uncertainties in the parameter estimation. This is a common issue for all anisotropic constitutive laws. Xi et al. [2] found that, with their synthetic LV models, a unique solution for estimating the seven parameters of Costa's law [33] cannot be achieved using a reduced-order unscented Kalman filter. Therefore, they used the transversely isotropic Guccione's law with a lower level of complexity.
Instead of changing the constitutive law, we reduce the complexity of the problem by estimating a total of five parameters only, C a , C b , a f , b f , C 3 . Even with the reduced set of parameters, it is not possible to establish the global minimization or the uniqueness of the solution given the ill-posed nature of the inverse problem. However, we have tested various cases to show that more or less the same mechanical responses in the physiological range could be achieved even though the parameters are somewhat different. One possible explanation for the same stress-strain relationship from different parameter values of the same constitutive law is that the law is based on coupled strain attributes [43].
In this study, only cine CMR imaging is used to provide data for the in vivo LV model. Although dedicated strain CMR imaging is able to provide 3D LV deformation, this requires complex image processing and additional scanning time, which may not be possible for some patients. Cine CMR images are widely available from routine scans. Hence the optimization procedure proposed based on cine CMR imaging can be readily used for clinical applications. The downside of the approach is that accurately estimating LV motion/strain from cine CMR images is more difficult because cine images are usually 2D images; therefore, the out-plane motion cannot be easily estimated. In addition, due to the lack of patterns or features for motion tracking, large uncertainties exist when estimating the pixel-wise strain. Our previous study [41] showed that regional circumferential strains estimated from cine images using a deformable image registration method compared well with those from dedicated strain CMR imaging for both healthy volunteers and patients with myocardial infarction, but greater discrepancies existed in the estimated regional radial strains. Thus, in the in vivo LV model, only regional circumferential strains are used for the objective function.
The dilemma is that, while fewer data make the inverse problem more ill-posed, demanding more data means more clinical measurements with longer acquisition times, which is not always possible. In a patient CMR study, only necessary measurements are performed routinely. Furthermore, CMR image at late-diastole is difficult to achieve, while the ED frame is always recorded with high quality; therefore, it is desirable to use information from the ED frame for material parameter estimation. However, issues associated with inversely estimating material properties from an ED frame have been reported in several studies [26]. Xi et al. [42] found that with an ED measurement, the parameters cannot be uniquely constrained, although such a measurement can provide a potentially robust indicator  [25] showed that with five frames of CMR data, the parameters from Guccione's law could be inversely estimated within a 5% margin of error. In our study, because of the lack of an in vivo EDP recording and strain CMR, we must rely on the limited ED frames for the parameter estimation. Therefore, rather than aiming to obtain the unique parameters, which is not possible, we try to extract the myocardial stiffness. It goes without saying that with more data available from routine CMR measurements, the parameters could be estimated with higher accuracy. This has been demonstrated by the synthetic model: the inversely estimated parameters are fairly close to the so-called true parameters.
In the in vivo LV model, an original set of parameters is required. We used a set obtained by fitting the simple shear tests on healthy swine myocardial samples [15,35]. These original parameters might affect the accuracy of inverse estimation. Parameters based on human myocardial samples may be a better candidate, but these are not yet available.
Validating the inversely estimated parameters in the in vivo LV model is difficult since it is impossible to perform mechanical tests on in vivo hearts. Previous studies tried to compare results with those of other studies [2,28]. Our estimated parameters are shown to be in the same range as in [12,29,42].
An important issue in LV modelling is the choice of a suitable constitutive law for the myocardium. In this paper, we choose the well-established Holzapfel-Ogden model. We are aware that the Holzapfel-Ogden law makes no use of the deformation invariants I 5 (= f 0 C 2 f 0 )) and I 7 (= s 0 C 2 s 0 ), which has recently been shown to be incompatible with linear elasticity [44,45]. However, it is well known that if you reduce the number of invariants, then you cannot fully capture linearly elastic responses, but the consequences of this should be viewed with extreme caution since experiments on soft tissue in the small-strain regime, by the very nature of the material, are not very accurate or reliable. Indeed, many of the current invariant-based laws do a very good job of fitting the data of a wide range of soft tissues to a wide range of deformations -and the fact that they do not fully capture linear elasticity is, under the present state of experimental knowledge of the mechanics of soft tissues, largely irrelevant. The Holzapfel-Ogden law, in particular, is the model that can fit all the data of the simple shear experiments well [15], and the computed LV dynamics based on the Holzapfel-Ogden law seem to predict the measured strains and volumes well within the physiological range. In addition, including more invariants will greatly increase the complexity of the inverse problem and make it extremely hard, if not impossible, to estimate all the parameters from in vivo measurements.
Other limitations of the work in common with many published LV models are as follows: (1) the in vivo LV geometry is reconstructed at ED, which is not stress-free, and residual stress is not considered [39]; (2) the heterogeneous distribution of material properties is not considered in the LV models; and (3) the proposed method can be potentially extended to diseased heart tissue with some changes in the optimization procedure to account for a remodelled micro-structure, such as myocardial infarction [46].

Conclusion
In this study, we have investigated, for the first time, the feasibility of inversely estimating parameters in the orthotropic Holzapfel-Ogden constitutive law for passive myocardium by proposing a multi-step optimization procedure using both strain and pressure-volume data. When applied to a synthetic LV model, the estimated parameters are very close to the known parameters, although some uncertainties exist in estimating parameters along the sheet direction due to the low sensitivity. For parameter estimation of in vivo models, a more comprehensive sensitivity study is performed due to the limited measurement data. The material parameters are scaled from the original parameters based on ex vivo experimental tests. A study of the sensitivity is also used to reduce the complexity of the problem. By matching the regional circumferential strains and pressure-volume at the ED frame, we have demonstrated that the parameters of the Holzapfel-Ogden law, and in particular the myofibre stress-strain relationship, can be estimated successfully when suitable constraints are introduced for the in vivo model.

Appendix
The proposed optimization procedure is further applied to two more healthy volunteers [volunteer 2 (male, age 22) and volunteer 3 (male, age 31)]. The LV geometry is reconstructed from CMR studies similarly to the volunteer study in the main text (volunteer 1, male, age 28), and 8 mmHg EDP is assumed. Table 8 contains the inversely estimated parameters of the Holzapfel-Ogden law. Parameters are similar for all three volunteers, though minor differences exist. Table 9 compares the ED regional circumferential strain in the three volunteers with strains from cine images, which shows that with optimized parameters, the simulated strains agree well with the measurements. Figure 9 in Appendix shows the corresponding myofibre stress-strain relationships for all three volunteers under uni-axial tension. Again, the mechanical responses for these healthy volunteers are very similar, particularly at the small strain (Fig. 9).