A transverse isotropic viscoelastic constitutive model for aortic valve tissue

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 (λ˙). Model (material) parameters were determined from uniaxial tensile deformation tests of porcine AV specimens at various deformation rates, by fitting the model to each experimental dataset. It is shown that the model provides an excellent fit to the experimental data across all different rates and satisfies the condition of strict local convexity. Based on the fitting results, a nonlinear relationship between η and λ˙ is established, highlighting a ‘shear-thinning’ behaviour for the AV with increase in the deformation rate. Using the model and these outcomes, the stress–deformation curves of the AV tissue under physiological deformation rates in both the circumferential and radial directions are predicted and presented. To verify the predictive capabilities of the model, the stress–deformation curves of AV specimens at an intermediate deformation rate were estimated and validated against the experimental data at that rate, showing an excellent agreement. While the model is primarily developed for application to the AV, it may be applied without the loss of generality to other collagenous soft tissues possessing a similar structure, with a single preferred direction of embedded collagen fibres.


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][5][6], whereas distinctions also exist in the load-bearing capacity and the distensibility of the AV tissue in these two loading directions [5][6][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][10][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][13][14][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][19][20][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  The fibre structure assumes a preferred direction along the circumferential axis (adapted from [3]).
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 ratė λ is established, characterizing a 'shear-thinning' behaviour for the AV tissue. Using this relationship, the viscous effects and subsequently the stress-deformation curves of the AV at a physiological deformation rate in both circumferential and radial directions are predicted and presented.

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 S = S(C,Ċ), (2.1) where C is the right Cauchy-Green tensor, which is related to the deformation gradient tensor F by We note thatĊ is the time derivative of C.
In the presence of viscous effects, the stress tensor S in equation (2.1) may be derived [22] as where 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 where p is the arbitrary Lagrange multiplier, enforcing the constraint of incompressibility.

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].

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,  square or rectangular specimens are prepared from the central region of the AV leaflet for biaxial or uniaxial tensile tests, respectively. Note that, for uniaxial tests, samples are obtained from both circumferential and radial directions; (c) for circumferentially cut samples, the fibre family is aligned with the principal direction. For radially cut samples, the fibre family makes an angle of 90°with the principal direction.
In order to develop equation (2.9) further, one needs to establish the expressions for ∂I i /∂C and ∂J i /∂Ċ, as follows where I denotes the identity matrix and ⊗ is the tensor product.
where, for simplicity, notations (W e ) i and (W v ) i have been adopted to represent ∂W e /∂I i and ∂W v /∂J i , respectively. Under the assumption of incompressibility, I 3 = det C=1, and therefore, equation (2.11) can be slightly simplified to S = −pC −1 + 2(W e ) 1 I + 2(W e ) 2 (tr CI − C T ) + 2(W e ) 4 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.

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  11 cos 2 ϕ + 2λ 6 1 (W v ) 12 cos 2 ϕ, σ 12 = λ 1 λ 2 (W e ) 4 sin 2ϕ + λ 1 λ 2 (λ 2 1 + λ 2 2 )(W e ) 5 sin 2ϕ 12 sin 2ϕ, 12 sin 2ϕ, with σ 13 = σ 23 = 0.

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 v ) 12 = −(W v ) 10 /(λ 2 1 + λ 2 2 ). Therefore, in order for the Cauchy stress tensor of a transversely isotropic viscoelastic material constructed from 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 11 sin 2ϕ, 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).

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 2b). 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 11 sin 2ϕ and σ 22 = −p + 2λ 2 2 (W e ) 1 + 2λ 2 2 (λ 2 1 + λ 2 3 )(W e ) 2 + 2λ 2 2 (W e ) 4 sin 2 ϕ + 4λ 4 2 (W e ) 5 sin 2 ϕ 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 and where, again, for a thin sheet membrane, the through-thickness (principal) Cauchy stress can be approximated to zero σ 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 ) i terms. Therefore, four equations are required to characterize the 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 inplane 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

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 2b. For a circumferential strip under uniaxial tensile deformation along that direction, σ 22 = 0 and ϕ = 0°, as shown in figure 2c. Therefore, the principal Cauchy stress in circumferential direction is established from equation (2.19) as where λ 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°a s shown in figure 2c, 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 ) i 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 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.  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 stressdeformation 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Ċ [22]. In addition, the value of W v must be zero whenĊ is equal to zero [22].

Model formulation
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]: where W iso e and W fibres e represent the influence of the isotropic matrix and the mean preferred fibre direction, respectively, on the overall elastic behaviour of the AV. For the purpose of our model, we employ an exponential-type elastic energy function for the isotropic matrix W iso e (as advocated in [15]), and a 'Holzapfel-type' elastic energy function for the contribution of the embedded fibre family W fibres e [19,21,31]: where α 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 W matrix v and the dissipative kinetics of the fibres W fibres v . Following Pioletti et al. [22], we choose W matrix v to depend on the invariants I 1 and J 2 , and assume W fibres v to depend on the invariants I 1 and J 5 , for the viscous energy function W v to have the following form: where η 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 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 stressdeformation 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 .

