Stochastic isotropic hyperelastic materials: constitutive calibration and model selection

1. Introduction

Mathematical models of solid materials often necessitate approaches that capture the randomness due to the uncertainties in the mechanical responses. As predictions depend on constitutive models, it may not be adequate for a mathematical model to depend on a single set of constant parameters, regardless of how well they seem to agree with certain experimental measurements. A statistical theory of heterogeneous linearly elastic solids is introduced by McCoy (1973) [1]. Stochastic strategies for the investigation of mesoscopic mechanical effects in random materials were proposed by Huet (1990) [2]. Further developments in the stochastic modelling of heterogeneous solids were reviewed in [3]. Recently, there has been a growing interest in probability and statistical techniques for engineering and biomedical applications, where the calibration of models using available data and the quantification of uncertainties in model parameters are of utmost importance [4–6]. There are, however, many challenges introduced by the consideration and quantification of uncertainties in mathematical models, and their use in making predictions, some of which are discussed in [7–12].

For natural and engineered materials, uncertainties in the experimental observations typically arise from the inherent micro-structural inhomogeneity [13,14], sample-to-sample intrinsic variability, or when elastic data are extracted from viscoelastic mechanical tests [15–19]. For these materials, hyperelastic models based on mean data values constitute a starting point for the development of more complex models. Stochastic models accounting also for data dispersion give additional insight about uncertainty and can provide useful bounds on a model’s prediction. A review of statistical approaches applied to the mechanical analysis of rubber-like networks is presented in [20]. Constitutive equations for soft tissues, including those based on statistical modelling for the evolution of the collagen network, are reviewed in [21].

A non-deterministic approach to model the stiffness variations of porcine liver tissue under compression was first proposed in [22]. In this case, experimental strain values at a given stress were assumed to vary according to a normal distribution, for which the mean value and standard deviation are independent, and a five-parameter Mooney–Rivlin hyperelastic model was calibrated numerically to the mean stress–strain curve.

Recently, stochastic strategies based on information theory, which aids with the calibration of hyperelastic models for isotropic elastic solids, were proposed in [23–25]. Specifically, stochastic–hyperelastic models were identified from experimental data consisting of the mean values and standard deviations of elastic stresses under finite strain deformations. Prior to this, in [26,27], a similar strategy was applied to the stochastic representation of tensor-valued random variables and random fields in linear elasticity. These strategies rely on the maximum entropy principle for a discrete probability distribution introduced by Jaynes (1957) [28–30]. The measure of entropy (or uncertainty) of a discrete probability distribution was first defined by Shannon (1948) [31,32] in the context of information theory. In [25], the stochastic approach was employed for the calibration of Ogden-type models to brain, liver and spinal cord data representing mean values and standard deviation of the first Piola–Kirchhoff stress under finite compression tests. Ogden-type strain–energy functions and their extension to compressible materials in a stochastic framework were formulated originally in [23,24]. Within this framework, strain–energy functions with random field parameters were obtained under a combination of physically realistic and theoretical restrictions, namely:

• (i) Material objectivity and symmetry. The principle of material objectivity (frame indifference) states that constitutive equations must be invariant under changes of frame of reference [33, p. 44]. In the case of isotropic materials, which have the same mechanical properties in all directions, material symmetry is taken into account by expressing the strain–energy function, equivalently, as a symmetric function of the principal stretches [33–36].

• (ii) Hadamard’s well-posedness property. Well-posedness is enforced by restricting strain–energy functions to satisfy the polyconvexity and coercivity conditions [37–39]. In particular, for Ogden-type stochastic–hyperelastic models, positive random field model parameters were assumed.

• (iii) Linear limit consistency. Mechanical consistency with the linear elasticity theory requires that the classical shear modulus μ>0 is recovered under small elastic strains.

• (iv) Finite mean and variance for the random linear shear modulus. For the random shear modulus, μ, and its inverse, μ⁻¹, the expectation is that they are second-order random variables, i.e. they have finite moments of order two (finite mean and variance). Under these constraints, the maximum entropy principle implies that μ follows a Gamma probability distribution [40,41].

