## Abstract

A new anisotropic viscoelastic model is developed for application to the aortic valve (AV). The directional dependency in the mechanical properties of the valve, arising from the predominantly circumferential alignment of collagen fibres, is accounted for in the form of transverse isotropy. The rate dependency of the valve's mechanical behaviour is considered to stem from the viscous (*η*) dissipative effects of the AV matrix, and is incorporated as an explicit function of the deformation rate (*η* and

## 1. Introduction

The prevalent structural component of aortic valve (AV) tissue is collagen. It comprises approximately 55% of an intact AV leaflet by dry weight [1], and is present within the tissue in the form of a network of fibres. The fibres are embedded within a viscous ‘gel-like’ matrix of glycosaminoglycans (GAGs) [2], assuming a preferred direction along the circumferential axis of each leaflet, as shown in figure 1. This preferred alignment of the collagen fibres along the circumferential direction endows the AV with strong directional dependency in its mechanical behaviour and material properties; uniaxial and biaxial tensile tests have demonstrated a significant distinction in the elastic properties of the tissue in the circumferential compared with the transverse (radial) direction [4–6], whereas distinctions also exist in the load-bearing capacity and the distensibility of the AV tissue in these two loading directions [5–7].

The viscoelastic characteristics of AV tissue have also been well documented. Tensile deformation tests on AV tissue specimens under various deformation rates have demonstrated a marked rate dependency of the ensuing stress–strain curves [6]. Moreover, studies investigating the time-dependent behaviour of the AV have reported stress relaxation under both uniaxial and biaxial tests [8,9], and creep under uniaxial conditions [9–11].

Based on the abovementioned attributes, the mechanical behaviour of AV tissue may be broadly classified as ‘anisotropic viscoelastic’, to reflect both directional and rate dependency of the mechanical properties of the tissue. Therefore, for mathematical continuum-based models to adequately and appropriately characterize the mechanical behaviour of the AV, it is imperative that both directional- and rate-dependency features are suitably incorporated and accounted for in such models. However, most mathematical continuum-based AV models developed to date have been derived under the assumption of hyperelasticity [7,12–15], and hence, in spite of providing a good fit to experimental stress–strain data obtained at any specific deformation rate, such models cannot, by definition, account for rate effects, nor model AV mechanics over a range of rate-dependent loading conditions. While discounting the rate effects may not produce significant discrepancies between experimental data and model predictions at lower deformation rates, as achievable in typical experiments, the strain rate experienced by the native AV *in vivo* is in the range of 15 000%** **min^{−1} [16]. Such high rates are not achievable by conventional material testing devices *in vitro*, and are likely to affect the material properties and the behaviour of the tissue. Therefore, models that incorporate deformation rate effects are required for accurate description of the *in vivo* mechanical behaviour of the AV.

In addition to discounting the rate effects, most of the currently developed continuum-based models also pose theoretical limitations and misperceptions in the way they incorporate directional dependency into the mechanical behaviour of the AV. To characterize a suitable class of anisotropy for the AV, an appropriate set of experimental stress–strain data (where the components of stress or strain can be independently controlled from one another) is needed. Because the thickness of the AV is much smaller than the other two in-plane dimensions (by two orders of magnitude), the AV is considered as a planar tissue [7,17]. This inherent geometrical characteristic of the AV implies that, from an experimental point of view, it would be very difficult, if not entirely impractical, to achieve stress–deformation data along the third dimension of the tissue, i.e. through its thickness, hence one is confined to only the in-plane dataset. Because the currently available uniaxial and biaxial material testing machines do not allow the independent control of in-plane shear from tensile deformation, both biaxial and uniaxial tensile tests only facilitate two components of strain/stress to be independently varied, which in turn will only sanction characterization of two independent components of the strain energy function *W*, in the form of two partial derivatives of *W* with respect to strain invariants (∂*W*/∂*I*). A rigorous analysis on this point, and the extent of suitability of the in-plane datasets, has been carried out and presented by Ogden & Holzapfel [18,19]. In this specific context, biaxial loading conditions, while more closely resembling the in-plane deformation experienced by the AV *in vivo*, do not offer any specific advantages over uniaxial datasets, which can also deliver stress–strain data on each of the two principal loading directions and provide enough data to characterize two independent ∂*W*/∂*I* terms [20,21]. However, whether the experimental datasets are achieved through uniaxial or biaxial tests, characterization of two independent ∂*W*/∂*I* terms will at best permit formulating a ‘transversely isotropic’ behaviour for inclusion into a mathematical continuum-based model [18–21]. This important point has often been overlooked in the literature and requires revisiting in formulating a new mathematical model for AV tissue. We further note that transverse isotropy is also structurally motivated in the case of AV tissue, as collagen fibres primarily assume a preferred direction along the circumferential axis (figure 1).

