An improved parameter estimation and comparison for soft tissue constitutive models containing an exponential function

Motivated by the well-known result that stiffness of soft tissue is proportional to the stress, many of the constitutive laws for soft tissues contain an exponential function. In this work, we analyze properties of the exponential function and how it affects the estimation and comparison of elastic parameters for soft tissues. In particular, we find that as a consequence of the exponential function there are lines of high covariance in the elastic parameter space. As a result, one can have widely varying mechanical parameters defining the tissue stiffness but similar effective stress–strain responses. Drawing from elementary algebra, we propose simple changes in the norm and the parameter space, which significantly improve the convergence of parameter estimation and robustness in the presence of noise. More importantly, we demonstrate that these changes improve the conditioning of the problem and provide a more robust solution in the case of heterogeneous material by reducing the chances of getting trapped in a local minima. Based upon the new insight, we also propose a transformed parameter space which will allow for rational parameter comparison and avoid misleading conclusions regarding soft tissue mechanics.


Introduction
Formulating an accurate constitutive law for soft tissues has been a contentious research topic for several decades (Maurel et al. 1998, chap. 4). Significant advances have been made since the seminal work by Fung and others, and a multitude of hyperelastic constitutive laws have been proposed for describing the stress-strain behavior of different soft tissues. These include myocardium (Humphrey and Yin 1987), arteries , ligaments (Natali et al. 2003;Weiss and Gardiner 2001), heart valves (May-Newman and Yin 1998).
A common feature of many of the proposed constitutive laws is the presence of an exponential function. This stems from the classic study by Fung et al. (1972), which demonstrated that stiffness of the soft tissues is proportional to stress. The exponential nature has been shown to be a result of collagen fiber recruitment and rotation that happens at the microstructural level (Lanir 1983;Billiar and Sacks 2000). However, fiber-level mechanics entails mesoscale calculations leading to high computational cost. Therefore, the phenomenological models containing an exponential are seen as better suited for tissue-or organ-scale biomechanical studies.
With a wealth of insight available on the suitability of constitutive laws, focus has been increasing on using them within biomechanical models to predict the mechanical behavior of soft tissues. Hence, accurate determination of the elastic parameters involved in the stress-strain relationships is a critical step in such predictive modeling. Various ex vivo and in vitro testing methods have been developed, and recently, methods applicable to in vivo dataset are gaining more attention as they allow the elastic properties to be characterized in tissues' native environment.
Exponential function results in a highly nonlinear stressstrain relationship, even if we discount the geometric non-linearity due to large deformation. In spite of this marked difference between soft tissue constitutive laws and other commonly used strain energy functions, e.g., for linear and rubber-like elastic material, standard techniques are used for parameter estimation and comparison of soft tissues. We aim to study the effect of nonlinearity of an exponential on soft tissue constitutive laws and related elastic parameters, and seek to improve upon the existing standard methods. This work is motivated by results from our recent study (Aggarwal and Sacks 2016), where an inverse model for bioprosthetic valve was developed. The study was designed to determine mechanical properties of a bioprosthetic valve leaflet by matching its deformed shape. The constitutive law contained an exponential function σ ∼ Ae B (although within a weighted integral and with neo-Hookean term), and the proposed framework included estimating two elastic parameters A and B (in the original work c 0 and c 1 were used, instead we use A and B for consistency in this manuscript). It was observed that the objective function contained a long and narrow valley in the parameter space ( Fig. 1a), which resulted in a slow convergence of the inverse model, especially once the iterative solution entered the valley.
The objective function with a long narrow valley is reminiscent of Rosenbrock function used in optimization textbooks (Fig. 1b), which is a challenging minimization problem because of its shape (Rosenbrock 1960). Furthermore, the flat shape of the valley means that parameters along the valley generate closely similar stress-strain response. Thus, the two parameters A and B are highly covariant along this valley. On the other hand, changing parameters transverse to the valley dramatically affects the stress-strain response (Fig. 2). Therefore, a simple comparison of elastic parameters for tissue samples might present a grossly wrong picture. For example, parameters corresponding to points 1 and 5 are much farther apart compared to points 8 and 9 (Fig. 2). However, the former produce a much closer stressstrain response compared to the latter.
These observations raise multiple questions: did the valley shape result from that one particular problem or is it a general feature of soft tissues? Can we improve the convergence of parameter estimation? Is there a more rational way  Aggarwal and Sacks 2016) to compare the elastic parameters of two tissue samples? If only two parameters lead to a challenging parameter estimation, how will the problem behave in case of heterogeneous media (four or more parameters), where one cannot visualize the objective function? Here, we aim to answer these questions by using ideas from elementary algebra and analyzing multiple cases by extending those ideas.

Methods and cases considered
Before starting the analysis, we present implementation details of the various functions and biomechanical models studied here. The problems were divided into two categories: displacement-controlled (DC) and force-controlled (FC). In DC cases, the input variable was the deformation or strain, and stress or forces were fitted to estimate parameters. On the other hand, in FC cases, the input variable was the force or stress, and strain or deformation were fitted.

One-dimensional curve fitting
Starting in one dimension, the simplest function with an exponential is Here represents strain and is the input, while σ is the output representing stress. Thus, (1) is a DC case, and the inverse of this relationship is an FC case. Since the stress function (1) is not zero at = 0, we also considered a more realistic representation of stress σ (A, B, ) = A e B − 1 and its inverse (3) We considered a model with higher nonlinearity: Lastly, the exponential function is sometimes truncated and linearized beyond an upper bound on strain ub (Fan and Sacks 2014). This may represent the condition when all the fibers have been recruited, thus producing a linear stressstrain response thereafter. Accordingly, we considered the following model The stress function (6) was defined such that both stress σ and stiffness ∂σ/∂ remain continuous at = ub . In all one-dimensional problems, the input was used to calculate output which was then matched to observed data in order to obtain the optimum parameters A and B. We call this process "curve fitting."