Here, we devise an explicit strategy for the calibration of homogeneous isotropic hyperelastic models with the random field parameters, following probability laws, to discrete mean values and standard deviation of either the stress–strain function or the nonlinear shear modulus, which is a function of the deformation under large strain and coincides with the classical shear modulus under small strain. For isotropic hyperelastic materials, the formal derivation of key nonlinear elastic parameters and their application to model calibration is reviewed in [42]. In practice, these quantities can meaningfully take on different values, corresponding to possible outcomes of experiments, and in general, more than one parametrized model will be available to explain their behaviour. Our modelling framework combines finite elasticity and information theory, as follows. At the level of finite elasticity, we consider the following conditions:

• (i) Material objectivity and symmetry. Material objectivity is guaranteed by considering strain–energy functions defined in terms of invariants. As usual, for isotropic materials, we assume the existence of a symmetric strain–energy function $\mathcal{W}({\mathrm{\lambda }}_1,{\mathrm{\lambda }}_2,{\mathrm{\lambda }}_3)$ of the principal stretches {λ_i}_i=1,2,3.

• (ii) Baker–Ericksen inequalities. In addition to the fundamental principle of objectivity and material symmetry, in order for the behaviour of a hyperelastic material to be physically realistic, there are some universally accepted constraints on the constitutive equations [33]. Specifically, for a hyperelastic body, the Baker–Ericksen (BE) inequalities state that the greater principal stress occurs in the direction of the greater principal stretch. In particular, under uniaxial tension, the deformation is a simple extension in the direction of the tensile force if and only if the BE inequalities hold [43,44]. Under these mechanical constraints, the nonlinear shear modulus, which varies with the deformation and is equal to the linear shear modulus in the small strain, is always positive [42].

• (iii) Non-polyconvexity. For a general theoretical framework that is nevertheless consistent with the observed mechanical behaviour of many materials operating in large strain deformation [18,19,45–48], restriction to the class of polyconvex strain-energy functions is not required [49–51]. The non-polyconvexity allows for more general a priori bounds on the random hyperelastic parameters to be chosen during the calibration process.

Our approach to stochastic elasticity and model selection further relies on the following assumptions:

• (iv) Finite mean and variance for the random nonlinear shear modulus. At any given deformation, the nonlinear shear modulus and its inverse are second-order random variables, i.e. they have finite mean and finite variance.

• (v) Model selection and Occam’s criterion. As alternative models that differ in form or number of parameters can be derived that reasonably approximate the data and its variability, we apply Bayesian statistics [52] and Occam’s principle [53–56] to select the best possible model from a given family of models [57].

In the next sections, we summarize the finite elasticity setting (§2) and develop our stochastic–deterministic strategy (§3), which we employ to construct explicit models from experimental data for rubber and brain tissue under uniaxial or multiaxial loads. We then apply Occam’s criterion (§4) to select a model among competing models calibrated to the available data (§5).

2. Finite elasticity prerequisites

We consider a unit cube of homogeneous isotropic incompressible hyperelastic material, subject to the following homogeneous deformation consisting of a simple shear superposed on a finite axial stretch [18,19,42,58–61],

where (X₁,X₂,X₃) and (x₁,x₂,x₃) are the Cartesian coordinates for the reference (Lagrangian) and the current (Eulerian) configuration, respectively, and k>0 and a>0 are positive constants representing the shear parameter and the axial stretch (0<a<1 for axial compression and a>1 for axial tension), respectively. We recall that homogeneous deformations are the same regardless of the geometry of the body, and are universal, i.e. they can be maintained in every homogeneous isotropic hyperelastic body by application of suitable traction [62–65]. For the homogeneous deformation (2.1), the constant gradient tensor and the corresponding left Cauchy–Green tensor are, respectively, where the superscript ‘T’ denotes the transpose. The principal stretches {λ_i}_i=1,2,3, such that ${\{{\mathrm{\lambda }}_i^2\}}_{i=1,2,3}$ are the eigenvalues of the left Cauchy–Green tensor, B, satisfy When the homogeneous isotropic incompressible hyperelastic material is described by a strain-energy function $\mathcal{W}({\mathrm{\lambda }}_1,{\mathrm{\lambda }}_2,{\mathrm{\lambda }}_3)$, where {λ_i}_i=1,2,3 are the principal stretches [34, p. 94], the Cauchy stress tensor, representing the force per unit area in the current configuration, is equal to where I is the identity tensor, p is the Lagrange multiplier associated with the incompressibility constraint ($det\mathbf{\text{F}}=1$), and are the constitutive coefficients. Then, the non-zero components of the Cauchy stress tensor in Cartesian coordinates are The principal components of the Cauchy stress tensor (i.e. its principal eigenvalues) are [33, p. 143] The associated first Piola–Kirchhoff stress tensor, representing the force per unit area in the reference configuration, is defined as [33, pp. 124–125]

