Estimating cardiac active tension from wall motion—An inverse problem of cardiac biomechanics

The contraction of the human heart is a complex process as a consequence of the interaction of internal and external forces. In current clinical routine, the resulting deformation can be imaged during an entire heart beat. However, the active tension development cannot be measured in vivo but may provide valuable diagnostic information. In this work, we present a novel numerical method for solving an inverse problem of cardiac biomechanics—estimating the dynamic active tension field, provided the motion of the myocardial wall is known. This ill‐posed non‐linear problem is solved using second order Tikhonov regularization in space and time. We conducted a sensitivity analysis by varying the fiber orientation in the range of measurement accuracy. To achieve RMSE <20% of the maximal tension, the fiber orientation needs to be provided with an accuracy of 10°. Also, variation was added to the deformation data in the range of segmentation accuracy. Here, imposing temporal regularization led to an eightfold decrease in the error down to 12%. Furthermore, non‐contracting regions representing myocardial infarct scars were introduced in the left ventricle and could be identified accurately in the inverse solution (sensitivity >0.95). The results obtained with non‐matching input data are promising and indicate directions for further improvement of the method. In future, this method will be extended to estimate the active tension field based on motion data from clinical images, which could provide important insights in terms of a new diagnostic tool for the identification and treatment of diseased heart tissue.


| INTRODUCTION
In the last years, computational cardiac modeling has advanced to provide more realistic and personalized approaches to simulate the human heart beat in silico. As a result of a numerical simulation, all quantities included in the model can be observed and evaluated to provide additional insight into cardiac function. In particular, during the simulation of the heart beat, passive and active forces interplay, which leads to the deformation of the myocardial wall and thus, the pumping of the heart. Despite progress in the field of cardiovascular imaging in the last years, there is no clinical modality which measures the passive and active forces. Still, the resulting deformation can be observed in various imaging techniques: echocardiography, computed tomography or magnetic resonance imaging (MRI). This leads to an inverse problem of cardiac biomechanics: an estimation of the active tension field when the deformation of the heart muscle is given. This corresponds to a general definition of an inverse problem, which is the calculation of the causal factors that produced the input set of observations.
In the last decade, solving such a problem has been the objective of some research groups. Different measures of deformation of the heart have been used as an input data. In particular, Finsberg et al. 1 and Balaban et al. 2 described an estimation method based on strain measurements in the AHA segments obtained by echocardiographic speckle tracking. Dabiri et al., 3 Asner et al. 4 and Hu et al. 5 used data from tagged MRI. Such input data contains information about the twist and rotation of the heart muscle contrary to displacement data of the endocardial surface, which we aim at using as an input for our method. Similarly, Otani et al. 6 used a volumetric displacement field, but they estimated the active stress on a cubically shaped tissue geometry. We use a truncated ellipsoid and a left ventricular geometry derived from MRI data. Different levels of spatial resolution of the active tension have been estimated. Asner et al. 4 and Hu et al. 5 estimated one tension value for the entire ventricle per time step. Dabiri et al. 3 estimated the maximal active tension for every ventricle segment (out of 16 segments). We aim at reconstructing a time dependent tension field with spatial resolution equal to the number of finite elements in the geometrical domain. To explore the accuracy of our method, we conducted a sensitivity study with synthetic data (generated in a forward simulation), where all the parameters are known. Asner et al. 4 and Otani et al. 6 also investigated the accuracy of their methods by adding noise with distribution of zero mean to the displacement data and to the fiber distribution in the tissue. In our study, we decided to add non-mean-free noise (constant shift) also to the displacement data and the fiber distribution to obtain more rigorous settings to explore the accuracy of our method.
The aim of this work is to present a numerical method that reconstructs the dynamic active tension field when the endocardial deformation is known. This is an ill-posed problem and therefore, spatial and temporal Tikhonov regularization is implemented. In this proof-of-concept work, the deformation information is extracted from in silico data obtained using the forward model. In contrast to data extracted from 3D tagged MRI data and echocardiographic speckle tracking, in the synthetic data used in the presented study, no information about the twist and the rotation of the tissue is available. In future, the deformation could be derived from any imaging technique which provides the wall motion (echocardiography, computed tomography or MRI). In particular, we are aiming at clinical cine MRI acquisitions. The rule-based assumption of cardiac fiber orientation could in future be personalized by diffusion tensor MRI 7 or quantified by echocardiography. 8,9 Since in case of using clinical data, there is no ground truth available to validate the method, we decided to first determine the accuracy of the reconstruction of active tension and its sensitivity to variation of model parameters using synthetic data.
We study the influence of the fiber orientation on the reconstruction because changes in fiber orientation affect the distribution of the active tension in the heart. 10 We approximate the fiber direction for the numerical simulations with a rule based method developed by Bayer et al. 11 We also explore the effect of variability in the endocardial deformation, considering the potential inter-and intra-observer variability in MRI segmentation. 12 We provide a correlation between the error in the reconstruction of the active tension field and the error in the input data. Finally, we assume that the myocardial regions which do not develop active tension are related to infarct scars and therefore, by reconstructing the dynamic active tension distribution, we can detect the size and the location of infarct regions. Additionally, the active stress distribution and dynamics could be complementary information. In clinical routine, scar tissue can be identified by abnormal regional wall motion using echocardiography or by injecting gadolinium to the patient in a late gadolinium enhancement (LGE) MRI study. Both ways, passive material properties and active tension in and around the scar cannot be depicted. Also, it has been shown that LGE has negative side effects and should be approached with caution. 13 Additionally, areas not showing hyper-enhancement might still have undergone remodeling and develop more or less or dynamically different active stress. In the sensitivity analysis, we defined different sizes of infarct areas.
In this work, first, the inverse method is outlined in its mathematical framework. Second, we provide the results of the sensitivity analysis based on synthetic input data. Third, we provide a study on infarct size detection and localization. Furthermore, different regularization parameters are investigated.

