## Abstract

Tendon exhibits anisotropic, inhomogeneous and viscoelastic mechanical properties that are determined by its complicated hierarchical structure and varying amounts/organization of different tissue constituents. Although extensive research has been conducted to use modelling approaches to interpret tendon structure–function relationships in combination with experimental data, many issues remain unclear (i.e. the role of minor components such as decorin, aggrecan and elastin), and the integration of mechanical analysis across different length scales has not been well applied to explore stress or strain transfer from macro- to microscale. This review outlines mathematical and computational models that have been used to understand tendon mechanics at different scales of the hierarchical organization. Model representations at the molecular, fibril and tissue levels are discussed, including formulations that follow phenomenological and microstructural approaches (which include evaluations of crimp, helical structure and the interaction between collagen fibrils and proteoglycans). Multiscale modelling approaches incorporating tendon features are suggested to be an advantageous methodology to understand further the physiological mechanical response of tendon and corresponding adaptation of properties owing to unique *in vivo* loading environments.

## 1. Introduction

### 1.1. Overview of tendon structure and composition

Tendons are stiff fibrous connective tissues that facilitate the transmission of muscular forces to bone [1,2]. Tendons have a complex hierarchical structure that spans several length scales: *tropocollagen* (approx. 300 nm long), formed helically by three intracellular peptide alpha chains, is assembled into fibrils; *fibrils* (approx. 1 µm long) are postulated to be cross-linked by glycosaminoglycan (GAG) side chains to form fibres; *fibres* (ranging from 1 to 300 µm) are bundled together to build fascicles, with spindle-shaped cells (i.e. tenocytes) sparsely situated in rows between fibres; and *fascicles* (on centimetre scale) are enclosed by a connective layer called the epitenon (or interfascicular matrix) and form whole tendon tissue (figure 1) [3–11]. The unique structural organization of tendon provides characteristics that enable proper functional properties. For example, waviness (or crimp) of fibres along the long axis of tendon is proposed to contribute to the distinct toe region of the nonlinear stress–strain curve observed in tensile testing [1,12]. The crimp pattern is also superposed by a three-dimensional helical superstructure of collagen fibres, detected by polarized light microscopy and interference microscopy with three-dimensional modelling [13,14]. This helical pattern is demonstrated to contribute to greater recovery and less hysteresis loss of energy-storing tendons compared with positional tendons [15]. Besides unique substructural features, tendon mechanical behaviour is also dependent on its composition such as collagen fibres embedded in non-fibrillar matrix [10].

Tendon is composed of approximately 55–70% water and predominantly consists of 60–85% collagen of dry weight [16,17]. Type I collagen comprises approximately 95% of total collagen, but tendon also has small amounts of types III, V, XII and XIV [18]. Tendons also contain proteoglycans (PGs, i.e. decorin, fibromodulin, biglycan, lumican, aggrecan) in relatively small quantities, with types and amounts varying according to different species, ages and anatomical locations [11,19–21]. In tendon, PGs have been found to play an important role in facilitating fibrillogenesis, regulating extracellular matrix assembly and modulating cell–matrix interaction [20,22–24]. Properties of, and interactions between, structural components at different length scales of the hierarchical tendon organization are significant for providing adequate mechanical function in a variety of *in vivo* loading environments. Elucidation and characterization of these multiscale properties and relationships are necessary to understand tendon in healthy conditions as well as in degeneration, damage and injury.

### 1.2. Motivation for modelling tendon mechanics

A number of recent studies have characterized the anisotropic, inhomogeneous and viscoelastic properties of tendons in different species and in different anatomical locations. For example, rat patellar tendon exhibits less ultimate stress than rabbit patellar tendon, and human supraspinatus tendon has smaller elastic moduli than human patellar tendon [25–29]. Several other factors, such as age, genetic modification, injury and especially *in vivo* mechanical environment, can also influence tendon mechanics [29]. Demonstrating the multiscale nature of tendon mechanics, applied strains were found to be attenuated from tissue to local matrix scale in both equine superficial digital flexor and common digital extensor tendons [15,30]. Illustrative of the multiaxial mechanical properties of tendon, our work has shown that distinct regions of bovine deep digital flexor tendon exhibit different stresses under shear and compression loading, and similar to tensile results, strains were attenuated from the tissue to nuclei scales under non-tensile deformation [6,31]. The modifications of composition and structural properties that occur across different length scales as a result of these factors certainly contribute to overall tendon mechanics and function [32–34]. However, it is time-consuming and costly to experimentally evaluate all tendons under all types of conditions. An integrated-throughput modelling system could provide the capacity to characterize such situation-dependent cases and perform multilevel and multiparameter evaluation. Towards this end, experimental studies that evaluate distinct tendons using consistent testing methods and protocols [26,31,35–39] are beginning to provide data necessary to develop appropriate model formulations. Furthermore, the utilization of multiscale modelling could help explore physical mechanisms regulating tendon behaviour, including controversial issues such as whether PGs mediate load share between fibrils [22,40,41].

