# Adjoint Multi-Start Based Estimation of Cardiac Hyperelastic Material Parameters using Shear Data

## Abstract

Cardiac muscle tissue during relaxation is commonly modelled as a hyperelastic material with strongly nonlinear and anisotropic stress response. Adapting the behavior of such a model to experimental or patient data gives rise to a parameter estimation problem which involves a significant number of parameters. Gradient-based optimization algorithms provide a way to solve such nonlinear parameter estimation problems with relatively few iterations, but require the gradient of the objective functional with respect to the model parameters. This gradient has traditionally been obtained using finite differences, the calculation of which scales linearly with the number of model parameters, and introduces a differencing error. By using an automatically derived adjoint equation, we are able to calculate this gradient more efficiently, and with minimal implementation effort. We test this adjoint framework on a least squares fitting problem involving data from simple shear tests on cardiac tissue samples. A second challenge which arises in gradient-based optimization is the dependency of the algorithm on a suitable initial guess. We show how a multi-start procedure can alleviate this dependency. Finally, we provide estimates for the material parameters of the Holzapfel and Ogden strain energy law using finite element models together with experimental shear data.

## 1Introduction

The personalization of computational models in cardiology is a key step towards making models useful in clinical practice and cardiac surgery. A computational model, once properly calibrated, has the potential to forecast cardiac function and disease, and can aid in planning treatments and therapies. To describe the mechanical function of the heart, the passive elasticity of the muscle tissue needs to be represented. Personalizing the effects of this elasticity in a computational model is typically accomplished by tuning a set of material parameters so that the output of the model fits observed data. Gradient-based optimization algorithms have successfully been used in the past to automatically perform the parameter tuning at an organ scale [2]. In these studies, the gradient of the objective functional is approximated using one-sided finite differences.

Compared to using a global optimization method, local gradient-based methods have the advantage of using relatively few optimization iterations. This is an important consideration when optimizing organ scale finite element models, for which running a single forward model can take hours or days. On the other hand, a disadvantage of using local optimization methods is the fact that they can converge to local, globally suboptimal, minima. One way to combine the speed of a local optimization with the robustness of a global optimization is to use the multi-start method. In this method, many local optimizations are run starting from various points in parameter space and the best fitting solution of the group is taken to be the global optimum.

Another popular approach to parameter fitting is the reduced order unscented Kalman filter. This approach was successfully used to fit a transversely isotropic passive mechanics model to synthetic data [31], to partially calibrate a multi-physics model [21], and to estimate regional contractility parameters [5]. Note however that the use of both unscented Kalman filtering and finite differences carries a computational cost that increases with the number of model parameters.

Assuming there are parameters to be estimated, an unscented Kalman filter with a minimal sigma-point configuration requires model evaluations at a single time level for each assimilated data point. An evaluation of a finite difference derivative on the other hand requires runs of the model throughout the full span of model configurations considered.

In contrast to these two techniques, the adjoint approach computes the objective functional gradient via the solution to an adjoint equation, which involves only a single solve of a linearized system for any number of model parameters. Thus, for models involving many parameters, either due to model complexity or spatiotemporal parameter variation, the adjoint approach offers a computationally attractive approach for parameter estimation.

There are some previous results involving adjoint equations and cardiac elasticity. Sundar et al. (2009) developed a framework for the estimation of wall motion based on cine-MRI images and adjoint inversion [26], and Delingette et al. (2012) used an adjoint equation to estimate contractility parameters [7]. However, both of these studies involve linear and isotropic elasticity models, which represent a significant simplification of the orthotropic and highly nonlinear behavior reported in the contemporary cardiac mechanics literature [6].

One reason why it is difficult to use an adjoint equation with modern nonlinear anisotropic models is the complexity required in deriving and implementing code for the solution of the adjoint problem. In order to resolve this issue, we make use of an automatic framework for generating adjoint code [9]. Here, we use this adjoint framework to estimate the material parameters of an invariant-based orthotropic myocardial strain energy law (the Holzapfel-Ogden model) [15]. This law is embedded here in an incompressible finite element framework, and we use the raw data from a simple shearing experiment [8] as a target for optimization. These data have previously been used to estimate material parameters for a variety of other strain energy functions using a finite element framework, but with a gradient obtained using finite differences [24]. The material parameters of the particular strain energy density that we are using have also been previously estimated using digitized data based on Figure 6 of [8], and a homogeneous deformation model [15]. Our study is however the first to use the adjoint approach for the estimation of cardiac hyperelasticity parameters and the first to provide optimized material parameters for the incompressible Holzapfel-Ogden model for non-homogeneous shear deformations.

The rest of this paper is organized as follows. In Section 2 we describe the variational formulation of the elasticity model, the optimization problem for identifying the material parameters, and how the adjoint gradient formula can be used to calculate a functional gradient. In Section 3 we describe the verification of the forward and inverse solvers, present timings to show the efficiency of the adjoint method, and show the results of parameter estimations. Finally, we test a multi-start optimization method in order to reduce the dependence of the gradient-based algorithm on the choice of initial parameter set. We conclude by discussing our findings in Section 4, and drawing some conclusions in Section 5.