| Method for solving the inverse problem to estimate active tension
While solving the forward problem of the cardiac biomechanics delivers the deformation of the myocardium, the deformation of the endocardial surface (target surface) over time is input data for the inverse problem to estimate active tension. The aim of solving this inverse problem is to reconstruct the active tension distribution over time which leads to the provided target deformation. In particular, we aim to estimate for each volumetric finite element e l ℰ (l = 1, …, M, M = j ℰj, where ℰ is the set of volume elements) for each time step t i (i = 1, …, T, where T = T CL /Δt, Δt = 10 ms and T CL = 800 ms is the heart cycle duration) a distinct value for the active tension s l i . Since the dimension of the input data is lower than the one of the output data (for each time step), the problem is underdetermined. To overcome this difficulty, we developed a model-based approach that combines models describing heart-specific phenomena and the finite element method. The numerical approach is described in the following.
For the inverse method, the distance between the measured endocardial surface (target surface) and the endocardial surface of the inverse model, called the reconstructed surface, is calculated. The distance vector between these two surfaces is denoted by g ℝ N . For every triangle on the reconstructed surface, each node is matched with the closest triangle on the target surface. The distance in x-, y-and z-direction is calculated between the reconstructed node and a point p T k ℝ 3 on the target triangle. This point p T k is the intersection point of the reconstructed triangle normal (starting at the reconstructed node) and the matched target triangle. Then, the component-wise distance g i from the j-th reconstructed node to the target surface is obtained as the mean over the distance (belonging to the j-th node) of the set of reconstructed triangles including this node (K j ): where where r = 0, 1 or 2 for the x-, y-and z-component, respectively and j is the index of the reconstructed node observed (j = 1,…,N ),N is the number of nodes on the reconstructed surface, N = 3N and ι R k T R is a triangle on the reconstructed surface with nodes n R j . The index k = 1,…, j T R j, where T R is the set of reconstructed triangles.
The target surfaces are available every 20 ms, but the distance is calculated for every simulation time step t i . Therefore, the coordinates of the nodes of the target triangles are interpolated linearly over time.
In order to obtain the change of the active tension contributing to the overall tension, the distance g i is minimized iteratively using a Newton scheme for every time step t i . The change of the active tension in the n-th iteration is denoted by τ i,n ℝ M and the gap by g i,n (the index i corresponds to the i-th time step). The overall active tension vector s i at the i-th time step adds up to s i = s i − 1 + τ i , where τ i = P n k = 1 τ i,k . Here, the number of Newton iterations is denoted by n and determined by the following stop criteria: where both the absolute tolerance ϵ a and the relative tolerance ϵ r are set to 10 −5 mm and the maximal number of iteration N max is 5. In every Newton step, the gap g i,n is recalculated to regard the deformation of the reconstructed surface, caused by the estimated active tension of the previous iteration. To calculate the deformation, the forward problem is solved for one time step to assert the balance of all forces (equilibrium equation). In particular, the active tension is balanced with the passive force and the pressure applied on the surface. After the forward solver is executed, the distance between target and reconstructed triangles is calculated once again.
In every Newton iteration step, the change of active tension is estimated. For this purpose, the relation of the gap g i, n and τ i,n is linearly approximated by the matrix A as follows: The matrix A ℝ N × M is the product of the inverse tangential stiffness matrix K − 1 d ℝ N × S and a matrix K T ℝ S × M , where S is the number of nodes of all volume elements. In general, the stiffness matrix describes the linearized relation of a small change of the nodal displacement Δd and the resulting change of the nodal forces Δf. In order to reduce the dimension of the stiffness matrix, only the rows related to the reconstructed nodes are used for the matrix K d . In particular, the reconstructed nodes are displaced in positive and negative direction and the resulting nodal forces in the entire volume are computed. A further dimension reduction is possible by removing some reconstructed nodes prior to the simulation. The latter reduction is only necessary for large number of endocardial nodes, which was not the case in the presented study. The second matrix K T describes the change of all nodal forces Δf due to the change of the active tension Δs in all volume elements. Both matrices K d and K T are calculated using central finite differences with ϵ d = 10 −12 m and ϵ t = 1 kPa, respectively. Altogether, the matrix A is the following: Due to the construction of the matrix A, the active tension τ i,n in each volume element is mapped to the displacement of the endocardial surface. Computing τ i,n from Equation (3) is an ill-posed problem since there is no unique solution, as the matrix A is not symmetric (N < M). Furthermore, this matrix has a high condition number (10 9 ). Thus, it requires regularization and we employed second order Tikhonov approach with spatial and temporal regularization. Therefore, Equation (3) was extended by a regularization term yielding the following minimization problem: argmin Aτ i,n −g i,n 2 + λ 1 1 The solution is obtained by solving the regularized normal equation: Here, I is the identity matrix and Δ c is the discrete Laplace operator, which contains the connectivity information of the volume elements. The parameter λ 1 controls the temporal regularization and λ 2 the spatial one. From Equation (5), each regularization parameter can be expressed as a ratio between the squared residual norm and the corresponding squared norm of the regularization term. The residual norm depends on the element discretization of the computational domain, as the volume of each element is involved in the construction of the matrix A. Thus, the choice of regularization parameters depends on the spatial discretization of the computational domain and is provided in Chapter 3. Equation (6) is solved in every iteration using LU decomposition to provide the change of the active tension τ i,n of the current iteration.
In Figure 1, the method for solving the inverse problem to estimate active tension is graphically represented.