(a) Nonlinear shear modulus

We note the following: (i) stresses are constant given a homogeneous deformation of the form (2.1); and (ii) the shear component of the first Piola–Kirchhoff stress tensor (2.8), P₁₂=σ₁₂/a, is proportional to the shear strain, ka. These observations justify the introduction of the nonlinear shear modulus [42]

This modulus is a function of the deformation, it is independent of the Lagrange multiplier, p, and it can be estimated directly from experimental observations if the shear force is known. Equivalently, by the representation (2.7) of the principal Cauchy stresses, the nonlinear shear modulus (2.9) can be expressed as [42] Hence, this modulus is always positive assuming that the following Baker–Ericksen (BE) inequalities hold [33, p. 158]:

When a→1, in the deformation (2.1), simple shear is superposed on infinitesimal axial stretch. Then, the nonlinear shear modulus given by (2.9) converges to the nonlinear shear modulus for simple shear [42],

where ${\hat{\beta }}_1=\underset{a\to 1}{lim}{\beta }_1$ and ${\hat{\beta }}_{-1}=\underset{a\to 1}{lim}{\beta }_{-1}$, and the corresponding principal stretches are ${\hat{\mathrm{\lambda }}}_i=\underset{a\to 1}{lim}{\mathrm{\lambda }}_i$, i=1,2,3.

Similarly, when k→0, the deformation (2.1) becomes an infinitesimal shear superposed on a finite axial stretch. In this case, the nonlinear shear modulus (2.9) converges to [42]

where ${\tilde{\beta }}_1=\underset{k\to 0}{lim}{\beta }_1$ and ${\tilde{\beta }}_{-1}=\underset{k\to 0}{lim}{\beta }_{-1}$, and the principal stretches are ${\tilde{\mathrm{\lambda }}}_i=\underset{k\to 0}{lim}{\mathrm{\lambda }}_i$, i=1,2,3.

In the linear elastic limit, i.e. when k→0 and a→1, the nonlinear shear moduli defined by (2.9), (2.12) and (2.13) converge to the classical shear modulus from the infinitesimal theory [[33, p. 179], [42]],

where ${\overline{\beta }}_1=\underset{a\to 1}{lim}\underset{k\to 0}{lim}{\beta }_1$ and ${\overline{\beta }}_{-1}=\underset{a\to 1}{lim}\underset{k\to 0}{lim}{\beta }_{-1}$, and the principal stretches are ${\overline{\mathrm{\lambda }}}_i=\underset{k\to 0}{lim}{\mathrm{\lambda }}_i$, i=1,2,3.

Examples of the shear moduli $\hat{\mu }(k)$ of (2.12), $\tilde{\mu }(a)$ of (2.13) and $\overline{\mu }$ of (2.14) for particular strain–energy functions are given in table 1. For these models, we can write the nonlinear shear modulus at small shear superposed on finite axial stretch as

where g_p(a), p=1,…,n, are functions of the stretch parameter a>0. For instance, the Yeoh model given in table 1 has g₁(a)=1 and g_p(a)=(a²+2/a−3)^p−1 for p>1. Similarly, under simple shear, the nonlinear shear modulus takes on the form

where h_p(k), p=1,…,n, are functions of the shear parameter k>0.

3. Stochastic–hyperelastic modelling

