## Abstract

Exogenous cross-linking of soft collagenous tissues is a common method for biomaterial development and medical therapies. To enable improved applications through computational methods, physically realistic constitutive models are required. Yet, despite decades of research, development and clinical use, no such model exists. In this study, we develop the first rigorous full structural model (i.e. explicitly incorporating various features of the collagen fibre architecture) for exogenously cross-linked soft tissues. This was made possible, in-part, with the use of native to cross-linked matched experimental datasets and an extension to the collagenous structural constitutive model so that the uncross-linked collagen fibre responses could be mapped to the cross-linked configuration. This allowed us to separate the effects of cross-linking from kinematic changes induced in the cross-linking process, which in turn allowed the non-fibrous tissue matrix component and the interaction effects to be identified. It was determined that the matrix could be modelled as an isotropic material using a modified Yeoh model. The most novel findings of this study were that: (i) the effective collagen fibre modulus was unaffected by cross-linking and (ii) fibre-ensemble interactions played a large role in stress development, often dominating the total tissue response (depending on the stress component and loading path considered). An important utility of the present model is its ability to separate the effects of exogenous cross-linking on the fibres from changes due to the matrix. Applications of this approach include the utilization in the design of novel chemical treatments to produce specific mechanical responses and the study of fatigue damage in bioprosthetic heart valve biomaterials.

## 1. Introduction

The application of exogenous-cross-links (EXLs) to native or biologically derived soft collagenous tissues finds its way into a wide range of medical therapies and device applications, such as surgical biomaterials, modification of corneal tissues (using riboflavin/UVA) and vascular grafts. Perhaps the most mechanically demanding application is the so-called bioprosthetic heart valve (BHV), which is fabricated from several types of biologically derived soft collagenous tissue membranes. From a clinical perspective, BHVs have important advantages in that they do not require permanent anticoagulation therapy, operate noiselessly, and have blood flow characteristics similar to the native valve, and thus have become the dominant heart valve therapy worldwide [1–3]. However, BHV durability continues to remain limited to the range of 10–15 years, resulting from leaflet structural deterioration mediated by fatigue and/or tissue mineralization [4,5]. In general, structural damage is a critical factor in BHV degeneration, and clearly implicates coupled material and design factors as major limiters to long-term durability [6,7]. However, a major reason why advances in the use of cross-linked tissues in BHVs and other biomedical applications is a dearth of knowledge on how EXLs affect the underlying tissue structure and macroscopic mechanical behaviour. Moreover, integration of such information into truly predictive modelling cannot proceed without accurate constitutive models of the EXL tissues and the subsequent fatigue processes [8,9].

The most common medical applications of cross-linked biologically-derived tissues are for dense collagenous tissues. Such tissues are typically composed of a dense, highly-organized network of type I collagen fibres, along with elastin, proteoglycans, glycosaminoglycans, cellular materials and a small amount of other fibrillar proteins. Type I collagen is the major determinant of its mechanical behaviour [10,11] and is the major tissue component affected by EXLs. At the molecular level, tropocollagen molecules are composed of a triple helix of three alpha chains [10,11] that arrange themselves into a quarter stacking array to form the collagen fibril [10,11]. Collagen fibrils form the functional subunits of the collagen fibres, as described by the Hodge–Petruska model [10,12], and then exhibit distinct large-scale structures (figure 1). Like the fibrils from which their functional properties are derived, collagen fibres exhibit high tensile but low flexural stiffness [14]. We note that there is no standard definition for a fibre and its relation to the fibril. They are very dependent on the specific tissues involved, and thus caution should always be exercised in the terminology used. In the mainstream tissue biomechanics literature, the fibre/fibril definition is often taken from the well-known work of Kastelic *et al*. [13], which focused on tendon. The pericardial tissues considered in this study clearly show fibre and fibril structures, including the fibre undulations commonly observed optically [15] (figure 1).

The fibre longitudinal/axial direction is the primary determinant of the stiffness of the tissue composite. Interestingly, collagen fibres typically have a stiffness of 1 GPa [15–18] and extend by no more than 4–5%. To increase the tissue-level compliance, collagen fibres at the macroscopic scale are sinusoidally crimped [10]. Tissue-level stress will not occur until the fibre-level crimp has been straightened. Moreover, the distribution of fibre straightening strains is the mechanism of tissue nonlinearity at large strains [19,20]. Thus, as in many other fields, the connection to the underlying structure can greatly inform our understanding of how collagenous tissues work and guide the development of mathematical models of their mechanical function.