| Geometry for the forward solver
For this proof-of-concept study, we used two different geometries. The first one was a hollow, truncated ellipsoid representing an idealized left ventricle as described in Land et al. 14 Gmsh 15 was used to obtain M = 2074 volumetric elements and 4237 nodes. The endocardial surface was extracted, resulting in j T R j = 544 triangles andN = 285 nodes. The long axis of the ellipsoid measured 25 mm and the diameter of the base plane was 20 mm. The mean edge length of the surface elements was 1.9 mm and the wall thickness was 3 mm. The ellipsoid clipped along its long axis is visualized in Figure 2B. The endocardial surface nodes were a subset of the volumetric nodes. This surface was considered as the reconstructed surface. All nodes on the open end of the truncated ellipsoid were fixed in all three directions as a boundary condition for the forward and inverse model. The fixed nodes are marked by red color in Figure 2B. The fiber direction was determined by the algorithm of Bayer et al. 11 In the reference case, the fiber orientation was transmurally changing from 90 on the endocardium to −90 on the epicardium. 16 In the following, this combination will be denoted by (90 , −90 ). The transmural spatial resolution of the ellipsoid was one element in most parts of the wall and two F I G U R E 1 A graphical representation of the method for solving the inverse problem to estimate active tension F I G U R E 2 A histogram of the fiber orientation assigned to the quadrature points in the finite elements of (A) the ellipsoid geometry and (C) the image-based geometry. The geometry clipped along its long axis with the fixed nodes at the base are marked by the red color for (B) the ellipsoid geometry and (D) the image-based geometry elements at the base. Every element had four quadrature points (QPs), distributed inside the element. The fiber orientation was assigned to the finite element at the QPs and the range (90 , −90 ) was defined at the surfaces. Thus, the range of the fiber orientation at the QPs is smaller than the range of the fiber orientation at the surfaces. Figure 2A shows the frequency of the assigned fiber orientation at the QPs. The most frequently assigned fiber orientations were about 65 and −65 .
The second geometry in this study was an image-based left ventricular geometry. The images were taken from a healthy volunteer at the University Hospital in Heidelberg with a 1.5 Tesla MR tomography (Ingenia CX, Philips Medical Systems). The voxel spacing of the images was 0.62 × 0.62 × 0.8 mm. The trigger time was 560 ms at a heart rate of 77 bpm. The volunteer gave informed consent and the study was approved by the institutional review board. The endocardial surface of the left ventricle was automatically segmented with the Philips IntelliSpace Portal. The epicardial surface was modeled as a smoothed and inflated version of the endocardial surface. The endocardial surface was discretized in j T R j = 432 triangles andN = 218 nodes. Gmsh 15 was used to obtain M = 1218 volumetric elements and 2475 nodes. The long axis of the left ventricle measured 100 mm and the maximal short axis length was 60 mm. The mean edge length of the surface elements was 4.3 mm and the wall thickness was 6 mm. The volume of the blood cavity was 133 mL. All nodes on the base were fixed in all three directions as a boundary condition for the forward and inverse model. In Figure 2D, the image-based geometry in end-diastolic state clipped along its long axis is visualized and the fixed nodes are marked by red color. The fiber direction was determined by the algorithm of Bayer et al. 11 with angles 80 on the endocardium to −70 on the epicardium as suggested by Wong et al. 17 Figure 2C shows a histogram of the fiber distribution. The most frequently assigned fiber orientations were about 60 and −50 .

| Geometry for the inverse solver
For the described inverse problem, the same geometry as for the forward problem was used. Furthermore, an additional surface was needed to prescribe the motion over time of the endocardial boundary, the target surface. It was defined by the same mesh topology as the reconstructed surface. However, this is not a prerequisite, a different topology can also be used. The node coordinates of the target surface were provided every 20 ms, because in the long term, we aim to apply the inverse method on clinical data (cine MRI), which are also available about every 20 ms. In this work, we used only synthetic data as this is the only case where the ground truth active tension distribution is known. Thus, the target surface and its deformation was extracted from the endocardial surface of the forward simulation.

| Forward model
The aim of the forward calculation is to compute the deformation of the geometry, which is caused by the interplay of the forces applied in the numerical model. The underlying physics of the cardiac biomechanics are described by the governing equation for the balance of linear momentum. If mass inertia is neglected, the governing equation reduces to the equilibrium equation ensuring that the external nodal forces match the internal nodal forces. 18 To calculate the deformation, the finite element method 18 was implemented as described in Fritz et al. 19 The computational domain was spatially discretized by tetrahedral elements with quadratic shape functions and triangles with linear shape functions on the surface. 19 The resulting non-linear algebraic equations are solved in a numerical framework called CardioMechanics, which was verified in the benchmark by Land et al. 14 In this framework, PETSc 20 software libraries were employed.
In this study, we solved the equilibrium equation. The internal forces are the sum of the passive and active forces and the external force is arising from the pressure applied to the endocardial surface. The forces are represented by models closer described in the following.

| Passive material model
The passive forces arise from a constitutive law, which relates the stress and the strain of the heart. We used the model of Guccione et al. 21 describing a hyperelastic, transversely isotropic material. Furthermore, incompressibility was imposed by adding a penalty term to the strain energy function: where C, b f , b t and b ft are the parameters of the Guccione model, E ij (i, j [1,2,3]) are elements of the Green strain tensor, det(F) is the determinant of the deformation tensor and K = 10 6 scales the penalty term.
To identify the passive material parameters, we applied a simple and robust method based on the work of Klotz et al. 22 A similar approach to determine passive parameters was also applied by Dabiri et al. 3 and Genet et al. 23 Klotz et al. established a relation between pressure and normalized volume of the ventricle based on measured data. Furthermore, they provided 0.42 as the ratio of the unloaded ventricle volume (at pressure zero) and the volume at a pressure of 30 mmHg. To estimate the passive parameters, we minimized the sum of the deviation of a simulated pressure volume curve from the Klotz curve and the error of the volume ratio. This resulted in the following parameters:

| Active tension model
The active tension model describes the development of the force, which leads to a shortening of the cardiac muscle cells and thus to the contraction of the heart. It was modeled to act purely in fiber direction. The course of the active tension curve is prescribed to be the same for each volume element without temporal offset. We implemented the Hill curve, which represents the elastance of the ventricle as it relates the ventricular pressure to its volume as described by Stergiopulos et al. 24 The maximal active tension was set to 40 kPa for all simulations with ellipsoid geometry. For the image-based geometry, the maximal active tension was increased to 180 kPa in order to counterbalance the force arising from the pressure, applied on the endocardial surface. The course of the active tension in the forward model, which is the ground truth for the inverse solution, is visualized by the green dotted line in each plot of active tension versus time (e.g., Figure 4). The value obtained from the active tension model s i is added to the first diagonal element of the second Piola-Kirchhoff stress tensor S for every time step (active stress approach): where W is the energy function of the passive material model and, E is the Green strain tensor obtained from the deformation tensor.