Multi-dimensional curve fitting
Hyperelastic constitutive laws for soft tissues are generally defined in two or three dimensions, based upon deformation gradient F (Holzapfel 2000). We studied two commonly used constitutive laws here; however, the results are expected to be extensible to most others. First is the Gasser-Ogden-Holzapfel (GOH) model (Gasser et al. 2006): where Q = (κ I 1 + (1 − 3κ)I 4 − 1) 2 , I 1 = tr(C) is the first invariant of C = F T F, I 4 = C : N ⊗ N is the fourth invariant with material direction vector N, and κ controls the "degree of anisotropy." matrix represents the contribution from the ground matrix in the tissue, and it is modeled as a compressible neo-Hookean solid: where J = √ det (C) represents the volume change. Second is the simplified structural model (SM) (Fan and Sacks 2014): where E = 1 2 (C − I) and (θ) is the fiber orientation function. (θ) is modeled as a truncated Gaussian in angle θ ∈ [−π/2, π/2] with standard deviation (SD) τ and peak at θ = ω: where P = dθ normalizes the distribution.
If we assume that microstructural parameters, such as κ, N, d e , ω and τ , and ground matrix properties are known, both constitutive models (7) and (9) have only two elastic parameters A and B to be determined. Even an approximation of the ground matrix elastic properties provides a good estimate of the overall mechanical behavior for many cases, such as bioprosthetic valves (Aggarwal and Sacks 2016). For both models, second Piola-Kirchhoff stress can be derived by the standard relation S = 2 ∂ ∂C , and it is easy to see that S ∼ Ae B(·) . In (7), the exponent is BQ, whereas in (9), we have an integral of exponentials. Thus, the two models represent two distinct classes of constitutive laws. GOH (7) and SM (9) models were used in multidimensional curve fitting, where known deformation gradients were used to calculate the stresses. These stresses were then matched to the observed stresses in order to determine the elastic parameters. Since the known input was deformation and stress was the output, these curve fitting problems belong to the DC category. We note that for these constitutive laws, there is no closed form solution of the inverse relation, i.e., strain or deformation gradient as a function of stress. Therefore, multi-dimensional curve fitting could not be performed for the FC case.

Inverse models
For cases where an explicit relation from input to output is not available (such as the FC case for multi-dimensional problems) or where the input parameters vary spatially and cannot be represented using a single set of values, curve fitting cannot be performed. In such situations, it is more appropriate to solve an inverse model, where the observations are matched to the outcome of a finite element simulation. This finite element model uses a predetermined constitutive law, and we tested both GOH (7) and SM (9) models.
We considered two inverse modeling problems, and all of the finite element simulations were performed using FEBio (Maas et al. 2012). First problem is that of a biaxial testing of a thin planar tissue sample (Fig. 3), which was studied for both DC and FC cases. In the DC case, known uniform displacement boundary conditions were applied on the sample edge, and total reaction forces on the edges were matched to the observed values. On the other hand, in the FC case, known uniform forces were applied on the sample edge, and average edge displacements were matched to observed values. A detailed description of the setup, such as the sample size, boundary conditions, is specified in "Details of biaxial simulation." The second problem studied here is the shape matching of a semilunar tissue sample under static pressure loading ( Fig. 4), which represents the closing of a bioprosthetic valve leaflet. Since the input for this problem is pressure traction, only FC case was possible. Details of the shape matching procedure were described in our previous work (Aggarwal and Sacks 2016) and are summarized in "Details of tissue pressurization simulation." All of the problems considered in this study are summarized in Table 1.

Generic notation
For all the cases considered, both DC and FC, we define a generic notation for the analysis.f (x i ) denotes the observed data at independent input variable x i for i = 1 . . . N , and f (A, B, x i ) denotes the model. A and B are the parameters to be determined by fitting f (A, B, x i ) tof (x i ). In the present context of soft tissue mechanics x denotes the input applied during testing, e.g., applied strain, displacement or load,f represents the observed output, e.g., stress, force, deformed shape or strain, and f (A, B, x) is the same quantity computed using our model. Since we are interested in functions of the form f (A, B, x) ∼ Ae h (B,x) or its inverse, A has units of stress and B is dimensionless. To fit our model to the observed data, we define a functional F = f −f . The minimum point of this functional corresponds to the optimum parametersȂ,B = arg min A,B F, where f andf are "closest." We use the norm · in a general sense, and two options were explored: The 2-norm is the standard Euclidean norm, which = 2 dx for continuous input x ∈ (0, X ). The "lognorm" uses logarithm of f andf in the standard Euclidean norm. The logarithm function is denoted as log to clarify its modified form In general, there could be constraints on A and B in our minimization problem. These constraints are physically motivated, for example from thermodynamics of strain energy density and convexity requirements for stress function. Here, we considered two constraints that both A and B must be positive.
In order to exclude the effect of noise, the observed data were generated synthetically for known values of the param-etersĀ andB. Hence, we denote the observed data as f (x) = f (Ā,B, x), whereĀ andB are known a priori. Clearly, in this case, the global minimum of F should occur atȂ =Ā andB =B. All the parameters to be determined are collectively denoted as c, and the model f and dataf at all inputs x i combined into a vector are denoted as f and f, respectively. Lastly, we define two curves in the (A, B) parameter space where one of the first derivatives of the functional F vanishes: We name them A-partial minima (APM) and B-partial minima (BPM) curves, respectively, or PM curves collectively. The two PM curves intersect at the global minimum (Ȃ,B), and the shape and adjacency of the two curves help develop an intuitive and qualitative understanding of the analysis. Whenever possible, closed form expressions were evaluated for PM curves. For other problems, the PM curves were determined numerically. For APM, this was done by fixing B at various values and minimizing F with respect to A, and vice versa for BPM.