The full evaluation of tendon mechanical properties and function *in vivo* is difficult owing to limitations in current technology. Although recent studies have demonstrated that ultrasonography and magnetic resonance imaging can be used to examine *in vivo* tissue abnormalities during the treatment of tendon injuries [42–44], more information, not just tissue stiffness and strain, needs to be extracted from experimental studies to improve understanding of mechanisms governing tendon responses in complex loading environments [45]. Modelling and simulation can function as complementary modules to mimic the physiological interactions of tendon with other structures, evaluate its mechanical properties under normal function, and predict disease-related perturbations or efficiency of injury treatments [46]. Importantly, most previous experimental studies of tendon structure–function relationships examined properties only at a single length scale (i.e. molecular, fibril/fibre or tissue). In contrast, a computational model using a multiscale approach could determine how forces/deformations applied at the tissue level are transduced to the cell level, and quantify the alternation and modulation of local extracellular environment for tenocytes as a function of applied physiologic loads.

The objective of this review is to summarize mathematical and computational modelling approaches that have been used to evaluate mechanical properties of tendon at different length scales. We summarize challenges for different spatial scale modelling and propose a multilevel modelling approach that could be used to better couple physics and mechanics between different structural scales of tendon.

## 2. Mathematical modelling

A number of mathematical models of tendon mechanics at the tissue-, fibre- or fibril-level can be categorized into two general types: (i) the *phenomenological approach*, which uses a mathematical formulation that can closely match the mechanical response of tendon, but often with parameters that lack a clear physical interpretation, and (ii) the *microstructural approach*, which describes gross tendon mechanical behaviour by combining or generalizing the mechanical behaviour of different components in tendon (e.g. collagen fibre distribution). Many models, of both types, are developed based on a continuum theory, wherein the mathematical formulation used to represent tendon behaviour is based on a continuum assumption on the tissue scale. In addition, several different modelling schemes have been developed to evaluate mechanical properties of individual tendon constituents, or their interactions, on the molecular level. Details and examples of these modelling formulations applied to different loading types (figure 2) are outlined in this section.

### 2.1. Mathematical modelling on the tissue scale using phenomenological approaches

On the tissue scale, phenomenological approaches have been widely used to evaluate tendon mechanics. A viscoelastic constitutive law, satisfying thermodynamic principles, works for large deformation of tendon [47], and one of the most commonly used is the quasi-linear viscoelastic (QLV) theory developed by Fung [48]. Based on time- and strain-dependent mechanics of soft tissues, the QLV model separates tissue mechanical properties into a time-independent stress component and a time-varying relaxation function [49–52]. Owing to its reduced mathematical interpretation relative to more complicated nonlinear formulations, this theory has been widely applied to cartilage, ligament and tendon [50,53–55]. However, because the QLV theory is developed on the assumption of an instantaneous step increase in strain and hold, the application of the OLV theory to more complicated *in vivo* loading patterns is limited and fits to experimental data can be limited by issues with convergence and non-unique solutions [56]. Another example of a viscoelastic model that has been used is the single integral finite strain viscoelastic model, which can yield classical linear viscoelasticity and reduce to be an appropriate finite elasticity for time zero [57]. Although this model has the potential to be extended and generalized to describe a variety of nonlinear anisotropic materials in three dimensions, its limitation is to assume an instantaneous material response after strain is applied [57]. More accurate mathematical models should be developed to fully interpret the complex nonlinear experimental results of the strain-dependent stress relaxation and stress-dependent creep behaviour of ligaments [58,59] and tendons [60].