| Pressure model
The pressure model describes the pressure (force per area) acting on the endocardial surface. For each triangle k with vertices a k , b k and c k , the contribution f k to the nodal force is obtained by: where P is the cavity pressure. For this study, we prescribed the pressure defined by analytical functions. For the image-based left ventricular geometry, we applied both the positive part of a sine function and a Hill curve, 24 scaled between 0 and 120 mmHg. Both functions are visualized in Figure 3.

| Sensitivity analysis
To explore the sensitivity of our novel inverse method regarding uncertainty of input parameters, we systematically varied them and evaluated the reconstructed active tension. In particular, we varied the fiber orientation and displacement of the target surfaces. We also defined infarct areas of different size. Except stated otherwise, we applied the parameters of the reference case, described in the following.

| Reference case
Each case consists of two simulations, one forward and one inverse. The forward simulation is executed with the ground truth of the active tension and the resulting endocardial surface deformation is extracted. For both simulations, respectively, a set of parameters is defined: related to the geometry (as described in Section 2.2)-the fiber direction, the number of nodes and elements; related to the solution method-time step size, tolerances, length of one heart beat (cycle length); related to the different models-passive parameters, active tension curve and its maximum. In the reference case, the parameters of the forward and inverse simulations were the same, except for sure that in the inverse simulation no active tension model was provided. The active tension curve was reconstructed instead. In particular, for the ellipsoid geometry, the reference fiber direction was (90 , −90 ) and for the image-based left ventricular geometry (80 , −70 ) as described in Section 2.2. A sine and a Hill pressure curves were applied to the endocardial surface of the image-based left ventricular geometry (respectively, Case 11 and Case 12 in Table 1). The cycle length was T CL = 800 ms and the time step size Δt = 10 ms. The target surfaces were extracted from the endocardial deformation of the forward simulation every 20 ms. The passive material parameters were the same in the forward and inverse model. Also, they were the same for the healthy tissue and the potential infarct area.

| Fiber variation
The fiber direction was varied only in the forward simulation, which led to a different target surface deformation compared to the reference case. Three additional fiber combinations were investigated: (80 , −90 ), (90 , −80 ) and (80 , −80 ). In Table 1, they correspond to Case 1, Case 2 and Case 3, respectively. For the model used in the inverse solver, (90 , −90 ) was used. In these cases, an infarct area was not present in the forward simulation. The fiber variation study was conducted only for the ellipsoid geometry.

| Target surface variation
The endocardial surfaces were extracted from the forward simulation and the position of a target node was varied in the direction of the node normal. The node normal was calculated as the mean of the normals of the adjacent triangles, which all point towards the barycenter of the ellipsoid. Only the nodes which were not fixed due to the boundary conditions were varied.
Here we defined nine experiments in which the surface node coordinates were varied with different percentage of wall thickness. The first three experiments were executed with the ellipsoid geometry: Case 4, Case 5 and Case 6 in Table 1. In Case 4, 20% of wall thickness were added to the node coordinates in positive normal direction, that is, Pressure curves applied to the endocardial surface: (A) Sine curve, (B) Hill curve towards the center. In Case 5, 20% of wall thickness were added to the node coordinates in negative normal direction. In Case 6, also 20% of wall thickness were added to the node coordinates, but the normal direction was alternating between positive and negative each 20 ms. The other six experiments were executed with the image-based geometry: Cases 13-18 in Table 1. In Case 13, Case 15 and Case 17, respectively, 10%, 20% and 30% of wall thickness were added to the node coordinates in positive normal direction, i.e., towards the center. In Case 14, Case 16 and Case 18, respectively, 10%, 20% and 30% of wall thickness were added to the node coordinates in negative normal direction. The variation was applied to the target surface over all time steps except the target surface at 0 ms. We assume there is no active tension present at the initial state of the simulation and therefore, the distance vector is the zero vector as the target surface and the reconstructed surface match at 0 ms. In these cases, an infarct area was not present in the forward simulation.

| Infarct size variation
We modeled simplified infarct areas as regions which developed no active tension in the ellipsoid geometry to see whether the proposed inverse method, which reconstructs the active tension, can detect these infarct areas. In order to explore the limits of the reconstructed infarct size, connected transmural infarct areas of different size were defined in forward simulations: 1%, 5%, 10% and 22% of the number of volume elements. In those elements, the active tension was constantly set to zero in the forward simulation. The passive material parameters were the same in the healthy and the infarct cells. The presence of an infarct area led to a different target surface deformation compared to the reference case.
In Table 1, the described variations correspond to Case 7, Case 8, Case 9 and Case 10, respectively.

| Evaluation of the inverse reconstruction
We defined four criteria to quantify the quality of the reconstructed active tension in the inverse simulation. The peak tension time (PTT) t l p (l = 1…M) was defined as the time step in which the maximum (over time) active tension of the l-th finite element is reached in the inverse reconstruction.
The peak tension error ε l a of the l-th finite element is the difference between the reconstructed active tension at PTT s l t l p and the ground truth active tension s l GT t l p , also at PTT: The mean tension error (MTE) of the inverse reconstruction is then defined as the mean over the absolute value of the peak tension error of all volume elements: The relative tension error ε l r of the l-th finite element is the peak tension error divided by the active tension at PTT s l GT t l p . The mean relative tension error (MRTE) ϵ r of the inverse reconstruction is the mean of relative tension errors of all elements: To capture the reconstruction error during the entire time course, we calculated the root mean square error (RMSE) for every element ε l m . It is the second-norm of the differences of the reconstructed active tension s l (t i ) and the ground truth active tension s l GT (t i ) divided by the root of the number of time steps T: To separate infarct tissue from healthy tissue, a relative threshold  for the reconstructed active tension was defined.
Each element e l ℰ was classified as an infarct element in ℰ i or as a healthy element in ℰ h . An element e l was considered to be infarcted if the reconstructed peak tension was less than the percentage threshold of the maximum over all reconstructed active tension for all elements: All other elements were considered healthy: The cardinalities of intersections of these four sets were used to calculate the sensitivity and specificity for each experiment. Since the elements in the sets are the IDs of the discrete finite elements of the geometry, the metrics inherently consider the locations of the healthy and unhealthy elements. The elements satisfying the true positive condition are the ones defined as infarct elements in the forward simulation as well as reconstructed as infarct elements in the inverse result. Thus, the sensitivity is calculated only for the experiments with defined infarct areas in the forward settings. Based on those cases, an overall receiver operating characteristic (ROC) curve is calculated, in order to determine one optimal threshold  for all cases with different infarct sizes. To determine the overall ROC curve, for each threshold value the mean over the corresponding values of all ROC curves (for each infarct case separately determined) is calculated.
For the cases of infarct size variation, MTE was calculated for both sets ℰ corr i and ℰ corr h . MRTE was calculated only for healthy elements since the active tension for infarct elements was set to be zero in the forward simulation. The infarct classification was applied even when no infarct area was defined in the forward simulation.