Parameter estimation
In order to calculate the elastic parameters for soft tissues, a numerical optimization has to be performed to minimize the function F(c). We focus on gradient-based line search methods, which proceed iteratively in two steps: (1) determine a direction along which F will decrease and (2) determine the step size to move along that direction (also known as line search). For calculating the direction in step 1, the simplest choice is the steepest descent direction −∇F. However, it leads to extremely slow convergence for functionals with a narrow valley, such as the Rosenbrock function (Nocedal and Wright 2006). On the other hand, using the full Newton's method gives second-order convergence. However, it requires calculation of second derivatives for the Hessian ∇ 2 F and solving a linear system of equations (15), making it computationally expensive. For least-square functionals, such as the 2-and log-norms (11,12), their form allows a simpler approximation of the Hessian. If J = ∂f/∂c, then the functional gradient is ∇F = J T (f − f) and the first-order approximation of the Hessian is ∇ 2 F ≈ J T J. This approximation leads to the Gauss-Newton algorithm which gives approximately second-order convergence and requires only the first derivative to be computed.
Comparing the different gradient-based methods for fitting the one-dimensional exponential function (1), we found that the Gauss-Newton algorithm performed the best. This is consistent with observations about the Rosenbrock function (Nocedal and Wright 2006). Henceforth, we used the Gauss-Newton algorithm (Algorithm 1) for all problems. For the line search in step 2, we used a simple backtracking algorithm, which took into account any constraints on the parameters and situations of failed f calculation (e.g., due to non-convergence of the finite element solver).

Algorithm 1: Parameter estimation using Gauss-Newton method with backtracking line search
Data: Observed data f and initial guess c 0 Result: Parameters that fit our model f to observed data f in chosen norm and thus minimize the functional Calculate the fitting model f(c) and its derivatives J = ∂f(c)/∂c using central finite difference f(c(1 ± δ)) ;

5
Calculate the search direction and step c by solving (16)