While there is a wide range of application-specific chemical agents (e.g. riboflavin/UVA cross-linking for corneas), for mechanically-demanding and blood-contacting applications (e.g. BHVs) using collagenous tissues, an aqueous solution of glutaraldehyde (GLUT) is used. GLUT application is necessary to both biochemically and mechanically stabilize the tissue for *in vivo* use. During the cross-linking process GLUT rapidly permeates the tissue, with the cross-linking process largely complete within an hour and essentially stabilized by 24 h. Much of what we know about how GLUT EXLs alter the structures of collagenous tissues was reported by Nimni, Cheung and co-workers [21–27]. Briefly, GLUT reacts primarily with e-amino groups of lysyl residues in proteins, with Michael-addition reaction products of Schiff bases usually the final stable products. Based on the spectral characteristics and the molecular weights of the reaction products, it has been predicted that GLUT reacts with free amines to form an intermediate with a molecular weight of about 200 Da. The GLUT–polymer amine complex is self-limiting in size and can undergo internal rearrangement to become chemically inert. An increased molecular length of GLUT polymers from the initial glutaraldehyde and lysyl-residue reaction is more likely than an increased number of cross-linked sites. Following free GLUT depletion by binding to reactive groups, additional GLUT molecules attached to already reacted molecules can give rise to larger GLUT polymers that are able to generate ‘long-range cross-links’ between further removed reactive sites (figure 2). It is apparent from a mechanical behaviour perspective that GLUT-associated chemical cross-linking of the collagen structure and biochemistry can produce complex changes from the native state at the molecular, fibril, fibre and tissue levels.

To enable improved application of the use of EXL tissues in *in situ* treatment and prosthesis design, the development of physically realistic constitutive models is clearly required. In a previous work, a structural approach was used that incorporated experimentally measured angular distribution of collagen fibres and an assumed isotropic form for the EXL matrix [28]. Good agreement with the experimental data was observed, supporting the basic approach. An important utility of that early model was its ability to separate the effects of the fibres and matrix. However, it was only a first step; other factors such as bending rigidity of EXL fibres, fibre–fibre interactions, and fibre–matrix interactions were not considered. Moreover, the experimental data were reduced under the assumption of an isotropic Fung model; no rigorous investigation of the most appropriate form was undertaken.

The focus of the present work is to more fully investigate the underlying characteristics of the effects of EXLs on soft collagenous tissues, and to use this information to develop a meso-scale (i.e. at the level of the fibre) structural constitutive model. In particular, we explored the following questions of the effects of EXLs on native collagenous tissues: (i) what are the effects on individual collagen fibres, (ii) what are the effects on single collagen fibre ensembles, (iii) are there interactions between fibre ensembles, and (iv) what is the functional form of the effective matrix response? This was done by exploiting experimental data from [29], wherein structurally controlled pericardial specimens were tested in the native state and then the EXL state. From these results, a comprehensive structural constitutive model was developed for EXL collagenous tissues and its predictive capability was evaluated. We note that, while we ultimately seek the micro-mechanical basis for macro-scale function, the present work is focused on a fibre-ensemble level approach.

## 2. Experimental methods and data post-processing

### 2.1. Tissue sources and experimental methods

Details of the tissue source, preparation and mechanical evaluation have been previously presented [29]. Briefly, large sections of native bovine pericardium were stored in phosphate-buffered saline (pH 7.4) at 4°C, then optically cleared using a hyperosmotic solution and the collagen fibre architecture (CFA) quantified. From the resulting CFA information, 25 × 25 mm test specimens exhibiting a high degree of structural uniformity suitable for biaxial testing were selected. The collagen fibre preferred and cross-preferred directions were aligned to the X_{1}–X_{2} axes (figure 3). A total of five specimens were prepared in the native state.