Our aim is to construct a stochastic–hyperelastic model $\mathcal{W}({\mathrm{\lambda }}_1,{\mathrm{\lambda }}_2,{\mathrm{\lambda }}_3)$ from a given dataset comprising the mean values and standard deviations for either the random shear stress P₁₂ or the nonlinear shear moduli $\hat{\mu }(a)$ or $\tilde{\mu }(a)$, defined by (2.12) or (2.13), respectively. In general, more than one parametrized model will be available that reasonably approximates the data and its variability. These models may differ in form or number of parameters. Here, we consider the constitutive models listed in table 1, where explicit forms for the shear moduli $\hat{\mu }(k)$ of (2.12) and $\tilde{\mu }(a)$ of (2.13), and their linear elastic limit $\overline{\mu }$ of (2.14) are provided. For these models, we focus on the implications of the variable data for the coefficients C_p, p=1,…,n, which are random constant parameters, independent of the other parameters, which are treated as deterministic constants. Other material models could then be treated in a similar manner. Here, we explain in detail the calibration of models from table 1 to experimental mean values and standard deviations of the nonlinear shear modulus $\tilde{\mu }(a)$ of (2.13), for small shear superposed on finite axial stretch. The calibration to data values of the elastic (shear) stress [25] or of the nonlinear shear modulus $\hat{\mu }(a)$ of (2.12), for simple shear, can then be performed analogously. Henceforth, the following notation is used: a quantity with an overbar denotes a value appearing in the theory of linear elasticity (e.g. $\overline{\mu }$); an underlined quantity denotes the mean value of that quantity (e.g. $\underset{\_}{\hat{\mu }}$, $\underset{\_}{\tilde{\mu }}$, $\underset{\_}{\overline{\mu }}$).

(a) Calibration of random field parameters

Whereas, in the deterministic models, we only require one mean value of the modulus provided for each of the m stretches, in the stochastic models, we also consider the measured standard deviation. Therefore, we assume that the given data consist of the mean values ${\{{\underset{\_}{\tilde{\mu }}}_s\}}_{s=1,\dots ,m}$ and the associated standard deviations {d_s}_s=1,…,m of the nonlinear shear modulus (2.13) at the prescribed stretches {a_s}_s=1,…,m. We employ the following two-step procedure:

Step 1. First, we carry on the traditional method used in the deterministic approach [17–19,71]. That is, we determine the mean value of the nonlinear shear modulus (2.13) by minimizing the residual

between the mean nonlinear shear modulus $\underset{\_}{\tilde{\mu }}$ and the mean data values ${\{{\underset{\_}{\tilde{\mu }}}_s\}}_{s=1,\dots ,m}$ at the prescribed stretches {a_s}_s=1,…,m. Doing so, we obtain the mean values {c_p}_p=1,…,n of the random constant coefficients {C_p}_p=1,…,n. If the exponents are not fixed a priori, we also identify the exponents {α_p}_p=1,…,n in the same process.

For the models listed in table 1, by (2.15), the mean value of the nonlinear shear modulus and its linear elastic limit (2.14) take on the respective forms

and

Step 2. Based on the mean values derived at the first step, the goal of the second step is to identify the probability distributions that the random model parameters follow. For the nonlinear shear modulus (2.15), we define the variance