## 2Mathematical models and methods

We shall use the notion of the directional derivative frequently throughout. For a functional for some vector space , we define the directional derivative of with respect to the argument named in the direction

Furthermore we denote the total derivative by the usual notation to mean the derivative of with respect to all arguments depending on .

### 2.1Hyperelasticity model

Let be an open and bounded domain with coordinates and boundary , occupied by an incompressible hyperelastic body. We consider the quasi-static regime of a body undergoing a large deformation and are interested in finding the displacement and the hydrostatic pressure that minimize the incompressible strain energy :

over the space of admissible displacements and pressures satisfying any given Dirichlet boundary conditions. In , is a set of material parameters, , where denotes the deformation gradient, is the identity tensor in , denotes a volume-preserving right Cauchy-Green strain tensor, and denotes an isochoric strain energy density.

The incompressible Holzapfel and Ogden hyperelasticity model [15] describes large deformations and stresses in cardiac tissue via the following energy density :

Here denote fiber and sheet directions, respectively; is a Heaviside function with a jump at , and the material parameters are

Moreover, are rotation invariant functions given by

where denotes the tensor trace and denote unit vectors pointing in the local myocardial fiber and sheet directions [15]. The strain energy density is rotation-invariant, and polyconvex if [15].

The Euler-Lagrange equations for the minimizing displacement and pressure of read: for given , find such that

for all admissible virtual variations . Inserting the total potential energy from and taking the directional derivatives, we obtain

### 2.2Parameter estimation as a PDE-constrained optimization problem

In the general case, the passive material parameters entering the constitutive relationship are not known. In order to estimate these parameters from data, we propose to use a numerical approximation in combination with a gradient-based optimization algorithm in which the gradients are computed via an adjoint model. The optimization algorithm seeks to minimize the misfit between model output and observations. Denoting the misfit functional by , the optimization problem reads:

together with suitable Dirichlet boundary conditions on . We also require that to ensure the functional is polyconvex [15]. For notational convenience we will sometimes use the reduced formulation of the misfit functional and its gradient with respect to the material parameters . In particular, we introduce the reduced functional

In our numerical experiments we use Sequential Least Squares Programming (SLSQP) as implemented in [18] and wrapped in the package SciPy [17] in order to solve .

### 2.3Multi-start Optimization

A common challenge with gradient-based algorithms is that the solution obtained depends on the choice of initialization point for the algorithm. Moreover, the optimized solution may be a local minimum only and not necessarily a global minimum. One way to attack these issues is to run many optimizations from randomly chosen initial parameter points, and to chose the resulting optimized material parameter set that gives the best fit. This method is often referred to as multi-start optimization [4] and is an example of combining global and local optimization.

Due to the presence of exponential functions in the strain energy , it is possible for calculated stresses to become very large, which may result in convergence issues for the numerical solution of the Euler-Lagrange equation . This can easily occur if several material parameters have large values. In order to minimize this problem we have designed a procedure to generate random initial guesses which limits the number of large material parameter values while still allowing for a large range of initial possible values for each parameter. The procedure works as follows: first set a maximum parameter value . Then choose (with in our case) points , from a uniform distribution defined over the interval and let . The parameter values are then set to be the distances between successive randomly drawn points, that is .

### 2.4Computing the functional gradient via the adjoint solution

Gradient-based optimization algorithms in general, and the SLSQP algorithm in particular, rely on the total derivative of the objective functional . By introducing an *adjoint* state variable, this derivative may be computed efficiently. We summarize this result below. Our presentation is based on [14], and is adapted here to the solid mechanics setting.

We define three abstract spaces , , and , where is the space of all possible solutions to the variational equation which also satisfy any given Dirichlet boundary conditions, is the material parameter vector space, and is the space of virtual variations. The Lagrangian is defined as:

For all , solving the state equation , we have

such that the total derivatives of and coincide,

If we choose such that

for all , which in particular includes , the total derivative of with respect to in the direction simplifies as follows using the chain rule:

Then, for any infinitesimal variation in the material parameters , combining , , and yields an efficient evaluation formula, not requiring derivatives of the state variable with respect to the material parameters , for the total derivative of :

We still need to compute . By defining the form and its adjoint ,

we can rewrite as

and thus recognize the adjoint equation: given , , find such that

for all .

In summary, the adjoint-based gradient evaluation formula is: given , first compute by solving the state equation , next compute by solving , and finally evaluate .

### 2.5Description of shearing experiments

We aim to optimize the material parameters of the Holzapfel-Ogden model with respect to target experimental data, in particular data resulting from an earlier set of simple shearing experiments [8]. In these experiments, pig hearts were extracted. From each heart, three adjacent cubic blocks were cut in such a way that the sides of the cubes were aligned with the local myocardial fiber and sheet directions. A device held two opposing faces of each cube between two plates using an adhesive. The top plate was displaced in order to put each specimen in simple shear. For each specimen 6 different modes of shear were tested. These modes are described using the coordinate system, which refer to the myocardial fiber, sheet and sheet normal directions, respectively. Each mode is denoted by two letters, where the first defines the normal of the face of the cube that is being displaced, and the second refers to the direction of displacement. These 6 modes are .