Biaxial mechanical testing methods have been previously described in detail [30,31]. Briefly, testing was performed with the specimen immersed in phosphate-buffered normal saline (pH 7.4) at room temperature. First the Piola–Kirchhoff stress **P** controlled test protocol was used, wherein the ratio of the normal stress components P_{11} : P_{22} was kept constant, with P_{12} = P_{21} = 0 and a maximum stress level of 1 MPa was used. Tissue deformations were quantified from the motion of four markers placed in the central third of the specimen, from which the deformation gradient tensor **F** was determined. For the first testing phase, an equi-biaxial stress protocol (i.e. P_{11} : P_{22} = 1 : 1) was used for both preconditioning and data acquisition. A total of 15 contiguous cycles were run with an approximate strain rate of 0.01 s^{−1}. Next, seven successive protocols were performed using ratios P_{11} : P_{22} = 1 : 0.1, 1 : 0.5, 1 : 0.75, 1 : 1, 0.75 : 1, 0.5 : 1 and 0.1 : 1. This range was chosen for extensive coverage of in-plane strain state. After testing, each native specimen was allowed to mechanically re-equilibrate by storing them in a stress-free state at 4°C for 24 h. Next, each specimen was chemically treated with 0.625% GLUT for a minimum of 72 h, with the tissue marker dimensions monitored throughout the cross-linking procedure, and then stored in phosphate-buffered normal saline at 4°C. As a final step, the above biaxial testing sequence was repeated. Data post-processing included computation of the second Piola-Kirchhof tensor **S** and deformation gradient tensor **F** using established methods [32]. This test design allowed a comprehensive planar mechanical behaviour dataset to be collected on matched native and EXL specimens, compensating for inter-specimen variations.

### 2.2. Kinematic considerations and mechanical data post-processing