To better capture complex nonlinear properties, significant advancements on viscoelastic models have been made (table 1). Because water, as a primary component of tendon, contributes to time-dependent properties by interacting with fibrillar and amorphous elements such as PGs, fluid-flow-dependent viscoelasticity of tendons has been considered in recent models [66–69]. Furthermore, tendon mechanical behaviour should be characterized in three dimensions, but simplification to two-dimensional models is beneficial to both successfully characterize tendon behaviour under uniaxial loading and eliminate material constants during the fitting process. Based on these assumptions, Duenwald *et al*. evaluated the accuracy of a modified QLV model with nonlinear superposition and a Schapery nonlinear viscoelastic model to interpret strain-dependent relaxation behaviour of biological tissues by including rise time and ramp strain history during loading [61,70]. Nonlinear superposition of this model failed to address recovery behaviour, but Schapery's model matched experimental curves and was able to predict both recovery and reloading behaviour [61–63,71]. Troyer & Shetye modified the single integral nonlinear superposition (SINS) formulation with the nonlinear relaxation modulus using a Prony series representation to develop a novel formulation, which can include nonlinear viscoelastic behaviour and model physiologically static and dynamic loading regimes [55]. A novel constitutive model formulated within the nonlinear integral representation of transversely isotropic tendon was also found to capture the nonlinear relaxation response of rat tail tendon fascicles [72,73]. In addition, Kahn *et al*. [64] proposed a model based on irreversible theory to characterize tendon behaviour under various conditions such as loading, unloading and reloading or successive relaxations at different strain levels. Finally, compared with QLV or exponential stress–strain models, DeFrate & Li found that the Mooney–Rivlin model more accurately predicted elastic stress at strain levels higher than the levels applied during testing for fitting the model [65]. These models have the ability to capture tendon response during quasi-static and dynamic tensile loading processes (i.e. loading, unloading, cyclical loading), but they are not evaluated under other loading scenarios (e.g. shear and compression). Additionally, comparisons among the many possible continuum models are lacking, making it difficult to determine which model to adopt in certain conditions.

Regarding the phenomenological approach, there have also been several models developed using various combinations of dashpots and springs (table 2) to address tendon mechanics observed in different *in vitro* loading scenarios [31,55,74–76]. For example, a component mechanical model composed of an elastic (steady-state) spring in parallel with four Maxwell components was proposed to describe tendon strain-dependent viscoelastic behaviour [55]. In our previous work, we used a two-relaxation-time solid linear model formed by two Maxwell elements and one spring in parallel to describe step strain relaxation of bovine flexor tendon [31]. In another study, a simple structural model incorporated two Kelvin models in series to represent non-fibrillar matrix and collagen fibrils, respectively, and the strain-rate dependent mechanical properties demonstrated by this model were consistent with experimental data of both normal and cross-link-deficient tendons [74]. Finally, a generalized Maxwell model, also called a Maxwell–Wiechert model, was shown to satisfactorily describe the viscoelastic behaviour of human subscapularis tendon under uniaxial tensile stress relaxation owing to its ability to include many spring–dashpot elements [75]. These viscoelastic models are relatively straightforward and can capture many aspects of tissue-level behaviour, yet this simplified approach does not consider the interface between collagen fibres and interfibrillar matrix and is unable to explore mechanisms of deformation at smaller length scales. Models taking structural organization and compositional constituents into consideration are necessary to address the influence of tendon's unique structural features (i.e. crimping, helical pattern) and specific components on normal function.

### 2.2. Mathematical modelling on the fibril/fibre scale using microstructural approaches

#### 2.2.1. Contribution of collagen fibril to tendon mechanics

It is well known that tendon under quasi-static tensile loading exhibits three stages of the stress–strain curve: an initial low-strain toe region, a higher-strain linear region, followed by failure (figure 2*c*) [40,77]. Because waviness (or crimp) of collagen fibrils and a helical superstructure of tendon are thought to contribute to the toe region [78,79], a few models have sought to address the influence of these structural features on tendon mechanics [12,73,80–84]. Hurschler *et al*. proposed that the straightening stretch ratio *λ*_{S} (i.e. the tissue stretch at which the fibre straightens and begins to bear load) of fibres in tendon was distributed based on a Weibull three-dimensional fibril orientation probability distribution function (PDF),
2.1
where *β* (*β* > 0), *δ* (*δ* > 0) and *γ* (*γ* > 0) are shape, scale and location parameters, respectively [73,80–83,85]. In this approach, the longitudinal tissue-level stress can be computed by integrating the contribution of all fibres over *λ*_{S},
2.2
where *σ*_{33}(*λ _{t}*/

*λ*) is the normal stress of the fibre and