In order to remove the effects of strain softening, preliminary displacements were applied to the tissue samples until no further softening was observed. After that, displacements were once again applied, and the forces in the shear direction were measured on the top plate. These measurements were taken for circa various states of shear per mode.

In Figure 1 we display the stress-strain relations for positive displacements that were obtained from the shearing experiments [8]. As can be seen in Figures 4 and 6 of [8] the experimentally obtained curves contain a high degree of symmetry through the line . We can expect the same symmetry in the stresses computed by finite element models which use the strain energy since changing the sign of the displacement map will change the sign of the resulting stresses but preserve their magnitude. In the previous studies [15], [12], and [28], only the data for positive shear displacements were used. For the sake of comparability, we restrict our data in the same way.

In our numerical experiments we use two data sets with reference to the numbering of [8]. The first is Data Set 6, and the second data is Data Set 2 with the and curves swapped. This swap and the choice of data sets are discussed further in Section 4. For clarity, we shall refer to Data Set 6 as “transversely isotropic” and Data Set 2 with the swap as “orthotropic”, as the respective stress-strain curves are typical of materials of these types. For each mode, the prescribed shear displacement is modelled as a Dirichlet boundary condition for the displacement on the respective top and bottom faces in the respective direction.

### 2.6Choice of objective functional

In order to estimate the passive material parameters of the Holzapfel-Ogden model, we make use of a least squares objective functional. This functional defines a distance from the model output to the data points of the shearing experiment, and we seek the material parameter set that minimizes this. Before introducing our objective functional, we define the set of directions , referring to fiber, sheet and sheet normal directions. We also use the notation to refer to a mode, with the index referring to the normal of the face that is shifted, and to the direction in which the shift occurs.

Our fit function is similar to that used in [23], and is given by

In , is the force measured during the experiment, and is the force generated by the finite element model at each prescribed shear displacement , where is the maximal prescribed displacement of the mode in the experiment. Each is chosen to be a Gauss point of a -point Gauss integration rule defined over , and is the value of the Gauss weight related to . Explicitly, for mode with top face , is given by

where is a shear component of the deformation gradient. Evaluating the inner loop of requires solving once for each given shear displacement . The motion given by the calculated displacements is then a quasi-static approximation of the motion undergone by the corresponding tissue in the shearing experiment.

Following [23], we evaluate the least squares fit at Gauss integration points, rather than for all 250 recorded points for each shear mode, in order to greatly reduce the computational expense of evaluating . At each Gauss point we obtain the corresponding shear stress by linearly interpolating between the two neighbouring stresses which were recorded in the experiments of Dokos et al. [8].

The use of Gauss integration is based on the observation that is an approximation to the following expression

By setting and approximating the integral by the midpoint rule applied to the full dataset we can determine the quality of the Gauss approximation. In order to do this we define the relative error

where is the midpoint rule approximation of evaluated over the full data, and , given by , is evaluated at a reduced set of Gauss points. We noticed that 9 Gauss points are sufficient to reduce to less than . However, in our numerical experiments we use Gauss points as this guaranteed small enough changes in the solution of the Euler-Lagrange equation from one Gauss point to the next, so that our Newton’s method solution of always converged.

### 2.7Finite element discretization of the hyperelasticity equations

We represent each tissue sample of the shearing experiments by a three-dimensional cube (). An mesh of this cube was constructed by uniformly dividing the mesh into boxes and then subdividing the boxes into tetrahedra. The local myocardial fiber and sheet orientations were represented as spatially constant vectors aligned with the coordinate axes.

On these geometries, we solve and its adjoint, using a Galerkin finite element method with the Taylor-Hood finite element pair [16]; e.g. a continuous piecewise quadratic vector field for the displacement and a continuous piecewise linear scalar field for the pressure. For the solution of the nonlinear system of equations, we use a Newton trust region method. The absolute tolerance of the nonlinear solver was set to in the numerical experiments below. Linear systems are solved by LU factorization.

Additionally, we model the case of a homogeneous deformation which corresponds to a linear displacement with a constant shear angle throughout the domain. Such a model can be represented by discretizing the cubes with a single layer of linear finite elements: the resulting displacement is completely determined by the prescribed boundary conditions. Figure 2 illustrates the two kinds of deformations on cube meshes.

The discrete variational formulation of the Euler-Lagrange equations is implemented using the FEniCS Project software [1] and dolfin-adjoint [9]. From a FEniCS forward model, dolfin-adjoint automatically generates the symbolic adjoint system of equations and computes the functional gradient using the adjoint solution. The FEniCS framework automatically generates and compiles efficient C++ code for the assembly of the relevant linear systems from the symbolic representations of both forward and adjoint equations, and solves the nonlinear and linear systems using e.g. PETSc [3]. With this setup, we observed that a typical solution of the Euler-Lagrange equation takes Newton iterations.

## 3Numerical Results

### 3.1Verification

Each of the finite element, adjoint, and optimization solvers have been carefully verified, separately and combined, as follows:

(i) The finite element solver was verified by the method of manufactured solutions [22]. Following this method we chose an analytic expression for the displacement and pressure fields