O L) and (I T E R < M AX I T E R) and
(max | c| > δ) ; Values of TOL = 10 −10 and δ = 10 −5 were used for all calculations. Convergence of nonlinear minimization can strongly depend on the initial guess or the starting point (denoted using subscript 0, c 0 , etc.). Therefore, for each problem, multiple minimizations were performed with starting points spanning the parameter space. Hence, 10 × 10 starting guesses were chosen uniformly distributed in the span of A 0 ∈ [0.005, 0.1] and B 0 ∈ [20, 100]. The resulting convergence statistics-number of iterations (#Iter), number of f(c) evaluations (#Eval) and number of parallel evaluations (#Paral)-were reported as mean ± SD. The number of parallel evaluations was important since the central difference calculations for determining J were done in parallel. However, during the line search, evaluations had to be performed sequentially, adding to the computational cost.

Noise
In order to study the effect of noise on the accuracy of the estimated parameters, we added random noise of varying magnitude to the target vector: Here ν represents the noise level, which was varied from 0.01 to 0.04, and the random number between −1 and 1 had uniform probability. The effect of the noise was quantified by an error in the estimated optimum parameters: Clearly, e(0) = 0 since without noise the global minimum coincides with the true minimum.

Heterogeneous model
In all of the problems considered so far, we assumed that the elastic parameters did not vary over the tissue sample. However, heterogeneity is a common feature of biomechanical systems. In order to study the parameter estimation properties for a heterogeneous system, we considered the simplest problem of two materials in a biaxial testing setup. That is, the tissue sample is made of two types of tissues with two sets of elastic parameters-(A 1 , B 1 ) and (A 2 , B 2 ) (Fig. 5). The boundary and loading conditions remained the same as in the biaxial inverse model ("Details of biaxial simulation"). We considered only the SM constitutive law for both DC and FC cases. If both tissues in the sample have the same microstructural properties, then the system is symmetric, i.e., both give exactly the same response. This leads to a singular Hessian whenever A 1 = A 2 and B 1 = B 2 . Therefore, in order to break this symmetry and make the Hessian non-singular, we used τ 1 = π/6 and τ 2 = π/7. In total, there are four unknown elastic parameters-A 1 , B 1 , A 2 and B 2 , which leads to a four-dimensional (4D) parameter space. As a consequence, the functional cannot be visualized and scanning the entire parameter space for starting points becomes prohibitively expensive. Therefore, we only scanned a diagonal plane in that 4D space, where the starting parameters were the same for the two materials: A 1,0 = A 2,0 and B 1,0 = B 2,0 . Hence, 8 × 8 starting guess points were chosen which were uniformly distributed in the span of log(A 1,0 ) ∈ [−4, −2] and B 1,0 ∈ [20, 100]. The parameter estimation was performed without adding any noise to the system.

Analysis
Before carrying out parameter estimation, we analyze the properties of various models described in the previous section and the functional F constructed for them. The analysis is divided into DC and FC cases.

Displacement-controlled cases
We start by looking at the simplest curve fitting involving an exponential function (1). Assuming that the "observed" data has no noise and that the number of observations is infinite, Here, E > 0 defines the maximum value of strain at which the stress has been observed. Our model (1) needs to be fit to σ and, thereby, determine the optimum parameters A and B.
As a first step, we simply take our functional as the 2norm of the difference between the observed data and model: (27) (details of all analytical derivations are in "Details of the analytical derivations"). Clearly, F = 0 at A =Ā and B =B, and F > 0 everywhere else. Thus, it has a global minimum at (Ā,B), which corresponds to the true solution. Also, it can be easily verified that there are no other local minima in this functional. That is, (Ā,B) is the unique global minimum and F is convex everywhere. Furthermore, this is a highly nonlinear functional, and its contour plot contains a long narrow valley similar to that observed in the inverse model of a bioprosthetic valve (Figs. 2, 6a). Hence, interestingly, curve fitting of a simple exponential function reproduces the behavior seen for a complex inverse modeling problem.
For this function, the APM and BPM curves can also be determined in closed form (29). These curves lie very close together in the valley of the functional (Fig. 6a). Since a point in the parameter space where both derivatives are zero represents the global minimum, the proximity of the two PM curves implies that both derivatives are approximately zero in the whole valley region. That is why, even though the functional is convex, parameters along the valley produce a similar response σ ( ), as observed previously (Fig. 2).
Based upon elementary algebra, for fitting an exponential function, it is significantly easier if the function is linearized by taking a logarithm before the 2-norm. That is, F = log(σ ) − log(σ ) 2 . This leads to a quadratic functional in log(A) and B (31), where both APM and BPM curves are straight lines (33) (Fig. 6c). Clearly, this norm also satisfies the condition of a unique global minimum at (Ā,B). More importantly, the two PM curves are not close anymore and are visually distinct (Fig. 6c). This is reflected in the functional shape as an absence of a valley.
Since σ > 0 and σ > 0 for this model for all strain values, the norm F = log(σ ) − log(σ ) 2 is equivalent to the log-norm (12). There is one subtle difference between the 2-norm and log-norm functionals: the latter is defined in a transformed parameter space of (log(A), B) instead of (A, B). Therefore, in order to objectively compare the two norms, we look at the 2-norm in the (log(A), B) parameter space (Fig. 6b). In this case, the APM and BPM curves are still close together; however, the valley becomes relatively straight.
We extend these ideas to function (3), which is a better representation of stress-strain relation. With 2-norm, we obtain a similar curved valley in the (A, B) space (Fig. 6d), where the PM curves are even closer together. In the (log(A), B) space using 2-norm, the PM curves and the valley are more straight (Fig. 6e). For the log-norm, modified log function (13) has to be used, since σ goes to zero for = 0 making its log undefined. For (3) with log-norm, it is not possible to obtain analytical expressions for the functional and PM curves. Thus, these were determined numerically using representative values ofĀ = 0.02 andB = 45. Again, we observe disappearance of the valley, and the two PM curves become distinct for B 10 (Fig. 6f).
We perform similar analysis for multi-dimensional constitutive laws using DC curve fitting. In this case, the functional and PM curves cannot be determined analytically for either of the norms, so numerical calculations were performed using SM and GOH models for the same representative parameter values (Ā = 0.02 andB = 45). The resulting functionals behave very similar to the one-dimensional problems. The 2-norm has a similar valley which is curved in the (A, B) space and practically straight in the (log(A), B) space, and the PM curves are extremely close in both cases ( Fig. 6g-l) 1 . Using the log-norm, the PM curves become distinct and no valley is observed (Fig. 6i, l). Lastly, a small difference can be noticed between the 2-norm functionals for two constitutive laws: the PM curves for SM model are closer together as compared to the GOH model.

Force-controlled cases
Similar to the analysis of DC cases, we start with the simplest FC problem: the inverse of an exponential function (2). We assume that the observed data is of the same form without any noise, i.e., (σ ) = log(σ/Ā) /B, and that the observations were made continuously in the range σ ∈ [1, ]. Here, is Fig. 6 Functional for displacement-controlled (DC) cases plotted as a contour; global minimum is indicated using a green circle, and the PM curves are plotted using colored lines-APM (blue) and BPM (red).
Left and center column plots are using 2-norm in (A, B) and (log(A), B) space, respectively, while the right column plots are using log-norm the maximum stress at which strain was observed. In order to fit our model to the observed data, as a first step, we take the 2-norm functional: F = − 2 . This can be evaluated analytically and results in a rational functional (35). Similar to the DC case, a long and narrow valley is observed (Fig. 7a).
Here also, we find that the PM curves (37) and (39) essentially overlap. Furthermore, the valley and PM curves become relatively straight in the (log(A), B) space (Fig. 7b)-another similarity to the DC case. However, unlike the DC case, changing the norm to lognorm, or any other norm, does not make the model linear, and the valley shape of the functional persists. Instead, if the function (2) is rewritten as: it is easy to see that defining a new pair of parameters makes the strain (σ ) = β log(σ ) − α linear in parameters α and β. Hence, the parameter estimation problem using 2-norm becomes a linear least-square problem. In other words, the transformation of parameters from A, B to α, β makes the 2-norm functional F(α, β) = (α, β) − 2 quadratic (Fig. 7c). Concomitantly, the PM curves (36) and (38) become straight lines that are distinct from each other.
This idea is extended to function (4), where, even though the transformation does not make the function exactly linear, it exhibits a similar behavior in the functional shape and PM curves (Fig. 7d-f). That is, in the original parameter space (A, B), the 2-norm functional has a narrow and curved valley. Transforming the parameter space to (log(A), B), the valley becomes straight, but the PM curves remain close together. Lastly, using the transformation (log(A)/B, 1/B), the valley ceases to exist, and the PM curves become distinctly different.
For the FC cases, curve fitting cannot be performed for multi-dimensional constitutive models. Instead, we calculate the functional and PM curves for biaxial inverse model using SM. Even this multi-dimensional problem, which involves a highly complex stress-strain function, behaves very similar to the 1D problems. With 2-norm, we obtain a narrow and curved valley in the (A, B) space (Fig. 7g), which becomes straight in the (log(A), B) space (Fig. 7h). In both of these cases the PM curves remain close together. However, in the (log(A)/B, 1/B) space we do not see any valley and the PM curves become distinct (Fig. 7i). The functional shape and PM curves can be useful indicators for the parameter estimation, as we see next.

Parameter estimation
Based upon the analysis, we compare the following four formulations for the DC cases: On the other hand, for FC cases, we compare the following three formulations-all with 2-norm:

Displacement-controlled cases
All four DC cases tested using the algorithm with four formulations were summarized (Fig. 8). The 2-norm functional was plotted versus iterations for one representative starting point, and average statistics for all starting points were tabulated. There was a consistent pattern that the number of iterations needed to converge decreased as we changed from formulation 1-4, with log-norm in (log(A), B) space producing the best results. The decrease was modest in some cases ( Fig. 8c), while it was as high as 60% in others (Fig. 8a,  d). Furthermore, convergence was slower with SM model (Fig. 8b, d) compared to the GOH model (Fig. 8a, c). In fact, the iterations for SM curve fitting with 2-norm diverged for a few initial guesses. Secondly, we also found a clear decrease in the SD of number of iterations as we go from case 1-4. This shows that by using log-norm and (log(A), B) space, the problem becomes uniformly converging.
A couple of worthwhile remarks about the implementation of the proposed algorithm: it was observed that when the minimization routine uses a finite element solver for model calculation (FEBio in this study), care must be taken with the floating point precision of input and output. If input and output are in ASCII format, severe truncation errors can be introduced into the minimization routine. These truncation errors can adversely affect the proposed minimization algorithm, especially with small TOL and δ values. This can be Parameter estimation results for four force-controlled (FC) problems with one representative example of convergence each (left), and statistics with starting points spanning the parameter space (right) summarized as mean ± SD alleviated to some extent by increasing the precision of floating point numbers in the input and output file format. Effect of this truncation error can be seen in convergence plots of inverse problems (Fig. 8c, d), where the functional did not decrease beyond a value. This minimum value did not change with the formulation and thus was only dependent on the problem setup. Curve fitting problems did not show such a minimum (Fig. 8a, b). Secondly, the Gauss-Newton equation (16) sometimes gives large step sizes ( c), such that the current guess c may go outside the physiological range. This can adversely affect the line search and calculation of f may fail occasionally. This is why the condition of function calculation being successful was added in the algorithm (Algorithm 1, line 10).

Force-controlled cases
For all FC cases, the parameter estimation improved as we changed from formulation 1-3 (Fig. 9). The formulation with 2-norm and (log(A)/B, 1/B) space performed the best in all problems. However, the improvements in convergence speed were smaller than those in DC cases, with the maximum reduction in iterations required being 35% (Fig. 9d). The SD also decreased as we went from formulation 1-3, showing that problem converged more uniformly. Lastly, the GOH model (Fig. 9a, c) performed marginally better than the SM (Fig. 9b, d) in the (A, B) space. However, the improvements due to new formulations were smaller for GOH model than those for the SM. Since all of the problems were inverse models that involved interfacing with FEBio, functionals did not decrease beyond a minimum value (similar to the DC inverse models).

Noise
When noise was introduced into the observed data, the global minimum of the functional moved away from the true solution for multi-dimensional curve fitting with both SM (9) and (a) (b) Fig. 10 Change in error as a function of random noise introduced into the target vector for DC cases. a structural model and b GOH model GOH (7) models. As the level of noise (ν) was increased, the error in estimated parameters also increased. However, for both models, the increase in error was significantly lower using the log-norm as compared to the 2-norm (Fig. 10). Change of parameter space between (A, B) and (log(A), B) had no effect on the error. Since norm was not changed in the FC cases, the variation of error as a function of noise remained the same irrespective of the parameter space used (results skipped for brevity).

Heterogeneous model
The true parameters (Ā 1 ,B 1 ,Ā 2 ,B 2 ) were chosen from the valley region of the 2-norm functional (i.e., on either APM or BPM curve-the two curves essentially overlapping). This led to a system where the two materials produced a similar response and were difficult to differentiate. Only the convergence/non-convergence 2 was recorded and compared between the two formulations that performed best and worst in the homogeneous parameters estimation. That is, for the DC case, (a) 2-norm with (A, B) space and (b) log-norm with (log (A), B) space were used. The results from various starting points for the DC case showed a clear difference between the two formulations (Fig. 11). Even though the second formulation did not give perfect results (50 out of the 64 initial guesses led to convergence), it performed significantly better compared to the first formulation (only 25 out of the 64 initial guesses led to convergence). In the FC case, two formulations with (a) (A, B) space and (b) (log(A)/B, 1/B) space were used, both with 2norm functional. The improvement in this case was much more significant (Fig. 12). Only 4 out of the 64 starting guesses resulted in convergence using the original parameters, whereas 63 out of 64 starting guesses resulted in convergence using the transformed parameters. Furthermore, the DC and FC cases had closely equivalent boundary conditions and true parameters (Ā 1 ,B 1 ,Ā 2 ,B 2 ). Therefore, it is worthwhile comparing their performance. In the standard  and (A, B) space), the DC case was better conditioned than the FC case and led to convergence more often (Figs. 11a, 12a). However, in the modified formulation, FC case significantly outperformed the DC case (Figs. 11b, 12b).