In this study, we derive a transversely isotropic viscoelastic model for application to the AV, incorporating the deformation rate as an explicit variable. The considered rate effects are reflected in the form of viscous damping *η* and are motivated by the dissipative effects of the valve's matrix, encompassing the viscous-like behaviour of GAGs and fibre kinematics. We start by demonstrating the general three-dimensional theory, and then apply it to uniaxial loading. Uniaxial tests were performed in circumferential and radial directions, under various deformation rates, to allow characterization of rate dependency as well as anisotropy. Appropriate mathematical constraints to comply with the condition of convexity are derived and verified against, to ensure appropriate parameter estimation. Based on the modelling results, a nonlinear relationship between the viscous damping effects and the deformation rate

## 2. Continuum mechanics framework

### 2.1. Preliminaries

Following Pioletti *et al*. [22], the second Piola–Kirchhoff stress tensor **S** for a viscoelastic material undergoing large deformations, with strain rate as an explicit variable, may be expressed as
**C** is the right Cauchy–Green tensor, which is related to the deformation gradient tensor **F** by
**C**.

In the presence of viscous effects, the stress tensor **S** in equation (2.1) may be derived [22] as
*W*_{e} and *W*_{v} are referred to as the elastic strain and the viscous dissipation energy functions, respectively.

For an incompressible material, equation (2.3) is replaced by
*p* is the arbitrary Lagrange multiplier, enforcing the constraint of incompressibility.

### 2.2. Material symmetry

Collagen fibres are the main load-bearing component of the AV extracellular matrix, providing reinforcement to the tissue. The fibres are embedded within a gel-like viscous ground substance formed of GAGs. This structure renders the AV tissue analogous to fibre composite materials [2]. The mean fibre orientation within this composite is identified by a single direction [3,23], referred to as the ‘circumferential’ direction. The transverse direction is referred to as the ‘radial’ direction. Together, these two directions are known as the principal loading directions of the AV, shown in figure 1 in relation to a single valve leaflet. This single preferred fibre direction endows the valve with pronounced directional-dependent mechanical properties between the two principal transverse directions, i.e. transverse isotropy. Accordingly, for formulating an appropriate continuum-based model, we shall specialize the class of anisotropy to transverse isotropy, with the preferred direction of the fibres aligned along the circumferential direction. We note that AV tissue is morphologically composed of three layers, namely the fibrosa, spongiosa and ventricularis, where collagen fibres are localized within the fibrosa and ventricularis layers. While the distribution of collagen fibre orientation within those layers reflects a certain degree of dispersity, small angle light scattering studies have established that the mean preferred fibre direction within the AV tissue is predominantly along the circumferential direction [14,23]. Therefore, in our approach, we treat the tissue macroscopically as a monolayer, with the global preferred direction of fibres along the circumferential loading direction of the tissue. For a detailed analysis on how to incorporate fibre orientation dispersion into mathematical models, the interested reader may wish to read the contributions made by Freed *et al*. [24], Gasser *et al*. [25] and Holzapfel *et al*. [26].

### 2.3. Energy functions *W*

Elastic and viscous potential functions appearing in equation (2.4) may be described by (**C**) and

In the case of transverse isotropy, *W*_{e} may be expressed as a function of 5 invariants [18], and *W*_{v} as a function of 17 invariants [28]:
*J*_{1}–*J*_{12} are the invariants of **M** in equations (2.7) and (2.8) denotes the preferred fibre direction given by **M **=** **[cos** ***φ*,sin** ***φ*,0]^{T} and is depicted schematically in figure 2*a*.

### 2.4. Viscoelastic stress tensors **S** and *σ*

The second Piola–Kirchhoff stress tensor **S** in equation (2.4) may be rewritten as

We note that, from matrix calculus,

In order to develop equation (2.9) further, one needs to establish the expressions for ∂*I _{i}*/∂

**C**and

**I**denotes the identity matrix and ⊗ is the tensor product.

Substituting the expressions from equation (2.10) into equation (2.9) yields
*W*_{e})* _{i}* and (

*W*

_{v})

*have been adopted to represent ∂*

_{i}*W*

_{e}/∂

*I*and ∂

_{i}*W*

_{v}/∂

*J*, respectively.

_{i}Under the assumption of incompressibility, *I*_{3 }= det **C**=1, and therefore, equation (2.11) can be slightly simplified to