_{s}*λ*is an integration variable [80]. After fitting parameters to experimental data via least-squares minimization, this model could predict the response of tendon in the nonlinear toe region. In a similar way, Szczesny

_{t}*et al*. applied a gamma-probability distribution to describe the distribution of fibre crimp [84,86,87], 2.3 where

*λ*is the uncrimping stretch, and

_{c}*α*and

*β*are positive constants that determine the shape of the distribution. After incorporating this PDF, the modelling results better predicted the initial nonlinearity of macroscale tissue behaviour. Likewise, Sverdlik & Lanir [81] proposed a beta PDF to illustrate the variability of the initial crimp pattern

*e*

_{s0}, 2.4 where the subscript ‘0’ refers to the initial state,

*α*and

*β*are parameters which determine the shape of the beta function,

*e*

_{s0 max}is the highest level of

*e*

_{s0}and

*Γ*is the gamma function. Additionally, the three-dimensional axisymmetric fibril orientation PDF can be formed by generalizing the two-dimensional von Mises PDF into a sphere [80], 2.5 where

*κ*is the concentration parameter, and

*Φ*and

*θ*are location parameters in the sphere. According to the studies described here, these fibril orientation PDFs could demonstrate the contribution of crimp pattern to mechanical behaviour of sheep digital extensor, rat tail and bovine Achilles tendon under

*in vitro*quasi-static relaxation or failure loading by adjusting the shape and location parameters. It remains to be seen if these functions can be applied to all different tendons under a range of unique loading protocols.

Collagen fibres also show time- and history-dependent stress–strain behaviour, similar to tendon on the tissue scale. With the assumption of fibres being independent of each other, collagen fibres have recently been characterized by linear elastic, bilinear elastic, linear viscoelastic and QLV formulations under tension, in studies which both consider the interfibrillar matrix or neglect it [57,80,81,88–90]. The various fibril PDFs (described above) can be incorporated into these formulations to examine the role of collagen fibres in tendon behaviour.

#### 2.2.2. Contribution of proteoglycans to tendon mechanics and their interaction with collagen fibrils

On the fibril/fibre scale, the microstructural approach is frequently applied to evaluate tendon mechanics by combining mechanical behaviours of collagen fibrils/fibres with PGs. In this formulation, tendon is modelled by embedding straight or wavy collagen fibrils/fibres in a surrounding non-fibrillar matrix, mainly composed of PGs [91–93]. The strain energy function for the composite tendon can then be simply decomposed into the elastic strain energy stored within reinforcing fibres and the isochoric energy of the deformation in the matrix [80,94]. The collagen strain energy can be further divided into three parts accounting for a slack length without bearing stress, the exponential and the linear forms [94,95]. In previous models, PGs and collagen fibrils have been proposed to be bridged by reversible noncovalent cross-links (figure 1*c*), where anionic GAGs of PGs extend between collagen fibrils at an interval of approximately 60 nm [96]. For example, thin 6–8 nm filaments arranged orthogonally to collagen fibrils distribute in a regular distance corresponding to the collagen fibril D-period, and connect collagen fibrils laterally to one another [97]. Hydrogen bonds, hydrophobic connections within carbohydrate backbone and van der Waals forces may also contribute to interactions between GAGs and collagen [98]. Decorin, a member of the small leucine-rich PG family, has an arch-shaped solenoid-like structure formed by leucine-rich amino acid repeats. Decorin regulates cellular growth, growth factor expression and collagen fibrillogenesis through binding ends of collagen helices to maintain interfibrillar spacing via dermatan sulfate side chains [99]. However, depletion of GAGs with chondroitinase ABC in fascicles from human patellar tendon did not significantly change tendon modulus, relative energy dissipation, peak stress or peak strain [100]. Additional studies provide further evidence that GAG chains of small leucine-rich PGs do not play a role in regulating dynamic elastic or viscoelastic behaviour of rat tail tendon fascicles, and that GAGs do not facilitate collagen fibril sliding under tensile loading [5,41,101]. Currently, experimental results evaluating the interactions between different tendon constituents are conflicting, and clear characterization of the role of many of these structures is still unknown.

