Auxetic behaviour , the unphysical transverse expansion during uniaxial tension, is a com- mon and undesirable feature of classical anisotropic hyperelastic constitutive models for soft tissue. In this study we uncover the underlying mechanism of such behaviour; high levels of in-plane compaction occurs due to increasing tension in strain-stiffening ﬁbres, leading to unphysical out-of-plane expansion. We demonstrate that auxetic behaviour is primarily inﬂuenced by the ratio of ﬁbre to matrix stiffness, and is accentuated by strain- stiffening ﬁbres in a constant stiffness matrix (e.g., the widely used neo-Hookean matrix with exponentially stiffening ﬁbres). We propose a new bilinear strain stiffening ﬁbre and matrix (BLFM) model which allows close control of the ﬁbre-matrix stiffness ratio, thereby robustly eliminating auxetic behaviour. We demonstrate that our model provides accurate prediction of experimentally observed out-of-plane compaction, in addition to stress-stretch anisotropy, for arterial tissue subjected to uniaxial tension testing. Keywords: Hyperelasticity, Constitutive modeling, Anisotropy, Soft tissue, Auxetic behaviour, Strain-stiffening ﬁbre-matrix model 1. Introduction The hyperelastic strain energy density function used to model ﬁbrous soft tissues is typically decomposed into an anisotropic part (representing the contribution of aligned families of ﬁbres), and an isotropic part (representing the contribution of non-ﬁbrous ground matrix). Common examples include models of arterial tissue (Gasser et al., 2006; Nolan and McGarry, 2016a; Fereidoonnezhad et al., 2016; Holzapfel et al., 2015), carti- lage (Dowling et al., 2013; Brown et al., 2009; Pierce et al., 2013; Deneweth et al., 2015), and the annulus ﬁbrosis (Guo et al., 2006; Peng et al., 2006). The most commonly imple- Corresponding author Email address: patrick.mcgarry@nuigalway.ie (J. P. McGarry) Preprint submitted to Journal of Biomechanics July 13, 2020 mented form for the strain energy density function entails exponentially strain-stiffening ﬁbres in parallel with a near-linear neo-Hookean ground matrix, following from the work of Holzapfel, Gasser and Ogden (2000). While such an additive decomposition of soft tis- sue stress into a ﬁbre contribution and a ground matrix contribution is attractive in terms of physical interpretation of the constitutive law, auxetic behaviour (i.e. orthogonal expan- sion during uniaxial tension) is a common and undesired feature of such models (Skacel and Bursa, 2016, 2019; V olokh, 2017). To better illustrate such auxetic behaviour , we re- produce a simulation from the ‘ Abaqus benchmarks Guide ’ (input ﬁle adventitia_circ.inp ) of uniaxial tension of an adventitial strip cut along circumferential direction (Abaqus, 2017). The material is modelled using the well-known HGO constitutive model (Gasser et al., 2006), with ﬁbre families prescribed at 49:98, as shown in Figure 1A. As shown in Figure 1 thickness in the out-of-plane 3-direction increases as the specimen is stretched. While auxetic behaviour has been reported for certain classes of synthetic materials (e.g. foams and lattice materials (Lim, 2015)), such auxetic behaviour is not observed in ex- perimental tests of artery tissue (Skacel and Bursa, 2016). To the best of the authors knowledge, auxetic behaviour has not been reported in any naturally occurring soft tissue or tissue engineered collagen constructs. Therefore, in the context of soft tissue material models, we refer to the computation of auxetic behaviour as unphysical. Auxetic behaviour of soft tissue hyperlastic models has been the subject of three recent studies. Skacel and Bursa (2016) demonstrated the computation of auxetic behaviour (which they refer to as a negative Poisson’s ratio) for the anisotropic hyperelastic soft tissue models of Holzapfel et al. (2000) and Gasser et al. (2006; 2012). The authors incorrectly suggest that such unphysical auxetic response might be the result of model incompressibility. In a follow-on study the authors later demonstrate that introduction of volumetric compressibility into constitutive model in fact increases auxetic behaviour of the models (Skacel and Bursa (2019)). The recent study by V olokh (2017) considers an anisotropic hyperelastic model with exponential stiffening and ﬁbre dispersion. The author suggested that a more accurate consideration of the ﬁbre dispersion could be used to reduce the auxetic effect. Fibre dispersion is also taken into account in the models investigated by Skacel and Bursa (Skacel and Bursa, 2016, 2019). However, in Section 7 of the current study, we demonstrate that dispersion can merely reduce, but not eliminate, auxetic behaviour. Furthermore, ﬁbre disperion is a physically motivated model feature 2 Figure 1: Results of the Abaqus Benchmark Guide example Abaqus (2017) for uniaxial tension of an adventitial strip cut along the circumferential direction. (A) Two ﬁbre families oriented as=49:98are considered. (B) Increase in the thickness of the strip in the case of fully aligned ﬁbres ( = 0) atL L0= 0:33. (C) An increase in the thickness of the strip is computed forL L0>0:2in the case of fully aligned ﬁbres ( = 0) and for 0:09<L L0<0:3in the case of dispersed ﬁbres ( = 0:226) . The reader is refer to Abaqus (2017) for details of material parameters and dimensions. that must be calibrated to reproduce observed material anisotropy. It is therefore not an independent parameter that can be used to control auxetic behaviour. The three previous studies (Skacel and Bursa, 2016, 2019; V olokh, 2017) demon- strate that incorrect computation of auxetic behaviour is a general problem for established anisotropic hyperelastic soft tissue models. However, these studies do not uncover the underlying mechanistic cause of such auxetic behaviour, but merely focus on secondary factors such as compressiblity and ﬁbre dispersion. No study to date has uncovered the root cause of computed auxetic behaviour in anisotropic soft tissue models, and no study has proposed a robust solution for such incorrect behaviour. The current paper (i) provides, for the ﬁrst time, an explanation for the underlying 3 mechanism of auxetic behaviour in anisotropic soft tissue models, and (ii) proposes a new constitutive law for anisotropic soft tissue that robustly controls and eliminates auxetic behaviour. The paper is structured as follows. In Section 2 we provide a mechanistic explanation for the computation of auxetic behaviour in ﬁbre reinforced soft tissue mod- els. In Section 3 we demonstrate the key role of the ﬁbre-matrix stiffness ratio in auxetic behaviour. In Section 4 we demonstrate that strain stiffening ﬁbres in a constant stiffness matrix promotes auxetic behaviour at high strains. In Section 5 we propose a novel bi- linear strain stiffening ﬁbre and matrix model (BLFM model) and we demonstrate that this new formulation allows robust control of out-of-plane deformation and elimination of auxetic behaviour, resulting in accurate prediction of experimentally measured artery anisotropic strain stiffening and out-of-plane contraction. 2. Why does auxetic behaviour occur? In order to provide a clear explanation of the origin of auxetic behaviour for consti- tutive models of ﬁbrous tissue, we consider the simple case of a cube that consists of a matrix material and two families of embedded collagen ﬁbres. Following the typical ap- proach for artery models, we assume that the ﬁbres are conﬁned to the (1-2) plane and are symmetric about the 1-direction (see Figure 1A). Following the established approach for the modelling of hyperelastic soft tissue, the strain energy density function is additively decomposed into an isotropic matrix contribution and an anisotropic ﬁbre contribution: = m+ f (1) This leads to an additive decomposition of the Cauchy stress tensor: i=m i+f i (2) Therefore, when the cube is subjected to uniaxial tension in the 1-direction, the equilib- rium equations in the 2- and 3-directions can be written as: 2=m 2+f 2= 0; 3=m 3+f 3=m 3= 0(3) wherem 2andm 3are the 2- and 3-components of the matrix Cauchy stress tensor, re- spectively.f 2is the 2-component of the ﬁbre Cauchy stress tensor. Because we assume 4 that ﬁbres are conﬁned to the (1-2) plane without any out-of-plane contribution, f 3is equal to zero. Additionally, due to symmetry of the two ﬁbre families about the loading direction, the 1-, 2-, and 3-direction are the principal directions. Typically the constitutive law for ﬁbres assume signiﬁcant strain stiffening (e.g. exponential stiffening of collagen ﬁbres) whereas matrix constitutive laws typically assume much lower strain stiffening (e.g. neo-Hookean model without progressive stiffening). Therefore, when a high level of stretching is imposed on the cube in the 1-direction, the ﬁbres exert signiﬁcant tension in the 2-direction (in addition to a larger contribution in the 1-direction). Therefore, to achieve stress equilibrium, the near constant stiffness matrix must develop large negative strains in the 2-direction such that m 2=f 2. As a result of extreme matrix compaction in the 2-direction, the matrix may undergo expansion in the 3-direction in order to satisfy stress equilibrium. Such expansion in a lateral direction during uniaxial stretching is re- ferred to as auxetic behaviour. Such behaviour is not observed experimentally for arterial tissue (Skacel and Bursa, 2016), or, to our knowledge, any other soft tissue. Therefore, new modelling strategies to control and eliminate the computation of non-physical auxetic behaviour must be developed for anisotropic ﬁbrous soft tissue. 3. The simple case of linear ﬁbres in a linear matrix We begin our investigation with a simple parametric analysis of the sensitivity of auxetic behaviour to (i) ﬁbre to matrix stiffness ratio, (ii) ﬁbres orientation and (iii) matrix compressibility. In this Section we initially limit our analysis to linear ﬁbre and linear matrix behaviour in order to elucidate the factors underlying auxetic behaviour. Once again, a cube with two ﬁbre families conﬁned to the (1-2) plane is considered, with applied uniaxial extension in the 1-direction (Figure 2). 5 Figure 2: Schematic of the uniaxial tension of a ﬁbre-reinforced cube with two ﬁbre families. (a) undeformed cube (b) deformed cube without auxetic behaviour , and (c) deformed cube with auxetic behaviour . The ﬁbres are distributed in (1-2) plane. Hyperelastic constitutive models for arterial tissue and other ﬁbrous soft tissues are deﬁned through the strain energy density function as follows: (C;a0i) = m(C) +X i=4;6 fi(C;a0i) (4) m(C)and fi(C;a0i)are the strain energy density contributions of the ground matrix and collagen ﬁbres, respectively. Two ﬁbre families ( i= 4;6) are represented in equa- tion (4), as commonly assumed for arterial tissue. C=FTFis the right Cauchy–Green tensor, where Fis the deformation gradient tensor. Typically for soft tissue (e.g. arter- ies) two ﬁbre families are represented, and a0iare two ﬁbre directions in the reference conﬁguration. Typically m(C)is decomposed into isochoric and volumetric matrix con- tributions, such that m(C) = vol(J) + m(C): (5) C=J2=3Cis the isochoric modiﬁed right Cauchy–Green tensor, while J= det F>0is the volumetric jacobian. The simple neo-Hookean formulation is generally used to model the ground matrix of soft tissues, such that: m(C) =K 2(J1)2+ 2(I13) (6) 6 whereKis the effective bulk modulus and is the effective shear modulus. I1= trCis the ﬁrst invariant of C. This formulation provides near-linear stress-strain behaviour. It is convenient to deﬁne Em, the effective matrix stiffness under conditions of uniaxial stress, such thatEm= 9K=(3K+). Note that the volumetric contribution in equation (6) is appropriate for soft tissues that do not undergo large deformations (e.g. arteries (Nolan and McGarry, 2016b)). However, for materials that exhibit large volume changes the reader is referred to Moerman et al. (2020) for more complex volumetric strain energy density functions. For the initial investigation in this section we consider a simpliﬁed constant stiffness ﬁbre strain energy density function is given as fi=E1f(2 33 fi2 fi+1 3) (7) which leads to a linear relationship between ﬁbre Cauchy-stress and ﬁbre stretch, where E1fis the ﬁbre stiffness and fiis the ﬁbre stretch obtained as: 2 fi=Ii=a0i:Ca0i. 1 2 30.20.40.60.81 1 2 311.521 2 30.40.60.81 1 2 30.40.60.811.21.41 2 30.40.60.81 1 2 311.21.4B A D EC F Figure 3: Inﬂuence of the ﬁbre to matrix stiffness ratio (A, D), ﬁbre angle (B, E) and compress- ibility (C, F) on the in-plane (ﬁrst row) and out-of-plane stretch 3(Second row) for the uniaxial extension of a cubic element of ﬁbrous materials. A neo-Hookean matrix with a linear ﬁbre model with the reference parameters value of =45,Em= 0:5,Ef=Em= 10 andK= = 100 have been used in all plots (blue curves). 1is the stretch in the loading direction. The results of our initial parameter study for linear ﬁbres in a linear matrix are shown in Figure 3. Computed in-plane contractions 2(Figure 3A, C, E) and the corresponding 7 out-of-plane deformations 3(Figure 3B, D, F) are plotted as a function of the applied stretch in the x-direction 1.3>1indicates that unphysical auxetic behaviour has been computed. As shown in Figure 3 (A, D), an increase in the ﬁbre to matrix stiffness ratio results in in-plane contraction in the 2-direction, and a corresponding increase in out-of- plane expansion in the 3-direction. Signiﬁcant auxetic behaviour is computed for all cases in Figure 3, except for the case where the ﬁbre and matrix stiffnesses are equal (i.e. the material is essentially isotropic). In summary, stretching of stiff ﬁbres results in a high component of ﬁbre tension in the 2-direction. Therefore, mechanical equilibrium requires a signiﬁcant matrix compaction in the 2-direction, with a corresponding matrix expansion in the 3-direction. Figures 3 B and E show the effect of initial ﬁbre orientation (with respect to the loading direction (1-direction)). As expected, no auxetic behaviour is observed for the case of = 0; ﬁbres contribute only to the 11 component of the stress tensor, and the deformation in the 2- and 3-directions is dependent only on the matrix properties. Similarly for the case of= 10, ﬁbres are closely aligned to the loading direction and contribute weakly to the 22 component of the stress tensor. Therefore deformations in the 2- and 3- directions are primarily governed by the matrix properties and, again, no auxetic behaviour is computed. For= 30auxetic behaviour is computed from the onset of stretching. This occurs due to a high component of ﬁbre tension in the 2-direction. Stress equilibrium therefore results in high matrix compaction in the 2-direction, and consequently an auxetic expansion in the 3-direction. 3is observed to peak at an applied stretch of 12. Further stretching results in further rotation of the ﬁbres towards the loading direction and a consequent slight reduction in 3. A similar trend is observed for a ﬁbre angle of = 45, with3 peaking at a higher stretch of 12:8. For ﬁbre angles greater than = 45, the ﬁbres are in compression at the start of the deformation. Since ﬁbres are not capable of bearing compressive loads (McEvoy et al., 2018; Holzapfel and Ogden, 2017; Holzapfel and Fereidoonnezhad, 2017), during the ini- tial stages of stretching the material behaves like an isotropic matrix with no ﬁbre contri- bution and consequently no auxetic behaviour . However, as the deformation increases the ﬁbres rotate towards the stretching (1-) direction and develop a tensile 22-Cauchy stress component. Eventually auxetic behaviour begins to appear and progressively gets worse. The effect of material compressibility on auxetic behaviour is also investigated ( Figure 3 8 C, F). It is observed that auxetic behaviour increases with increasing material compress- ibility. Figure 3 presents a limited parametric investigation to illustrate the mechanisms underlying auxetic behaviour. It should be noted that the full range of interactions be- tween the governing parameters is not presented. For example, auxetic behaviour will be computed for = 10if a ﬁbre to matrix stiffness ratio ofEf Em= 30 is considered. 4. Auxetic behaviour is increased by strain–stiffening ﬁbres in a linear matrix The isotropic neo-Hookean model with a linear ﬁbre model has been used for the illus- tration of the auxetic behaviour and parameter studies in Section 2. However, experimen- tal data of ﬁbrous soft tissue exhibits anisotropic strain-stiffening behaviour . Typically this is captured using a strain-stiffening constitutive model for ﬁbres in parallel with a neo-Hookean matrix. In order to systematically analyse the role of ﬁbre strain stiffening on auxetic be- haviour we propose a new bilinear strain-stiffening ﬁbre (BLF) constitutive law. Fibre deformation is considered to consist of three distinct regimes: (i) an initial compliant quasi-linear region to model the stretching of coiled collagen ﬁbres; (ii) a non-linear tran- sition region to model the uncoiling of ﬁbres; (iii) a stiff quasi-linear region to model fully uncoiled ﬁbres. The proposed bilinear stress-strain relation for the collagen ﬁbres is given as: aniso=X i=4;6fiai ai; fi=8 >>>< >>>:E1f(fi1) fi1D1f pf(fi1)2+qf(fi1) +rf D1f<fi1<D 2f E2f(fi1D2f) + (pfD2 2f+qfD2f+rf)fi1D2f (8) where,E1fandE2fare slopes of the linear regimes, D1fandD2fare the values of the nominal ﬁbre strain at the end of the ﬁrst linear regime and at the beginning of the second linear regime, ai;i=4;6are the directions of the ﬁbres in deformed conﬁguration and  stands for the Cauchy stress. It is noted that the material parameters of the BLF model in equation (8) have a clear physical meaning that can be related to a soft tissue experimental stress-stretch curve. During the initial regime the collagen ﬁbres are tightly coiled. As the initial regime reaches D1fand enters the transition regime the collagen ﬁbres uncoil 9 which leads to the rapid stiffening of the stress-stretch curve. At D2f, the point where the second linear regime is entered, the collagen ﬁbres are completely uncoiled and behave as straight stiff elastic ﬁbres. Moreover, pf,qfandrfare not independent parameters; in order to maintain C0andC1continuity the following relations must be enforced: pf=E1fE2f 2(D1fD2f); qf=E1f2D1fpf; rf= (E1fqf)D1fpfD2 1f;(9) The corresponding ﬁbre strain energy density function is given as: aniso=X i=1;2 fi; fi=8 >>>< >>>:E1f(2 33 fi2 fi+1 3) fi1D1f 2 33 fi(qf2pf) +pf 24 fi+2 fi(pfqf+rf) + 01 D1f<fi1<D 2f 2E2f 33 fi+2 fi(pfD2 2f+qfD2f+rfE2fE2fD2f) + 02fi1D2f (10) where, 01and 02are two constants which ensure the continuity of strain energy. Figure 4 contains a schematic illustration of the material parameters of the proposed BLF model and functionality of stress and strain energy. An alternative form of the BLF model that uses a single function, rather than the piece-wise function, is presented in Appendix A. 10 Figure 4: A schematic illustration of a typical f-f(A) and f-f(B) curves for the proposed BLF model. The nature of the parameters E1f,D1f,E2fandD2f, and the uncoiling of the collagen ﬁbres during each deformation regime are illustrated. D1fis the nominal strain at the end of the ﬁrst linear regime which is controlled by the initial slope E1f.D2fis the nominal strain at the beginning of the second linear regime which is controlled by the slope E2f. The low- stiffness and the high-stiffness linear regime are connected by a transition regime between D1f andD2f. 1 2 3020406080 1 2 30.40.60.81 1 2 311.11.21.3 1 2 302040 1 2 30.40.60.81 1 2 311.11.2B C E D FA Figure 5: The inﬂuence of ﬁbre strain-stiffening and transition strain on the stress (A, D), in- plane transverse stretch (B, E) and out-of-plane stretch (C, F). Reference parameters E1f=Em= 1:5; = 45o; E2f=Em= 6; D1f= 0:2; D2f=D1f+ 0:2have been used. The BLF model is placed in parallel with a linear neo-Hookean matrix and subjected 11 to uniaxial tension boundary conditions. The inﬂuence of the ﬁbre strain stiffening in regimes (ii) and (iii) is investigated in Figure 5. In the low-stiffness regime no aux- etic behaviour is observed. The onset of ﬁbre strain stiffening in the transition regime leads to an increased rate of compaction in the 2-direction and the onset of auxetic be- haviour (3>1). The degree of auxetic deformation increases with increasing ﬁbre strain-stiffening. However, even at low ratios of ﬁbre stiffness to matrix stiffness (e.g. E2f=Em= 4) auxetic deformation occurs upon the onset of ﬁbre strain stiffening. In order to accurately model the anisotropic strain-stiffening behaviour observed in exper- imental testing of arteries, a high ﬁbre stiffness is required at high strains, so that aux- etic behaviour becomes unavoidable if a linear neo-Hookean (constant stiffness) matrix is used. Moreover, increasing the transition strain D1fleads to postponing the auxetic behaviour as it postpones the ﬁbre strain-stiffening (Figure 5 D, E, F). However, it is im- portant to note that the parameters explored in Figure 5 (i.e. the ﬁbre-matrix stiffness ratio (E2f=Em) and the values of transition strain (e.g. D1forD2f) are calibrated to ﬁt experimental stress-strain curves, and cannot be arbitrarily chosen to eliminate auxetic behaviour. 5. A strain–stiffening matrix with a strain–stiffening ﬁbre reduces auxetic behaviour As demonstrated in Figures 3 and 5, auxetic deformation is strongly inﬂuenced by the ratio of ﬁbre to matrix stiffness. Auxetic deformation is accentuated by strain-stiffening ﬁ- bres in a constant stiffness matrix. Therefore, in this section we propose a strain-stiffening isotropic matrix model with near bilinear behaviour (BLM model), placed in parallel with the BLF model described in Section 4. This approach allows precise control of the ﬁbre-matrix stiffness ratio at small and large deformations, facilitating the simulation of anisotropic strain- stiffening behaviour without unphysical auxetic behaviour. Assuming the isotropic strain energy function m(1;2;3)for matrix, where i;i= 1;2;3are the modiﬁed principal stretches, the Cauchy stress components, iare deter- mined from the following relationship (Holzapfel, 2000): Ji=i@ m @i1 33X j=1j@ m @j; i = 1;2;3 (11) Using the Valanis-Landel hypothesise the strain energy function m(1;2;3)is 12 presented as m(1;2;3) = iso(1) + iso(2) + iso(3) (12) In order to simulate an isotropic matrix material that exhibits bi-linear type strain-stiffening we propose a isochoric strain energy density function such that: Si=i@ iso(i) @i =8 >>>< >>>:E1m(i1) ji1jjD1mj pm(i1)2+qm(i1) +rm jD1mj<ji1j<jD2mj E2m(i1D2m) + (pmD2 2m+qmD2m+rm)ji1jjD2mj (13) where the parameters E1m,E2m,D1mandD2mhave the same meaning as the corresponding parameters in the ﬁbre model (equation (8)) and pm,qm, andrmare obtained in a similar manner as equation (9) by using the corresponding matrix parameters. After a straightforward manipulation the strain energy function associated with equation (13) is obtained as: iso(i) = 8 >>>< >>>:E1m(ilni1) ji1jjD1mj pm(2 i 22i+ lni) +qm(ilni) +rmlni+ 01mjD1mj<ji1j<jD2mj E2m(i(1 +D2m) lni) + (pmD2 2m+qmD2m+rm) lni+ 02mji1jjD2mj (14) where, 01mand 02mare two constants which ensure the continuity of strain energy. As an alternative formulation to the piece-wise strain-stiffening constitutive law above, we propose a non-piece-wise equivalent formulation in Appendix A (based on error function strain stiffening). It is shown that identical material behaviour can be obtained, both in terms of stress and strain energy density. This proposed BLM model is used in parallel with our BLF (equation (10)) to inves- tigate the control/elimination of auxetic deformation under uniaxial tension. Henceforth, this combination of matrix and ﬁbre models is called BLFM model for convenience. In Figure 6 we explore the inﬂuence of matrix strain-stiffening on auxetic behaviour. Ma- trix strain-stiffening is characterised by the ratio of the matrix tangent modulus at high strains (E2m) to that at low strains ( E1m). Increasing E2m=E1m>1results in reduced compaction in the 2-direction at higher applied stretch, and consequently a signiﬁcant reduction in auxetic behaviour. For the parameter sets considered in Figure 6, values of 13 3>1are computed only for the case of a non-strain stiffening matrix ( E2m=E1m= 1). We also explore the inﬂuence of the matrix transition strain parameter D1mon auxetic behaviour (with E2m=E1m= 3). If the value of D1mis increased, the matrix strain- stiffening is delayed, and may be preceded by signiﬁcant ﬁbre strain-stiffening. Auxetic behaviour occurs if D1mis signiﬁcantly higher than D1f. However, if the ﬁbre and the matrix undergo strain-stiffening at similar levels of deformation, auxetic behaviour may be eliminated. 1 2 30204060 1 2 30.40.60.81 1 2 30.850.90.9511.05 1 2 30204060 1 2 30.40.60.81 1 2 30.850.90.9511.05A C DB E F Figure 6: The inﬂuence of matrix stiffening and transition strain on the stress (A, D), in- plane transverse stretch (B, E), and auxetic behaviour (C, F). Reference parameters D1f= 0:2; D2f= 0:4; E1f=E1m= 1:5; = 45o; E2f=E1m= 8; D1m= 0:2; D2m=D1m+ 0:2;E2m=E1m= 3; K=E 1m= 30 have been used. To demonstrate the practical importance of our ﬁndings, we use our new ﬁbre/matrix strain-stiffening formulation to simulate experimental test data for aorta tissue (FitzGib- bon and McGarry, 2020). The results, as shown in Figure 7, reveal that the BLFM model is not only able to acccurately capture experimentally measured anisotropic stress-stretch data, but also correctly predicts the non-auxetic ( 3<1) out-of-plane stretch with rea- sonable accuracy. In contrast, if a non-strain stiffening linear matrix (L-BLF) model is used, signiﬁcant auxetic behaviour ( 3>1) is incorrectly computed. 14 Figure 7: The capability of the BLFM model to reproduce the stress-stretch experimental data of the uniaxial tension of ascending aorta FitzGibbon and McGarry (2020) (A) and out-of-plane stretch in circumferential test (B). FE results of the uniaxial tension tests for BLFM model and L-BLF models are also shown (C). ”L-BLF” stands for a linear matrix model in parallel with the BLF model. 6. Auxetic behaviour using an exponentially stiffening ﬁbre model The strain energy density function for the classical exponential form for collagen ﬁbre strain-stiffening (Holzapfel et al. (2000)), is given as: aniso(C;a04;a06) =k1 2k2X i=4;6 exp k2(Ii1)2 (15) wherek1andk2are material parameters. Note that equation (15) is the modiﬁed anisotropic (MA) form for slightly compressible materials ((Nolan et al., 2014)) requiring the use of the full anisotropic invariants I4;6instead of the isochoric anisotropic invariants I4;6. This exponential ﬁbre component is placed in parallel with a neo-Hookean matrix (equa- tion (6)). This exponential-ﬁbre formulation combined with neo-Hookean matrix has been widely used following implementation in leading commercial FE software programs such 15 as Abaqus (Abaqus, 2017), FEAP (Taylor, 2013), Ansys (Ansys, 2017) and ADINA (Ad- ina, 2017). It worth to note that although the contribution of the ﬁbres in compression is usually excluded, ANSYS implementation does not include the tension-compression switch of the ﬁbres. The response to a unixial tension test of this exponential-ﬁbre formulation combined with neo-Hookean matrix is demonstrated in Figure 8 using material parameters cali- brated for human coronary arteries (Holzapfel, 2002). Strong strain stiffening and anisotropy is evident in Figure 8A. However, signiﬁcant auxetic behaviour is evident (Figure 8B), with3>1when1>1:15. 1 1.2 1.4 1.6 1.800.050.10.150.2 1 1.2 1.411.11.21.31.41.51.61.7 Circumferential AxialA B Figure 8: Comparison of Ogden and neo-Hookean matrix in parallel with the exponential ﬁbre model. The combination of Ogden with exponential ﬁbre model postpones the development of auxetic behaviour in comparison to the neo-Hookean matrix with exponential ﬁbre model (B) while the stress-stretch curves are the same for both models (A). The parameters = 0:0014 MPa , k1= 0:0051 MPa andk2= 15:4have been used for NH-Exp model and two parameter sets (i.e., set 1:1=2= 0:0014 MPa , 1= 2= 1:5,k1= 0:0051 MPa ,k2= 15:4and set 2: 1=2= 0:0024 MPa , 1= 2= 4,k1= 0:0015 MPa ,k2= 15:4) have been used for Ogden-Exp model. Based on our proposed concept of using a non-linear matrix to control auxetic be- haviour (Section 5), here we next explore the combination of the exponential ﬁbre model (equation (15)) with the established Ogden strain-stiffening isotripic hyperelastic model 16 for the matrix. The Ogden hyperelastic strain energy function is given as: m(1;2;3) =NX i=12i 2 i i 1+ i 2+ i 33 (16) where, 1,2and3are isochoric principal stretches and iand iare material constants and we assumed N= 2 in this paper. Additionally, we retain the isotropic volumetric strain energy density function (ﬁrst term on RHS of equation (6)) and enforce near in- compressibility with K= 1:4 MPa . The material parameters of the ‘Ogden-Exponential’ model was chosen to reproduce the anisotropic strain-stiffening behaviour of the afore- mentioned ‘neo-Hookean-Exponential’ model (see Figure 8A). Two parameter sets have been found for ‘Ogden-Exponential’ model which ﬁt the corresponding stress-stain curve of the ‘neo-Hookean-Exponential’ curve, as shown in Figure 8A. In parameter set 1 a low rate of matrix strain-stiffening is considered. In parameter set 2 a higher rate of matrix strain-stiffening is considered. As shown in Figure 8B, using the Ogden model in place of a neo-Hookean model for the matrix does not eliminate auxetic behaviour, but merely postpones it (in this speciﬁc case 3>1when1>1:4.) The stiffness ratio between ex- ponentially stiffening ﬁbres and the power-law stiffening matrix will changes rapidly with increasing applied strain, resulting in uncontrolled auxetic behaviour. In contrast, our pro- posed BLFM model facilitates control of stiffness ratios and, consequently, elimination of auxetic behaviour. 7. Does including ﬁbre dispersion reduce auxetic behaviour? In this section the inﬂuence of ﬁbre dispersion on auxetic behaviour is investigated. To this end, we have implemented a discrete ﬁbre dispersion method around two mean ﬁbre directions for the exponential (MA) model and the BLF model via the user-deﬁned subroutine UMAT in Abaqus where the mechanical contribution of all ﬁbres under com- pression are excluded. We consider that ﬁbres can exist in mdiscrete directions in a 3D sphere at each integration point of a ﬁnite element model. A non-uniform segmentation of the unit sphere has been done such that higher density of segmentation around the mean ﬁbre direction is obtained. A sensitivity analysis of the number of ﬁbre directions mrevealed that convergence is achieved for m= 50 . For each ﬁbre family, the transversely isotropic and -periodic von-Mises probability distribution function (PDF) i, as proposed by Gasser et al., (2006), is used to describe 17 ﬁbre dispersion, whereby i=1 r b 2exp 2b(a0;i:a0;mean) er (p 2b)(17) where a0;iis a unit vector indicating one of the mdiscrete directions, a0;mean is a unit vector indicating the mean ﬁbre direction, bis a constant dispersion parameter, and er (x)denotes the imaginary error function. The introduction of ﬁbre dispersion results in a ﬁbre contribution in the 3-direction. Therefore, in equation (3), the 33-component of the ﬁbre stress tensor f 3will not be zero, and may therefore inﬂuence the computation of auxetic behaviour. The results for the BLF model and exponential (MA) ﬁbre model, both with a neo- Hookean matrix, are presented in Figure 9. 18 Figure 9: Effects of ﬁbre dispersion on the stress-stretch curve and auxetic behaviour, for bilinear (B,C) and exponential (D,E) ﬁbre models with neo-Hookean model for the matrix. The probability density function, , for the ﬁbres dispersion around the mean ﬁbre directions ( =45) is also shown in (A). Figure 9 suggests that the inclusion of ﬁbre dispersion does decrease the auxetic de- formation in both BLF and exponential MA models. As ﬁbre dispersion is increased the level of auxetic behaviour decreases. However, even for highly dispersed ﬁbres ( b= 10 ) auxetic behaviour occurs in both the BLF and exponential MA model. Only in the case of a nearly isotropic ﬁbre distribution (i.e. an isotropic material with nearly equal ﬁbre contribution in all directions) auxetic behaviour is eliminated. 19 8. Summary and conclusions In this study we uncover the underlying mechanism of auxetic behaviour in anisotropic hyperelastic soft tissue models that comprise of strain-stiffening ﬁbres in parallel with an isotropic matrix. High levels of in-plane matrix compaction due to increasing tension in strain-stiffening ﬁbres leads to unphysical out-of-plane expansion, i.e. auxetic be- haviour. Previous papers have reported the secondary inﬂuence of compressibility (Ska- cel and Bursa, 2016, 2019) and ﬁbre dispersion (V olokh, 2017) on the degree of auxetic behaviour, but no previous study has uncovered this primary underlying mechanism. Fur- thermore, we propose a new bilinear strain-stiffening ﬁbre and matrix (BLFM) model which facilitates close control of the ﬁbre-matrix stiffness ratio, thereby robustly elimi- nating auxetic behaviour. We demonstrate that our model provides accurate prediction of experimentally observed out-of-plane compaction, in addition to stress-stretch anisotropy, for arterial tissue subjected to uniaxial tension testing. The key ﬁndings of our study are summarised as follows: Auxetic deformation is strongly inﬂuenced by the ratio of ﬁbre to matrix stiffness. Auxetic deformation is accentuated by strain-stiffening ﬁbres in a constant stiffness matrix (e.g., the widely used neo-Hookean matrix with exponentially stiffening ﬁ- bres). Our proposed bilinear strain-stiffening ﬁbre-matrix constitutive model (BLFM model) facilitates robust control of the ﬁbre to matrix stiffness ratio in the low and high de- formation regime, consequently eliminating the computation of unphysical auxetic behaviour. Our BLFM model is shown to accurately capture experimentally reported anisotropic stress-stretch data for aortic tissue, in addition to correctly predicting the degree of (non-auxetic) compaction in the out-of-plane direction. In contrast, if a non-strain stiffening linear matrix model is used, signiﬁcant auxetic behaviour is incorrectly computed. Several studies demonstrate that removal of active cellular contractility, using cytochalasin- D, exposes non-linear strain-stiffening passive behaviour, e.g. (Reynolds et al., 2014; 20 Reynolds and McGarry, 2015). Therefore, the concept of a strain stiffening matrix intro- duced in our model may be motivated by the presence of the non-linear passive contri- bution of endothelial cells, smooth muscle cells, and ﬁbroblasts in arterial tissue, or the passive contribution of cardiac myocytes in cardiac tissue. Of course, the complete de- scription of the non-linear anisotropic behaviour of such tissues requires the addition of our anisotropic ﬁbre component. The new modelling framework will provide improved predictions of artery deforma- tion during physiological loading and stent deployment without incorrectly predicting auxetic behaviour under such large-strain deformation. Therefore, the new BLFM model can potentially improve current practices in the ﬁeld of design and regulatory approval of medical devices. Future extensions of the model should consider viscoelasticity (Jacopo and Paola, 2020), viscoplasticity (Weickenmeier and Jabareen, 2014), damage (Ferei- doonnezhad et al., 2016; Holzapfel and Fereidoonnezhad, 2017), and non-afﬁne matrix- ﬁbre interaction. The modelling framework should also be considered for anisotropic contractility of smooth muscle cells and ﬁbroblasts (McEvoy et al., 2019), in addition to passive non-linear behaviour of the cell cytoplasm (Reynolds and McGarry, 2015; Weafer et al., 2013). While we speciﬁcally calibrate and validate our novel modelling approach using experimental data for aortic tissue, we also explore the ability of our model to ro- bustly eliminate computation unphysical auxetic behaviour over a wide range of applied deformations (including very high deformations that may be encountered for skin (Lanir and Fung, 1974) or urinary bladder samples (Farhat et al., 2008). Future work should precisely calibrate the model using experimental data for several types of soft tissues, including myocardium (McEvoy et al., 2018), cartilage (Dowling et al., 2013) and skin (Annaidh et al., 2012b,a; Limbert, 2019). Acknowledgements This work was partially funded by a European Union Horizon 2020 Research and Innovation Program, under grant agreement No. 777072 and Science Foundation Ireland grant number SFI-12/IP/1723. 21