where Var[C_p] is the variance of C_p, and Cov[C_p1,C_p2] is the covariance of C_p1 and C_p2. The standard deviation of the nonlinear shear modulus is the square root of the variance, and, similarly, for every random constant coefficient, C_p, p=1,…,n, the standard deviation is $\parallel C_p\parallel =\sqrt{\text{Var}[C_p]}$. To find ∥C_p∥, we then need to minimize the residual between the standard deviation (3.5) and the associated data {d_s}_s=1,…,m at the prescribed stretches {a_s}_s=1,…,m. Before we do so, we fix the value of the stretch parameter to a particular value a₀>0 that is used for calibration. The corresponding random shear modulus (2.15) is Note that, when a₀=1, $\tilde{\mu }(a_0)$ is simply $\overline{\mu }$. Assuming C_p>b, for all p=1,…,n, where $b>-\mathrm{\infty }$ is chosen a priori, such that the mean values of the random coefficients, which are known from step 1, are bounded away from b, c_p>b, p=1,…,n, we define the auxiliary random parametersThese parameters are such that R_p>0 and, by (3.7), satisfy i.e. they form a complete probability distribution. Then, by (3.8), the random coefficients take on the form Next, following [23–27], for the random nonlinear shear modulus $\tilde{\mu }(a_0)$, defined by (3.7), we set the mathematical expectations: and where, by the constraint (3.11), the mean value $\underset{\_}{\tilde{\mu }}(a_0)$ is fixed and greater than zero, and the logarithmic constraint (3.12) implies that both $\tilde{\mu }(a_0)$ and $\tilde{\mu }{(a_0)}^{-1}$ are second-order random variables (i.e. they have finite mean and finite variance). Critically, equations (3.11) and (3.12) imply that $\tilde{\mu }(a_0)$ follows a Gamma distribution (the maximum entropy distribution) [72,73], Γ(ρ₁(a₀),ρ₂(a₀)), with ρ₁(a₀)>0 and ρ₂(a₀)>0 satisfying where the mean value $\underset{\_}{\tilde{\mu }}(a_0)$ is obtained at step 1. For the random vector (R₁(a₀),…,R_n(a₀)), applying the constraints [23,25] and this vector follows a Dirichlet distribution [41,74], D(ξ₁(a₀),…,ξ_n(a₀)). Then, every random variable R_p(a₀), p=1,…,n, follows a standard Beta distribution [40,41], B(ξ_p(a₀),ψ_p(a₀)), with ξ_p(a₀)>0 and ${\psi }_p(a_0)=\sum _{q=1,q\ne p}^n{\xi }_q(a_0)>0$ satisfying where r_p(a₀) is the mean value and Var[R_p(a₀)] is the variance of R_p(a₀), with the standard deviation $\parallel R_p(a_0)\parallel =\sqrt{\text{Var}[R_p(a_0)]}$. By (3.8), the mean value is equal to and is calculated from the mean values obtained at step 1. Finally, the optimal hyperparameter vectors (ρ₁(a₀),ρ₂(a₀)) and (ξ₁(a₀),…,ξ_n(a₀)) are identified by minimizing the residual for the standard deviation (3.6), and taking into account relations (3.10), (3.13) and (3.16).

(b) The particular case of one-term models

We now specialize the above approach to one-term models, which are of particular interest because, for these models, there is only one random coefficient that needs to be determined, C₁, and only one random auxiliary parameter, R₁=1. For the one-term model, at any stretch a=a₀, the random shear modulus (2.13) is equal to

where $\overline{\mu }$ is its linear elastic limit (2.14). In this case, we apply the two-step procedure, as follows:

Step 1. We determine the mean coefficientc₁, and any other unknown constant parameter appearing in the expression of the strain–energy function, by minimizing the residual function for the mean values (3.1). The mean value of the random shear modulus (3.18) is equal to $\underset{\_}{\tilde{\mu }}(a)=c_1{g}_1(a)$, and its linear elastic limit is $\underset{\_}{\overline{\mu }}=c_1$.

Step 2. The variance defined by (3.4) simplifies to $\text{Var}[\tilde{\mu }(a)]=\text{Var}[C_1]g_1{(a)}^2$, and the corresponding standard deviation, given by (3.5), is equal to $\parallel \tilde{\mu }(a)\parallel =\parallel C_1\parallel g_1(a)$, where $\parallel C_1\parallel =\sqrt{\text{Var}[C_1]}$ is the standard deviation of C₁. By (3.18), assuming (3.11) and (3.12) for $\tilde{\mu }(a_0)$ is equivalent to assuming

and Then, $\overline{\mu }$ follows a Gamma distribution, Γ(ρ₁,ρ₂), with ρ₁>0 and ρ₂>0 satisfying where c₁ is obtained at step 1. After the optimal value of ∥C₁∥ is computed by minimizing the residual (3.6) for the standard deviation, the hyperparameters (ρ₁,ρ₂) are obtained from (3.21).

4. Bayesian model selection and Occam’s principle

Here, we show how Bayesian inference can be employed to select a model among competing models calibrated to the same given data. We denote by P(M) the prior probability of the model M before the data values D are taken into account, and by P(D|M) the likelihood of the data D given the model M, describing the probability of obtaining the data values D from the model M. Then, Bayes’ theorem [11,30,52] is used to update the probability of the model M in the light of the data D. This theorem states that the posterior probability of the model M, denoted by P(M|D), is proportional to the product of the prior and the likelihood, i.e.