The Cauchy stress **σ**, also known as the true stress, is obtained from **σ** = **FSF**^{T}, and, in view of equation (2.12), may be expressed as

Equations (2.12) and (2.13) present the generic relationships for the second Piola–Kirchhoff stress tensor **S** and the Cauchy stress tensor **σ**, respectively, as a function of the right Cauchy–Green tensor **C** for a transversely isotropic incompressible viscoelastic continuum, based on equation (2.4). It must be noted that equation (2.13) contains similar expressions to that provided by Ogden for transversely isotropic elastic materials, formulated in relation to the left Cauchy–Green tensor **B** [18].

In the following, we shall consider appropriate assumptions and conditions that best describe the deformation of the AV specimens in mechanical tensile tests, in order to specialize equation (2.13) for suitable application to experimental data and estimation of the material parameters.

### 2.5. Pure homogenous deformation

When the principal axes of deformation coincide with the Cartesian coordinate directions, and the principal stretches *λ*_{1}, *λ*_{2} and *λ*_{3} are independent of the coordinates (say *x*, *y* and *z*), the deformation is said to be a pure homogeneous deformation [18]. This is often the case in biaxial and uniaxial tensile deformation tests of the AV, as the specimens are prepared such that the circumferential and radial directions are often in line with the *x-* and *y*-directions of the Cartesian coordinate system. In this case, the components of the deformation gradient tensor **F** have a diagonal form diag[*λ*_{1},*λ*_{2},*λ*_{3}]. It therefore follows

Substituting the expressions given in (2.14) into (2.13), and noting that
*σ*_{13} = *σ*_{23} = 0.

#### 2.5.1. Point of caution

The principle of conservation of angular momentum for a continuum in static equilibrium enforces symmetry upon the Cauchy stress tensor; that is, *σ*_{12} = *σ*_{21}. We note that in equation (2.15) this is achieved only if *W*_{e} and *W*_{v} functions, defined by the invariants in equations (2.7) and (2.8), to comply with the principle of conservation of angular momentum and hence be symmetrical, the function *W*_{v} must be such that the above relationship between (*W*_{v})_{10} and (*W*_{v})_{12} holds. This point has been overlooked in the literature concerning anisotropic viscoelastic constitutive models developed for application to soft tissues. Models that are derived based on theoretical criteria that do not ascertain the symmetry of the Cauchy stress tensor inevitably describe unrealistic and infeasible stress–deformation relationships and subsequently result in erroneous parameter estimations.

In the light of the interrelationship between (*W*_{v})_{10} and (*W*_{v})_{12}, the components of the Cauchy stress tensor given in equation (2.15) may be presented as
*σ*_{13} = *σ*_{23} = 0. Additionally, owing to incompressibility, *λ*_{3} = 1/*λ*_{1}*λ*_{2}; however, for simplicity, we leave *λ*_{3} as it is shown in the expressions (2.16). The expressions in equation (2.16) describe the Cauchy stresses in a general case. We note that, considering only the elastic contribution, these expressions are similar to those presented by Ogden in [18], equations (65)–(68).

### 2.6. Application to biaxial tensile deformation

Biaxial tensile tests characterizing the mechanical behaviour of the AV tissue overwhelmingly use square specimens cut from the central region of the cusp (e.g. [4,5], as shown in figure 2*b*). Subsequently, the specimens have been considered as thin sheet ‘membranes’, and therefore appropriate ensuing assumptions are applied to model the experimental data using mathematical expressions. One such assumption is that, for a thin sheet membrane, the through-thickness (principal) Cauchy stress can be approximated to zero, *σ*_{33} = 0. The expressions in equation (2.16) therefore may be reduced to

We note that the hydrostatic pressure *p* can now be determined from *σ*_{33} = 0. Alternatively, following Ogden [18], the expressions in equation (2.16) may be rewritten as
*σ*_{33} = 0.

Notwithstanding the viscous terms in the expressions in equations (2.17) and (2.18), (i.e. terms that include (*W*_{v})* _{i}*), equations (2.18) render three independent components of deformation and stress, while containing four (

*W*

_{e})

*terms. Therefore, four equations are required to characterize the*

_{i}*W*

_{e}function, and thereby the properties of the AV tissue, while biaxial tensile tests at best could provide information regarding three independent deformation–stress sets. This problem has been discussed and analysed at length by Holzapfel & Ogden [18,19]. Equations (2.17) and (2.18) suggest that this problem is further exacerbated by inclusion of viscous terms. From this perspective, therefore, biaxial tests do not have much advantage over other loading protocols that may render lower ranks of datasets than the number of unknowns [20], particularly in characterizing the anisotropic viscoelastic behaviour of soft tissues such as the AV.

It must be further noted that experimental systems that could facilitate independent control of in-plane shear have not yet been introduced and employed in performing the tensile deformation tests on AV tissue specimens, as reflected in recent reviews of the state of the art [7,29,30]. Moreover, in the light of equation (2.18), if the preferred fibre direction is along one of the coordinate axes, i.e. *φ* = 0 or *φ* = *π*/2, then *σ*_{12} = *σ*_{21} = 0. This is often the case in square and/or rectangular specimens used in biaxial and uniaxial deformation tests of the AV, prepared from the valve cusps as shown in figure 2. Therefore, expressions in equation (2.18) may be reduced to