To further examine how applied force is distributed in collagen fibrils or transferred between different adjacent fibrils, two basic theories have been proposed. The *global load-sharing model* hypothesizes a perfect, no-slip bond of the fibre/matrix interface and negligible contact among fibres. In this formulation, the loading modules are unbranched fibres and non-fibrillar matrix with a uniform strain [102,103], and the representative volume of material is developed to guarantee transverse homogeneity [92]. Different mechanical parameters (e.g. Young's modulus and Poisson's ratio) were used to describe collagen fibres and PG matrix as ground substance, and structural features such as crimp or helical organization were incorporated using PDFs as mentioned above. In these models, the mass ratio of collagen fibres to ground matrix tended to correlate with tendon strength. In one study, the load redistribution in a partially torn or lacerated tendon was evaluated using the global load-sharing model and demonstrated that applied load was equally supported by the surviving fibres within the tissue [102–104]. Compared with the shear lag model discussed below, the global load-sharing model overestimates the strength of remaining intact fibres owing to ignoring the local effect of laceration.

In contrast, the *shear lag model* hypothesizes that the mechanism regulating tendon strain transfer across different lengths scales is interfibrillar shear load transfer between discontinuous fibres [86,105,106]. The shear lag model has been shown to provide good fits to experimental strain–stress curves in the toe, linear and even failure regions [22,86]. Compared with finite-element (FE) models (discussed more below), this model more efficiently showed the deformation around fibres because it took cylindrical fibres as one-dimensional spring elements (with constant axial stress and no radial stress perpendicular to fibre alignment) and treated the non-fibrillar matrix as a component bearing no tension [107]. Additionally, the perfect side bonding between the non-fibrillar matrix and collagen fibres facilitates continuous strain, but without any bonds between the fibre ends and non-fibrillar matrix [102]. In these types of models, several different formulations have been used to represent the interfibrillar matrix [86]. For example, the interfibrillar matrix was assumed to be linearly elastic, such that axial displacements of the matrix relative to the bonding region between the matrix and fibres were linearly related to the matrix shear stress [22,102]. In contrast, plasticity of the interfibrillar matrix suggested that a constant interfibrillar shear stress should be applied if relative sliding occurs between adjacent fibrils, but no interfibrillar stress if sliding between fibrils does not occur [86,87]. The combined elastoplastic shear lag model showed initial elasticity of the interfibrillar matrix, but a plastic response after the interfibrillar shear stress gradually increased beyond a threshold value with increasing sliding [86]. In another study, a shear lag model was developed wherein the interfibrillar matrix was modelled as piece-wise linear springs, which stretched with zero stiffness under the critical strain, then further uncoiled with substantially increased stiffness beyond the critical strain [22]. This model incorporated three-dimensional effects such that fibrils interacted with different numbers of adjacent fibrils and altered force transfer following GAG depletion [22]. One conclusion of this work was that the mechanical properties of tendon do not change with the deletion of GAGs if fibril lengths are significantly larger than the characteristic length scale relative to GAG spacing [22].

Similar to examples of normal (intact) tendons described above, several shear lag models have been developed to examine load redistribution in partially lacerated tendons [102,108–110]. Hedgepeth's damage model assumed that transverse displacement was independent of the longitudinal equilibrium equation, that fibres broke transversely, and that fibrils and interfibrillar matrix could be approximated by small deformation elasticity [102,110]. Compared with this model, Wagner and Eitan's model, which was upgraded from Cox's shear lag model by linearly superposing the stress increases caused by many broken fibres, better fitted the experimental properties of partially lacerated porcine flexor tendons [102,103,105]. Instead of a shear lag damage model, Natali *et al*. [111,112] also developed a hyperelastic constitutive model to analyse mechanical behaviour of healthy tendons during physiological loading and define the damage phenomena by a Helmholtz free energy function.

Besides shear lag models, fibre–matrix interactions in tendon have been incorporated into microstructural models in different ways. A nonlinear fibre-reinforced strain energy model was formulated to include the effects of isotropic interfibrillar matrix, interactions between fibre stretch and bulk matrix strain, fibre stretch and interactions between deviatoric fibre and matrix deformations [113,114]. The incorporation of fibre–matrix interactions in this model led to better fits to experimental data, and suggested that fibre–matrix interactions contribute to longitudinal tendon stress–strain curve in both the toe and linear regions [113,114]. Using a different approach, a dissipative theory was applied to account for temporary reversible interfibrillar bridges between collagenous and non-collagenous components in interpreting tendon softening and non-recoverable strains in cyclic mechanical testing [76,115]. In this case, the time-dependent strain energy function was taken as the sum of the viscoelastic strain energy functions for the anisotropic non-fibrillar matrix and collagen fibrils [91,115]. The breaking and reformation of active regions in both matrix and fibrillar components were represented by the transient network theory [76] and contributed to the model's constitutive equation [76,114]. This dissipative model was able to reproduce the complex time-dependent response of tendon and was used to explore microscopic mechanisms of how temporary realignment in extracellular matrix can contribute to viscoelastic behaviour [116–118]. These fibril scale models using the microstructural approach have attempted to explore important roles of different tissue constituents on tendon mechanical properties. Although their assumptions/conclusions need further validation by experimentation, they bridge studies at the levels of tissue and collagen molecules, and facilitate evaluation of collagen behaviour under different loading regimes.