Here refer to Cartesian coordinates and is a scaling parameter which we set to . Using this analytic expression we derived Dirichlet boundary conditions over a unit cube, and a loading term which satisfied a pointwise form of equation

Note that the chosen displacement field satisfies the incompressibility constraint . We then computed finite element approximations to and observed the expected second-order convergence of the displacement gradient to the analytical displacement gradient [16].

(ii) We verified the computation of stresses in the finite element model by prescribing a homogeneous deformation and comparing the resulting numerically integrated top face shear stress values to analytically computed values. The analytic values were based on the calculations found in [15] and the numerical values were observed to match closely.

(iii) We confirmed the correctness of the adjoint gradients by considering the linearization of the functional around with perturbation and using Taylor’s theorem: the expression

converged to at a rate of as , which can only be expected if is computed accurately.

### 3.2Parameter estimation with synthetic data

Additionally, we verified the optimization solver by performing a synthetic data test. In this test we chose a target set of material parameters, Table ?, 2nd line, and used them to compute synthetic integrated stress values for all 6 shear modes of the tissue experiment [8]. These synthetic stresses were then matched by an optimization starting from material parameter values higher than the target.

We performed this test using our two models for deformation. The first model assumed a homogeneous shear angle through the material and the second model was a finite element model with a mesh. Since the displacement field of the finite element model was element-wise quadratic, it allowed for more flexibility in the deformation field. The results of this synthetic data test are presented in Table ? and show that the optimization algorithm was able to closely match the target material parameters.

(kPa) | (kPa) | (kPa) | (kPa) | (mN) | |||||

Initial | 0.059 | 8.023 | 18.472 | 16.026 | 2.481 | 11.120 | 0.216 | 11.436 | |

Target (80%) | 0.047 | 6.418 | 14.778 | 12.821 | 1.985 | 8.896 | 0.173 | 9.149 | |

Homogeneous | 0.047 | 6.418 | 14.778 | 12.821 | 1.985 | 8.896 | 0.173 | 9.149 | 4.611 |

Finite Element | 0.047 | 6.406 | 14.778 | 12.821 | 1.983 | 8.938 | 0.173 | 9.155 | 0.00082 |

### 3.3Parameter estimation with experimental stress data

In the following, we present the results of fitting the Holzapfel-Ogden strain energy law using the objective function and a SLSQP optimizer with bound constraints. The SLSQP algorithm makes use of the gradient of the objective functional which we obtain using the adjoint gradient formula .

As the numerical solution of the nonlinear Euler-Lagrange equation easily fails to converge when a material parameter becomes too small, we set a lower bound of on the components of while optimizing finite element models. This bound was not necessary for the homogeneous deformation models as no Euler-Lagrange equation is solved. All optimizations were carried out until the optimizer was unable to further reduce the objective functional or an absolute tolerance of in the 2-norm of the functional gradient was reached.

#### Material parameter estimation using *a priori* knowledge

The material parameters of the Holzapfel-Ogden model have previously been estimated using a homogeneous deformation model (Table 1, 2nd row in [15]). We first used these values as the initial values for optimization of our homogeneous model targeting the transversely isotropic and orthotropic data sets. The optimized results are listed in Table 1 with the label Homogeneous.

We next consider finite element models that allow for heterogeneous shear displacements. Beginning with a cube and the optimal material parameters from the homogeneous model as initial values, we computed optimal values for the case. This procedure was repeated for cubes with , using the results of the previous optimization as the initial condition for the next case. The resulting parameter values are presented in Table 1, and the corresponding optimal stress-strain curves are shown in Figure 3.

We note that going from to using both the transversely isotropic and the orthotropic data does not change the material parameters rounded to two 2 significant digits, and therefore consider our finite element models to be sufficiently refined at this resolution. We also note that the fit values, , decreased with mesh refinement up to about 2 digits accuracy. We expect this decrease since increased mesh refinement gives more flexibility in the deformation field of the finite element model.

Ev. | Grad | ||||||||||

(kPa) | (kPa) | (kPa) | (kPa) | (mN) | Ev. | ||||||

Transversely Isotropic |
|||||||||||

Homogeneous | 0.544 | 6.869 | 23.220 | 39.029 | 0.0001 | 0.172 | 0.248 | 5.310 | 3.291 | 41 | 21 |

N = 1 | 0.593 | 6.841 | 23.209 | 38.826 | 0.010 | 0.010 | 0.243 | 9.531 | 3.173 | 44 | 37 |

N = 2 | 0.732 | 6.818 | 22.110 | 39.946 | 0.010 | 0.010 | 0.183 | 13.614 | 3.010 | 24 | 18 |

N = 4 | 0.807 | 6.737 | 21.349 | 40.468 | 0.010 | 0.010 | 0.122 | 17.936 | 2.819 | 25 | 18 |

N = 6 | 0.794 | 6.859 | 21.212 | 40.537 | 0.010 | 0.010 | 0.129 | 17.462 | 2.802 | 22 | 15 |

N = 8 | 0.784 | 6.973 | 21.149 | 40.584 | 0.010 | 0.010 | 0.145 | 16.401 | 2.815 | 21 | 14 |