### 2.7. Application to uniaxial tensile deformation

In uniaxial tests, rectangular strips are cut from the central region of the valve cusp along the preferred fibre direction (circumferential), and the transverse direction (radial), as shown schematically in figure 2*b*. For a circumferential strip under uniaxial tensile deformation along that direction, *σ*_{22} = 0 and *φ* = 0°, as shown in figure 2*c*. Therefore, the principal Cauchy stress in circumferential direction is established from equation (2.19) as
*λ*_{1}·*λ*_{2}·*λ*_{3} = 1 and *λ*_{2} = *λ*_{3}. Substituting these into the above, equation (2.20) may be rewritten as

Similarly, for a radial strip under uniaxial tensile deformation along that direction, *σ*_{11} = 0, *φ* = 90° as shown in figure 2*c*, and *λ*_{1} = *λ*_{3}. Furthermore, because the fibres are aligned in the circumferential direction, we expect negligible contribution from (*W*_{e})_{4} and (*W*_{e})_{5} when the strip is stretched in the radial direction, as the fibres do not support compression [18], and the same may be assumed for the contribution of (*W*_{v})_{4}, (*W*_{v})_{5}, (*W*_{v})_{10} and (*W*_{v})_{11}. The principal Cauchy stress in the radial direction is therefore established from equation (2.19) as

Equations (2.21) and (2.22) express the principal Cauchy stress components in the circumferential and radial directions, respectively, under uniaxial tensile deformation. However, we note that uniaxial tensile tests provide two independent stress–deformation datasets, whereas equations (2.21) and (2.22) include 15 unknowns, four (*W*_{e})* _{i}*s and 11 (

*W*

_{v})

*s. Therefore, reasonable assumptions and appropriate forms of energy functions shall be considered for specialization of equations (2.21) and (2.22), to enable formulation of admissible models that can suitably describe the stress–deformation data. Such considerations include*

_{i}*a priori*assumptions regarding the appropriate number of invariants in energy functions, as well as certain physical and mathematical conditions which ensure model validity and material stability. In the following, we shall invoke these considerations and introduce the energy functions in mathematical form.

## 3. Model formulation

### 3.1. *W*_{e} and *W*_{v} functions

As equations (2.19), (2.21) and (2.22) indicate, biaxial tensile tests in which only two strain components are varied independently, and uniaxial tensile tests in transverse directions, inherently do not provide enough independent datasets for a complete characterization of energy functions *W* in transversely isotropic viscoelastic tissues. This is mathematically inferred, as the number of constitutive functions (*W*_{e})* _{i}* and (

*W*

_{v})

*supersedes the number of existing relationships between stress and deformation. In such cases, it is admissible to assume*

_{i}*a priori*that

*W*

_{e}and

*W*

_{v}are a function of only certain invariants, i.e. some of the invariants are considered absent from the general form of the energy functions [19,22,31,32]. Therefore, one may be faced with the task of choosing an appropriate form of

*W*

_{e}and

*W*

_{v}.

Standardized theoretical frameworks that facilitate axiomatic choices of elastic and/or viscous energy functions have not been articulated in the literature concerning soft tissues, if, indeed they are possible to develop. For elastic energy functions, Ogden [18] advocates three baseline factors that provide a sound reference for a valid starting point. First, *W*_{e} must be chosen, so that the ensuing stress–deformation relationships are consistent with the experimentally observed behaviour of the subject tissue. For example, most collagenous soft tissues exhibit an initial ‘soft’ stress–deformation phase, followed by a stiffening phase at higher deformations. An appropriate *W*_{e} should therefore accommodate this nonlinear stress–deformation behaviour. Second, *W*_{e} must reflect the relevant material symmetry of the subject tissue. For example, if a tissue is transversely isotropic, an appropriate *W*_{e} function for that tissue must reflect this characteristic material symmetry. Third, *W*_{e} must satisfy the condition of ellipticity and convexity, in order to furnish well-posed boundary-value problems and material stability. For the viscous energy function, thermodynamic requirements enforce *W*_{v} to be continuous, positive and convex with respect to *W*_{v} must be zero when