where P(D) is the normalization value, also known as the marginal likelihood.

The Bayesian formula (4.1) then provides a methodology for estimating the odds for the model M⁽ⁱ⁾ to the model M^(j) in light of the data D,

where is the Bayes factor. In other words, the posterior odds O_ij for the model M⁽ⁱ⁾ against the model M^(j), given the data D, are equal to the prior odds multiplied by the Bayes factor. We note that, when the two models have equal prior probabilities, i.e. P(M⁽ⁱ⁾)=P(M^(j)), meaning that there is no prior favourite model, the prior odds are 1 and, by (4.2), the posterior odds are equal to the Bayes factor. If the Bayes factor is 1, then Occam’s razor [53–56] implies that one should assign a larger prior probability to the simpler model than to the more complex one for reasons of parsimony. The Bayesian approach that we develop here does not depend on the choice of the prior probabilities.

To maintain a general framework, we assume P(D|M) to be an arbitrary probability that is symmetric about the mean value D=0 and decreasing in the absolute value of D. In this case, the Bayes factor B_ij satisfies the inequality [57]

where ∥D⁽ⁱ⁾∥ and ∥D^(j)∥ designate the standard deviation that the predicted quantity of interest computed with the model M⁽ⁱ⁾ and M^(j), respectively, deviates from the observed data value D. The lower bound on the Bayes factor B_ij given by (4.4) represents an estimate of the amount of evidence against the model M⁽ⁱ⁾, i.e. the maximum support for the model M^(j) provided by the data. The expression on the right of the inequality (4.4) gives an explicit lower bound on the Bayes factor B_ij. Then, by (4.3), B_ij=1/B_ji, and the inverse of the lower bound on the Bayes factor B_ji, provided by exchanging i and j in (4.4), constitutes an upper bound on the Bayes factor B_ij, i.e. Thus, assuming equal prior probabilities, i.e. prior odds are equal for all models, the explicit lower and upper bounds on the Bayes factor provided by (4.4) and (4.5), respectively, represent bounds on the posterior odds. In the following examples, these bounds will be estimated and applied to select a model among competing models calibrated to experimental data.

5. Examples and applications

In this section, we construct explicit stochastic–deterministic models with the random hyperelastic parameters calibrated to discrete mean values and standard deviations of either the stress–strain function or the nonlinear shear modulus estimated from experimental data for rubber and brain tissue, respectively. In all cases, multiple models that differ in form or number of parameters are obtained that reasonably approximate the data and its variability. We further employ the Bayesian model selection, whereby we calculate explicitly the lower and upper bounds on the Bayes factor given by (4.4) and (4.5), then rely on these bounds to select a model among competing models calibrated to the available data. If the bounds are similar, then the simplest model is chosen, as simpler models are more likely to be used even if their approximation of the experimental data is not the best, as advocated in [21].

(a) Rubber-like material

Example 5.1

First, we calibrate the random Piola–Kirchhoff shear stress P₁₂ of three different models from table 1 to the experimental data for rubber material under simple shear reported in [75]. The stochastic–hyperelastic models are as follows:

For these models, the hyperparameters of the Gamma distribution for the random nonlinear shear modulus $\hat{\mu }(k_0)$ of (2.12) at k₀=0.1 were identified. In addition, for the Yeoh model (5.3), the hyperparameters of the Dirichlet probability distribution for the random coefficients {C_p}_p=1,2,3 were also obtained. The corresponding calibrated parameters and prior distribution hyperparameters are listed in table 2. We note that, for the multiple-term model, the choice of b in (3.8) is not unique. The random Piola–Kirchhoff shear stress is plotted in figure 1a. Our results show that the three stochastic models perform very similarly when calibrated to the given data, and the mean value parameters of the Lopez–Pamies model are similar to the values μ≈0.52±0.03 and α≈0.69±0.05 identified in [75]. However, as the mean relative errors (not shown here) are less than 5% when k>0.7 and nearly 43% when k=0.1, it is instructive to also compare the nonlinear shear modulus $\hat{\mu }(k)=P_{12}/k$ of (2.12) for the models and for the data. In this case, we find a sharp increase in the mean values and standard deviation of the nonlinear shear modulus computed from the shear stress data when k<0.7 (figure 1b). Further numerical tests (not shown) reveal that attempting to improve the model approximation when k<0.7 by increasing the number of terms, in any of the three models, causes a steep increase in the mean relative error when k>0.7. This suggests that significant noise may be present in the data when k<0.7.