Parameter comparison
The valley in the functional is a region of high covariance between the two parameters A and B, and parameters along this valley produce similar stress-strain response while those across it do not (Fig. 2). Therefore, knowing the equation of the valley will greatly facilitate parameter comparison. Using  (log(A), B) parameter space with 2-norm, a straight valley was observed in the functional for all problems (Figs. 6, 7). Therefore, as a generalization, we approximate the equation  (3), (5) and (6), respectively, and d and e for multi-dimensional GOH (7) and structural (9) models, respectively each specific setup. Thus, we determined the relation between the slope of valley and various parameters for several 1D and multi-dimensional problems, with the aim of establishing a general criterion on how to approximate the slope without explicitly calculating it. In all problems, we numerically verified that the slope of the valley did not depend on the true parameters (Ā,B) (within the physiological range). In 1D, for exponential function (3), the slope was well approximated by the maximum value of strain: m ≈ E (Fig. 13a). This can also be confirmed by manipulating the analytical expression of PM curves (29). For exponential function with higher nonlinearity (5), m ≈ E 2 (Fig. 13b). For truncated exponential function (6), we found that m ≈ if ≤ ub and m ≈ ub otherwise (Fig. 13c). Thus, in general, the slope of the valley is well approximated by the highest coefficient of B in the exponential.
Similarly, for multi-dimensional GOH model (7), slope was approximately equal to the maximum value of the exponent Q: m ≈ Q max (Fig. 13d). For the structural model (9), stress is an integral of the exponential function with varying exponent B (N · E · N), where ρ(N) = N · E · N is the extension in direction N. Since there is not a single value of the exponent at the maximum load, we denote the maximum and minimum coefficients of B as ρ max and ρ min , respectively. Next, we assume that the slope m lies between these two values, i.e., ρ min ≤ m ≤ ρ max . In order to verify this assumption, we define an indicator We found 0 ≤ χ ≤ 1 (Fig. 13e) for all cases, thus proving our assumption. Moreover, for most of the cases, the slope was well approximated by the maximum extension ρ max (i.e., χ ≈ 1). Only in the case when the principal stretches differed significantly and the maximum stretch was not aligned with the material direction, χ ≈ 1. Thus, in most of the situations, it is possible to approximate the valley's slope simply from the maximum applied strain and tissue microstructure.

