Bayesian inference of constitutive model parameters from uncertain uniaxial experiments on murine tendons

Journal of the Mechanical Behavior of Biomedical Materials 96 (2019) 285–300

Contents lists available at ScienceDirect

Journal of the Mechanical Behavior of Biomedical Materials journal homepage: www.elsevier.com/locate/jmbbm

Bayesian inference of constitutive model parameters from uncertain uniaxial experiments on murine tendons

T

Akinjide R. Akintundea, Kristin S. Millera, Daniele E. Schiavazzib,∗ a b

Department of Biomedical Engineering, Tulane University, New Orleans, USA Department of Applied and Computational Mathematics and Statistics, University of NotreDame, Notre Dame, USA

A R T I C LE I N FO

A B S T R A C T

Keywords: Uncertainty quantification Bayesian inference Constitutive models for murine patellar tendons Parameter sensitivity and identifiability

Constitutive models for biological tissue are typically formulated as a mixture of constituents and the overall response is then assembled by superposition or compatibility. This ensures the stress response of the biological tissue to be in the range of a given constitutive relationship, guaranteeing that at least one parameter combination exists so that an experimental response can be sufficiently well captured. Another, perhaps more challenging, problem is to use constitutive models as a proxy to infer the structure/function of a biological tissue from experiments. In other words, we determine the optimal set of parameters by solving an inverse problem and use these parameters to infer the integrity of the tissue constituents. In previous studies, we focused on the mechanical stress-stretch response of the murine patellar tendon at various age and healing timepoints and solved the inverse problem using three constitutive models, i.e., the Freed-Rajagopal, Gasser-Ogden-Holzapfel and Shearer in order of increasing microstructural detail. Herein, we extend this work by adopting a Bayesian perspective on parameter estimation and implement the constitutive relations in the tulip library for uncertainty analysis, critically analyzing parameter marginals, correlations, identifiability and sensitivity. Our results show the importance of investigating the variability of parameter estimates and that results from optimization may be misleading, particularly for models with many parameters inferred from limited experimental evidence. In our study, we show that different age and healing conditions do not correspond to statistically significant separation among the Gasser-Ogden-Holzapfel and Shearer model parameters, while the phenomenological FreedRajagopal model is instead characterized by better indentifiability and parameter learning. Use of the complete experimental observations rather than averaged stress-stretch responses appears to positively constrain inference and results appear to be invariant with respect to the scaling of the experimental uncertainty.

1. Introduction Mechanical tests, reproducing relevant physiologic loading conditions are commonly employed to understand structure-function relationships in soft tissues. In this context, some commonly used testing methods are compression, shear, uniaxial and biaxial tensile tests. However, experimental data from these tests is invariably collected with some degree of uncertainty stemming from various sources, e.g., error in estimating the unloaded or reference configuration of the tissue, accuracy of the chosen deformation tracking method (i.e. optical vs. grip-to-grip tracking), the precision level of the load cell as well as biological variability. Constitutive models describe the response of materials to loads under specific conditions (Humphrey and DeLange, 2013), and include parameters that are trained from experimental data, in order to investigate structure-function relationships in the tissue of

∗

interest. These models are then used in mathematical modeling and computational simulations of growth and remodeling (Miller et al., 2014), material failure phenomena (Ionescu et al., 2005) and for patient-specific simulations with the aim of guiding clinical decisions (Schiavazzi et al., 2017). Often, particularly for complex constitutive models, parameter training is achieved through optimization, i.e., without explicitly characterizing the variability associated with parameter estimates, and also a thorough study of the parameter identifiability and sensitivity is not performed, limiting the translational relevance of these approaches and delaying their clinical adoption. In this context, transition to robust approaches for parameter estimation under uncertainty will lead to more dependable analyses and simulations, thereby enhancing clinical relevance. Recent studies in the literature have also used robust Bayesian approaches to increase physiological and clinical relevance of

Corresponding author. 101G Crowley Hall, Notre Dame, IN 46556, USA. E-mail address: [email protected] (D.E. Schiavazzi).

https://doi.org/10.1016/j.jmbbm.2019.04.037 Received 4 February 2019; Received in revised form 20 March 2019; Accepted 18 April 2019 Available online 30 April 2019 1751-6161/ © 2019 Elsevier Ltd. All rights reserved.

Journal of the Mechanical Behavior of Biomedical Materials 96 (2019) 285–300

A.R. Akintunde, et al.

soft tissue models (Schiavazzi et al., 2017; Seyedsalehi et al., 2015; Sankaran et al., 2013; Madireddy et al., 2016; Doraiswamy and Srinivasa, 2013; Ritto and Nunes, 2015). In this study, we discuss the problem of identifying the parameters for constitutive models of the murine patellar tendon, under uncertainty due to random experimental and modeling errors (Madireddy et al., 2016). In prior work, three constitutive models with varied levels of microstructural detail and a varying number of parameters were considered (Akintunde and Miller, 2018). Herein, we extend this work using a statistical model with Gaussian additive errors, thus leveraging adaptive Markov chain Monte Carlo (Vrugt et al., 2009) to sample from the posterior distribution of the constitutive parameters for each model. We investigate the conditions facilitating learning for these parameters, examine how age and injury affect their distributions, inspect their correlations and study their identifiability and sensitivities. It may be argued that an additional form of uncertainty relates to partial information provided by the nature of the collected experimental data. In other words, the availability of a complete set of multiaxial experiments would significantly improve the identifiability of even very complex microstructurally motivated constitutive models. On the other hand, we don't have data other than experimental results from uniaxial tests (as typically the case), and therefore we ask ourselves about the best parameter estimates we can provide under these specific settings. Thus, the conclusions we draw from our study are to be understood in light of the constitutive models we selected and experimental evidence discussed in Section 2.1. For patient-agnostic applications, it is of interest to determine some sort of average behavior of the soft tissues. However, the wide interspecimen variability commonly observed leads to a lack of consensus on the optimal approach to use for computing this average and a debate exists between fitting model to each data sample and then computing average of the model parameters versus fitting model to the average experimental data (Madireddy et al., 2016; Bellini and Di Martino, 2012; Ferruzzi et al., 2013; Robertson and Cook, 2014). In our study, we address this issue and compare model training by assembling all experimental data in the likelihood or just the average stress levels. The paper is organized as follows. Section 2.1 discusses the experimental dataset, while the mathematical formulation of the three constitutive models selected for this study is explained in Section 2.2. The prior ranges for each parameter are justified in Section 2.3. Section 2.4 explains how the parameter estimation problem is formulated based on a Bayesian paradigm, while Sections 2.5 and 2.6 focus on the parameter identifiability and global sensitivity, respectively. Inference results are presented in Section 3 and critically analyzed in Section 4.

Table 1 Number of uniaxial experiments repetitions by age group and injury time points. Age

0 week (uninjured) 3 weeks 6 weeks

Mature (120 days)

Aging (270 days)

Aged (540 days)

5 17 14

14 15 15

20 9 11

physiologic loading, and up to 8% under acute stress. 2.2. Constitutive models Collagen in tendons has a unique hierarchical arrangement (Kannus, 2000). Intact tendons mostly consist of uniaxially organized collagen fascicles, i.e., bundles of collagen fibers formed by aggregation of wavy collagen fibrils (Fratzl, 2008). Historically, most models employed to characterize the response of tendons and ligaments to mechanical loading are phenomenological, i.e., providing a purely mathematical description of the nonlinear stress-stretch response (Holzapfel et al., 2000; Weiss et al., 1996), but offering little to no insight on tissue level mechanisms governing the response to external loads. In addition, the parameters of such models vary widely, and are therefore prone to a subjective interpretation. In our prior work, three constitutive models with varied level of microstructural motivation were selected for study (Akintunde and Miller, 2018). The Shearer (SHR) constitutive model has the greatest microstructural detail, describing tendons at both the fibrillar and fascicular levels (Shearer, 2015). The Gasser-Ogden-Holzapfel (GOH) model contains some phenomenological parameters, but uniquely captures the dispersion of collagen fibers and their alignment about an identified direction (Gasser et al., 2006). Finally, the simplest model is the Freed-Rajagopal (FR) model for tissue fibers, which contains three interpretable parameters (Freed and Rajagopal, 2016). 2.2.1. Hyperelastic formulation for the SHR and GOH models The patellar tendon was assumed to be incompressible and with transversely isotropic solid cylindrical geometry (Vergari et al., 2011; Dourte et al., 2013). Hence, for uniaxial tensile test, the deformation gradient tensor F and corresponding right Cauchy-Green tensor C are: −1/2 0 0⎤ ⎡λ F=⎢ 0 λ−1/2 0 ⎥, ⎢ 0 0 λ⎥ ⎣ ⎦

2. Materials and methods

−1

⎡λ C=⎢ 0 ⎢ 0 ⎣

0 0⎤ λ−1 0 ⎥, 0 λ2 ⎥ ⎦

(1)