N = 10 | 0.778 | 7.048 | 21.112 | 40.585 | 0.010 | 0.010 | 0.150 | 16.036 | 2.819 | 24 | 17 |

Orthotropic |
|||||||||||

Homogeneous | 0.556 | 7.940 | 33.366 | 14.224 | 2.804 | 0.0001 | 0.588 | 8.216 | 6.804 | 31 | 20 |

N = 1 | 0.766 | 6.857 | 31.640 | 15.210 | 2.069 | 0.010 | 0.352 | 15.243 | 5.880 | 29 | 19 |

N = 2 | 1.040 | 6.557 | 29.375 | 15.979 | 1.742 | 0.010 | 0.118 | 23.296 | 4.565 | 39 | 24 |

N = 4 | 0.979 | 7.364 | 28.882 | 15.813 | 2.058 | 0.010 | 0.107 | 24.039 | 3.952 | 28 | 16 |

N = 6 | 0.961 | 7.495 | 28.762 | 15.783 | 2.088 | 0.010 | 0.114 | 23.549 | 3.899 | 21 | 13 |

N = 8 | 0.962 | 7.510 | 28.649 | 15.806 | 2.044 | 0.010 | 0.122 | 23.027 | 3.899 | 20 | 11 |

N = 10 | 0.959 | 7.542 | 28.565 | 15.813 | 2.017 | 0.010 | 0.123 | 22.750 | 3.981 | 25 | 12 |

#### Material parameter estimation using multi-start optimization

In this section, we present the results of using the multi-start method to estimate the optimal material parameters, rather than relying on a good initial guess. For the calculation of random initial guesses we set , cf. Section Section 2.3. This value is close to the largest material parameter found in Table 1. Note that this choice gives a conservative set of initial parameters for the optimization algorithm (low initial values) which in turn enhances the robustness of the procedure. We also set as an upper bound for each material parameter value during the optimization. Without this upper bound we observed that many optimizations crashed or converged to suboptimal local minima.

In each multi-start experiment, random starting points were used. The mesh fineness was set to the level of , which was sufficient to give converged material parameter sets when using a priori knowledge in Section ?. In Table ? we present the best fitting results of the multi-start experiments and note that they are very close to those obtained with a priori knowledge in Table 1.

I | |||||||||

(kPa) | (kPa) | (kPa) | (kPa) | (mN) | |||||

Transversely Isotropic |
|||||||||

N = 8 |
0.784 |
6.973 |
21.149 |
40.584 |
0.010 |
0.010 |
0.145 |
16.401 |
2.815 |

Multistart Best Fit | 0.795 | 6.855 | 21.207 | 40.545 | 0.010 | 0.010 | 0.130 | 17.446 | 2.802 |

Orthotropic |
|||||||||

N = 8 |
0.962 |
7.510 |
28.649 |
15.806 |
2.044 |
0.010 |
0.122 |
23.027 |
3.899 |

Multistart Best Fit | 0.964 | 7.510 | 28.654 | 15.791 | 2.051 | 0.010 | 0.118 | 23.230 | 3.959 |

#### Objective functional values for alternative material parameters

Several other studies [15] have used of the Dokos et al. 2002 shear data [8] to calibrate the Holzapfel and Ogden strain energy . These studies used homogenized deformation models for the optimization. In Table ? we list the computed objective functional value of parameter sets originating from previous studies using the orthotropic dataset and finite element model . The results indicate that our parameter set fits these data better than the previously computed ones.

We also note that our finite element parameter set with finite element model has a better fit value than the homogeneous parameter set with the homogeneous model. Indeed we expect the finite element fit to be at least as good as the homogeneous fit, as the finite element model allows for greater flexibility in the the deformation field, above and beyond that of the homogeneous model.

Source | ||||||||||

(kPa) | (kPa) | (kPa) | (kPa) | (mN) | (mN) | |||||

Holzapfel et al 2009 | 0.059 | 8.023 | 18.472 | 16.026 | 2.481 | 11.120 | 0.216 | 11.436 | 36.143 | 36.825 |

Goektepe et al 2011 | 0.496 | 7.209 | 15.193 | 20.417 | 3.283 | 11.176 | 0.662 | 9.466 | 28.583 | 29.480 |

Wang et al 2013 | 0.2362 | 0.810 | 20.037 | 14.154 | 3.7245 | 5.1645 | 0.4108 | 11.300 | 33.271 | 34.195 |

Current (hom) | 0.556 | 7.940 | 33.366 | 14.224 | 2.804 | 0.0001 | 0.587 | 8.216 | 6.804 |
9.653 |

Current (fem) | 0.962 | 7.510 | 28.649 | 15.806 | 2.044 | 0.010 | 0.122 | 23.027 | 41.622 | 3.899 |

### 3.4Computational efficiency of the adjoint-based functional gradient

Adjoint solver efficiency may be measured by comparing the runtime of the adjoint and forward solves. Here, we examine the overall gradient efficiency in a similar manner. We consider the evaluation of the gradient of the objective functional , though in a reduced case with only a single shear mode included in the sum and a reduced forward solve consisting of a single nonlinear solver iteration. In this case, the forward and adjoint models each consist of a single linear solve in addition to a number of residual evaluations. For larger linear system sizes, the runtime of a linear solve is expected to dominate the runtime of assembly, and thus these forward and adjoint models are of roughly the same computational expense.