• Download figure
• Open in new tab
• Download PowerPoint

Next, for each model recorded in table 2, we calculate the standard deviation that the mean shear modulus $\underset{\_}{\hat{\mu }}(0.1)$ deviates from the known mean data value D=0.9, and obtain:

By (4.4) and (4.5), the corresponding Bayes factors satisfy Assuming equal prior probabilities, i.e. prior odds 1, the posterior odds for each model against the other is equal to their respective Bayes’ factors. Then, the bounds on the Bayes factors estimated above suggest that the posterior odds for any of the models against another are also close to 1. Thus, based on which model provides a more accurate approximation for the mean data at k₀=0.1, any of the three models is equally acceptable.

(b) Brain mechanics

Example 5.2 mouse brainmouse brain

In this example, we calibrate the random nonlinear shear modulus $\tilde{\mu }(a)$ of the stochastic one-term Ogden model (5.1) and of the stochastic multiple-term Ogden models of the form

with n=3,4,5, respectively, to mean values and standard deviation data for mouse brain tissue tested under 2% shear superposed on up to 40% tension or compression, in 10% increments. The nonlinear shear modulus of deterministic Ogden models with multiple terms was previously calibrated to similar mean values in [18]. Experimental tests on mouse brain tissue under small shear combined with axial deformations were reported on in [16]. The calibrated parameters and prior distribution hyperparameters of the stochastic models are recorded in table 3. For these models, the hyperparameters of the Gamma distribution for the random linear shear modulus $\overline{\mu }$ of (2.14) were obtained. Again, for multiple-term models, the choice of b is not unique. For a clear illustration, the random nonlinear shear modulus $\tilde{\mu }(a)$ of the one-term and four-term models are plotted in figure 2.

• Download figure
• Open in new tab
• Download PowerPoint

Then, for each model listed in table 3, we estimate the standard deviation that the mean shear modulus $\underset{\_}{\overline{\mu }}$ deviates from the mean data value D=0.1915 at 2% simple shear. Taking the models in the order of their complexity, from the simplest, one-term model, to the most complex, five-term model, we obtain:

When we compare each model to the next more complex one, by (4.4) and (4.5), the Bayes factors satisfy: Assuming the prior odds of one model against another are all equal, the estimated bounds on the Bayes factors represent bounds on the posterior odds. From the above estimates, we first notice that the lower bounds appear to decrease, while the upper bounds tend to increase as the complexity of the models increases. However, as these changes are slow, we can only infer that the posterior odds for each model against another are approximately better than 1/2 (one to two) and worse than 2/1 (two to one), i.e. the posterior odds satisfy In this case, it is instructive to calculate also the approximate bounds on the corresponding posterior probabilities arising from these bounds. Taking P(D|M^(j))=1−P(D|M⁽ⁱ⁾), from the above double inequalities, we obtain 1/3≤P(D|M⁽ⁱ⁾)≤2/3, i=1,2,3,4. As these probabilities are close to 1/2, we conclude that the given data are equally probable according to the computations performed by the simplest or the more complicated models. Thus, given the accuracy of the approximations attained by different models, the explicit Bayesian approach cannot favour any particular model listed in table 3. Occam’s razor then implies that the simplest model should be chosen. In general, within a large set of possible models with no Bayesian favourite, it may not always be clear what ’the simplest model’ means. How do we quantify the notion of simplicity for different functions? However, in our case, we have a well-defined family of models with increasing number of terms. Therefore, we understand ’the simplest model’ as the one with the fewer number of terms, i.e. the one-term (two-parameter) stochastic Ogden model.

Example 5.3