Effect of exponential function
Soft tissue's mechanical response is characterized by high nonlinearity, anisotropy and non-homogeneity. The nonlinearity feature is represented as an exponential function in many of the stress-strain relationships (Maurel et al. 1998, chap. 4). Motivated by some observations from our previous work (Aggarwal and Sacks 2016), we aimed to analyze the effect of this high nonlinearity on elastic parameters. We found that it is a fundamental feature of the exponential function that leads to a narrow valley in the functional, irrespective of the complexity of the problem setup or constitutive law details. Furthermore, we categorized all the problems into DC and FC cases and found a similar functional shape for both problems. That is, the functionals contained a characteristic valley shape which was related to the proximity of two PM curves (Figs. 6, 7). This functional valley resulted in an ill-conditioned system, high covariance between the two elastic parameters, slow convergence during param-eter estimation, and misleading conclusions in parameter comparison.

Significance of presented results
In order to solve these challenges, we drastically simplified the problem by looking at an exponential function (1). Although this may seem like an oversimplification at first, the function reproduced the behavior observed for soft tissues. Thus, extending the ideas from elementary algebra, we found two different solutions for the two categories of problems. For DC cases, we found that using log-norm and (log (A), B) space made the parameter estimation close to linear and improved the condition of the problem. This was graphically depicted using two PM curves and separation between them. On the other hand, for FC cases, we found that by keeping the 2-norm but using (log(A)/B, 1/B) space, the problem's condition was significantly improved. These qualitative changes were consistently observed irrespective of the specific details of the problem (Figs. 6, 7). The new formulations led to an increase in the convergence speed for all problems considered (Figs. 8,9). The increase in convergence was significant in some problems, while limited in others. In general, the improvements were smaller for the FC cases compared to those in DC cases. Also, the closeness of two PM curves was found to be related to the convergence speed, thus proving to be a useful indicator. Secondly, we found that using 2-norm and (log(A), B) space, the functional contained a straight valley for all problems studied here. The slope of the valley was found to be well approximated by the highest coefficient of B for most of the problems (Fig. 13). Therefore, it is possible to estimate the slope of the valley in the 2-norm functional for a specific problem based upon the applied strains and microstructural properties, without performing any minimization or parameter estimation. It was established that changing parameters along the valley minimally affects the overall stress-strain response, whereas changing the parameters across the valley dramatically varies the response (Fig. 2).

Proposed modifications
Based upon these results we propose to use log-norm for DC cases and 2-norm for FC cases. Furthermore, we propose the use of (log(A), B) space and (log(A)/B, 1/B) space for DC and FC cases, respectively. We note that for the second proposition, since A has the units of stress, taking its log would be physically nonstandard. Therefore, we propose a modified form of the constitutive laws. Instead of σ (A, B, ) ∼ Ae B , we propose to use a form where η is a known and specified number with units of stress (e.g., stiffness of collagen fibrils). This modification, in addition to making the parameterÂ dimensionless, can also be used to scaleÂ so that its value is in the same range as B, thus improving the numerical conditioning of the problem. Another benefit of using this modified form is that noŵ A ∈ R is not restricted to be positive as A was, thus reducing the constraints on the parameter estimation problem.
Using the valley's slope in the functional space, which can be approximated for a given problem setup, we define a new set of parameters (Fig. 14 These parameters may be used to rationally compare elastic parameters of tissue samples, so as to truly represent the difference between their resulting stress-strain relationship. We call ζ the "major parameter" as it varies the stress-strain relation significantly, and η the "minor parameter" as its effect on the stress-strain relation is minimal.