| Reference case
In the reference case, the reconstructed active tension was matching the ground truth active tension with an MRTE of 0.7% and a maximal peak tension error <0.5 kPa ( Figure 4A). All elements were correctly identified as non-scar elements. This result was obtained with spatial regularization λ 2 = 10 −14 . Applying stronger spatial regularization (e.g., λ 2 = 10 −12 and λ 2 = 10 −10 ) slightly improved the result as the MRTE decreased to 0.6%. Less spatial regularization led to higher MRTE, for λ 2 = 10 −16 it was 2%. The peak tension errors for all elements for different spatial regularization  Figure 5. The following results were obtained with spatial regularization λ 2 = 10 −14 , except stated otherwise.

| Fiber variation
Next, the fiber direction in the forward model was varied. In the inverse model, the fibers of the reference case were used. For both cases of the endocardial or epicardial angles varied by 10 (Case1, (80 , −90 ) or Case 2, (90 , −80 )) in the forward simulation, the reconstructed active tension had an MRTE of about 26%. Furthermore, a large spread of the active tension at the PTT was observed ( Figure 4B,C), ranging from 15 to 55 kPa, whereas the maximum of the forward tension was 40 kPa for all elements. In both cases, the tension in the apex region was underestimated and led to an incorrect detection of an infarct region (specificity 0.92 for threshold  = 0:52 ). Nevertheless, the shape of the time course of the reconstructed active tension curve resembled the shape of the Hill curve used in the forward model as visualized in Figure 4B,C. The variation of both endocardial and epicardial angles by 10 (Case 3, (80 , −80 )) in the forward simulation led to an MRTE of 68% for the reconstructed active tension. The spread of the tension ranged from 0.7 to 76 kPa ( Figure 4D) and thus, led to a misclassification of elements with specificity 0.85.

| Target surface variation
Here, the endocardial surfaces extracted from the forward simulation were varied and provided as input target surfaces for the inverse solver. In Case 4 in Table 1, the nodes were shifted towards the inside of the ellipsoid by 20% of the initial wall thickness. This led to a severe overestimation of the active tension, ranging from 66 to 140 kPa at the PTT, as visualized in Figure 6A. Thus, a misclassification of some elements was observed with specificity of 0.95.
In Case 5 in Table 1, the nodes were shifted towards the outside of the ellipsoid by 20% of the initial wall thickness. In contrast to the case described before, an underestimation of the active tension was observed ( Figure 6B). The reconstructed tension ranged from 3.5 to 30 kPa, also leading to a misclassification of elements (specificity 0.42).
These two experiments explored the maximal variation of the target surface regarding the potential uncertainty in segmentation of clinical data. In Case 6 in Table 1, variation of the target surface was added in alternating directions, changing between positive and negative displacement every 20 ms. The reconstruction of the active tension had an MRTE of 95% and values at the PTT ranging from 46 to 86 kPa ( Figure 6C). All elements were classified correctly (specificity 1.0).
In a further experiment for Case 6, additional to the spatial regularization, temporal regularization was imposed. For temporal regularization λ 1 = 10 −20 , the MRTE decreased to 12% and the active tension values at the PTT were between 40 and 48 kPa ( Figure 6D). Less temporal regularization λ 1 < 10 −20 did not improve the reconstruction compared to the result with no temporal regularization. More temporal regularization λ 1 = 10 −18 led to an MRTE over 200% F I G U R E 5 The peak tension error of all volume elements and their distribution for different values of the spatial regularization parameter λ 2 . The black dot represents the median of the peak tension error and the gray bar represents the error between the first and third quartile as the PTT for all elements was shifted to 600 ms with a range for the reconstructed active tension between 62 and 100 kPa.

| Infarct size variation
In order to estimate the infarct size based on the inverse solution, all finite volume elements were classified as healthy or infarct elements by a relative threshold applied on the reconstructed active tension. The optimal threshold value was  = 0:52 as determined by the overall receiver operating characteristic (ROC) curve (Figure 7). The area under the curve (AUC) was 0.968. The calculated threshold value is close to the optimal threshold for the case with 1% infarct area, which is 0.53 (AUC 0.999) and further away from the optimal threshold for the case with 22% infarct area, which is 0.28 (AUC 0.975).
For an infarct with a size of 1% (Case 7 in Table 1) of the whole left ventricular wall volume in the forward simulation ( Figure 9A), the sensitivity was 0.95 and the specificity was 0.99. For cases of infarct sizes of both 5% and 10% (Case 8 and 9 in Table 1, respectively), the sensitivity was 1.0 and the specificity was 0.89. For the infarct size of 22% (Case 10 in Table 1), the sensitivity was 0.98 and the specificity was 0.71. The threshold value 0.52 was applied for all experiments in this work, in order to calculate the specificity of the classification.
In the forward simulation, the elements of the infarct area experienced some deformation due to the surrounding healthy elements. For all infarct cases, the mean reconstructed active tension of the infarct elements varied up to 22 kPa ( Figure 8A-D), which is an overestimation compared to the zero tension in the forward simulation. On the other hand, as the inactive (only passive) behavior of the infarct elements restricts the overall contraction, the mean reconstructed active tension of the healthy elements was slightly underestimated and varied between 34 and 40 kPa ( Figure 8A-D). Additionally, this effect was amplified by the spatial regularization (λ 2 = 10 −14 ), which led to blurring of the active tension distribution across the elements. In particular, in the border zone of the detected infarct, the active tension gradually rose from the inside of the infarct area in direction of the surrounding tissue ( Figure 9A-D).
In general, the reconstruction of the infarct size was improved by applying less spatial regularization (λ 2 < 10 −14 ). For example, for the 5% infarct case, the specificity rose to 0.96 and the MRTE decreased from 14% to 11% for F I G U R E 7 The overall ROC curve is the multi-color line. The single-color curves are the ROC curves calculated for each experiment with an infarct scar of specific size: 1% (blue line), 5% (orange line), 10% (yellow line) and 22% (purple line) Reconstruction of active tension for infarct size 1%, Case 7; (B) for 5%, Case 8; (C) for 10%, Case 9; (D) for 22%, Case 10. The reconstructed mean active tension time course of healthy and infarct elements is visualized by the black and gray line, respectively. The time-varying distribution percentiles (from 5% to 95%) are shown as shaded bands (5%) around the median (not shown) of the reconstructed active tension time course for all volume elements. For every time step, each band represents the active tension range in 5% of the finite elements. The blue bands represent the healthy elements and the red ones the infarct elements (as defined in the forward simulation). The green dotted line is the ground truth active tension of the healthy elements λ 2 = 10 −16 . In case more spatial regularization was applied, the peak tension error of the reconstructed active tension for elements in the scar region rose compared to the case with λ 2 = 10 −14 . The distribution of the peak tension error for all infarcted elements for different values of the spatial regularization parameter is visualized in Figure 10.