[human brain] Next, we calibrate the random nonlinear shear modulus $\tilde{\mu }(a)$ of the stochastic Ogden models (5.1) and (5.4), where n=3,4,5, to experimental data for human brain tissue under infinitesimal shear superposed on up to 25% tension or compression, in 5% increments. The nonlinear shear modulus of deterministic Ogden models with one or multiple terms was previously calibrated to the same mean values in [19]. Extensive experimental results for human brain tissue under combined shear and axial deformations were reported in [17], where the mean elastic responses were calculated as the average between the viscoelastic loading and unloading paths. The standard deviation considered here represents the range of viscoelastic responses. For the stochastic models, the calibrated parameters and prior distribution hyperparameters are given in table 4. For these models, the hyperparameters of the Gamma distribution for the random linear shear modulus $\overline{\mu }$ of (2.14) were identified. The random nonlinear shear modulus $\tilde{\mu }(a)$ of the one-term and four-term models are illustrated in figure 3.

• Download figure
• Open in new tab
• Download PowerPoint

As in the previous example, for each model recorded in table 4, we estimate the standard deviation that the mean linear shear modulus $\underset{\_}{\overline{\mu }}$ deviates from the known mean data value D=0.3379. We find:

By (4.4) and (4.5), the corresponding Bayes factors satisfy Taking prior odds 1, the bounds on the Bayes factors represent bounds on posterior odds. Again, the above bounds are approximately at most 1/2 against the simpler model and at least 2/1 for the more complex one. Then, comparing one model with another, the posterior probability for each model satisfies 1/3<P(D|M⁽ⁱ⁾)<2/3, i=1,2,3,4, and is close to 1/2. Hence, by taking into account how well each model approximates the mean data, the explicit Bayesian approach does not favour any particular model listed in table 4. Applying Occam’s razor, the simplest model, i.e. the one-term (two-parameter) stochastic Ogden model, should be selected.

Remark 5.4

It is important to remark here that the estimated bounds on the Bayes factors in the two examples involving brain tissue data [18,19] are very similar. Further numerical tests that are not included here also show that similar bounds, i.e. approximately 1/2 and 2, respectively, are found when simple shear superposed on the maximum compression or tension is considered, or when models with more terms are calibrated to the given datasets. By the likelihood principle, we infer that the two datasets have the same likelihood [76]. This is very striking, given that the respective data arise from different brain tissue types tested by different experimental procedures. In this case, the likelihood principle seems to play a unifying role for the experimental testing. A natural question is then: could this principle be further applied to guide experiments?

6. Conclusion

Homogeneous isotropic hyperelastic models can capture characteristic mechanical behaviours of many deformable solids and underpin their analyses and computer simulation. However, natural and bioinspired materials exhibit inherent variations in their elastic properties, which play important roles in their functional performance, and are not represented by a nonlinear elastic constitutive law. For these materials, stochastic representations accounting also for data dispersion contain additional information about the variability of material properties. We combine finite elasticity and information theories to construct homogeneous isotropic hyperelastic models with random field parameters calibrated to discrete mean values and standard deviations of either the stress–strain function or the nonlinear shear modulus, which is a function of the deformation, estimated from experimental tests. These quantities can take on different values, corresponding to possible outcomes of the experiments. In summary, we cast the model parameters as random variables and use the maximum entropy probability distribution to express the uncertainty of the data variability. In our approach, the mean values and standard deviations of the model parameters and the hyperparameters of the underlying probability distribution are calculated formally, although these quantities are not unique in general. As multiple models, which differ in form or number of parameters, can be derived that adequately represent the observed phenomena, we apply Occam’s razor by providing an explicit criterion for model selection, based on Bayesian statistics. We then employ this criterion to select a model among competing models calibrated to the available data for rubber and brain tissues under single or multiaxial loads. Our modelling strategy can further be used to study the variation in the elastic behaviour of solid materials in different applications. In medicine, this research enhances the current solid mechanics research as it will enable better predictions from ensemble data.

Data accessibility

The datasets for this article have been uploaded as supplementary material.

Authors' contributions

All authors contributed equally to all aspects of this article and gave their final approval for publication.

Competing interests

We declare we have no competing interests.

Funding

The support for L.A.M. by the Engineering and Physical Sciences Research Council of Great Britain under research grant no. EP/M011992/1 is gratefully acknowledged.

Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.4020850