As observed in our other studies [29,32], the chemical fixation process will affect the specimen dimensions, and any analysis must carefully account for these effects on the collagen fibre kinematics. We thus defined the following configurations: *β*_{0}-native, *β*_{1}-EXL (figure 3*b*), used as the referential configurations for the native and EXL states, respectively. We represented all deformations using the notation for the deformation gradient tensor where *i* and *j* represent the initial and final configurations, respectively (table 1). Values for the components of were determined using the same method from §2.1 for the displacements of the four markers pre- and post-cross-linking for each specimen. Next, as first described by Lanir [19], we defined a fibre ensemble as a group of fibres with a common orientation. It has been shown that the ensemble stress–strain relation can be obtained from the interpolated equi-biaxial strain path, where **F** = diag[*λ*, *λ*, 1/*λ*^{2}] using [8]. To derive the equi-biaxial strain path (*λ*_{1} = *λ*_{2}) from the stress-controlled experimental data, all mechanical data were combined and interpolated using cubic Hermite patches [33]. A strain path with *λ*_{1} = *λ*_{2} was interpolated within the range of the data as defined by the convex hull of {*λ*_{1}, *λ*_{2}}, and was implemented separately for each stress component (figure 4). To reliably overcome regions of sparse data we enforced the surface to be strictly convex everywhere. Finally, since the S_{12} component was negligible in all specimens, it was ignored in the subsequent analyses.

### 2.3. Establishing and modelling the mechanical behaviour of the native collagen fibre ensemble

Based on our previous tissue model findings [8,33–35], the dominant cause of the nonlinearity of the tissue-level mechanical behaviour of collagenous tissues is the gradual recruitment of collagen fibres [19]. The collagen fibres themselves behave linearly under typical fibre strains experienced under physiological stresses in tissues (2–5%). Once all fibres are fully straightened, the summed response should appear linear in the ensemble response. When applied to the equi-biaxial strain derived fibre-ensemble data, the upper bound can be directly determined as the transition point between the nonlinear and linear regions, with the slope of the linear region establishing the maximum tangent modulus (MTM) of the collagen fibre ensemble (e.g. [33]). The matrix response can also be determined from the pre-recruitment region where no collagen fibre have contributed. To determine the recruitment upper bound in the native tissue, we started at the largest measured strain and decreased the strain level until the region above was no longer linear. Linearity was defined from the mean squared error (MSE) of the linear regression to be less than 0.005% of the total MSE of all data, where Similarly, we determined the lower bound by starting at zero strain and increasing the strain until a deviation from linearity was determined.

Next, we note that in some previous structural models a fibre stress–strain relation linear in the second Piola-Kirchhoff stress and Green Lagrange strain has been used [8,34,36]. However, SAXS studies have demonstrated a linear force–displacement relation for collagen fibrils in the tendon [37,38] and MV tissue [39,40]. This is further corroborated by the atomistic modelling results by Buehler [41], where the force–displacement relation is essentially linear at strains lower than 0.35. We have recently determined that for the mitral valve leaflet the tissue level-derived collagen fibre mechanical behaviour is actually quite linear, with an effective modulus of approximately 160 MPa. Based on these considerations, we assumed for the native collagen fibres that

(1) they exhibit a linear P–

*λ*response(2) slack stretch of the collagen fibres does not vary with orientation.

From these two basic considerations, we used the following effective native collagen fibre model [33,34]. We start by defining the native collagen fibre strain energy as
2.1
where is the true stretch of the fibre. This leads to the following P–*λ* form using
2.2
where P_{f} is the first Piola-Kirchhoff stress of the fibre, *η*_{c} is the modulus of the fibre, *λ*_{f} is the fibre stretch and *λ*_{s} is the fibre slack stretch. Next, we use this fibre model in the expression for the native collagen fibre ensemble using
2.3
where *ϕ*_{c} is the collagen fibre mass fraction, *λ*_{θ} is the fibre-ensemble stretch along the direction defined by *θ* (computed from the tissue-level deformation using ), and *D*(*λ*_{s}) is the probability distribution function describing the distribution of collagen fibre slack length within the ensemble. We assumed *D*(*λ*_{s}) is Beta distributed, so that
2.4
where *λ*_{lb} and *λ*_{ub} are the lower and upper bound stretch of the collagen fibre recruitment, respectively. Note that in preliminary examinations of the data we found that all specimens exhibited distinct pre- and post-transition locations (figure 5), allowing *λ*_{lb}, *λ*_{ub} to be determined directly from the collagen fibre ensemble data. Thus, the complete initial ensemble model (equation (2.3)) has three parameters () to fit to the data using standard techniques [33].

## 3. Delineation and modelling of the tissue-level mechanical effects of exogenous cross-links

### 3.1. Rationale

While extensive work has been done on the characterization of the biomechanical effects of EXL formation on soft tissues [29,42–54], there has been surprisingly little work done known to the authors on the development of formal constitutive models (other than [28]). Related work on proteoglycan and related collagen fibril sub-forms have revealed complex micromechanical interactions (e.g. [55,56]), but micromechanical interactions modified by EXLs on the macroscale tissue responses remain largely unknown. Thus, prior to developing the constitutive model form, we first carefully examined the effects of EXL formation on the measured tissue-level biomechanical behaviours in the present dataset.

### 3.2. Effects of exogenous cross-link formation on collagen fibre ensembles

All specimens exhibited anisotropic dimensional changes due to preconditioning and cross-linking (figure 3*b*). Interestingly, we found that about 6% shrinkage occurred in the preferred direction, and approximately 7% expansion in the cross-preferred directions. Such changes can alter both the angular dependence on collagen recruitment and the collagen fibre orientation distribution. We determined basic characteristics of the collagen fibre stress–strain relations directly from the data (no modelling), including the lower and upper bound-associated stresses, the initial tangent modulus, and the MTM from a running 15-point window. From this analysis, we were able to determine a number of important mechanical characteristics (table 2 and figure 5). These include:

(1) All specimens exhibited an approximately 6.5-fold increase in the initial tangent modulus (table 2), very similar to values and native/EXL ratios reported in [54], which were conducted under flex conditions.

(2) The upper bound stress and MTM were found to be unaffected by EXL formation (table 2).

(3) EXL formation induced a reduction in achieved strain levels compared with the native state (figure 5

*a*).(4) The effective collagen modulus was unaffected by EXL formation (figure 5

*c*).

Collectively, these results reveal some important features of the effects of EXLs on collagen tissues. First, as noted in our previous studies [28,54] EXLs produce a substantial increase in the low-strain modulus. Next, equation (2.3) indicates that the MTM is proportional to the collagen fibre modulus and the recruitment function D. Use of equation (2.3) compensates for the effect of the changes in tissue dimensions due to cross-linking on the fibre recruitment, allowing separation of changes in fibre architecture from the modulus on the ensemble stress–strain curve. Thus, the lack of changes in effective modulus are independent of any effects of changes resulting from tissue dimensions and represent an accurate modulus estimate.

## 4. Initial model formulation

### 4.1. General approach

The above findings provide sufficient information to develop a new model. In the present study, we assume EXLs induce fibre–fibre and fibre–matrix interactions that are mechanically significant. We ignore any time-dependent effects, as we have found that native and cross-linked valvular tissues exhibit minimal time-dependent effects [57–60]. Next, we assume that the pericardial tissues considered are only composed of collagen fibres and a matrix constituent that represents non-cross-linked and cross-linked components, and water. The contributions from elastin or other tissue components are ignored, because they have either negligible mass or stiffness. In all previous structural models of soft tissues, interactions between components (fibres, matrix) have been ignored. As we cannot assume this in the present investigation, we use the following hyperelastic general form:
4.1
where *ϕ*_{c} is the mass fraction of the collagen fibres, *Ψ*_{c}, *Ψ*_{m} and *Ψ*_{int} are the strain energy density functions of the collagen, matrix and interaction terms, respectively, *J* = det(**F**), and *p* is the Lagrange multiplier to enforce incompressibility. The resulting tissue-level response in terms of the second Piola-Kirchhoff stress tensor **S** is given by
4.2

### 4.2. Accounting for changes in tissue dimensions for the collagen phase

Results from §3 suggest that the native collagen fibre modulus is unaffected by cross-linking. However, the observed changes in tissue dimensions can also induce changes in tissue-level mechanical behaviour by altering the structure due to tissue shrinkage. This essentially results in a different reference configuration. There is thus a need to reformulate structural models to account for these effects directly. The formulation described in the following allows handling of changes in tissue reference state geometry. The key assumptions are:

(1) Changes are due to alterations in the initial geometric configuration only, so that

(a) Mass fractions of each phase remain unchanged.

(b) The internal mechanical energy remains zero. Thus, all changes in internal component configurations are not associated with any change in internal energy (which remains zero)—just initial configuration (e.g. fibre orientation, degree of undulation, thickness and length).

(2) Tissue dimensions and internal architecture change under the

*affine*kinematic assumption. Thus, the configuration of all constituent fibres in the altered reference state (after all changes in initial specimen geometry have taken place) can be predicted. Moreover, the configurational change is homogeneous, and can be thus described by a deformation gradient tensor with constant components.(3) To be consistent with the fibre recruitment mechanisms (e.g. [8,33]), all fibres remain undulated in the new reference state.

(4) The matrix phase is unaffected by the geometric configuration changes and is referenced to

*β*_{1}(figure 3*b*) for all subsequent stress calculations.

As a first step, we recast the recruitment function parameters determined in *β*_{0} but mapped to *β*_{1} using
4.3
where the left subscript denotes the configuration state. Note carefully that will in general be a function of *θ*_{1}, so that even though the collagen fibre recruitment is orientation-independent in *β*_{0}, it will have angular dependence in *β*_{1} unless The resulting expression for the ensemble stress is
4.4

To complete the tissue-level formulation, we use the affine transformation assumption and the formulation described in [34] to obtain the collagen fibre orientation distribution function *Γ*_{0} in the native state to that in *β*_{1} using
4.5
Note that the angle *θ*_{1} of a fibre originally oriented at *θ*_{0} can be determined using
4.6

The final form of the native collagen fibre phase expressed in the EXL state is thus 4.7

The above equation has a total of five fitted model parameters *η*_{c}, *μ*_{m}, *μ _{Γ}*,

*σ*,

_{Γ}*σ*0 (

*λ*

_{lb},

_{0}

*λ*

_{ub}are determined from the experimental data), all with a physical meaning and all referred to

*β*

_{0}.

### 4.3. The matrix phase

The matrix response can be estimated from the low stress region where collagen does not contribute any stress (figures 5 and 6). A careful examination revealed that the toe region is a convex curve in S–*λ*, which is inconsistent with a neo-Hookean material model which is concave in S–*λ* (figure 6). While we have used an exponential isotropic function for the matrix before [28], we considered the Yeoh model but found it unable to fit the data. Therefore, we developed a modified Yeoh model as
4.8

When applied to all pre-transition mechanical data (i.e. from all protocols and where there is no collagen fibre contribution), this form was found to fit the low-stress data quite well (figure 6).

### 4.4. Interactions

Our key working assumption is that all fibre–fibre and fibre–matrix interactions can be represented at the *fibre-ensemble level*. This is done to simplify the model formulation, and further since the exact micromechanical mechanisms of cross-linking have yet to be determined. Based on the form assumed in equation (4.2), we now have the ability to estimate the form and magnitude of the fibre-ensemble interactions from the fibre-ensemble data using where **S**_{ens} is the ensemble stress (figure 5). Note here that the collagen stress and is thus the contribution of the collagen fibres expressed in the EXL configuration *β*_{1} using equation (4.7). This approach allowed us to exploit the matched native–EXL mechanical data by fitting the native responses then mapping them to the EXL state, so that they are a known quantity rather than one that required data fitting. Results of this analysis indicated some intriguing results. First, while the collagen phase contributed substantially to the total ensemble stress, it was only about 50%, and the matrix only about 20%. This revealed that the remaining approximately 30% portion of the total ensemble stress must be a result of the interaction mechanisms.

To model the interactions, we first consider two fibre ensembles with orientation vectors **n**_{0}(*α*) and **m**_{0}(*β*) in the reference configuration (figure 7). These two ensembles can mechanically interact by elongation and relative rotation. Kinematically, these mechanisms can be captured using the pseudo-invariant I_{8} [61,62]
4.9
Note that we can also use [62] to account for the relative change in fibre rotations if necessary. Thus, I_{8} can be considered the product of an extensional term and a rotational term cos(*θ*). We consider two sub-aspects of ensemble-level effects: intra- and inter-ensemble levels. The intra-ensemble incorporates all fibre–fibre interactions that occur within a single ensemble and are limited to extensional effects only. By contrast, inter-ensemble effects can include both extensional and rotational effects.

To best capture these phenomena, we do not use I_{8} directly but rather the following forms. The results for the ensemble stress suggest that the interaction terms are exponential in character (figure 8). Thus, for the extensional intra-ensemble effects we use
4.10
where and *c*_{0}, *c*_{1} are constants, and the associated single ensemble stress extensional interaction is
4.11
Next, for the *extensional inter-ensemble* interactions, we use a similar form
4.12
where *d*_{0} and *d*_{1} are constants, with associated stresses
4.13

Note that the exponential term was set to zero if In this approach, we integrate over all ensembles, weighted by their respective orientation distribution functions, to obtain the total contribution.

Next, we developed a rotational pseudo-invariant defined as the change in the cosine between the two ensemble fibre directions, which is simply 4.14 Using it can be shown that for planar distributed fibre ensembles the stress tensor is 4.15 Combining equations (4.7), (4.8), (4.11), (4.13) and (4.15), we obtain the final form of the full model stress as 4.16

## 5. Model simplifications and parameter estimation

### 5.1. Final form

While in principle equation (4.16) can be implemented within a robust parameter estimation procedure, simplifications are clearly in order given its complexity (14 parameters). We first consider an equi-biaxial test wherein all fibre rotations are zero, so that The total interaction stress is thus given by just the following extensional contributions: 5.1

In practice, we found that while the intra-ensemble form (first r.h.s. term) was used alone it was able to capture the equi-biaxial strain behaviour (figure 5), it was unable to capture the response from all test protocols. Moreover, given the similarity in form, the two components on the r.h.s. of equation (5.1) can capture similar responses. We thus chose to ignore the intra-ensemble stress contribution removing two parameters. Next, while it is intuitive that fibre-ensemble rotations produce important contributions to the total tissue stress, closer analysis of equation (4.15) indicated that it will produce compressive stresses in the direction of lesser stretch. These characteristics were not consistent with any of the observed experimental data. Even when choosing various forms of we could not match the experimentally observed responses. Interestingly, only the contribution of equation (4.16) was found to model the interaction stresses well.

Thus, we are left with the following interaction stresses:
5.2
leading to the following *final form* of the constitutive model:

5.3

It is understood that **n**_{0} and **m**_{0} are referred to *β*_{1} and that we merged the Lagrange multiplier with the matrix by assuming a planar tissue to simplify the formulation. This final model has 11 independent fitted parameters and three directly determined parameters all with a physical meaning.

### 5.2. Parameter estimation procedures

While at first glance this appears to be a major nonlinear optimization undertaking with all the usual pitfalls, we can use the following sequence to make actual parameter estimation quite tractable:

(1) From the native tissue mechanical data, predict the collagen phase parameters using standard procedures [33,63].

(2) From the pre-transition collagen recruitment portion of all of the EXL tissue mechanical data, determine the matrix parameters .

(3) Using the and responses, determine the interaction stress responses for

*all test protocols*using(4) Using the results of step 3, determine the final two parameters by fitting equation (5.3) but only allowing them to vary while keeping the other terms to their above-fitted values.

We found that this basic sequence ensured a robust parameter is obtained, because the entire model is never fitted at once. Moreover, this approach allowed us to separate the contributions to the stress of each of these mechanisms. As in our previous studies [33,63], we employed the genetic based Differential Evolution algorithm in Mathematica to perform the optimization. All parameter estimation was performed using a custom program written in Mathematica (Wolfram Research Corp.).

## 6. Primary results

From the five specimens used, the model was able to successfully fit all data quite well (total fit *r*^{2} > 0.97). Moreover, the final parameter values were quite consistent, with generally low standard errors (table 3). The mean collagen fibre modulus (approx. 279 MPa) and fibre splay (approx. 38°) were comparable to previous studies [34]. Interestingly, the lower bound stretch was small (1.01 or approx. 1% strain), so it is likely that it could be set to 1 (i.e. zero strain). The native collagen fibre recruitment parameters were also consistent (table 3), and indicated a very rapid recruitment at stretch of approximately 1.18–1.2 (figure 9). This is a more complete picture of the entire fibre recruitment than in our previous work [8,34], and suggests that the collagen fibres are effectively well ordered with a small deviation in crimp amplitude and wavelength.

One advantage of our approach is that the various contributions to the total stress can be separated (figure 8). To better reveal the present findings, it is useful to examine the effects on the individual stress components under various loading paths. Following the trends of the ensemble results (figure 9), we noted that the interactions produced substantial contributions to the total stress (figure 10). Interestingly, for S_{11} the interactions actually produced the largest contributions, followed by the matrix and collagen fibres. By contrast, for S_{22} the contributions were much more dependent on the particular loading path, with the collagen phase dominating when *λ*_{2} > *λ*_{1}. When *λ*_{1} > *λ*_{2}, the matrix phase dominated S_{22}. We further note here that the contribution of the matrix was much less loading path sensitive, due to its near-linear, isotropic behaviour.

## 7. Discussion

### 7.1. Major findings

This study represents the first rigorous *full structural model* (i.e. explicitly incorporating various features of the CFA) for cross-linked soft tissues, and also includes a specific interaction term. An important utility of this model is its ability to separate the effects of EXLs on the fibres and matrix, so that the matrix, collagen and interaction effects could be clearly identified. This was made possible, in-part, with the use of the native–EXL matched experimental dataset and a modification to the structural model so that the uncross-linked collagen fibre responses could be mapped to the EXL configuration. As in our previous study, we found that the matrix could be well modelled as an isotropic material. However, we found that a much more linear-like response was necessary. Perhaps the most novel findings of this study were that (i) the effective collagen fibre modulus was unaffected by cross-linking and (ii) the interaction term played such a large role in stress development, often dominating the response (depending on the component and loading path being considered).

The lack of change in the effective collagen fibre modulus has been corroborated by experimental results from Gentleman *et al*. [16]. In that study, they found a modulus range of 269.7 ± 11.9 to 484.7 ± 76.3 for cross-linked collagen fibres in the bovine Achilles tendon, which corresponds to the same range as another study for native collagen fibres from various sources [17]. Yang *et al*. [18,64] determined that for the mechanical properties of *hydrated* native and cross-linked type I collagen fibrils that cross-linking the collagen fibrils with a water-soluble carbodiimide did not significantly affect the bending modulus. The work by Shen *et al*. [15] noted a modulus of 0.86 ± 0.45 GPa (range, 0.36–1.60 GPa; *n* = 13), in reasonable agreement with our results. Further, six of the 13 fibrils showed linear behaviour. At the tissue level, our findings are also consistent with the findings of Lee *et al*. [65–67], who found no change in MTM in cross-linked pericardial tissues as in our study (table 2). It is interesting to note that, when using the native tissue fibre-ensemble model (equation (4.4)) on both the native and cross-linked data (figure 5), one can increase the fibre modulus determined from the native state to match the post-EXL data (figure 11*a*). However, this will induce a parallel increase in the MTM of approximately 75% (figure 11*b*), which is inconsistent with the experimental findings (table 2). This is the case even when compensating for the effects of tissue contraction. This simple simulation lends support to the collagen modulus results (figure 5*c*).

Our model results suggest that the major effect of EXL formation was not the formation of a mechanically substantial matrix or stiffening of the collagen fibres, but rather a *dramatically enhanced bonding both within and between fibre ensembles.* Note that our specific interaction term (equation (4.13)) captured the effects of both individual ensemble stretch and relative stretching between ensembles. This is consistent with what is known about GLUT bond formations (figure 1*a*) [22,23,25–27]. Yet, we found that relative rotations between fibre ensembles as modelled by equation (4.15) could not capture the observed responses. The underlying structural mechanisms for this behaviour remain largely uncharacterized. One possibility is that the protein core of the proteoglycans that bind collagen fibrils become strongly cross-linked and are thus substantially stiffened, acting to more cohesively bind the collagen fibres. This is supported by findings of Liao & Vesely [68], who observed substantial deformations of proteoglycans in mitral valve chordae under stretch. Moreover, the present model suggests that such mechanisms are strongly associated with fibre-ensemble orientation distributions. While not the final word, our results suggest that EXLs produce a stiffening effect via an isotropic matrix, with the interaction effects being the dominant effect.

### 7.2. Modelling approach

The present approach was a direct extension of the stochastic, tissue-level meso-structural models first pioneered by Lanir [19] and used in various applications and extensions by our group [8,33,63]. By meso-scale, we refer to the fibre-ensemble scale, which is fibre features down to approximately 100 µm. Moreover, recent evidence has demonstrated that the underlying affine kinematic assumption is valid [34,35], so that the strain energy of each fibre ensemble can be kinematically connected with the macroscopic strain tensor. We took the approach to develop a more general model first, then show what components could be modified or removed entirely. Given the lack of knowledge and modelling efforts in this area, we felt this was appropriate and helped to illustrate what underlying mechanisms should be incorporated. We also considered a more extensive approach based on elastica-based theory for sinusoidal fibres based on [69] under the assumption that EXL formation dramatically increased the fibre modulus, which ultimately proved to be unnecessary. An inter-fibre sliding model was also developed based on the relative sliding between fibres due to differences between their respective slack lengths, and used as a means to model the intra-ensemble EXL effects. However, this approach was found to be unable to capture the individual fibre-ensemble responses, suggesting that relative sliding between fibres at the ensemble level was not a major mechanism. A final question that may impact physical plausibility of the current model is convexity and physical plausibility. Lanir [70] demonstrated that the native tissue structural model is compatible with a physically plausible response. Based on both this and the experimental observations, we focused on monotonically increasing functions of strain for all model terms to ensure physical plausibility and that convexity was maintained.

### 7.3. Limitations and future directions

The current model is limited in the fact that it is quasi-static and does not account for permanent set phenomena we have observed [71]. The homogenization method used in §3 is fairly standard and has been used by others and the author for some time (e.g. [8,33,34]). We have observed that for soft tissues, whose composition is dominated by a dense network of distinct fibres, especially collagen type I and elastin, these structures can be mechanically treated as fibre ensembles; that is, groups of fibres with a common orientation. The scale of the representative volume element is about 100 µm, which is sufficient to capture the salient mechanical features of the fibre ensemble. We emphasize that this is not meant to be a universal model of all types of soft tissue structures, which are very diverse for a single theoretical treatment. Rather, we focus on a sufficiently generalized approach for exogenously cross-linked dense collagenous tissues, such as pericardium, heart valves and sclera. These structures fit our approach well and also have important biomedical therapeutic applications.

A complete understanding of the current phenomena must be based on well-characterized micro-scale events. For example, Kojic *et al*. [72] developed a model for fibre–fibre kinetics that uses Coulomb friction, which results in a simple and robust approach for tissue simulations. However, our knowledge of even native tissues at the micro-level remains limited. The situation remains more complicated by the fact that our knowledge of the interrelationships between the physical chemistry of EXLs and the macro-scale mechanics of collagenous tissues remains limited at the present time for more sophisticated models to be reliably attempted. To date, no material model is able to fully account for such complex observed microstructural and biological behaviour. The next steps include incorporation of the permanent set effect commonly observed in cross-linked tissues into the present model, and exploration of how alternative cross-linking methods affect macro-scale tissue behaviours.

## Authors' contributions

M.S.S. developed the research approach, the analytic approaches and wrote the paper; W.Z. coded the parameter estimation program, reanalysed the experimental data, fitted all data and tabulated results; S.W. assisted with the extension of the native tissue model to handle changes in reference configuration.

## Competing interests

We declare we have no competing interests.

## Funding

This study was supported by NIH grant no. R01 HL108330.

## Footnotes

One contribution of 19 to a theme issue ‘Integrated multiscale biomaterials experiment and modelling: towards function and pathology’.

- © 2015 The Author(s)