Transversely isotropic viscoelastic model
and The expressions in equation (3.4) represent the final form of our transversely isotropic viscoelastic model, describing the stress-deformation behaviour of the AV leaflet, using uniaxial tension in transverse directions.

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λ of 0.001, 0.01, 0.1 and 0.5 s −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 2b). The tensile tests to failure were performed using two material testing machines, a Bionix 100 (MTS, Cirencester, UK) for tests underλ = 0.001 s −1 ,λ = 0.01 s −1 andλ = 0.1 s −1 , and a Bose Electroforce 3200 (Minnesota, USA) for tests aboveλ = 0.1 s −1 . The initial distance between the grips was set at 10 mm for all test protocols, after which a tare load of 0.01 N was applied to the specimens, to establish a consistent zero position. The adjusted distance between the grips was then used as the initial sample length. For each tensile test, three specimens were used. No preconditioning was applied to the specimens prior to the start of the tests.
The stress-deformation curves of the AV specimens obtained atλ = 0.001 s −1 ,λ = 0.01 s −1 andλ = 0.1 s −1 were reported in a previous study [6], and are reproduced here in conjunction with the new data collected specifically for this study, at a stretch rate ofλ = 0.5 s −1 . It must be noted that the experimental results obtained from the tensile tests provide data in terms of λ 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λ. However, the strain rate-associated stiffening appears to be rate-limited, especially in the radial direction, suggesting the data are approaching a thresholdλ whereby increasing the deformation rate would not significantly alter the deformation curves.
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.

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λ circum = 0.0011 s −1 andλ rad = 0.0055 s −1 , assuming a linear relationship between elongation and time. Therefore, an elongation rate ofλ = 0.001 s −1 would ensure achieving a stress-deformation curve intractable to additional reductive viscous effects, which we consider as the elastic curve.

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 We note that 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), which implies that H is a symmetric matrix. Therefore, in order for H to be positive definite, it is necessary that ∂ 2 W/∂λ 2 1 , ∂ 2 W/∂λ 2 2 , det (H) and the eigenvalues of H are all positive. Thus and eigenvalues(H) > 0 ⇒ both roots of the 'characteristic equation' must be positive.
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°o r 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].

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 (λ = 0.001 s −1 ) in each respective direction. As equation (3.2) indicates, parameters α 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λ. With this consideration, the best fit in both the circumferential and radial directions is sought, and the numerical values of η 1 and η 2 at eachλ are established. The convexity of the total energy function 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 stressdeformation data in transverse directions.

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λ are presented in figure 5. The continuous line represents the model predictions.
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 6a) and the radial (figure 6b) 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 6c illustrates

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λ as an explicit parameter. The rate effects were considered in the form of viscous damping η, formulated as a function of C andĊ invariants. We therefore note that the model is applicable to monotonic proportional loading conditions, i.e. tensile deformations, and may not be suitable for directly modelling relaxation or creep behaviours. While applied to uniaxial data, the model was developed within the general three-dimensional context, with the appropriate mathematical and mechanical conditions introduced at each step to tailor the model for application to uniaxial tensile data. 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λ. In rheological terms, this behaviour is referred to as a 'shearthinning' behaviour. We recall that η 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λ in the AV has been observed and reported previously [2,6], and associated with shear-thinning behaviour of the GAG constituents of the valve. Therefore, the reduction in η 1 values we have reported with an increase inλ may be interpreted as the reflection of the shear-thinning behaviour of the GAGs. The physical interpretation of the behaviour of η 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λ may stem from a decrease in these dissipative effects at higher deformation rates. At higherλ, less time is afforded to the fibres to slide against each other, or to return back to their initial orientation during the deformation, compared with the case at lowerλ. Further microstructural evidence concerning fibre kinematics within the AV tissue during deformation and relaxation/creep is required to reach a more concrete conclusion.
To estimate the stress-strain behaviour of the AV tissue at physiological rates, whereλ = 2.5 s −1 (corresponding to the physiological strain rate of 15 000% min −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 7a shows the variation of η 1 and η 2 withλ, as given in table 1, in a logarithmic scale. Using an allometric function to fit to the (η,λ) data points, the values of η 1 and η 2 atλ = 2.5 s −1 are calculated to be 0.15 and 5.54 MPa s −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 7b 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 7b). 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λ = 0.2 s −1 . Using these values, and the values for α, β, k 1 and k 2 in table 1, we used the model to predict the σ − λ curves atλ = 0.2 s −1 in both the circumferential and radial directions. We then performed tensile tests at that rate and obtained the corresponding σ − λ 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.