### 2.3. Mathematical modelling on the molecular scale

To fully evaluate the mechanical behaviour of tendon, including on tissue and fibre scales, it is necessary to understand the mechanical properties of individual tissue constituents, such as collagen, PGs and the interaction between them on the molecular scale. Viscoelastic properties of collagen molecules were evaluated from the level of amino acids using steered molecular dynamics simulations using *in silico* creep tests [119,120]. A Kelvin–Voigt model was fitted the to the experimental stress–strain response of singe collagen molecules [120]. In this molecular context, the physical interpretation of the Kelvin–Voigt model was that the elastic spring represented the dihedral angle and the bond deformation of the protein backbone, whereas the viscous dashpot described the breaking and reforming of hydrogen bonds between collagen chains [120,121]. Using such an approach, the values of both viscous parameters and the elastic modulus of the collagen molecules were found to be significantly greater than corresponding parameters on the fibril length scale [119,120,122]. To evaluate these molecular mechanisms further, the influence of representative divalent and trivalent cross-links in tropocollagen molecules on the mechanics of collagen fibrils was examined through a three-dimensional coarse-grained model [123–125]. It was found that fibrils with certain amounts of divalent cross-links had similar stiffness and stress levels to fibrils with half the amount of trivalent cross-links, and that mechanical properties of cross-links contributed to fibril behaviour in both the linear and failure regions [125].

Regarding non-collagenous components of tendon, the stiffness of chondroitin-6-sulfate, a typical GAG attached to decorin, has been analysed through the molecular mechanics approach [126–128]. The potential energy and physical properties of GAGs were obtained using analytical functions of potential energies from direction bonding and non-bonding interactions of atoms. The piece-wise linear stiffness of GAGs, which shows an initial plateau region at low strains followed by a high stiffness region, was the second-order derivative of the energy with respect to the GAG length, normalized by its molecular weight after interpolation of energy data [127]. A multi-fibril model, which introduced the calculated GAG stiffness through the molecular mechanics approach, reproduced experimental results of fibre incremental elastic modulus in mature tendon [127,128]. On this scale, single collagen molecules and GAGs were also shown to be viscoelastic via molecular mechanics approach, which is beneficial for conformational analysis to investigate the equilibrium configuration of molecular chains. However, unclear understanding of the collagen/fibril assembly architecture prevents accurate characterization of tendon properties with mathematical modelling. Computational modelling might provide an efficient method to simulate a variety of cross-links/connections between components on different length scales.

## 3. Computational modelling

As the previous sections described, in certain cases, an appropriate constitutive equation combined with assumed homogeneity and simplified boundary conditions (i.e. constant stress, constant strain) will enable a model to approximate the mechanical behaviour of tendon under relatively simple *in vitro* loading regimes. However, tendon has very complicated hierarchical structures from the molecule to tissue scales, which are often heterogeneous, and may function in complex, multiaxial physiological loading environments. Additionally, the compositional and structural properties of tendons vary between different species and different developmental stages. To address these more complex aspects, many studies have used computational modelling to analyse tendon behaviour and explore structure–function relationships [129]. The FE method is a popular, and powerful tool that can accommodate many of these complex issues related to tendon mechanics. Because most computational models of tendon to date have used the FE method, this section focuses on computational models that have been developed using FE approaches to evaluate tendon mechanics on different length scales.

### 3.1. Computational modelling on the tissue scale