For this reduced case, we evaluated the adjoint-based gradient for a range of linear system sizes. For each system size, we calculated the gradient runtime ratio; that is, the runtime used by the evaluation of the gradient divided by the runtime of the forward solve. The resulting ratios are plotted in Figure 4. The curve indicates that the gradient run-time ratio gets close to the theoretically optimal value of as we increase the system size.

## 4Discussion

### 4.1Choice of shearing experiment datasets

Of the six shearing experiment datasets, cf. Figure Figure 1, we have used two for parameter estimation. One of the reasons for this choice is an incompatibility of most of the datasets with assumptions made in the design of the strain energy functional . In particular, the strain energy dictates an ordering of the shear mode stiffnesses in the case of a homogeneous shear displacement. We can see this by adapting the analysis that leads to equations (5.23) – (5.28) of [15]. In this analysis a parameter is introduced to represent the amount of simple shear displacement present in a homogeneous deformation. For example for the FS mode

Using this deformation gradient, and the respective deformation gradients of the other modes, the shear component of the Cauchy stress in the shearing direction can be calculated for each mode. If we consider the same invariants as in , that is , and use the notation , we arrive at the following equations for shear stress as a function of shear displacement

For further details regarding the derivation of these equations we refer the reader to [15]. The simple shear stresses reveal two assumptions built into the design of , namely for homogeneous simple shear deformations

Out of the six datasets, only one is consistent with these orderings, namely the 6th one, which was used here under the label transversely isotropic. In this dataset the stress-strain relationship is typical of a transversely isotropic material with a stiffer fiber direction. In several other cardiac mechanics simulation studies [19], the Holzapfel and Ogden energy functional has been simplified to model transversely isotropic behavior by removing the terms involving the invariants . For such a simplified model one could use the parameter estimates for that we obtained from the Transversely Isotropic dataset.

However, the Holzapfel and Ogden model was originally proposed to model orthotropic behavior. This motivates also targeting a dataset displaying fully orthotropic properties. In particular, dataset 2 in Figure 1 is such and compares well with Figure 6 of [8] and Figure 2 of [15]. By switching the and curves of Dataset 2 we were able to reinterpret this data in a way that is consistent with the interpretation in [15], and the shear stiffness orderings .

### 4.2Discussion of optimal material parameter values

We have obtained two sets of material parameters: one corresponding to an orthotropic case and one corresponding to a transversely isotropic case. We observe that for both sets of material parameters, the parameter essentially vanishes. For the Transversely Isotropic case, both and essentially vanish, which is in excellent agreement with the transversely isotropic stress-strain pattern. Furthermore we note that the magnitude of both and parameters in the best fitting parameter sets presented in Table ? are very small. In light of the shear stress calculations we can see that the and parameters are related to the degree of extra stiffness in the sheet direction over the sheet normal direction. Indeed when we examine the shear data, Figure 3, we can see that the curves are only slightly stiffer than the curves, which explains why the optimal values of and are so small.

Comparing the orthotropic material parameter values to the previously published values in Table ?, we observe that the fit of our material parameters is significantly better, as expected. By using a finite element model we have been able to relax the homogeneous shearing angle assumption and more realistically model the motion of the cubes in the shearing experiment. We note that our material parameters differ from those previously published, and also that there is a significant variability in the parameter values previously reported. Some of this variability is most likely due to the differences in the selection of points during the digitization of [Figure 2 of [15]], which was done in the studies whose material parameter sets we compare in Table ?. By using original data from the shearing experiment, we were able to remove the uncertainty due to digitization in our parameter estimates. Finally we note that even after the SF-SN curves are swapped in Dataset 2 of Figure 1, there are still minor differences when compared to [Figure 7 of [15]] and [Figure 3 of [12]] and [Figure 4 of [28]]. This also explains why our parameter sets differ from those calculated in the previous studies.

### 4.3Computing functional gradients in cardiac mechanics

Figure 4 demonstrates that the computational cost of the adjoint gradient computation is comparable to that of a single iteration of the nonlinear solution algorithm of for larger system sizes. For smaller system sizes, the cost of symbolic computation and the cost of residual and Jacobian assembly contribute significantly yielding higher ratios – as expected. Wang et al.’s 2013 simulations of a human left ventricle in diastole use system sizes of approximately degrees of freedom [28]. Given the trend in Figure 4, we can expect that the adjoint method and solver implemented in this work will continue to be efficient at this scale and beyond.

Comparatively, assuming the use of Newton’s method for the solution of nonlinear systems, the evaluation of a finite difference gradient requires a linear system assembly and solve for each Newton iteration, and one nonlinear solve is required per component of the gradient. Counting the 8 parameters in the Holzapfel-Ogden model , and assuming a typical solution of the Euler-Lagrange equation takes 6 Newton iterations, we can expect the computational cost of finite difference gradient evaluation to be circa 48 times greater than that of the adjoint method.