Faster convergence
One of the primary effects of the new formulations was the change in the functional shape (absence of a valley), which led to faster convergence of the parameter estimation (Figs. 8,9). Although faster convergence may seem inconsequential for curve fitting procedures, it can save substantial computational time for inverse models where each iteration may take minutes or hours to compute. This is becoming more important as inverse models are increasingly being recognized as the most appropriate way of characterizing mechanics of soft tissues (Rausch et al. 2013;Zhang et al. 2015). Curve fitting assumes uniform microstructural properties and uniform strain, which usually do not hold true for soft tissues because of their heterogeneity. Therefore, inverse modeling is seen as a more accurate tool, and improved parameter estimation has the potential for making them more widely usable and, possibly, suitable for bedside computations.

More than improved speed
Although faster convergence is an important result, the changes in the functional shape with the proposed formulations had other implications as well. We found that the parameter estimation becomes more robust with respect to noise when using log-norm for DC cases (Fig. 10). However, the effect of noise only depends on the norm and not the parameter space used. Thus, no improvement was observed for FC cases with the new formulation. Interaction between parameters of multiple materials in a heterogeneous system may lead to the formation of local minima, even with simple constitutive models. Here, a problem with two sets of elastic properties was tested. We found that when using the modified formulations, chances of getting trapped in a local minima were significantly reduced (Figs. 11,12). Furthermore, we also found that FC case performed better than the DC case in heterogeneous inverse model. Thus, if there is an option to choose between FC and DC settings for a study, the former might provide more robust results using the new formulation.

Extension to multiple parameters
In this work, we used two forms of the constitutive laws; however, the results presented are expected to be general and applicable to other forms involving an exponential function. The analysis could also be extended to the case of multiple parameters, although with additional complexity. For instance, the model by May-Newman and Yin (1998) describes strain energy density (F) = A e Q − 1 , where 1 4 , and thus has three elastic parameters (A, B 1 , B 2 ). In this case, in the three-dimensional parameter space, we obtain a valley in the AB 1 -plane for each fixed B 2 , and vice versa in AB 2 -planes ( Fig. 15a, b). Despite these differences, the proposed modification leads to an improved convergence (Fig. 15c). Similarly, in models with multiple exponential terms Holzapfel and Ogden 2009), each term will lead to a valley while other terms are held fixed. Interaction of different valleys may lead to a complex behavior. For example, the heterogeneous model presented in this study is qualitatively similar to having multiple fiber families in the constitutive model , which led to the formation of local minima. Further analysis of such interactions among different terms is necessary, but is beyond the scope of present manuscript. Similar convexity analysis is also required for constitutive models with a large number of parameters, e.g., Fung's law (Labrosse et al. 2016).

General implications
We chose a minimization algorithm based upon Gauss-Newton and line search methods (Algorithm 1). However, the proposed modifications are expected to benefit all minimization routines for this class of problems. This is because of the fundamental changes in the functional shape that should positively affect the convergence properties irrespective of the algorithm used. In fact, many of the other methods, such as gradient-based Levenberg-Marquardt and evolutionbased genetic algorithm (Lee et al. 2014;Labrosse et al. 2016), perform poorly in the presence of narrow valleys in the functional (Nocedal and Wright 2006). Therefore, the improvements from proposed modifications may be even greater for those methods. Covariance between parameters plays an important role in the design of experiments. Here, we found that the slope of high covariance region is related to the maximum strain applied. That is, applying a higher strain leads to a more horizontal valley, which represents a lower covariance between the two parameters. The analysis presented here could be extended to find optimum points of observation so as to establish maximum confidence in the estimated parameters. This is extremely important for non-convex problems, since the solution could depend upon initial guess due to the presence of local minima (Abbasi et al. 2016). The proposed modifications are simpler for FC cases since it conserves the use of 2-norm. We also observed that the FC case performed better than the DC case in the heterogeneous parameter estimation, which indicates that FC setup may be better suited for inverse modeling. For all these reasons, we consider this work as the discovery of a new direction in soft tissue biomechanics that would lead to multiple future studies as discussed next.