On the tissue scale, computational modelling using the FE method frequently includes tendons and their surrounding anatomical structures (i.e. muscle, ligament, cartilage, meniscus, bone) to analyse functions of musculoskeletal joints [46,130,131]. For example, a three-dimensional model of the glenohumeral joint of the shoulder incorporated the long head of the biceps tendon as a transversely isotropic, hyperelastic material with a constant elastic modulus. This model showed increased displacement of the labrum caused by tension of the biceps tendon [131,132]. However, the mechanical behaviour of the tendon in the joint was not fully explored, and further studies focused on tendon properties under different physiological loads should be conducted. Wren *et al*. used the poroelastic constitutive model to describe rabbit flexor digitorum profundus tendon [133]. The simulation sought to examine changes in permeability and compressive stresses of an assumed homogeneous tendon in response to local cyclic hydrostatic pressures [133]. This model determined that fluid pressures contributed to the low permeability in the adaptive fibrocartilaginous regions and were beneficial for maintaining structure–function of wrap-around tendons under cyclic compressive loading [133]. Besides permeability, other aspects of tendon mechanics such as compressive elastic modulus, shear stresses, strains on different length scales and Poisson's ratio have also been evaluated with computational modelling [31,134–136]. For example, three-dimensional models with hyperelastic constitutive properties were created to evaluate heterogeneous compressive moduli over the surface of the supraspinatus tendon [135]. In another study, the humeral head and the attached supraspinatus tendon with a graded elastic modulus were included in a two-dimensional FE model, with partial-thickness tendon tears created on different surfaces [136]. Not surprisingly, tears were found to cause stress concentrations, which also occurred near the articular surface at the tendon-to-bone insertion [136,137]. In other work, strain-dependent Poisson's ratio was predicted in biphasic simulations using a strain energy equation that could characterize the nonlinear, transversely isotropic compressive behaviour of tendon [138]. Compared with generalized Maxwell mathematical models, computational models using three-parameter Mooney–Rivlin formulations and 64-parameter Prony series were found to be more efficient because they demonstrated it was not necessary to use so many Maxwell elements to provide high accuracy to experimental data [75,134].

### 3.2. Computational modelling on the fibre/fibril scale

Compared with the method by which fibril orientation is implemented into mathematical models using PDFs (described earlier), crimp and helical structures can be incorporated into computational models of tendon by modifying fibril subunits with sinusoidal and helical structural transformations [139]. Based on the observation that crimp period is constant along the fibre/fascicle direction and that the helical pattern follows a mean fibril pitch [93,140–142], a sinusoidal crimp waveform and a helical twist are capable of endowing relevant structural features to computational fibres [139]. In one study, an isotropic, neo-Hookean hyperelastic constitutive equation was proposed to model both fibrils and interfibrillar matrix [139]. After development of hexagonally shaped unit fibres and FE analysis of the unit fibres under tension, the model not only showed nonlinear stress–strain behaviour of fibres with different crimp angles, but also indicated that mean helical pitch, crimp angle and the ratio of fibre modulus to interfibrillar matrix determined Poisson's ratio [139]. The interactions and mechanisms of stress transfer between collagen fibrils and interfibrillar matrix in tendon are still under study in computational modelling. In one study, Goh *et al*. [143] constructed a computational model of connective tissues as a tapered fibre surrounded by a coaxial cylinder of PG gel. Stress transfer was incorporated by applying constraints on FEs and adjusting the Young's modulus of collagen relative to the surrounding PG gel. Significantly, more slender fibrils and large collagen/PG moduli differences were found to contribute to high stress in fibrils and low probability of interfacial rupture [143].

### 3.3. Computational modelling using multiscale approaches

#### 3.3.1. Modelling of tendon

Different from experimental methodology, computational modelling has the capacity to resolve the challenge of coupling the multiscale relationships among joint function, tendon macroscopic mechanical properties, microscopic structure of extracellular matrix and tenocyte mechanotransduction. Unfortunately, very few multiscale models have been developed and used to study tendon properties across different length scales. In one example, a computational cell composed of hyperelastic fibril bundles, interfibrillar matrix and embedded tenocytes was created, and FE analyses were conducted to investigate the cause of crimp in tendon [144,145]. In this study, shear stresses on the interface between collagen fibrils and interfibrillar matrix—attributed to cell contraction—caused crimp structure, and the ratio of fibrils to interfibrillar matrix determined crimp length [144]. Although there have been few multiscale models of tendon, a few models that have been developed to evaluate other tissues provide some insight into potential options for future work in computational modelling of tendon.

#### 3.3.2. Modelling of other biological materials