In the optimization results of Table 1, we observed iteration counts of up to for the optimization of parameters using our gradient-based method. This compares favorably with the circa iterations needed to estimate parameters using a global method in [Figure 5 of [30]].

### 4.4Implications for organ-scale image-based parameter estimation with spatially resolved material parameters

Although we have tested our adjoint-based multi-start optimization method on the 2002 shear data of Dokos et al [8], we believe our methods will provide the biggest advantage in the case of optimizing cardiac model parameters in high spatial resolution at the organ scale to MRI or echocardiographic image data. In this case the high spatial resolution would allow for detailed modelling of regional differences in tissue stiffness, which is for example present in patients with post-infarct fibrosis.

In such an application a model parameter could be represented as a finite element function similarly to the displacement or hydrostatic pressure fields (**u**, p). Doing this would increase the number of components of the gradient by the number of degrees of freedom needed to spatially represent the parameter of interest. Using a finite difference or reduced order Kalman filter approach in this case would require an additional evaluation of the Euler-Lagrange equation for each degree of freedom introduced, whereas the adjoint gradient formula only needs to be calculated once regardless of the number of additional degrees of freedom. In the current study the adjoint gradient is estimated to be

times faster than finite differencing. In the case of a spatially varying model parameter the speedup is potentially a lot more significant.

When fitting material parameters to the Dokos experiment data, we were able to generate good initial guesses for the local optimization by progressively refining the mesh and using the optimal results from the previous coarser refinement level as an initial guess in the successive finer level. It would be more challenging to apply this technique using image based ventricular geometries, due to the problem of accurately representing the geometry with few elements. As an alternative we propose the multi-start approach, which we have shown here to be accurate and viable using the Dokos experiment data.

One issue that would arise in using the multi-start approach with image based geometries would be the choice of the number of multi-start points; using less points is more computationally efficient, while using more is potentially more robust. Possible solutions are the use of optimal stopping criteria [4] or more sophisticated local-global searches [27].

## 5Conclusions

In this work, we have presented a new application of efficient gradient-based optimization methods in the context of estimating cardiac hyperelastic material parameters from experimental data. In particular, we have demonstrated how an adjoint solution can greatly speed up the evaluation of functional gradients. These methods have produced two new sets of material parameter values that yield simulated stress-strain curves that fit closely to orthotropic and transversely isotropic shear data. For future parameter estimation studies using image based geometries and a local search algorithm, multi-start or a similar method should be used in order to avoid suboptimal minima.

### References

**ACM Transactions on Mathematical Software (TOMS)****40**(2), 9 (2014)

Alnæs, M.S., Logg, A., Ølgaard, K.B., Rognes, M.E., Wells, G.N.: Unified form language: A domain-specific language for weak formulations of partial differential equations.**Journal of biomechanical engineering****127**(1), 148–157 (2005)

Augenstein, K.F., Cowan, B.R., LeGrice, I.J., Nielsen, P.M., Young, A.A.: Method and apparatus for soft tissue material parameter estimation using tissue tagged magnetic resonance imaging.**URL http://www.mcs.anl.gov/petsc**

Balay, S., Brown, J., Buschelman, K., Gropp, W., Kaushik, D., Knepley, M., McInnes, L.C., Smith, B., Zhang, H.: Petscc web page (2015).**Mathematical Programming****37**(1), 59–80 (1987)

Boender, C.G.E., Kan, A.R.: Bayesian stopping rules for multistart global optimization methods.**Biomechanics and modeling in mechanobiology****11**(5), 609–630 (2012)

Chabiniok, R., Moireau, P., Lesault, P.F., Rahmouni, A., Deux, J.F., Chapelle, D.: Estimation of tissue contractility from cardiac cine-mri using a biomechanical heart model.**Philosophical transactions of the Royal Society of London. Series A: Mathematical, physical and engineering sciences****359**(1783), 1233–1250 (2001)

Costa, K.D., Holmes, J.W., McCulloch, A.D.: Modelling cardiac mechanical properties in three dimensions.**Biomedical Engineering, IEEE Transactions on****59**(1), 20–24 (2012).

Delingette, H., Billet, F., Wong, K.C.L., Sermesant, M., Rhode, K., Ginks, M., Rinaldi, C.A., Razavi, R., Ayache, N.: Personalization of Cardiac Motion and Contractility From Images Using Variational Data Assimilation. DOI 10.1109/TBME.2011.2160347**American Journal of Physiology-Heart and Circulatory Physiology****283**(6), H2650–H2659 (2002)

Dokos, S., Smaill, B.H., Young, A.A., LeGrice, I.J.: Shear properties of passive ventricular myocardium.**SIAM Journal on Scientific Computing****35**(4), C369–C393 (2013)

Farrell, P.E., Ham, D.A., Funke, S.W., Rognes, M.E.: Automated derivation of the adjoint of high-level transient finite element programs.**URL http://www.malt-meeting.net/2015/abstracts/balaban.pdf**

Finsberg, H.N., Balaban, G., Sundnes, J., Odland, H.H., Rognes, M., Wall, S.T.: Mechanical imaging of dynamic patient stress patterns (2015).**In: Statistical Atlases and Computational Models of the Heart-Imaging and Modelling Challenges, pp. 93–104. Springer (2015)**