Taking the above considerations into account, a widely acceptable elastic energy function *W*_{e} for incompressible transversely isotropic tissues has been devised to depend only on invariants *I*_{1} and *I*_{4}, of the following form [19,21,31]:
*α* and *k*_{1} are positive stress-like material parameters, and *β* and *k*_{2} are positive dimensionless parameters. We note that there are now only two invariants incorporated in the *W*_{e} function, which, given the fact that biaxial and uniaxial tensile tests in transverse directions provide two independent stress–deformation equations, in principle should enable one to characterize the elastic behaviour of the valve if the elastic response of the tissue specimens is established from the experiments.

For devising an appropriate viscous energy function, we note that the viscous effects of the bulk AV tissue may stem from the gel-like viscous GAG matrix as outlined in §2.2, in addition to the dissipative kinematics of the fibre–matrix and the fibre–fibre sliding and interaction [6]. Therefore, we consider the overall viscous energy function *W*_{v} of the valve as the sum of the contribution of the valve's viscous matrix *et al*. [22], we choose *I*_{1} and *J*_{2}, and assume *I*_{1} and *J*_{5}, for the viscous energy function *W*_{v} to have the following form:
*η*_{1} and *η*_{2} are viscosity-like parameters reflecting the dissipative effects of the viscous matrix and the fibre kinematics, respectively, and are positive. We note that according to the definition of *J*_{2} and *J*_{5} given in equation (2.8), *W*_{v} is a quadratic function of *W*_{v} = 0 when *J*_{2} and *J*_{5}) to the stress–deformation equations. However, the only constitutive component of *W*_{v} in equation (3.3) that appears in the stress–deformation equation in the radial direction (equation (2.22)) is (*W*_{v})_{2}. Therefore, theoretically, (*W*_{v})_{2} may be characterized using an additional set of stress–deformation data obtained from tensile tests performed under a different strain rate compared with that of the elastic response, in the radial direction. Then, the stress–deformation equation in the circumferential direction (equation (2.21)) facilitates the characterization of (*W*_{v})_{5} using a set of stress–deformation data obtained from tensile tests performed in the same direction but under a different strain rate compared with that of the elastic response. Therefore, from a theoretical point of view, stress–deformation curves obtained from AV specimens under various deformation rates in transverse directions should, in principle, allow characterization of the viscous energy function *W*_{v}.

### 3.2. Transversely isotropic viscoelastic model

Equations (3.2) and (3.3) may now be inserted into equation (2.19) to develop a model to describe the stress–deformation behaviour of the AV under biaxial tension or similarly into equations (2.21) and (2.22) for a uniaxial model describing the deformation in transverse directions. For the purpose of this study, we have performed uniaxial tensile tests on AV specimens. We shall therefore employ equations (2.21) and (2.22) to develop a transversely isotropic viscoelastic model applicable to the uniaxial data:

## 4. Tensile deformation tests

For the purpose of this study, we used experimental stress–deformation data of AV specimens subjected to uniaxial tensile tests to failure in both circumferential and radial directions. The data were obtained under four stretch rates ^{−1}. Porcine hearts (*n* = 10 in total) were obtained from mature animals, ranging from 18 to 24 months, within 2 h of slaughter from a local abattoir. The three AV leaflets were dissected from the aortic root and maintained in Dulbecco's modified Eagle's medium (DMEM, Sigma, Poole, UK) at room temperature (25°C). From each leaflet, a 5 mm wide circumferential or radial strip was excised from the central (belly) region (figure 2*b*).

The tensile tests to failure were performed using two material testing machines, a Bionix 100 (MTS, Cirencester, UK) for tests under

The stress–deformation curves of the AV specimens obtained at *λ* and the nominal ‘engineering’ stress **P**. To enable the application of the experimental data to the developed model in equation (3.3), one must first convert the engineering stress **P** to Cauchy stress **σ** via **σ** = **P***λ* [33]. The resulting *σ* − *λ* curves under the corresponding stretch rates are shown in figure 3, for representative circumferentially and radially loaded samples. The data highlight increasing sample stiffness with increasing

## 5. Procedure for model application and parameter estimation

### 5.1. Rate dependency of the parameters