Multiscale constitutive models including structural, compositional and mechanical information have been used to address strain and stress transfer in different tissues, such as vascular tissue [146,147], meniscus [148], cartilage [149] and engineering collagen-based gels [150,151]. Instead of implementing a macroscopic constitutive law for the whole tissue, vascular tissue was modelled as PG-cross-linked collagen fibrils with a piece-wise analytical stress–stretch response [146]. Macroscopic volume-averaged Cauchy stresses for vascular tissue were computed by integrating corresponding collagen Cauchy stresses with volume averaging theory, which realized stress transfer from fibril level to tissue level [146,147]. To resolve the problem that different local micromechanical environments of fibrochondrocytes are influenced by meniscus tissue mechanics, a compressive normal load was applied macroscopically to nonlinearly elastic, anisotropic and biphasic axisymmetric FE elements to predict strains in different directions; macroscale strains were transferred to the microscale mesh as a boundary condition input [148]. To demonstrate the effects of collagen and PG distribution on knee joint mechanics, Halonen *et al*. [149] varied the value of mechanical parameters in different locations of a three-dimensional FE model, which could predict stresses or strains during dynamic joint loading. Based on the concept of evaluating continuum tissue-level responses as a function of properties of the collagen network in order to form a link between nanoscale collagen features and tissue mechanics [152], the interaction among physiological loading, microscopic fibre realignment and cellular responses can be coupled and examined using a multiscale modelling approach [151,153,154]. Although different biological materials have varying structures, components, and function in unique *in vivo* loading environments, the methodology applied to develop multiscale models for other biological materials could provide motivation and inspiration for developing computational studies of tendon across different length scales.

## 4. Challenges and remarks

As an essential member of the musculoskeletal system, tendons function in a complex multiscale physiological environment. Injuries are very common in many tendons (e.g. rotator cuff, Achilles and flexor/extensor tendons), and can lead to tendon degeneration, pain and altered joint kinematics [155]. Recent experimental studies have explored some important aspects of tendon mechanical properties, mechanotransduction and mechanoadaptability [11,84,102]. However, much about tendon mechanics and structure–composition–function relationships remains to be understood, in both healthy and injured tissues. Computational modelling offers an approach that can provide insights into tendon relationships that are impossible to measure *in vivo*. In addition, models can be beneficial in handling large amounts of data and complicated loading situations, including interactions with neighbouring anatomy. Appropriately defined models could be formulated for a specific tendon with its unique properties (e.g. multiaxial mechanical properties, location-specific microstructural organization, inhomogeneous distribution of non-fibrillar matrix, site-specific geometry, etc.) and, using a multiscale approach, could be constructed to enable mechanical evaluation across the full range of length scales relevant to the hierarchical organization of tendon. Importantly, well-validated models enable predictions of tissue behaviour under different loading conditions, including those that are difficult to explore experimentally, and can be used to study the effects of disease, degeneration and/or injury on functional mechanical properties.

Although some studies have developed mathematical or computational models to evaluate tendon behaviour, there are still some challenges to resolve for improving the efficiency and predictability of these modelling approaches. More advanced technology and methodology should be incorporated to accurately consider inhomogeneous material properties of tendon. Elastic, viscoelastic and poroelastic constitutive equations have been formulated to characterize mechanical properties of tendon on the tissue scale, collagen fibrils on the microscale or interfibrillar matrix, but a well-accepted standard constitutive material model is still lacking to more accurately address the unique multiaxial behaviour, remodelling and damage of tendons. To the best of our knowledge, most mathematical/computational models aim to reproduce or predict tendon stress or strain response under specific loading types and they fail to characterize tendon behaviour under different loading scenarios. Additionally, there have been very few studies focused on multiscale transfer and correlations of tendon structures and properties. For example, computational models of tendon could be used to resolve how strains and stresses are passed down from macro- to micro- and nanoscale, and to explore mechanisms regulating this unique hierarchical structure spanning multiple length scales. The lack of progress in this field can be attributed to the complexity and difficulty to incorporate biological, mechanical and chemical factors in different spatial, physical and temporal scales. Simulation-based interpretation of governing mechanisms could facilitate the determination of functional properties of microstructural components in tendon, many of which have not yet been addressed. Useful information extracted from multiscale modelling can assist future patient-oriented treatment and even soft tissue replacements in tissue engineering.

## Authors' contributions

F.F. and S.P.L. have made substantial contributions to drafting and critically revising the paper. Both authors have read and approved the final submitted manuscript.

## Competing interests

We declare we have no competing interests.

## Funding

We received no funding for this study.

## Footnotes

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

- © 2015 The Author(s)