| Reference case
The reference case of the image-based left ventricular geometry did not consider pressure in the forward and the inverse simulation. The reconstructed active tension was matching the ground truth active tension with an MRTE of 0.1% and a maximal peak tension error <0.05 kPa for λ 2 = 10 −12 . Applying stronger spatial regularization (e.g., λ 2 = 10 −10 and λ 2 = 10 −8 ) slightly improved the result as the MRTE decreased to 0.06%. Less spatial regularization led to higher MRTE: 3.8% for λ 2 = 10 −16 . The peak tension errors for all elements for different spatial regularization parameters λ 2 = 10 −16 , …, 10 −6 are visualized in Figure 11. The results with applied pressure and varied target surfaces were obtained with spatial regularization λ 2 = 10 −12 . This was the smallest value of λ 2 to lead to a small MRTE (<1%) for the reference case. No temporal regularization was imposed.
F I G U R E 9 The reconstructed active tension at 400 ms is color coded on the geometry for the different infarct sizes: (A) for 1%, Case 7; (B) for 5%, Case 8; (C) for 10%, Case 9; (D) for 22%, Case 10. The infarct area as defined in the forward simulation is indicated by green wireframe

| Pressure application
For Cases 11-18, we applied pressure on the endocardial surface in the forward and inverse simulation. As described in Chapter 2.3, two pressure curves were considered: a sine curve (Case 11 in Table 1) and a Hill curve (Case 12 in Table 1). In both cases, the reconstructed active tension was matching the ground truth active tension with an MRTE of 0.2% and a maximal peak tension error <0.05 kPa. Figure 12 visualizes the time course of the corresponding volume change. The volume change of the ventricular cavity with the Hill curve was only 2 mL between 320 and 420 ms. At 400 ms, which lies inside this quasi-isovolumetric time interval, the maximal tension is applied in the forward simulation and correctly reconstructed in the inverse reconstruction. For Case 11 and Case 12, the Newton solver reached the maximum number of iterations ( N max = 5) for 77% of the time steps of the whole simulation (800 ms) and for more than 90% of time steps for which the ground truth active tension was 0 kPa. In both cases, the median of the mean gap was 0.004 mm (gray horizontal line in Figure 14), the maximal value of the mean gap was 0.03 mm and the minimum was 0.0005 mm. The mean gap is the mean of the distances between the nodes on the target and the nodes on the reconstructed surface. The median of the mean gap is the median over all time steps for which the Newton solver reached the maximum number of iterations.

| Target surface variation
In the image-based left ventricular geometry, endocardial surfaces extracted from the forward simulation were systematically altered to study the influence of imaging and segmentation uncertainty. In both forward and inverse simulations, the Hill pressure curve was applied. In Case 13, the nodes were shifted towards the inside of the ventricle by 10% of the initial wall thickness. This led to an overestimation of the active tension, ranging from 209 to 211 kPa at the PTT, F I G U R E 1 0 The peak tension error and its distribution for all infarcted elements e l ℰ corr i of the ellipsoid geometry for different values of the spatial regularization parameter λ 2 . The black dot represents the median of the peak tension error and the gray bar represents the error between the first and third quartile F I G U R E 1 1 The peak tension error of all volume elements of the image-based left ventricle and their distribution for different values of the spatial regularization parameter λ 2 . The black dot represents the median of the peak tension error and the gray bar represents the error between the first and third quartile as visualized in Figure 13A. The median of the mean distance between the target and the reconstructed surface was 0.01 mm as visualized in Figure 14. For Case 15 and Case 17, the nodes were shifted by 20% and 30%, respectively, towards the inside of the ventricle. This led to an overestimation with a mean value of 251 kPa at the PTT for Case 15 and 306 kPa at the PTT for Case 17, as visualized in Figure 13C,E. For both cases, the median of the mean gap was above 0.02 mm (Figure 14).  (Table 1), the nodes were shifted towards the outside of the ellipsoid by 10%. The mean of the active tension across all elements (156 kPa) was slightly underestimated compared to the ground truth, visualized by the green dotted line in Figure 13B. The median of the mean gap was 0.007 mm, which is the closest to the one of Case 12 (gray horizontal line in Figure 14). In Case 16 and Case 18, the nodes were shifted towards the outside of the ellipsoid by 20% and 30%, respectively. Here, an underestimation of the active tension was observed ( Figure 13D,F), ranging from 138 to 142 kPa for Case 16 and from 126 to 131 kPa for Case 18. For Case 16, the median of the mean distance between the target and the reconstructed surface was above 0.014 mm and for Case 18 -0.026 mm ( Figure 14). For all cases of target surface variation with the image-based left ventricular geometry, all elements were classified correctly (specificity 1.0). Figure 14 shows the relation between the mean relative tension error (MRTE) and the median of the mean gap for all cases with varied target surface. The median was calculated only for the time steps when the Newton solver met the maximal number of iterations ( N max = 5). In these cases, the Newton solver could not provide an active tension field which leads to distances between the target and the reconstructed surface below the absolute tolerance ϵ a = 10 −5 mm. The gray line in Figure 14 is the median of the mean gap for Case 12 and it is below the median of the mean gap of all cases with varied target surfaces. As the MRTE rises, the median of the mean gap also increases.