where λ is the stretch measured along the z axis of the tendon during the uniaxial extension test. The first and fourth strain invariants are → → I1 = tr(C) = 2 λ−1 + λ2 and I4 = M C M . The forth invariant is also exSHR I4 = λ−1 sin2 ψ + λ2 cos2 ψ I4GOH = and pressed as − 1 2 2 2 λ sin ξ + λ cos ξ , for the SHR and GOH models, respectively. The → vector M represents the preferred collagen fascicle/fiber direction in the reference configuration, measured from the z axis. In cylindrical →SHR →GOH coordinates this is M = [sinψ, 0, cosψ] and M = [sinξ , 0, cosξ ], where ψ is the fascicle alignment angle, similar to ξ in the GOH model, i.e., the fiber alignment angle, both measured from z and determined through imaging, if available. The hyperelastic framework for transverse isotropy (Ogden, 2003) is employed, with the Cauchy stress tensor obtained as a function of the strain energy-density function (SEF), W = W (I1, I4 ) :

2.1. Experimental data Stress-stretch data for murine patellar tendon at various aging and healing timepoints before and following injury are gathered from studies in the literature (Dunkman et al., 2013, 2014a, 2014b; Mienaltowski et al., 2016). These studies followed an established excisional injury model to the mid-substance of the tendon (Beason et al., 2012; Lin et al., 2006), with their experimental protocols approved by the University of Pennsylvania Institutional Animal Care and Use Committee. Pre-(0 week) and post-(3 and 6 weeks) injury, real-time load-displacement and the undeformed cross-sectional area data were obtained from three age groups of female C57BL/6 mice: 120, 270, and 540 day-olds which were considered as mature, aging and aged, respectively (Table 1). The nominal stress is computed from the load and undeformed cross-sectional area and converted to the true stress via push forward (Taber, 2004; Humphrey, 2013). A quasi-static (strain rate = 0.1%/sec) uniaxial ramp-to-failure test was performed with a 10% stretch cut-off (Józsa, 1997; Kannus, 2000), justified considering 1–4% elongation experienced by healthy tendons under normal

σ = −p I + 2

→ ∂W ∂W → b+2 FM ⊗ FM , ∂I1 ∂I4

(2)

where b is the left Cauchy-Green strain tensor. The SHR SEF yields the following expression for the axial Cauchy stress, after the Lagrange 286

Journal of the Mechanical Behavior of Biomedical Materials 96 (2019) 285–300

A.R. Akintunde, et al.

multiplier (p) is found by applying a no-traction boundary condition at the tendon's outer radial surface within the equilibrium equation → ∇⋅σ = 0 ,

sin2ψ ⎞ 1 SHR , σzz = (1 − ϕ) μ ⎛λ2 − ⎞ + 2 χ SHR ⎛λ2cos2ψ − 2λ ⎠ λ ⎝ ⎠ ⎝ ⎜

⎟

(3)

where

χ SHR

0, ⎧ ⎪ ϕE ⎪ ⎡2 − = 6 I4 sin2θ0 ⎣ ⎨ ϕE ⎪ β− ⎪ 2 I4 ⎩ 3 1 − cos θ0 =2 , 3sin2θ0

(

I4 < 1

(

3 I4 − 1 I4 I4 1 I4

) ⎤⎦, 1 ≤ I

4

),

I4 >

≤ λ*

and β

λ*

(4) Fig. 1. Relation between the Shearer transition stretch λ* and a transition stretch based on the stress-stretch curve derivative. For an increasing alignment angle ψ, the large transition stretch affects the extension of the linear portion of the stress-stretch curve.

where ϕ ∈ (0,1] is the collagen fiber volume fraction, μ is the shear modulus of the ground matrix which consists of non-collagenous substances such as elastin, glycosamino-glycans, and proteoglycans, E is the fibril Young's modulus, θ0 is the crimp angle of the outermost fibrils and λ* =

1 cosα

1 cos2θ0

− sin2α ≈ 1/cosθ0 is the stretch along the fascicle

2.2.2. The Freed-Rajagopal (FR) model The FR model splits the total tissue strain in the contribution of the elastic filament network and collagen fibers, respectively. While the collagen fibers are treated as neo-Hookean, the elastin filament network serves as strain limiter, with mechanical response obtained from an implicit SEF (Freed and Rajagopal, 2016). The three-parameter stretch as a function of stress FR model is

direction that straightens the outer fibrils, i.e., the critical or transition stretch. Note also how the expression (4) assumes a zero fibril helix angle α for the murine patellar tendon, as justified by preliminary computational studies on the same experimental dataset, and also alluded to by Shearer for the patellar tendon in humans (Shearer, 2015). Likewise, for the GOH model, the following is obtained for the axial Cauchy stress GOH σzz

1 = (1 − ϕ) μ ⎛λ2 − ⎞ + χ GOH λ⎠ ⎝

⎧ 1 ⎪ 1 λ = λ e λ c , ε = ln λ, ε = εe + εc = 1− β ⎨ ⎪ ⎡1 + (β − 1) ⎣ ⎩ σzzexp , + Ec

⎡κ ⎛λ2 − 1 ⎞ ⎢ ⎝ λ⎠ ⎣

sin2ξ

⎞ ⎤, + (1 − 3 κ ) ⎛λ2cos2ξ − 2 λ ⎠⎥ ⎝ ⎦ ⎜

⎟

(5)

where

⎫ ⎪ β exp σzz ⎤β − 1 Ee

⎦

⎬ ⎪ ⎭ (6)

where Ee (> 0) and Ec (> Ee) are the elastic moduli of the toe and linear regions - dominated by elastin (e) and collagen fibers (c), respectively and β = 1/ε max (> 1) is the inverse of the true strain at which the linear e region starts, and the collagen fibers are straight.

2

χ GOH = 4 ϕ c1 ν e c2 ν , ν = κ I1 + (1 − 3 κ ) I4 − 1. In (5), κ ε [0, 1/3] quantifies the collagen fiber dispersion about →GOH M . It is zero when the fibers are perfectly aligned and 1/3 for an isotropic distribution. The quantity c1 > 0 is a modulus-like parameter, while c2 > 0 is a dimensionless parameter. These two parameters provide a phenomenological description of the structural response of the collagen fibers, but their physical significance is unclear and subject to personal interpretation as neither parameter can be measured experimentally.

2.3. Parameter bounds To ensure physical validity and motivate physiologic relevance, bounds were placed on all model parameters (Table 2). For structural model parameters, all angles were constrained between 0 and π/2 and collagen dispersion parameter κ from 0 to 1/3. Material parameters that

Remark. (Fascicle and fibril alignment angle parameters ψ and ξ). The → formulation discussed in this section was derived for direction vectors M , not necessarily aligned with the longitudinal tendon axis z. The alignment angle parameters thus introduced interact with the other model parameters, for example in quantifying the transition between an elastin-dominated lowstretch response and the collagen-dominated linear response. To see this, we consider the relation between the transition stretch λ* = 1/cosθ0 in the SHR model and the beginning of the linear response region evaluated through the change in the stress-stretch curve derivative. This is shown in Fig. 1. For small alignment angles, the transition stretch λ* fully correlates with the beginning of the linear region, as expected. For increasingly large alignment angles, the linear region appears to be confined close to the 10% stretch limit and to be independent from the crimp angle parameter θ0 . This demonstrates how the crimp and alignment angles interact in determining the transition stretch λ* , and how the identifiability of (say) θ0 is progressively reduced for an increasing alignment angle ψ. This also motivates an investigation on how the alignment angles affect the inferred GOH and SHR model parameters. As discussed in Section 3.1.2, we repeated Bayesian inference by including and excluding ξ and ψ.

Table 2 Theoretically and experimentally motivated bounds on GOH, SHR and FR model parameters and associated priors. Adapted from Table 1 in (Akintunde and Miller, 2018).

Material

Structural

287

Model

Parameter

Constraint

References

GOH & SHR GOH GOH SHR

(1 − ϕ) μ

>0

[0, 108]

c1 c2 ϕE

>0 >0 0–14 GPa

[0, 108] [0, 103] [0, 108]

FR FR GOH

Ec Ee κ

>0 >0 0 - 1/3

GOH SHR SHR FR

ξ θ0 ψ β

0 - π/2 0 - π/2 0 - π/2 10–1000

Andriotis et al. (2015)

Gasser et al. (2006) Shearer (2015)

Prior range

[10−3, 108] [10−3, 108] [0, 1/3] [0, π/2] [10−3, π/2] [0, π/2] [101, 103]

Journal of the Mechanical Behavior of Biomedical Materials 96 (2019) 285–300

A.R. Akintunde, et al.