Limitations and future work
In this study, we considered only elastic or hyperelastic problems, but similar analysis could be done for time-dependent models (such as viscoelastic), which are widely used for soft tissues. Furthermore, we restricted to only two parameters in order to facilitate the analysis and visualization of the functional shape. Extending the work to multiple parameters will be of significant value and will be completed in a future study. Although we only used planar tissues as examples here, the ideas are expected to be applicable to thick tissues, e.g., myocardium. Taking a log-norm assumes that all the stresses are positive, which is reasonable since exponential laws should not produce negative values. However, this may not hold true in general, especially for shear stresses. This issue was not addressed here and will be investigated in the future.
Our treatment of the heterogeneous media was not exhaustive. Examples were used to only demonstrate the advantages of the proposed modifications. As heterogeneous inverse problems become exponentially complex as the number of materials increase, we limited the current study to only two material types. Although the proposed modifications led to improved results, they were still not perfect as occasional non-convergence was observed. A more thorough analysis of these factors will be performed in the future, especially when noise is present in a heterogeneous system.
The minimization algorithm proposed here also could be improved, e.g., by using a better line search method, using forward difference instead of central difference for the gradient, carefully chosen δ value, etc. Hence, the focus here was not on designing the best minimization algorithm, but to use the insight on fundamental features that will benefit all minimization methods. Similarly, innovative ideas could be used to improve the convergence for specific problems. For example, using (ζ, η) parameter space and minimizing the functional with respect to only the major parameter ζ first, followed by complete minimization could provide significant speedup.
In addition to parameter estimation, the proposed modifications could also improve the material and shape optimization of bioprosthetics (Fan et al. 2013). These ideas could also be applied to other techniques, such as elasticity imaging (Oberai et al. 2009) and deformable image registration (Veress et al. 2005). Lastly, the non-elastic parameters were assumed to be known, which is true for several tissues, such as heart (Nielsen et al. 1991), aorta (Haskett et al. 2010) and aortic valve (Aggarwal et al. 2014). In the future, we will study the case when non-elastic parameters also need to be determined from inverse models, e.g., for pathological tissues.

Conclusion
In this study, the aim was to elucidate some of the issues related to soft tissue constitutive laws observed in our last study. We found that these features are fundamental to the exponential function, which is commonly used in soft tissue mechanics. By simplifying the problem and using elementary algebra, we proposed solutions that showed improved convergence and robustness with respect to noise. More importantly, the modified formulations were found to be less susceptible to being trapped in a local minima for heterogeneous problems. We also proposed modified parameters that will facilitate a rational comparison of elastic parameters. The insight obtained in this study will be used in the future to significantly improve inverse models and provide higher confidence in the elastic parameters used for soft tissue biomechanical studies. thin shell elements. A 40 × 40 quadrilateral element mesh was used with the tissue in x y-plane (Fig. 3). The fiber direction was predominantly in the x-direction (ω ≈ 0) and fiber dispersion τ was approximately 35 • , with small perturbations which made it more appropriate to be solved using inverse modeling rather than by curve fitting. The neo-Hookean matrix stiffness in (8) was set at λ = 172.84 kPa and µ = 74.1 kPa (corresponding to a linearized Young's modulus of 200 kPa and Poisson's ratio of 0.35). The lower and left edges were fixed in the yz and xz directions, respectively, and free to deform in the third direction (Fig. 3). For DC case, uniform displacement was applied on the right and upper edge of the tissue so as to generate a maximum strain of 13 and 8% in x-and y-direction, respectively. For FC case, a uniform force was applied on the right and upper edges in the normal direction: x-force of 0.4 N on the right edge and y-force of 0.26 N on the upper edge. The applied force or deformation was linearly increased at each quasi-time step, and system was solved for static equilibrium. The resulting normal reaction force or deformation was averaged at right and upper edges at 151 uniformly spaced steps between unloaded and maximally loaded states (both inclusive).

Details of tissue pressurization simulation
A semilunar-shaped single tissue representing a bioprosthetic valve leaflet with an area of 2.3 cm 2 and thickness of 0.38 mm was simulated using Reissner-Mindlin thin shell elements. The sample geometry was meshed using 2789 quadrilateral elements with constant thickness (Betsch et al. 1996). Displacement was interpolated using bilinear shape functions, and stresses were integrated through the thickness using three point Gauss quadrature rule to obtain bending moments. The fibers were aligned approximately in the diagonal direction with a fiber dispersion τ ≈ 35 • and relatively small perturbations over the sample. Contact of the sample with other leaflets was modeled using two idealized rigid planes placed at ±60 • with respect to the sample's plane of symmetry (Fig. 4).
The neo-Hookean matrix stiffness in (8) was set at λ = 172.84 kPa and µ = 74.1 kPa (corresponding to a linearized Young's modulus of 200 kPa and Poisson's ratio of 0.35). A uniform normal pressure was applied on the tissue with a maximum value of 120 mm of Hg. The contact was solved using augmented Lagrange method as pressure was linearly increased. At each load step, static equilibrium equations were solved to obtain the deformed shape of the tissue sample, which was recorded at four uniformly spaced points between unloaded and maximally loaded states (both inclusive). The leaflet shape was compared using 129 virtual markers, where L 2 -norm of the closest-point projection onto the deformed mesh was the objective function, as defined previously (Aggarwal and Sacks 2016).

For Ae B and A e B − 1 using 2-norm
For both functions (1) and (3), we use a common representation σ ( ) = Ah(B ). Thus, h(B ) = e B corresponds to (1) and h(B ) = e B − 1 corresponds to (3). Using 2-norm, we have If we denote Thus, the equations of two curves are similar except for the replacement of g(B) with g (B).

For Ae B using log-norm
For function (1)  If we denote α = log(A), α = log(Ā), β = B and β =B, we obtain The gradient of this functional w.r.t parameters α and β is The condition gradient ∇F = 0 gives us equations of two PM curves, which, in this case, are straight lines,
We note that solving ∂F/∂ A = 0, ∂F/∂ log(A) = 0, or ∂F/∂α = 0, all lead to the same equation. Or in other words, any of those equations could be used to determine APM curve. Therefore, for deriving the equation for APM curve, we use ∂F/∂α = 0 because of its simplicity: which is clearly a straight line in the (log(A)/B, 1/B) space with slope φ 1 /φ 3 . This curve can be expressed in the original (A, B) coordinates as: However, ∂F/∂ B = 0 and ∂F/∂β = 0 are different curves. For BPM, its equation for the (log(A)/B, 1/B) space is simply