| Overview of the results of the sensitivity analysis
For the ellipsoid geometry, the root mean square error (RMSE) for all elements and its distribution is visualized in Figure 15 for the reference case and Cases 1-10 as described in Table 1. For these experiments, spatial regularization with λ 2 = 10 −14 was applied. The errors correspond to the courses of the reconstructed active tension for all visualized cases. For example, the highest RMSE was calculated for Case 4, where the reconstructed tension was overestimated ( Figure 6A). The lowest RMSE (median of 1.2 kPa) was obtained for the reference case.
F I G U R E 1 4 Relation between the mean relative tension error (MRTE) and the median of the mean gap distance for Cases 13-18 in Table 1. The gray line is the median of the mean distance for Case 12 (no variation of the target surface) F I G U R E 1 5 The root mean square error and its distribution for all elements of the ellipsoid geometry for the reference case and Cases 1-10 as described in Table 1. The black dot represents the median of the peak tension error and the gray bar represents the error between the first and third quartile For the image-based left ventricular geometry, the RMSE for all elements and its distribution is visualized in Figure 16 for Cases 11-18 as described in Table 1. For those experiments, spatial regularization with λ 2 = 10 −12 was applied. Here, again, the errors correspond to the courses of the reconstructed active tension for all visualized cases. The highest RMSE was calculated for Case 17, where the reconstructed tension was overestimated ( Figure 13E). All cases with target surface variation towards the inside of the ventricle (Case 13, 15 and 17) led to higher RMSE compared to the cases (Case 14, 15 and 18, respectively) with the same amount of target surface variation, but pointing in the opposite direction (towards the outside of the ventricle). Low RMSE were obtained for Case 11 and 12 (median of 3.7 and 4.9 kPa, respectively), where no variation was included. Furthermore, the RMSE obtained for the cases of surface variations (+20% and −20%) for both geometries (Case 4, 15 and Case 5, 16, respectively) were in a similar range (around 40 and 20 kPa, respectively), despite the difference in the maximal tension value of the ground truth (40 kPa for the ellipsoid and 180 kPa for the image-based geometry).