Equation (2.3) associates the overall stress in a viscoelastic continuum to the superposition of the elastic and the viscous contributions of its constituents. The premise of elasticity requires the elastic response of the continuum to be independent of the deformation rate. The viscous effects, by contrast, are dependent on the rate. Therefore, when the model is fitted to the stress–deformation curves obtained at various deformation rates, the parameters related to the elastic behaviour are to remain unchanged, whereas the viscous-related parameters reflect the rate dependency and alter at each rate. To this end, it is important to experimentally establish the elastic response of the tissue, i.e. the elastic stress–deformation curve, from which the associated elastic parameters of the model may be derived. Those parameters are then set to remain unchanged, while fitting the whole model to the stress–deformation curves obtained at different rates, to characterize the viscous-related parameters. The flowchart in figure 4 illustrates this procedure.

### 5.2. Elastic response

It is perhaps impractical to characterize and obtain a ‘pure’ elastic response from tissue samples that are inherently viscoelastic. The stress–strain curves of viscoelastic tissues often exhibit a marked rate dependency, particularly in the case of AV [6]. Pioletti & Rakotomanana [34] postulate that in these circumstances the choice of elastic curve is a matter of definition and identify the curves obtained at lower rates as the elastic response. We qualify this definition further by countenancing the role of characteristic time. Stress–relaxation tests enable the quantification of the characteristic times *τ*, fast and slow, whereby 99% of the relaxation fades within a time *t* = 5*τ* [9]. Therefore, if the deformation rate of the tensile test of tissue specimens is chosen sufficiently low to allow enough time for the viscous processes to take effect and fade, the ensuing stress–deformation curve may be deemed intractable to further reductive viscous effects. Such a curve may therefore provide a baseline that, with a degree of tolerance, may be referred to as the elastic curve.

To apply this definition to our AV samples, we take the average slow relaxation time as reported previously [9] to be *τ*_{circum} = 81.14 s and *τ*_{rad} = 32.11 s, in the circumferential and radial directions respectively. We note that by allowing enough time for the slow relaxation to take effect, the fast relaxation process has taken effect and completed *a priori*. Given that the maximum failure elongation of the specimens is reported to be *λ*_{circum} = 1.46 and *λ*_{rad} = 1.89 [6], a sufficiently low elongation rate to allow for the slow relaxation may be approximated as

### 5.3. Convexity

The common approach to characterize the model parameters *α*, *β*, *k*_{1}, *k*_{2}, *η*_{1} and *η*_{2} is to fit the model in equation (3.4) to the experimental data obtained from the tensile deformation tests described above. In the process of fitting, however, due care must be observed to ensure that the achieved parameter values do not result in undesirable material instabilities or implausible physical behaviours. Therefore, appropriate mechanically and mathematically motivated constraints need to be derived and applied to restrict the solutions to a meaningful domain in the ‘parameter space’. These constraints are often expressed in the form of mathematical inequalities, imposed on the constitutive parameters during the process of fitting. Following Holzapfel & Ogden, we employ the condition of strict local convexity in order to obtain the relevant inequalities [18,21].

From a mathematical point of view, the condition of strict local convexity requires that the matrix containing the second derivatives of the energy function *W* with respect to the Green–Lagrange strain tensor (**E**), or alternatively with respect to the principal stretches (*λ*), to be positive definite [18,21]. This matrix, also known as the Hessian of *W*, in terms of *λ* is presented by
**H** in equation (5.1) represents the Hessian matrix and *W* represents the total energy function, i.e. *W*_{e} + *W*_{v}. In view of equations (3.2) and (3.3),
*H* is a symmetric matrix. Therefore, in order for **H** to be positive definite, it is necessary that **H**) and the eigenvalues of **H** are all positive. Thus

These inequalities indicate that the material parameters cannot be chosen arbitrarily, and ensure the strict local convexity of *W*. In the light of equation (2.7) and the constraint that *φ* can assume either 0° or 90°, the inequalities in equation (5.2) may be further simplified to elicit *k*_{2}, *β* > 0. No further explicit interrelationships between the model parameters or their numerical range may be directly elucidated from equation (5.2). However, this equation is used to check whether the parameters obtained by fitting the model in equation (3.4) to the experimental data satisfy the inequalities. Graphical representation of the condition of strict local convexity reflects itself in the convexity of the projections of the contours of constant *W* in the (*λ*_{1},*λ*_{2}) and (*E*_{11},*E*_{22}) planes, and will be presented in §6. For more in-depth analysis on conditions of convexity, the interested reader is referred to Holzapfel *et al*. [35] and Balzani *et al*. [36].