Gjerald, S., Hake, J., Pezzuto, S., Sundnes, J., Wall, S.T.: Patient–specific parameter estimation for a transversely isotropic active strain model of left ventricular mechanics.**International Journal for Numerical Methods in Biomedical Engineering****27**(1), 1–12 (2011)

Göktepe, S., Acharya, S., Wong, J., Kuhl, E.: Computational modeling of passive myocardium.**Urbana****51**, 61,801 (1999)

Goldberg, D.E., Voessner, S.: Optimizing global-local search hybrids.**Siam (2003)**

Gunzburger, M.D.: Perspectives in flow control and optimization, vol. 5.**Philosophical transactions. Series A, Mathematical, physical, and engineering sciences****367**(1902), 3445–75 (2009).

Holzapfel, G., Ogden, R.W.: Constitutive modelling of passive myocardium: a structurally based framework for material characterization. DOI 10.1098/rsta.2009.0091.**Finite element methods in flow problems pp. 121–132 (1974)**

Hood, P., Taylor, C.: Navier-stokes equations using mixed interpolation.**URL http://www.scipy.org/.**

Jones, E., Oliphant, T., Peterson, P., et al.: SciPy: Open source scientific tools for Python (2001–).**Tech. Rep. DFVLR-FB 88-28, DLR German Aerospace Center – Institute for Flight Mechanics, Koln, Germany (1988)**

Kraft, D.: A software package for sequential quadratic programming.**Journal of computational physics****244**, 4–21 (2013)

Krishnamurthy, A., Villongco, C.T., Chuang, J., Frank, L.R., Nigam, V., Belezzuoli, E., Stark, P., Krummen, D.E., Narayan, S., Omens, J.H., et al.: Patient-specific models of cardiac biomechanics.**Springer (2011)**

Logg, A., Mardal, K.A., Wells, G.N., et al.: Automated Solution of Differential Equations By the Finite Element Method.**Journal of the mechanical behavior of biomedical materials****20**, 259–71 (2013).

Marchesseau, S., Delingette, H., Sermesant, M., Sorine, M., Rhode, K., Duckett, S.G., Rinaldi, C.a., Razavi, R., Ayache, N.: Preliminary specificity study of the Bestel-Clément-Sorine electromechanical model of the heart using parameter calibration from medical images. DOI 10.1016/j.jmbbm.2012.11.021.**Tech. rep., Sandia National Labs., Albuquerque, NM (US); Sandia National Labs., Livermore, CA (US) (2000)**

Salari, K., Knupp, P.: Code verification by the method of manufactured solutions.**Journal of biomechanical engineering****129**(2), 279–283 (2007)

Schmid, H., Nash, M., Young, A., Röhrle, O., Hunter, P.: A computationally efficient optimization kernel for material parameter estimation procedures.**Biomechanics and modeling in mechanobiology****7**(3), 161–173 (2008)

Schmid, H., O’Callaghan, P., Nash, M., Lin, W., LeGrice, I., Smaill, B., Young, A., Hunter, P.: Myocardial material parameter estimation-a non–homogeneous finite element study from simple shear tests.**Computer methods in biomechanics and biomedical engineering****12**(3), 283–295 (2009)

Schmid, H., Wang, Y., Ashton, J., Ehret, A., Krittian, S., Nash, M., Hunter, P.: Myocardial material parameter estimation-a comparison of invariant based orthotropic constitutive equations.**In: Medical Image Computing and Computer-Assisted Intervention–MICCAI 2009, pp. 257–265. Springer (2009)**

Sundar, H., Davatzikos, C., Biros, G.: Biomechanically-constrained 4d estimation of myocardial motion.**Water resources research****39**(2) (2003)

Tsai, F.T.C., Sun, N.Z., Yeh, W.W.G.: Global-local optimization for parameter structure identification in three-dimensional groundwater modeling.**International journal for numerical methods in biomedical engineering****29**(1), 83–103 (2013)

Wang, H., Gao, H., Luo, X., Berry, C., Griffith, B., Ogden, R., Wang, T.: Structure-based finite strain modelling of the human left ventricle in diastole.**Medical image analysis****13**(5), 773–84 (2009).

Wang, V.Y., Lam, H.I., Ennis, D.B., Cowan, B.R., Young, A.a., Nash, M.P.: Modelling passive diastolic mechanics with quantitative MRI of cardiac structure and function. DOI 10.1016/j.media.2009.07.006.**Journal of the mechanical behavior of biomedical materials****43**, 35–52 (2015)

Wong, K.C., Sermesant, M., Rhode, K., Ginks, M., Rinaldi, C.A., Razavi, R., Delingette, H., Ayache, N.: Velocity-based cardiac contractility personalization from images using derivative-free optimization.**Journal of the mechanical behavior of biomedical materials****4**(7), 1090–1102 (2011)

Xi, J., Lamata, P., Lee, J., Moireau, P., Chapelle, D., Smith, N.: Myocardial transversely isotropic material parameter estimation from in-silico measurements based on a reduced-order unscented kalman filter.