et al., 2005) and justifies our adoption of a Bayesian approach to inference. Due to the non linearity of the three constitutive models in θ, the likelihood (8) is generally not normally distributed, and generation of samples from the stationary posterior (9) through random walk Metropolis may require a very large number of forward evaluation of the selected constitutive models. Differential evolution adaptive Metropolis (DREAM) (Vrugt et al., 2009; Vrugt, 2016) combines differential evolution (Storn and Price, 1997) and self adaptive randomized subspace sampling, leveraging parallel Markov chains to automatically adapt the scaling and orientation of the trial distribution. It also shows a substantial improvement with respect to other adaptive MCMC schemes when sampling from heavy tailed or multi-modal posteriors.

cannot be measured experimentally were simply required to be positive. The collagen fibril elastic modulus was constrained based on available experimental data. A detailed discussion on the chosen bounds can be found in (Akintunde and Miller, 2018). Additionally, uninformative uniform priors are assumed for all parameters, with ranges also reported in Table 2. 2.4. Bayesian parameter estimation We denote all model parameters with the vector θ = (θ1, θ2, …, θd ) ∈ d and embrace a Bayesian perspective on parameter estimation, treating θ as a vector of random variables. Specifically, we consider a complete probability triplet (Ω, , P ) equipped with a space of elementary events ω ∈ Ω, a Borel σ-algebra  of events defined over 2Ω and a probability measure P assigning values in [0,1] to elements in  . The arbitrarily distributed and correlated parameters θi (ω): Ω→ Θi ⊂  , i = 1, …, d are characterized in probability through their joint distribution ρ (θ ) .

2.4.3. Convergence assessment Convergence, i.e., the ability of the adaptive Markov chain to generate samples from the desired posterior distribution is assessed, in this study, through the Gelman-Rubin diagnostic metric (RGR ) (Gelman and Rubin, 1992), but several other approaches are available through the literature. GR evaluates the similarity of m parallel chains with n sample each based on an analysis of variance. Consider the j-th sample from the i-th parameter generated from the k-th Markov chain θi(j, k ) , and the quantities

2.4.1. Statistical model for uniaxial stress and stretch The experimental uniaxial stress (stretch) in the z-direction for the GOH, SHR and FR material models is considered as a random variable with distribution induced by the following statistical model

σzz (λ ) = σzzM (λ, θ ) + ε, or λ (σzzM ) = λ zzM (σzzM , θ ) + ε,

(7)

where M = {GOH,SHR,FR} and ε ∼ N (0,  [ε ]) is an error term assumed Gaussian, with variance  [ε] and standard deviation σε . The model (7) suggests σzz (λ ) ∼ N (σzzM (λ, θ ),  [ε]) with likelihood function expressed as

ℓθ [σzz (λ )] =

1 2π  [ε]

⎛ 1 exp − ⎜ 2 ⎝

S

N

∑∑

[σzz(i)

(λj ) −

σzzM (λ,

 [ε]

i=1 j=1

1 k θ¯i = n

θiinter =

(8)

where σzz(i) (λj ) represents the i-th experimental result available under the j-th measured axial stretch λj for a single age/healing group (see Table 1), S is the total number of samples per group and N is the total number of stretch levels. For convenience, the vector λ = [λ1, λ2, …, λN ] will contain all experimental stretch levels, while the tensor σzz (λ ) denotes all experimental stresses at all stretch levels within a single agehealing group. Note how the additive nature of the noise term in the statistical models (7) essentially determines the expression (8) for the likelihood. While we do not consider here alternative likelihood formulations, the results in Section 3.2 suggest σε to be an unimportant parameter that does not affect significantly the inference process. The effect of systematic errors and more complicated error distributions can be taken into account through the selection of an appropriate likelihood function.

=

P (σzz |θ, λ, M ) P (θ )

∫Θ P (σzz |θ, λ, M ) P (θ) dθ

m

∑

(j, k ) θ¯i .

(10)

k=1

n m−1 1 m

m

∑ k=1

m

∑

2 k (θ¯i − θ¯i ) , and θiintra

k=1

1 n−1

n

∑

k 2 (θi(j, k ) − θ¯i ) .

(11)

j=1

Based on (10) and (11), an estimate of the parameter marginal posterior variance and a convergence indicator RGR can be computed as

 [θi |σzz , M ] ≈

1 n − 1 intra θi + θiinter , RGR = n n

 [θi |σzz , M ] . θiintra

(12)

2.4.4. MAP estimate A maximum a posteriori (MAP) estimate of the model parameters θ is determined by optimization. Specifically, we employ a restarted Nelder-Mead (NM) algorithm (Nelder and Mead, 1965) and set the initial guess equal to the MCMC sample with maximum posterior. Additionally, we restart the NM iteration several times, according to references in the literature showing how this consistently improves convergence (McKinnon, 1998).

2.4.2. Adaptive Markov Chain Monte Carlo The joint posterior distribution for θ is obtained through Bayesian conjunction of likelihood and prior (see e.g., Tarantola, 2005) as:

P (θ|σzz , λ, M ) =

j=1

1 θi(j, k ) , and θ¯i = m

The inter-chain and intra-chain variance for parameter θi are expressed as

2

θ )] ⎞ , ⎟ ⎠

n

∑

, (9)

2.5. Local parameter identifiability