### 5.4. Fitting procedure

The model in equation (3.4) was fitted to the experimental data by the curve fitting toolbox in MATLAB^{®}, using the Levenberg–Marquardt algorithm. The flowchart in figure 4 illustrates the different steps and sequences in estimating the model parameters, using the uniaxial stress–deformation curves in the transverse directions, obtained under different deformation rates. The first phase includes the estimation of the elastic parameters of the model, namely *α*, *β*, *k*_{1} and *k*_{2}. In this phase, the elastic terms in *σ*_{circumferential} and *σ*_{radial} in equation (3.4) are fitted to the experimentally obtained elastic curves *α* and *β* are related to the ‘isotropic’ matrix, and shall therefore assume the same numerical values in both directions. With this consideration, the best fit in both the circumferential and radial directions is sought and the numerical values of the elastic parameters *α*, *β*, *k*_{1} and *k*_{2} are estimated. Once the parameters are verified to result in a convex *W*_{e}, their numerical values are taken as the output of the fitting procedure for the elastic curves.

In the next phase, the values of *α*, *β*, *k*_{1} and *k*_{2} are incorporated into the model in equation (3.4) and are set fixed. The model is then fitted to the stress–deformation curves in circumferential and radial directions obtained under each considered deformation rate. As equation (3.3) indicates, *η*_{1} is related to the viscous properties of the matrix and shall therefore assume the same numerical value in both directions at each *η*_{1} and *η*_{2} at each *W* is then verified graphically by plotting *W* and its contours in (*λ*_{1},*λ*_{2}) and (*E*_{11},*E*_{22}) planes (*E*_{11} and *E*_{22} represent the principal Green–Lagrange strains).

Taken together, this procedure enables quantification of the parameters of the transversely isotropic viscoelastic model in equation (3.4), describing the mechanical behaviour of the AV using uniaxial stress–deformation data in transverse directions.

## 6. Results

Following the procedure described in §5, the model in equation (3.4) was fitted to the experimentally obtained *σ* − *λ* data. The model adequately captured the deformation behaviour of the specimens at each corresponding rate, reporting *R*^{2} values more than 0.97. Representative curves for both loading directions at each

The characterized model parameters are summarized in table 1, presented as mean ± s.e. Parameters *α*, *β*, *k*_{1} and *k*_{2} are the ‘elastic’ parameters and by definition are independent of the deformation rate; *η*_{1} and *η*_{2} are the parameters related to the viscous behaviour of the continuum, and therefore are rate-dependent. Note that *α* and *β* are the elastic parameters associated with the ‘isotropic’ matrix, and therefore are also independent of the direction of loading. These characteristics are reflected in the numerical values of the parameters listed in table 1.

The plots for *W*_{e} and its contours in the (*λ*_{1},*λ*_{2}) and (*E*_{11},*E*_{22}) planes, constructed using the values given in table 1, are shown in figure 6, for the circumferential (figure 6*a*) and the radial (figure 6*b*) directions, confirming the convexity of the elastic energy function. The numerical values of the model parameters listed in table 1 all report convex total energy functions. Figure 6*c* illustrates those energy functions in the circumferential loading direction. Contours of *W* in the circumferential direction for

## 7. Discussion

A new transversely isotropic viscoelastic model was developed and presented in this paper to describe the behaviour of the AV tissue under uniaxial tensile deformation. The model accounts for the rate effects, by incorporating the rate of deformation *η*, formulated as a function of **C** and

Embedded within the final form of the model in equation (3.4) are the assumptions of ‘thin sheet membrane’ and ‘pure homogenous deformation’. It is therefore important to note that equation (3.4) may not be suitable for application to situations where through-thickness or deformation under shear analyses are required. For further discussions on the scope of application of such continuum-based models, the interested reader is referred to the contributions made by Holzapfel and co-workers [20,21]. We further note that uniaxial tensile data alone are not sufficient for complete characterization of the multidimensional behaviour of the AV tissue.

As customary in studies modelling the biomechanical behaviour of soft tissues, model parameters and material properties are obtained in relation to experimental data, which itself is affected by the specimen samples and the experimental set-up. The AV tissue, similar to other biological soft tissues, is subject to sample variability. The deformation tests themselves are also affected by the properties of the experimental set-up such as the gripping mechanism, shape and size of the samples, and alignment of the gripped specimens, to name but a few. In particular, it has been previously shown that for radially cut specimens the gripping effects may compound the observed stress–strain behaviour of the samples [37,38], as the characteristic decay length may be the same order of magnitude as the gauge length of the samples (10 mm). Therefore, a degree of tolerance and caution has to be afforded to the numerical values of the model parameters reported in this study, and indeed in any such studies. Nonetheless, the detailed mathematical basis upon which our model was developed, together with the rigorous experimental campaign employed in this study, provides a solid basis for better understanding of the material properties of the AV and modelling its biomechanical behaviour.