| DISCUSSION
In this work, we presented a novel numerical method for solving an inverse problem of cardiac biomechanics-a reconstruction of the dynamic active tension field. The algorithm delivers the estimated active tension for every time step, resolved on the scale of the individual finite element. In the work of Balaban et al. 2 and Finsberg et al., 1 a sequential quadratic programming algorithm was used to minimize the misfit between model and measurement data. To further constrain the optimization, they introduced a first order Tikhonov regularization, which smooths the amount of muscle shortening. In contrast, we implemented spatial and temporal second order Tikhonov regularization to smooth the change of active tension among the neighboring elements and in the time course of the active tension.
The first experiment was the reconstruction of the reference case with the ellipsoid geometry without any model parameter mismatch and it delivered promising results for this least error prone case. Already in this control case, spatial regularization ( λ 2 ) was necessary and was applied for all other cases as well. λ 2 was set to 10 −14 striking a reasonable trade off between the reconstruction in presence of infarct regions ( Figure 10) and the reconstruction accuracy for entirely healthy tissue ( Figure 5). Furthermore, temporal regularization improved the result in case the added variation to the displacement field of the target surface was changing over time. For the reference case with the image-based geometry, we set λ 2 to 10 −12 . As mentioned in Chapter 2.1, the regularization parameter depends on the squared value of the residual norm. The matrix involved in the residual norm depends on volume of each finite element, which we approximated by the third power of the mean surface edge length. For the used geometries, the third power of the ratio between edge mean length of the image-based left ventricle and the ellipsoid is (4.38/1.9) 3 = 12.25. This leads to a factor of about 12 2 between the regularization parameters for the image-based and ellipsoid geometries. In this study, we determined the regularization parameters with the knowledge of the ground truth. However, when using clinical data and a priori method such as the L-curve criterion is needed to determine the values of the regularization parameters.
In the case of fiber variation, an infarct area of varying size was always falsely detected (specificity was below 1.0). As the limits of the fiber variation lay just at the level of accuracy with which fiber orientation can be measured by diffusion tensor MRI, 25 the reconstruction result is not satisfying. Furthermore, the generation of the fiber orientation F I G U R E 1 6 The root mean square error and its distribution for all elements of the image-based left ventricular geometry for Cases 11-18 as described in Table 1. The black dot represents the median of the peak tension error and the gray bar represents the error between the first and third quartile using a rule based model differs from the real fiber orientation and thus introduces further fiber variation. 26 Additionally, our geometrical model consists of a rather coarse mesh, as it has only one transmural element and thus has only two distinct fiber directions assigned to the quadrature points across the wall. Since we were not using structured meshes, we obtained diverse fiber orientations (Figure 2A,C). Nevertheless, the fibers pointing in circumferential direction are underrepresented compared to the fiber distribution described by Streeter et al. 16 This can be improved by increasing the number of finite elements in transmural direction, which will also lead to a smoother transition between the fiber orientation defined on the endocardial and epicardial surfaces and more continuous transmural rotation of the fibers. In addition, it might be generally impossible to obtain certain deformations with different fiber orientations. In future, to overcome the difficulties due to the approximation of fiber direction, we will extend the inverse solver to first estimate the direction of the fibers and then to estimate the active tension. For the estimation of the fiber orientation, additional input data might be used, for example, including rotational information. In addition, a cross-fiber active tension can be implemented to alleviate the difficulty arising from uncertain fiber orientation. In contrast to our results, Otani et al. 6 stated that their inverse model performed reasonably successfully with variation of the fiber orientation. However, they perturbed the fiber direction with a random error from a distribution of zero mean and a standard deviation of 10 . We chose a constant shift in fiber directions, of 10 on the epicardial or endocardial surface and respectively on both surfaces, which results in 20 deviation. Also, the group of Asner et al. 4 performed the reconstruction of active tension with different fiber angle ranges that deviated by 5 on both the endocardial and epicardial surface and concluded that the result is thereby not significantly affected. The reconstruction method of Asner et al. 4 used tagged MRI data as input data, where information about the rotation and twist of the tissue is available, whereas such information is not provided by the surface deformation used as input for our inverse method. In conclusion, uncertainty in fiber orientation impairs the accuracy of the estimated tension field.
In case of target surface variation, the results were as expected. In Case 4, Case 13, Case 15 and Case 17 in Table 1, an overestimation of the active tension was expected, as a stronger deformation towards the inside can be achieved by higher tension. In Case 5, Case 14, Case 16 and Case 18 in Table 1, again, the obtained result was expected, as lower tension leads to less deformation. The Case 6 in Table 1 provided very promising results by adding temporal regularization. The results with the ellipsoid and the image-based left ventricular geometry show the same tendency of the method to overestimate (Case 4 and Case 15) or underestimate (Case 5 and Case 16) the tension when the input deformation is noisy. For both geometries, the root mean square error in the same range, although the absolute values differ due to the different maximal tension values assumed in the forward models. In the image-based geometry, pressure was applied, which counteracts the active tension. Therefore, the maximum of the active tension (180 kPa) is higher compared to the idealized geometry (40 kPa), where no blood pressure was considered. Additionally, the initial wall thickness differs between the two geometries as described in Chapter 2.2. Thus, 20% of the relative thickness change results in different absolute displacements of target surfaces. Our results are in line with the work of Asner et al., 4 where the displacement noise in the synthetic data was set to 10%-20% of the maximal displacement. The reconstruction of the global active tension parameter resulted thereby in minimal increase in error. Also, in the work of Otani et al., 6 the method performed reasonably well when white noise of relative amplitude of 10% was added to the displacement data. Case 6 is the one expected to occur also using clinical data, due to the intra-observer variability during the segmentation of MRI data. 12 In the reconstruction with the image-based geometry, the overestimation was severe, while the underestimation was rather minor. This suggests that the method is robust against an overly thin segmentation of the ventricular wall. In particular, in future we intend to segment cine MRI short axis slices and use their endocardial surface as input of the inverse method.
For all four experiments in presence of an infarct in the forward simulation (Case 7-10 in Table 1), the size and the location of the infarct were reconstructed well. Additionally, the active tension distribution inside and around the infarct area might provide valuable diagnostic information. Nevertheless, the correct identification of an infarct area strongly depends on the choice of the regularization parameters. Too much spatial smoothness imposed on the solution will blur the boundary of the infarct area and will lead to misclassification of the infarct area. Furthermore, in this work we have not investigated how noisy input data influence the reconstruction of the size and location of an infarct area.
In case pressure was included in the model (Case 11 and Case 12 in Table 1), the reconstruction of the tension matched the ground truth very well. We applied the Hill pressure curve and obtained quasi-isovolumetric phases around the peak tension. This suggests that the presented method will be able to reconstruct the active tension correctly also during isovolumetric phases. In Case 11, we applied the sine pressure curve since it is more similar to the shape of measured pressure curve compared to the Hill pressure curve (Case 12).
In this work, the material properties were assumed to be homogeneous in the entire geometry, although they are expected to vary in the infarct region compared to healthy tissue. 27 The wall thickness of the image-based ventricle and the fiber orientation can also be affected significantly after an ischemic injury. In future, the sensitivity analysis can be extended to investigate the variation of the material properties, the wall thickness and the fiber orientation in infarct regions. The presented method can be adapted to estimate the parameters of the material law for each finite element based on the deformation information in the relaxation phase. Whereas the active tension can be estimated in the systolic phase, as it dominates the overall tension in this phase.
For the presented inverse method, the deformation of the endocardial wall is required as input data, whereas other methods expect three-dimensional displacement data [4][5][6] or strain and volume data. 1,2 We observed that the deformation of the epicardial surface aligned very well with the expected deformation of the forward simulation, probably due to the combination of an incompressibility term added to the material law and the transmural resolution of one element. This advantage of our method can be especially valuable for the reconstruction of the active tension of the right ventricle, as in MRI data the right ventricular epicardial surface is hardly visible. To apply the presented method for a hypertrophic left ventricle, higher transmural resolution could be required.
In this work, all quality measures are defined based on a comparison with the ground truth. For applications with clinical data, this will not be available. Therefore, further quality assessment methods should be investigated. We observed that higher tension error occurs together with an increased median of the mean gap for the cases with varied displacement data ( Figure 14). Thus, the distance between target and reconstructed surface can be potentially a useful indicator for the quality of the solution.
One of the possible applications of the proposed inverse method is to locate and detect the size of infarct scars and to identify the elastomechanical consequences. Furthermore, a reconstruction of the temporal and spatial distribution of the active tension field might be useful for diagnostic and therapeutic assessment. For example, wall stress is a sensitive marker of hypertrophic cardiomyopathy. 28 In addition, the presented inverse method can be adapted to provide the temporal and spatial distribution of the material properties of the myocardium as they also influence the stress field in the tissue.

| CONCLUSION
In this work, we presented a novel numerical method including second order Tikhonov regularization for solving an inverse problem of cardiac biomechanics. By imposing spatial regularization in the reference case, the active tension can be reconstructed well (maximal peak tension error <0.5 kPa). To obtain a maximal RMSE <20% with the presented method, it was necessary to provide the fiber orientation with an accuracy of 10 in our experiments. In case temporal variation up to ±20% of wall thickness was present in the input deformation data, imposing temporal regularization in the inverse method could decrease the MRTE to 12% (from 95% without temporal regularization). Despite the reconstruction errors in the active tension, different sized infarct areas could be identified accurately in the inverse solution (sensitivity >0.95). The morphology of the reconstructed active tension curve matched the one of the ground truth in all cases, whereas the amplitude of the reconstructed tension differed compared to the ground truth. The results obtained through this novel inverse method are promising but further effort is needed to apply the method on motion data from clinical images and improve the performance considering reconstruction with uncertain fiber orientation.