and, in practice we use Markov chain Monte Carlo (MCMC) (Gilks and Roberts, 1996; Jackman, 2000) to generate independent samples from (9). Complex constitutive models with many parameters, inferred from limited experimental evidence (e.g., uniaxial stress-stretch tests for a multiaxial constitutive model) may present significant negative correlations among parameters or, in other words, a lack of identifiability, meaning that multiple combinations of parameters could generate outputs matching the experimental data equally well. We will see in the result sections how the GOH and SHR models are affected by this phenomenon, particularly when the alignment angles ξ and ψ are included in estimation. Specification of a prior for the parameters can substantially improve parameter learning and identifiability (Gustafson

In an identifiable system, the parameters can be uniquely determined from the available observations. If this happens in a neighborhood of a specific location θ the system is locally identifiable and globally identifiable if this happens over the whole parameter space Θ. A review of identifiability for nonlinear ODE models in viral dynamics is provided, for example, in Miao et al. (2011). In this study, we use several metrics to detect non identifiable combinations of parameters using both a frequentist and Bayesian perspectives. First, we investigate local identifiability using the rank of the Fisher information matrix I (θ ) (see e.g., Rothenberg, 1971). As a direct result of assuming the likelihood in equation (8), we have: 288

Journal of the Mechanical Behavior of Biomedical Materials 96 (2019) 285–300

A.R. Akintunde, et al. T

d

∂ σ M (λ, θ ) ⎤ ⎡ 1 ∂ σzzM (λ, θ ) ⎤ I (θ ) = ⎡ zz . ⋅I ⎤ ⎡ ⎢ ⎥ ⎢  [ε] ⎥ ⎢ ⎥ ∂θ ∂θ ⎦⎣ ⎣ ⎦ ⎣ ⎦

 [σzzM (θ˜)] =

(13)

 [θj |σzz (λ )]/  [θj] ,

Si =

∼ σzzM, i (θi ) +

i=1

∑

θi [ θ∼ i [σzzM |θi]] ≈

r

(15)

∼ σzzM, i1, i2, … , ir (θ˜i1, θ˜i2, …, θir ) dθ˜i1 dθ˜i2…dθ˜ir = 0, ∀ 1 ≤ i1 < i2 < …

[0,1]

(16)

σzzM (θ˜)

is decomposed into zero-mean contributions or, in other words, of increasing dimensionality (Sobol, 2001), and, as a consequence  [σzzM (θ )] = σzzM,0 . If we elevate (15) to the second power and integrate over the d-dimensional hypercube, we obtain 2

2

∫ [σzzM (θ˜)]2 dθ˜ − [σzzM,0] = ∫ ∑id= 1 [σzzM, i (θ˜i)] dθ˜ + ⋯ d

1 2N

N

∑ j = 1 [f (A) j − f (ABi) j]2 .

(21)

All computations discussed in the previous sections are carried out using the tulip library for uncertainty analysis. Tulip is a C++ library that provides a Python API layer and implements a complete separation between models, actions and data. Models are generic black box relations between inputs and outputs and can be used to wrap the evaluation of constitutive models, physics-based solver, particularly lumped parameter representations of circulation networks. Actions are generic steps in uncertainty quantification, such as forward propagation with Monte Carlo simulation or stochastic spectral expansion, optimization, local and global sensitivity analysis, local identifiability analysis, adaptive Markov chain Monte Carlo sampling and others. Finally, likelihood or cost functions can be evaluated by combining model outputs with data from various sources, e.g., CSV files, XML files or remote repositories. The library is designed to be extensible and new models, actions or data sources can be added by implementation of virtual functions.

σzzM,1,2, … , d

< ir ≤ d,

N

∑ j = 1 f (B) j [f (ABi) j − f (A) j],

2.7. The tulip software library

is a constant and and are function of one, where two and d parameters respectively, subject to the condition

∫

1 N

The interested reader is referred to the discussions in Saltelli (2002), Owen (2013) and Janon et al. (2014) for further details on Monte Carlo estimators for Sobol’ sensitivity indices.

∼ ∼ σzzM, i, j (θi , θj )

σzzM, i, j

(20)

 θ∼ i [θi [σzzM |θ∼ i]] ≈

∼ ∼ ∼ + ⋯+σzzM,1,2, … , d ( θ1, θ2, …, θd ),

σzzM, i ,

Ditot ,  [σzzM (θ˜)]

all contributions to the variance of terms having at summing in least one index equal to i. Note how, in our case, the parameter realizations are available through arbitrary correlated and distributed samples from MCMC. Computation of the sensitivity indices is accomplished by first creating a polynomial chaos surrogate of the quantity of interest using regression. Second, the direct and total sensitivities are determined through Monte Carlo estimates. Consider a sample matrix of size N × 2d , where each row is a sample of dimension 2d, generated from the probability distributions of the random inputs. Denote with A the first d columns of this sample matrix, and the remaining d columns as B ∈ N × d . This effectively gives two independent samples of the d-dimensional random inputs at N locations. Consider d additional N × d matrices ABi , i = 1,2, …, d , such that the i-th column of ABi is equal to the i-th column of B, and the remaining columns are from A . Now, each constitutive model is evaluated at the N ⋅(d + 2) parameter realizations stored in the matrices A , B, and ABi , i = 1, …, d and the sensitivity indices (direct and total, respectively) are computed from the estimators

i
σzzM,0

(19)

Ditot

Consider the stress σzzM (λj , θ˜) , i.e., the stress response of a generic constitutive model M at a given stretch level λj (the dependence on λj will be omitted below for simplicity) as a function of rescaled inputs θ˜ defined over the hypercube [0,1]d . Under these circumstances, the stress response can be decomposed as

∑

(18)

Di ,  [σzzM (θ˜)]

Sitot =

2.6. Variance based global sensitivity analysis

σzzM (θ˜) = σzzM,0 +

Di, j + ⋯+D1,2,3, … , d .

i
quantifying the contribution to the variance of the i-th random input θ˜i alone. Similarly, total sensitivity indices can be computed as

(14)

d

∑

Following this approach, direct, first order global sensitivity indices (i.e., effect of single random inputs) can be determined as

where  [θj |σzz (λ )] denotes the marginal posterior variance for parameter θj while  [θj] denotes its prior marginal variance.

d

d

Di +

i=1

Note that Equation (13) is easily obtained by combining the matrix ∂ σzzM (λ, θ )/ ∂θ of local derivatives (e.g., computed through finite differences) with the diagonal precision matrix 1/  [ε]⋅I of the experimental stress (GOH, SHR) or stretch (FR). As discussed in the literature (Rothenberg, 1971), rank deficiency in I (θ ) reveals the presence of non-identifiable combinations of parameters, and null eigenvectors can be computed to identify such combinations. In practice, the numerical eigenvalues could be small but not negligible and therefore it may be difficult to identify parameter combinations lacking identifiability, particularly if these involve a large number of parameters. On the other hand, it is rather easy to identify locally unimportant parameters, where the eigenvalue is very small (e.g., < 10−10 ) and the eigenvector is dominated by a single component. Also, changes in the rank deficiency of I (θ ) as a result of applying different scaling factors to the experimental precision matrix are useful to recognize true singularities. Other identifiability metrics are obtained as a byproduct of Bayesian inference, for example, the parameter correlations identified from the MCMC stationary posterior samples. Large negative correlations between two parameters may be indicative of unidentifiability, i.e., stress changes compensate in the model response. The shape of the marginal posterior distribution may also be indicative of the importance of a given parameter. Specifically, a flat posterior distribution typically indicates limited importance. Finally, reduction in marginal variance from prior to posterior is an indication of how much a certain parameter was learned upon conditioning on the available data. In this context, a learning factor γj for the j-th parameter was introduced in (Schiavazzi et al., 2017) as

γj = 1 −

∑

d

[0,1]

[0,1]

∼ ∼ 2 d dθ˜ + ∫ ∑i < j [σzzM, i, j (θi , θj )]

3. Results

[0,1]d

2 … + [σzzM,1,2, … , d (θ˜1, θ˜2, …, θ˜d )] ,

3.1. Bayesian inference

(17)

Results obtained from Bayesian estimation through adaptive MCMC are presented next. We investigate the convergence of the parallel

suggesting how the variance can be decomposed as a sum of contributions with increasing dimensionality, i.e. 289

Journal of the Mechanical Behavior of Biomedical Materials 96 (2019) 285–300

A.R. Akintunde, et al.

Fig. 2. Mean value and 5–95% percentiles for Gelman-Rubin diagnostic across all parameters and age-healing groups.

understand the variability associated to the parameter estimates, we report the marginal posterior median together with the 10–90% confidence intervals. Due to the higher learning ratios consistently provided, Figs. 5–9 were produced by training the models using the complete stress dataset for each age-injury group. They show how evaluation of parametric uncertainty from Bayesian estimation is essential to be able to identify changes in the tendon ligament microstructure from the uniaxial stress-stretch response. In this regard, the GOH and SHR models do provide a more complete description of the microstructure but at the expense of large uncertainties in the estimated parameters. Conversely, the FR model is characterized by only three parameters each associated to a well defined stretch range (i.e. lowstretch elastin response modulus Ee , transition parameter β and collagen response parameter Ec ) and considerably smaller uncertainties, as shown in Fig. 9. Figs. 5–9 confirm that microstructural parameters, whose median estimate may look different for various age/healing groups, cannot be distinguished under uncertainty. Note how classical parameter fitting through optimization might have led to the opposite conclusion, as a results of the inability to quantify the confidence in these estimates. This suggests how measures of confidence are essential for the interpretability of the results, e.g., when it is of interest to quantify function from the identified micro-structural model parameters. This also shows how the complexity of the underlying constitutive model should play a role in the interpretation of the results, where the presence of non identifiable parameter combinations may contribute substantially to the marginal parameter variance. For the GOH and SHR models, several factors contribute to the parameter marginal variance for each age/healing group. First, the aged group appears to be associated with the largest variability in the parameter estimates. This is expected due to the increasing difficulties to produce experimental data at 540 days. Additionally, the marginal variance appears to change for the same age groups and across healing conditions. As reported in Table 1, the number of uniaxial experiment repetitions may vary sensibly with the healing condition and surely affects the marginal variance. Finally, a significant contribution to the parameter marginal variance is due to non-identifability, as confirmed by the large negative correlations from MCMC, shown in Section 3.1.4. We conclude this section by reporting the MAP parameter estimates for the three models in Tables 3–5. The most interesting behavior relates to the GOH and SHR models, specifically how the MAP estimate changes by setting the alignment angle to zero. While the estimates for the (1 − ϕ) μ , ϕ c1, c2 and ϕ E parameters seem to be consistent, more noticeable changes are observed in κ and θ0 . In the GOH model, the estimated ξ angle is considerable (33.7–49.3°). When ξ = 0 the fiber dispersion (κ) become more isotropic to reduce the stress response that has increased due to a perfect alignment of the fibers. A similar trend is observed for the SHR model between the two parameters θ0 and ψ, but with smaller estimated fascicle alignment angle.

chains, the marginal distribution of each θi , i = 1, …, d , the uncertainty in the identified stress-stretch constitutive relations, the correlations among parameters. 3.1.1. Convergence from parallel inference We look at the convergence diagnostic discussed in Section 2.4.3 for all parameters, all models and all aging-healing groups. The number of parallel Markov chains was selected equal to the number of parameters and 30,000 generations (i.e., model evaluations) were performed for each chain. Fig. 2 shows how the 5–95% percentiles for the GelmanRubin statistic (RGR ) fall below an acceptance threshold of 1.1 in all cases. Moreover, the model complexity (i.e., the number of constitutive parameters) clearly affects how many posterior evaluations are needed for convergence, with the FR model showing the best convergence followed by SHR and GOH. Similar convergence plots (not shown) are obtained by introducing the alignment angles ξ and ψ in the GOH and SHR model, respectively. 3.1.2. Marginal parameter learning and distributions To investigate the ability to identify constitutive model parameters from experimental data with a sufficiently large confidence, we evaluated the learning factors (14) for all models and parameters. We also differentiate cases where the complete experimental dataset was included in the likelihood (denoted by ALL in Figs. 3 and 10) rather than just average stress levels (denoted by AVG). For the GOH model, the only parameter showing satisfactory variance reduction relates to the ground matrix response parameter (1 − ϕ) μ , which is the only load-carrying mechanism occurring at a sufficiently low stretch, i.e., before full engagement of collagen fibers. The maximum learning factors for parameters associated with the collagen fiber response remain instead limited to approximately 0.2, indicating correlation among parameters and a reduced identifiability. Additionally, larger learning factors are obtained for the GOH model when the complete experimental dataset is included. Finally, a zero alignment angle ξ significantly improves learning of the fiber dispersion parameter κ. The SHR model shows improved learning factors in all parameters due to a reduced model complexity, but also highlights the benefits of including the complete set of experimental stresses. The learning factor at 540 days age and 3 weeks healing appears associated with a particularly small learning factor even after including the complete experimental dataset during training (not shown). This can be attributed to two factors, first a relatively small number of experimental data points for this specific age-healing group (see Table 1) and second to the larger impact of injury for this age group. For the same stretch range, the stresses are lower but also the variability in the stress response is reduced as shown in Fig. 4. Therefore, for this experimental group there is little difference between training using average stresses or the complete experimental dataset. Finally, the parameters of the FR constitutive model appear to be equally well learned in all cases, with a slightly reduced learning of the collagen elastic modulus Ec from all stress values. Next, we look at the changes in the parameter marginal median for each model, following injury in various age groups. To better

3.1.3. Optimal and MAP parameter estimates and stress-stretch response Lack of identifiability has also a deleterious effects on the ability to gather microstructurally relevant parameters from experiments through optimization. To show this, we have determined the parameters using 290

Journal of the Mechanical Behavior of Biomedical Materials 96 (2019) 285–300

A.R. Akintunde, et al.

Fig. 3. Learning factors for GOH, SHR (arbitrary and zero alignment angles) and FR.

checked the uncertainty in the stress-stretch response. To do so, we have separated the stress-stretch curves based on age and showed how the response changes with the healing conditions. We also superimposed on the same figure the experimental data, the stress-stretch response corresponding to the MAP parameter estimate, the associated 5–95% percentile variation and the results obtained by straight optimization without taking experimental uncertainty into consideration. The experimental stresses are characterized by significant overlap following 6 weeks injury recovery, as shown in Figs. 11–13, leading to very similar constitutive responses. Stress response remain similar after

both the Nelder-Mead optimization algorithm (see e.g., Nelder and Mead, 1965) and by identifying the maximum of the posterior distribution from the MCMC samples. The results in Fig. 10 show how the parameters are inconsistent in the GOH and SHR models while perfect agreement is observed for FR. This suggests a lack of convexity in the posterior distribution, where either the selection of different initial guesses leads to convergence to a different local maxima or the optimizer encounters a region of flat posterior. This is indeed a manifestation of the lack of identifiability. After verifying the agreement between estimated parameters, we 291

Journal of the Mechanical Behavior of Biomedical Materials 96 (2019) 285–300

A.R. Akintunde, et al.

Fig. 4. Number of experimental data points at various stretch levels, for the same healing group (3 weeks) at age 120 days (a), 270 days (b) and 540 days (c).

Fig. 5. Changes in the parameter median value and 10–90% percentile due to injury and age for the GOH model with ξ ≠ 0 .

Fig. 6. Changes in the parameter median value and 10–90% percentile due to injury and age for the GOH model with ξ = 0 .

parameters. The strong negative correlation between the fiber dispersion parameter κ and the alignment parameter ξ is particularly evident. This simply suggests that unaltered experimental stresses can be obtained, for example, increasing the degree of random orientation of the fibers and reducing the alignment angle, and confirms the difficulty in discerning between these two phenomena from uniaxial experiments. Finally, we observe a recurrent block structure in the correlation matrix for the GOH and SHR models, where the two blocks correspond to the low and high-stretch response, dominated by the ground matrix and

setting the alignment angles ξ and ψ to zero. 3.1.4. Parameter Correlations The average parameter correlations over all age-healing groups are shown in Fig. 14, from models trained using both the complete set of experimental stresses or only their average levels. It appears clear from this plots how the availability of complete experimental information accentuates the negative correlations for the GOH and SHR models, while confirming the lack of correlation among FR constitutive 292

Journal of the Mechanical Behavior of Biomedical Materials 96 (2019) 285–300

A.R. Akintunde, et al.

Fig. 7. Changes in the parameter median value and 10–90% percentile due to injury and age for the SHR model with ψ ≠ 0 . Fig. 8. Changes in the parameter median value and 10–90% percentile due to injury and age for the SHR model with ψ = 0.

Fig. 9. Changes in the parameter median value and 10–90% percentile due to injury and age for the FR model. 293

Journal of the Mechanical Behavior of Biomedical Materials 96 (2019) 285–300

A.R. Akintunde, et al.

Fig. 10. GOH, SHR and FR representative model parameters computed through optimization and maximum a posteriori estimates from MCMC.

collagen fibers, respectively.

Due to the scaling provided by the inverse variance from the experiments in (13), definition of a sharp cutoff for zero eigenvalues may result in variability in the identifiability results. To overcome this difficulty, we amplified the experimental precisions with factors in the range [10−2, 103] and observed the corresponding change in the FIM eigenvalues. For the GOH and SHR models with zero alignment angle and FR, all eigenvalues are affected by the scaling suggesting these models to be locally identifiable (i.e. associated with a non singular FIM around the MAP parameter estimates), while zero eigenvalues remain for the GOH and SHR models with non-zero alignment angles. We then set a cutoff equal to 10−10 allowing for the selection of eigenvalues remaining small even after significantly amplifying the experimental precisions. As expected, Fig. 17a suggests a sensible increase in the eigenvalues as a result of setting the alignment angles ξ and ψ to zero. Additionally, the analysis confirms how the standard deviation parameter is locally unimportant as suggested by the dominant null-eigenvectors in Fig. 17b and c. Fig. 17b suggests instead how the GOH parameters ϕ c1 and (1 − ϕ) μ form an unidentifiable combination. This is consistent with the large negative correlation shown by these two parameters in Figs. 14 and 15, with two possible constituents able to withstand the low stretch response, i.e., elastin matrix or collagen fibers. Finally, Fig. 17c suggests the fibril elastic modulus ϕ E to be unimportant for the SHR model with non-zero alignment angles. This can be explained by looking at Fig. 1 and the relatively large MAP estimates for ψ. For large alignment angles the transition to the linear region of the stretchstress curve occurs at very large stretches and therefore the importance of the collagen modulus ϕ E remains limited.

3.2. Including data variability in inference Instead of simply assuming a deterministic value for σε , we performed numerical tests to understand how specification of this quantity affects the results from parameter estimation. Specifically, we assumed an experimental uncertainty in the GOH and SHR models equal to 106 σε (i.e., the experimental uncertainty is expressed in MPa). Stresses for these two models are measured in Pa and therefore the final uncertainty is assumed uniform over the whole stretch range and equal to 0.1–10.0 MPa. Since the FR model has stretch outputs, the uncertainty was assumed instead as 0.01 σε . The marginal posterior of these parameters was found completely flat as shown for the SHR model in Fig. 16, suggesting how this parameter can be set to any reasonable value and it does not affect parameter identification. The parameter correlations in Fig. 14 also show how the experimental uncertainty is completely uncorrelated with respect to all other parameters for the three models. 3.3. Local identifiability The identifiability of the three models has been formally analyzed using the Fisher information matrix rank discussed in Section 2.5. Local identifiability was determined for parameter configurations in a neighborhood of the realizations associated with maximum posterior distribution from MCMC, and using all experimental data rather than just the average experimental stresses at each stretch level. Table 3 Parameter estimates for the GOH constitutive model. Parameter

Units

(1 − ϕ) μ

MPa

ϕ c1

MPa

c2

–

κ

–

ξ

degrees

Age

Mature Aging Aged Mature Aging Aged Mature Aging Aged Mature Aging Aged Mature Aging Aged

Healing - ξ = 0

Healing - ξ ≠ 0 Uninjured

3 Weeks

6 Weeks

Uninjured

3 Weeks

6 Weeks

2.0 ( ± 1.7) 0.8 ( ± 0.7) 0.7 ( ± 0.7) 53.8 ( ± 28.0) 37.4 ( ± 23.2) 27.0 ( ± 20.5) 647.3 ( ± 213.4) 748.1 ( ± 187.0) 724.9 ( ± 224.8) 0.18 ( ± 0.09) 0.20 ( ± 0.08) 0.21 ( ± 0.08) 37.4 ( ± 20.4) 34.0 ( ± 16.8) 33.9 ( ± 22.4)

0.7 ( ± 0.7) 0.8 ( ± 0.7) 1.4 ( ± 0.8) 35.6 ( ± 32.0) 31.2 ( ± 25.9) 45.5 ( ± 29.9) 649.1 ( ± 244.0) 750.5 ( ± 190.4) 505.8 ( ± 287.4) 0.20 ( ± 0.09) 0.20 ( ± 0.08) 0.24 ( ± 0.09) 37.1 ( ± 18.4) 36.9 ( ± 21.0) 49.3 ( ± 19.3)

0.8 ( ± 0.7) 0.8 ( ± 0.7) 0.9 ( ± 0.8) 26.4 ( ± 23.5) 22.9 ( ± 23.6) 34.1 ( ± 26.8) 699.3 ( ± 220.6) 700.9 ( ± 234.6) 681.3 ( ± 218.5) 0.19 ( ± 0.09) 0.17 ( ± 0.09) 0.21 ( ± 0.09) 33.7 ( ± 19.3) 37.2 ( ± 22.0) 35.0 ( ± 20.4)

1.5 ( ± 1.3) 0.8 ( ± 0.7) 0.6 ( ± 0.6) 57.6 ( ± 23.9) 32.1 ( ± 20.8) 25.6 ( ± 19.7) 614.3 ( ± 208.7) 670.0 ( ± 215.5) 671.0 ( ± 238.6) 0.28 ( ± 0.02) 0.28 ( ± 0.02) 0.27 ( ± 0.02) –

0.6 ( ± 0.5) 0.8 ( ± 0.7) 1.0 ( ± 0.7) 29.3 ( ± 28.7) 27.7 ( ± 23.0) 39.1 ( ± 29.2) 605.6 ( ± 263.5) 720.2 ( ± 203.8) 501.7 ( ± 288.6) 0.27 ( ± 0.06) 0.28 ( ± 0.03) 0.31 ( ± 0.03) –

0.8 ( ± 0.7) 0.7 ( ± 0.6) 0.9 ( ± 0.8) 24.4 ( ± 22.6) 19.8 ( ± 17.5) 37.5 ( ± 30.0) 670.1 ( ± 229.0) 680.6 ( ± 231.1) 655.8 ( ± 239.3) 0.27 ( ± 0.03) 0.27 ( ± 0.03) 0.28 ( ± 0.03) –

294

Journal of the Mechanical Behavior of Biomedical Materials 96 (2019) 285–300

A.R. Akintunde, et al.

Table 4 Parameter estimates for the SHR constitutive model. Parameter

Units

Age

(1 − ϕ) μ

MPa

ϕE

MPa

θ0

degrees

ψ

degrees

Mature Aging Aged Mature Aging Aged Mature Aging Aged Mature Aging Aged

Healing - ψ = 0

Healing - ψ ≠ 0 Uninjured

3 Weeks

6 Weeks

Uninjured

3 Weeks

6 Weeks

1.6 ( ± 1.3) 0.5 ( ± 0.5) 0.4 ( ± 0.4) 95.3 ( ± 4.4) 84.5 ( ± 11.5) 84.6 ( ± 11.8) 21.5 ( ± 1.4) 24.3 ( ± 3.4) 24.14 ( ± 3.7) 5.4 ( ± 4.0) 13.6 ( ± 7.4) 17.1 ( ± 7.7)

0.9 ( ± 0.8) 0.6 ( ± 0.6) 1.6 ( ± 0.8) 66.9 ( ± 22.3) 78.0 ( ± 15.8) 50.2 ( ± 28.6) 37.9 ( ± 16.8) 26.6 ( ± 6.3) 49.64 ( ± 24.2) 18.5 ( ± 12.8) 17.1 ( ± 9.4) 42.1 ( ± 25.2)

0.6 ( ± 0.6) 0.6 ( ± 0.6) 1.0 ( ± 0.9) 75.8 ( ± 16.8) 77.1 ( ± 16.4) 73.0 ( ± 18.5) 26.4 ( ± 5.9) 26.4 ( ± 6.5) 27.88 ( ± 7.9) 16.6 ( ± 9.3) 18.0 ( ± 9.6) 16.1 ( ± 10.0)

1.4 ( ± 1.14) 0.4 ( ± 0.51) 0.3 ( ± 0.36) 95.1 ( ± 4.53) 79.3 ( ± 13.06) 76.7 ( ± 14.46) 21.9 ( ± 1.14) 27.3 ( ± 2.89) 28.6 ( ± 3.28) –

0.8 ( ± 0.66) 0.6 ( ± 0.55) 0.9 ( ± 0.68) 66.6 ( ± 23.56) 71.2 ( ± 18.41) 43.9 ( ± 26.18) 50.9 ( ± 15.76) 32.2 ( ± 5.50) 62.3 ( ± 18.85) –

0.6 ( ± 0.55) 0.6 ( ± 0.55) 0.9 ( ± 0.86) 71.5 ( ± 18.55) 70.9 ( ± 19.03) 69.6 ( ± 19.69) 32.3 ( ± 5.53) 33.1 ( ± 5.96) 34.2 ( ± 8.21) –

the interval [ [θi] − 3 σ [θi],  [θi] + 3 σ [θi]], where  [θi] is the average value of the i-th parameter computed from the MCMC traces, while σ [θi] is its standard deviation. We also have determined the polynomial chaos coefficients using numerical integration on a third-order Smolyak sparse grid and the Gauss-Legendre integration rule. We Remark how the polynomial chaos-based sensitivity indices do not account for the correlations among parameters, assumed independent. The average sensitivity indices are computed together with their 10%–90% confidence intervals over age/healing conditions, where the likelihood in MCMC is assembled from the complete set of experimental stretchstress data, and we have also compared non-zero and zero alignment angles for the GOH and SHR models. Local sensitivity indices (LSI) are instead computed using normalized derivatives from central finite differences (FD). Common to all models, the sensitivity indices at zero stretch are zero as expected, as a zero stress is obtained irrespective of the selected realization in parameter space.

Table 5 Parameter estimates for the FR constitutive model. Parameter

Units

Ec

MPa

Ee

MPa

β

–

Age

Mature Aging Aged Mature Aging Aged Mature Aging Aged

Healing Uninjured

3 Weeks

6 Weeks

99.3 ( ± 0.7) 99.3 ( ± 0.7) 99.5 ( ± 1.2) 0.9 ( ± 0.9) 0.5 ( ± 0.8) 0.6 ( ± 1.4) 30.4 ( ± 4.0) 19.0 ( ± 2.6) 17.8 ( ± 3.1)

96.2 ( ± 3.3) 87.4 ( ± 7.0) 75.9 ( ± 13.6) 0.6 ( ± 0.6) 0.9 ( ± 1.0) 0.5 ( ± 0.5) 14.1 ( ± 2.3) 16.6 ( ± 3.6) 13.0 ( ± 2.2)

98.0 ( ± 1.8) 96.2 ( ± 3.4) 98.3 ( ± 1.6) 0.8 ( ± 0.6) 0.7 ( ± 0.6) 0.4 ( ± 0.2) 16.2 ( ± 3.9) 16.0 ( ± 2.9) 16.9 ( ± 1.4)

3.4. Local and global sensitivity analysis The next sections report the global and local sensitivity indices for the three models. We have plotted the direct and total sensitivity indices (SI) using the Monte Carlo estimates mentioned in Section 2.6 using 50,000 samples from the posterior distribution, thus accounting for the correlations between model parameters from MCMC. For comparison, we have also computed the direct sensitivities from the polynomial chaos (PC) expansion coefficients. Specifically, independent model parameters were assumed with uniform marginal distribution in

3.4.1. GOH model For non-zero helix angles, the stress response of the GOH model is dominated by the fiber dispersion and alignment angle parameters κ and ξ, respectively and this is captured both using polynomial chaos expansion and from the local sensitivities. The contribution of these two parameters to the total stress variance is observed to vary significantly at larger stretches across age-healing groups (see Figs. 18 and 19). When ξ = 0 , the MC based sensitivities seem to capture a different

Fig. 11. Assimilated and experimental stress-stretch response for the GOH model. The experimental stress at various ages has been slightly shifted on the stretch axis to improve clarity. 295

Journal of the Mechanical Behavior of Biomedical Materials 96 (2019) 285–300

A.R. Akintunde, et al.

Fig. 12. Assimilated and experimental stress-stretch response for the SHR model. The experimental stress at various ages has been slightly shifted on the stretch axis to improve clarity.

large correlation from MCMC occurring between the helix angle ψ and the parameter θ0 , justifies the difference between the Monte Carlo and polynomial chaos sensitivity estimates for these two parameters. By setting ψ = 0 , the importance of the elastin modulus is amplified in the MC estimates, particularly when the underlying correlation of elastic modulus and collagen fiber modulus is considered (i.e., sensitivities from MCMC). Increased importance of the elastin modulus for zero alignment angles is also a feature observed in the GOH model sensitivities. In a neighborhood of the MAP parameter estimate, the response is predominantly affected by the collagen modulus and transition parameter θ0 , as confirmed through the local sensitivities. Finally, setting the alignment angle to zero seems to affect the SHR model less than the GOH model.

trend than polynomial chaos or local sensitivity indices. Specifically, local sensitivities assign a greater importance to the fiber response through the orientation parameter κ and, in second instance, to the ϕ c1 and c2 parameters both at small and large strain regimes. Conversely, the Monte Carlo estimates for the global sensitivities appear to assign a greater importance to the stiffness contribution of the elastin matrix rather than collagen fibers. Finally, the Monte Carlo estimates for the direct and total sensitivity indices appear to be very similar, with limited resulting interactions among the parameters. Note how the importance of parameters for the GOH model is significantly altered by setting ξ = 0 , consistent with the large values of the alignment angle estimates in Table 3. 3.4.2. SHR model The modulus of the elastin matrix appears to be an important parameter for the SHR model with no-zero alignment angle ψ, particularly in the small stretch regime, as expected (see Figs. 20 and 21). Additionally, the large sensitivities of the stress response to the alignment angle ψ and transition parameter θ0 appear to be consistent across global and local metrics, while the collagen elastic modulus shows an increasing effect on the stress response for increasing stretch levels. The

3.4.3. FR model The most important parameters are the elastin modulus Ee and the transition parameter β. This latter parameter affects how much of the stress-stretch response relates to a linear regime and so it is expected for β to have larger sensitivity than Ec . Additionally, the importance of the elastin modulus descreases for an increasing stretch, while the importance of the transition parameter increases. Finally, the importance

Fig. 13. Assimilated and experimental stress-stretch response for the FR model. The experimental stress at various ages has been slightly shifted on the stretch axis to improve clarity. 296

Journal of the Mechanical Behavior of Biomedical Materials 96 (2019) 285–300

A.R. Akintunde, et al.

Fig. 14. GOH parameter correlations from MCMC under average (a) and complete (d) experimental data. The parameter correlations for the SHR (b,e) and FR (c,f) models are also shown.

structure and function from these parameters. In this context, i.e., inverse problems, superposition of stress generated by different microstructural mechanisms typically reduces the identifiability of the model parameters from experiments, phenomenon that is even more apparent under limited experimental evidence. In this study, we employ Bayesian inference to compare three constitutive models for the murine patellar tendon representing three degrees of microstructural complexity and investigate how their parameters may be inferred from uncertain uniaxial stress-stretch experimental data. The experimental data was collected on mice under various conditions related to age and degree of injury (Dunkman et al., 2013, 2014a, 2014b; Mienaltowski et al., 2016). For two of these models, i.e., GOH and SHR, an additional parameter for the fiber and fascicle alignment angle, respectively was

of the collagen modulus Ec seem to increase with the stretch level (see Fig. 22). 4. Discussion Micro-structurally motivated constitutive models may be formulated as a superposition of the response provided by various constituents, possibly operating at overlapping stretch regimes. This approach is typically verified in terms of forward model evaluations and the ability to have a certain number of experimental results in the range of the model, i.e., the set of representable stress-stretch responses. A different but perhaps more important problem relates to identify the parameters of the models from uncertain experimental data and infer

Fig. 15. Parameter correlations under zero alignment angle. GOH parameter correlations from MCMC under average (a) and complete (b) experimental data. The parameter correlations for the SHR (c,d) are also shown. 297

Journal of the Mechanical Behavior of Biomedical Materials 96 (2019) 285–300

A.R. Akintunde, et al.

Fig. 16. Representative marginal distributions for the standard deviation of the experimental stress σε in the SHR model.

Fig. 17. Fisher Information matrix eigenvalues for GOH and SHR models. Models with zero and no-zero alignment angles are compared in (a) to emphasize the change in local identifiability promoted by suppressing the alignment angles. The null eigenvectors are also illustrated in (b,c) for the GOH and SHR models, respectively.

Fig. 18. Average global and local sensitivity indices for the GOH model and associated 10%–90% percentile confidence across age-healing groups. Model parameters are identified using the complete dataset of uniaxial tests.

Fig. 19. Average global and local sensitivity indices for the GOH model and associated 10%–90% percentile confidence across age-healing groups. Model parameters are identified using the complete dataset of uniaxial tests and with zero alignment angle ξ.

parameters, as confirmed by a flat marginal posterior if this parameter is added to the inference. Third, we observe a trade-off between the complexity of the underlying microstructural constitutive model and effective inference of tendon structure-function from the data. Supporting our prior work, the FR phenomenological model was the only model for which changes in function caused by age/healing could be directly observable from the model parameters (Akintunde and Miller, 2018). This model is however not perfect as the transition

added as a free parameter, and the effects on the MAP parameter estimates, identifiability and sensitivity observed. We report the following findings. First, from numerical experiments comparing inference with complete and average stress responses, we conclude that use of all available data has a positive effect on constraining inference and leads to highest learning factors, i.e., reduces the variability in the estimated parameters. Second, changes in the data standard deviation σε do not seems to affect the estimates of the 298

Journal of the Mechanical Behavior of Biomedical Materials 96 (2019) 285–300

A.R. Akintunde, et al.

Fig. 20. Average global and local sensitivity indices for the SHR model and associated 10%–90% percentile confidence. Model parameter are identified using the complete dataset of uniaxial tests.

Fig. 21. Average global and local sensitivity indices for the SHR model and associated 10%–90% percentile confidence. Model parameters are identified using the complete dataset of uniaxial tests and with zero helix angle.

Fig. 22. Average global and local sensitivity indices for the FR model and associated 10%–90% percentile confidence across age-healing groups. Model parameters are identified using the complete dataset of uniaxial tests.

simple statistical model (which can be further refined based on a more detailed analysis of the source of experimental errors and their correlation) and by entering the prior distribution for each parameter (based on histological and biomechanical admissibility). We have generated parameter samples from the joint posterior distribution through MCMC with a large number of model solutions, examined distributions of model results, parameter correlations and performed an analysis of the parameter sensitivity and identifiability. It is evident how this process requires more work than classical parameter fitting, particularly due to the large volume of generated results, and the additional time required to verify their consistency, where the comparisons need to be performed between distributions rather than point values. Acceptable results are evaluated based on various metrics. Essential is to verify the convergence of the MCMC sampler, which may be far from trivial for non-identifiable and overparameterized models. Additionally, results are not set in stone and maybe refined taking advantage of additional structure by specifying, e.g., a different prior thus leading to improved estimators characterized by reduced variance (see, e.g., the multi-level estimator discussed in Schiavazzi et al. (2017). As discussed in the introduction, additional experimental evidence (e.g., multi-axial testing) might also be provided to further reduce the variance of the estimates. Thus, parameter estimation under uncertainty is more a dynamic process subject to iterative refinements than a static, certain but possibly myopic answer.

parameter determines the extent of the linear response region, and therefore inference of β interacts with inference of Ec . Fourth, this study highlights the importance of investigating the amount of variability of the parameter estimates, particularly in case of point estimates obtained through optimization. While MAP and optimal parameter estimates may be possibly close in simpler phenomenological models, the same does not seem to occur with models that include a more complex microstructure. Fifth, introduction of structural parameters such as the fiber and fascicle alignment angle reduces the identifiability of the other model parameters due to a non trivial interaction, for example, with the transition parameter as shown in Fig. 1 for the SHR model. Finally, this paper demonstrates how Bayesian inference and the analysis of parameter sensitivity and identifiability are essential complements to the knowledge of the underlying microstructural mechanics to provide constitutive models that are informative on the structurefunction of a certain biological tissue and, in particular, able to distill the information obtained from experimental testing. We also showed how all analysis tasks in uncertainty quantification can be delegated to a library able to include a constitutive model as a black box and presented results using the tulip library. We hope this study may provide guidance to other researchers in the solution of inverse problems under uncertainty, for applications in biomechanics. In this spirit, we have organized the inference tasks in this study following a rational order, starting with the specification of a 299

Journal of the Mechanical Behavior of Biomedical Materials 96 (2019) 285–300

A.R. Akintunde, et al.

This work investigates the difficulties in inferring structural or functional changes from the solution of an inverse problem starting from uniaxial stress-stretch experimental data for various age/healing groups in the murine patellar tendon. Future work will instead be devoted to the development of improved Bayesian approaches for estimating the contribution of multiple constituents. Future research will also be devoted to drive data collection in order to estimate specific parameters.

Springer Science & Business Media. Humphrey, J.D., DeLange, S., 2013. An Introduction to Biomechanics: Solids and Fluids, Analysis and Design. Springer Science & Business Media. Ionescu, I., Guilkey, J., Berzins, M., Kirby, R.M., Weiss, J., 2005. Computational simulation of penetrating trauma in biological soft tissues using the material point method. Stud. Health Technol. Inf. 111, 213–218. Jackman, S., 2000. Estimation and inference via Bayesian simulation: an introduction to Markov chain Monte Carlo. Am. J. Pol. Sci. 375–404. Janon, A., Klein, T., Lagnoux, A., Nodet, M., Prieur, C., 2014. Asymptotic normality and efficiency of two Sobol index estimators. ESAIM P. S. 18, 342–364. Józsa, L.G., 1997. Human Tendons: Anatomy, Physiology, and Pathology. Human Kinetics, Champaign, IL. Kannus, P., 2000. Structure of the tendon connective tissue. Scand. J. Med. Sci. Sports 10 (6), 312–320. Lin, T.W., Cardenas, L., Glaser, D.L., Soslowsky, L.J., 2006. Tendon healing in interleukin4 and interleukin-6 knockout mice. J. Biomech. 39 (1), 61–69. Madireddy, S., Sista, B., Vemaganti, K., 2016. Bayesian calibration of hyperelastic constitutive models of soft tissue. J. Mech. Behav. Biomed. Mater. 59, 108–127. McKinnon, K.I.M., 1998. Convergence of the Nelder–Mead simplex method to a nonstationary point. SIAM J. Optim. 9 (1), 148–158. Miao, H., Xia, X., Perelson, A.S., Wu, H., 2011. On identifiability of nonlinear ODE models and applications in viral dynamics. SIAM Rev. 53 (1), 3–39. Mienaltowski, M.J., Dunkman, A.A., Buckley, M.R., Beason, D.P., Adams, S.M., Birk, D.E., Soslowsky, L.J., 2016. Injury response of geriatric mouse patellar tendons. J. Orthop. Res. 34 (7), 1256–1263. Miller, K.S., Lee, Y.U., Naito, Y., Breuer, C.K., Humphrey, J.D., 2014. Computational model of the in vivo development of a tissue engineered vein from an implanted polymeric construct. J. Biomech. 47 (9), 2080–2087. Nelder, J.A., Mead, R., 1965. A simplex method for function minimization. Comput. J. 7 (4), 308–313. Ogden, R.W., 2003. Nonlinear elasticity, anisotropy, material stability and residual stresses in soft tissue. In: Biomechanics of Soft Tissue in Cardiovascular Systems. Springer, pp. 65–108. Owen, A.B., 2013. Variance components and generalized Sobol’ indices. SIAM/ASA J. Uncertain. Quantification 1 (1), 19–41. Ritto, T.G., Nunes, L.C.S., 2015. Bayesian model selection of hyperelastic models for simple and pure shear at large deformations. Comput. Struct. 156, 101–109. Robertson, D., Cook, D., 2014. Unrealistic statistics: how average constitutive coefficients can produce non-physical results. J. Mech. Behav. Biomed. Mater. 40, 234–239. Rothenberg, T.J., 1971. Identification in parametric models. Econometrica. J. Econom. Soc. 577–591. Saltelli, A., 2002. Making best use of model evaluations to compute sensitivity indices. Comput. Phys. Commun. 145 (2), 280–297. S. Sankaran, J.D. Humphrey, and A.L. Marsden. An efficient framework for optimization and parameter sensitivity analysis in arterial growth and remodeling computations. Comput. Methods Appl. Mech. Eng., 256:200–210, 2013. Schiavazzi, D.E., Baretta, A., Pennati, G., Hsia, T.Y., Marsden, A.L., 2017. Patient-specific parameter estimation in single-ventricle lumped circulation models under uncertainty. Int. J. Num. Method. Biomed. Eng. 33 (3). Seyedsalehi, S., Zhang, L., Choi, J., Baek, S., 2015. Prior distributions of material parameters for bayesian calibration of growth and remodeling computational model of abdominal aortic wall. J. Biomech. Eng. 137 (10) 101001. Shearer, T., 2015. A new strain energy function for modelling ligaments and tendons whose fascicles have a helical arrangement of fibrils. J. Biomech. 48 (12), 3017–3025. Sobol, I.M., 2001. Global sensitivity indices for nonlinear mathematical models and their monte carlo estimates. Math. Comput. Simulat. 55 (1–3), 271–280. Storn, R., Price, K., 1997. Differential evolution–a simple and efficient heuristic for global optimization over continuous spaces. J. Glob. Optim. 11 (4), 341–359. Taber, L.A., 2004. Nonlinear Theory of Elasticity: Applications in Biomechanics. World Scientific. Tarantola, A., 2005. Inverse Problem Theory and Methods for Model Parameter Estimation, vol. 89 SIAM. Vergari, C., Pourcelot, P., Holden, L., Ravary-Plumioën, B., Gerard, G., Laugier, P., Mitton, D., Crevier-Denoix, N., 2011. True stress and Poisson's ratio of tendons during loading. J. Biomech. 44 (4), 719–724. Vrugt, J.A., 2016. Markov chain Monte Carlo simulation using the DREAM software package: theory, concepts, and MATLAB implementation. Environ. Model. Softw 75, 273–316. Vrugt, J.A., Ter Braak, C.J.F., Diks, C.G.H., Robinson, B.A., Hyman, J.M., Higdon, D., 2009. Accelerating Markov chain Monte Carlo simulation by differential evolution with self-adaptive randomized subspace sampling. Int. J. Nonlinear Sci. Numer. Stimul. 10 (3), 273–290. Weiss, J.A., Maker, B.N., Govindjee, S., 1996. Finite element implementation of incompressible, transversely isotropic hyperelasticity. Comput. Methods Appl. Mech. Eng. 135 (1–2), 107–128.

Acknowledgements and disclosures This work was supported by National Institute of Health through the awards P20GM103629 (KSM) and RO1-EB018302 (PI Alison Marsden, subcontractor DES). We wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome. This research used computational resources provided through the Center for Research Computing at the University of Notre Dame and high performance computing resources at Tulane University. References Akintunde, A.R., Miller, K.S., 2018. Evaluation of microstructurally motivated constitutive models to describe age-dependent tendon healing. Biomechanics Model. Mechanobiol. 17 (3), 793–814. Andriotis, O.G., Chang, S.W., Vanleene, M., Howarth, P.H., Davies, D.E., Shefelbine, S.J., Buehler, M.J., Thurner, P.J., 2015. Structure–mechanics relationships of collagen fibrils in the osteogenesis imperfecta mouse model. J. R. Soc. Interface 12 (111) 20150701. Beason, D.P., Kuntz, A.F., Hsu, J.E., Miller, K.S., Soslowsky, L.J., 2012. Development and evaluation of multiple tendon injury models in the mouse. J. Biomech. 45 (8), 1550–1553. Bellini, C., Di Martino, E.S., 2012. A mechanical characterization of the porcine atria at the healthy stage and after ventricular tachypacing. J. Biomech. Eng. 134 (2), 021008. Doraiswamy, S., Srinivasa, A.R., 2013. A Bayesian Approach to Accounting for Variability in Mechanical Properties in Biomaterials. 1312.2856. Dourte, L.M., Pathmanathan, L., Mienaltowski, M.J., Jawad, A.F., Birk, D.E., Soslowsky, L.J., 2013. Mechanical, compositional, and structural properties of the mouse patellar tendon with changes in biglycan gene expression. J. Orthop. Res. 31 (9), 1430–1437. Dunkman, A.A., Buckley, M.R., Mienaltowski, M.J., Adams, S.M., Thomas, S.J., Satchell, L., Kumar, A., Pathmanathan, L., Beason, D.P., Iozzo, R.V., et al., 2013. Decorin expression is important for age-related changes in tendon structure and mechanical properties. Matrix Biol. 32 (1), 3–13. Dunkman, A.A., Buckley, M.R., Mienaltowski, M.J., Adams, S.M., Thomas, S.J., Kumar, A., Beason, D.P., Iozzo, R.V., Birk, D.E., Soslowsky, L.J., 2014a. The injury response of aged tendons in the absence of biglycan and decorin. Matrix Biol. 35, 232–238. Dunkman, A.A., Buckley, M.R., Mienaltowski, M.J., Adams, S.M., Thomas, S.J., Satchell, L., Kumar, A., Pathmanathan, L., Beason, D.P., Iozzo, R.V., et al., 2014b. The tendon injury response is influenced by decorin and biglycan. Ann. Biomed. Eng. 42 (3), 619–630. Ferruzzi, J., Bersi, M.R., Humphrey, J.D., 2013. Biomechanical phenotyping of central arteries in health and disease: advantages of and methods for murine models. Ann. Biomed. Eng. 41 (7), 1311–1330. Fratzl, P., 2008. Collagen: structure and mechanics, an introduction. In: Collagen. Springer, pp. 1–13. Freed, A.D., Rajagopal, K.R., 2016. A promising approach for modeling biological fibers. Acta Mech. 227 (6), 1609–1619. Gasser, T.C., Ogden, R.W., Holzapfel, G.A., 2006. Hyperelastic modelling of arterial layers with distributed collagen fibre orientations. J. R. Soc. Interface 3 (6), 15–35. Gelman, A., Rubin, D.B., 1992. Inference from iterative simulation using multiple sequences. Stat. Sci. 457, 472. Gilks, W.R., Roberts, G.O., 1996. Strategies for improving MCMC. In: Markov Chain Monte Carlo in Practice. Springer, pp. 89–114. Gustafson, P., Gelfand, A.E., Sahu, S.K., Johnson, W.O., Hanson, T.E., Joseph, L., Lee, J., 2005. On model expansion, model contraction, identifiability and prior information: two illustrative scenarios involving mismeasured variables. Stat. Sci. 111–140. Holzapfel, G.A., Gasser, T.C., Ogden, R.W., 2000. A new constitutive framework for arterial wall mechanics and a comparative study of material models. J. Elast. Phys. Sci. Solids 61 (1–3), 1–48. Humphrey, J.D., 2013. Cardiovascular Solid Mechanics: Cells, Tissues, and Organs.

300