Based on the modelling results summarized in table 1, a reduction is observed in the numerical values of *η _{1}* and

*η*

_{2}with increase in

*η*

_{1}and

*η*

_{2}are parameters associated with the dissipative effects of the viscous matrix and the fibre kinematics, respectively. The reduction in the values of

*η*with

*η*

_{1}values we have reported with an increase in

*η*

_{2}, however, may be less obvious. As

*η*

_{2}is associated with the dissipative effects of fibre kinematics, e.g. fibre sliding or reorientation, the reduction in

*η*

_{2}values with increase in

To estimate the stress–strain behaviour of the AV tissue at physiological rates, where ^{−1} as reported in [16]), the values of *η*_{1} and *η*_{2} at that rate shall be extrapolated from the available data. The graph in figure 7*a* shows the variation of *η*_{1} and *η*_{2} with *η*,*η*_{1} and *η*_{2} at ^{−1}, respectively (figure 7). Incorporating these values into equation (3.4), the predicted *σ* − *λ* curves of the AV tissue at physiological loading rates are illustrated in figure 7*b* for both loading directions. To the best of the authors' knowledge, this study presents the first prediction of *σ* − *λ* curves for the AV in the principal loading directions at physiological loading rates using a continuum-based model incorporating the deformation rate as an explicit variable. We note that some experimental data exist in the literature in relation to the peak stretches experienced by the AV *in vivo*, measured at the systolic and diastolic phases of the cardiac cycle [39,40]. Owing to experimental limitations, complete stress–strain curves were not established. Nonetheless, the reported values for stress at maximum diastolic stretches of approximately 1.13 in the circumferential direction (approx. 2.9 MPa) [39] correspond well with that from our predicated circumferential *σ* − *λ* curve (approx. 2.6 MPa; figure 7*b*). There is also some literature reporting membrane tension versus stretch curves for porcine AV specimens, obtained at rates close to physiological rates under equi-biaxial loading conditions *in vitro* [8]. However, these loading conditions appear to result in far larger circumferential and radial stretches than those measured *in vivo* [39,40], and drawing a direct relevance between those data and the mechanical behaviour of the AV at physiological rates, or our reported *σ* − *λ* curves, may be problematic.

To verify the reliability of the estimated values of *η*_{1} and *η*_{2}, we used our newly established equation of the line of best fit to calculate the values of *η*_{1} and *η*_{2} at *α*, *β*, *k*_{1} and *k*_{2} in table 1, we used the model to predict the *σ* − *λ* curves at *σ* − *λ* curves in both loading directions. The graphs in figure 8 illustrate the degree of conformity between the model predictions and the experimental data, reporting *R*^{2} values more than 0.99. Based on this result, the model predictions for physiological *σ* − *λ* curves may be treated with a high degree of confidence.

While the model presented in this study was primarily developed for application to the AV, the mechanical and mathematical criteria within which the model was derived are general and universal. Therefore, the model and the modelling approach presented here can be applied to other collagenous soft tissues with similar structural building blocks and a single preferred direction of the embedded collagen fibres, without the loss of generality.

## Ethics

AV tissue samples were porcine in origin and all collected from a local abattoir. No ethics approval was therefore required.

## Data accessibility

Data for the reported *σ* − *λ* curves from a previous study [6] and the new datasets obtained for this study are available at University of Portsmouth research repository at: http://dx.doi.org/10.17029/754df175-2aa1-4859-993f-4cb775c7122d.

## Authors' contributions

A.A.-B. and H.R.C.S. were involved in conception and design of the work. A.A.-B. and A.Bu. carried out the experiments, and together with S.L.E. were involved in developing the mathematical model. A.A.-B. and A.Bu. drafted the article. H.R.C.S. and S.L.E. were involved with data analysis and interpretation, as well as with critically reviewing the article. All authors gave their final approval for publication.

## Competing interests

The authors have no competing interests.

## Funding

Funding for this study was provided by the School of Engineering, University of Portsmouth.

## Acknowledgement

The authors are grateful to Prof. Gerhard Holzapfel for stimulating discussions and his valuable suggestions on revising the manuscript.

- Received August 12, 2016.
- Accepted November 14, 2016.

- © 2017 The Authors.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.