Exploiting Viscoelastic Experimental Observations and Numerical Simulations to Infer Biomimetic Artificial Tendon Fiber Designs

Designing biomimetic artificial tendons requires a thorough, data-based understanding of the tendon's inner material properties. The current work exploits viscoelastic experimental observations at the tendon fascicle scale, making use of mechanical and data analysis methods. More specifically, based on reported elastic, volumetric and relaxation fascicle scale properties, we infer most probable, mechanically compatible material attributes at the fiber scale. In particular, the work provides pairs of elastic and viscous fiber-scale moduli, which can reproduce the upper scale tendon mechanics. The computed range of values for the fiber-scale tendon viscosity attest to the substantial stress relaxation capabilities of tendons. More importantly, the reported mechanical parameters constitute a basis for the design of tendon-specific restoration materials, such as fiber-based, engineering scaffolds.


INTRODUCTION
Tendons are natural fibrous tissues that transfer mechanical loads from the muscles to the bones. They are structured in a highly hierarchical manner (Maceri et al., 2012), consisting of a series of inner fibrillar scales that are immersed in a matrix substance (Shen, 2010;Zhang et al., 2014). As for their structural arrangement, both human and animal tendon specimens have been commonly described as a multiscale composite materials. The tendon unit consists of fascicles (at the scale of hundreds of micrometers), which are in turn composed of matrix-immersed fibers (micrometer). Figure 1 provides a schematic representation of the tendon's inner multiscale architecture (Figure 1) (Goh et al., 2008;Thorpe et al., 2014).
The tendons' inner fibrillar components are not parallel to the length evolution of the tendon. On the contrary, they are structured in helical patterns (Orgel et al., 2006), forming an undulated inner structure; an observation reported in different microscopy based studies (Yahia and Drouin, 1989;de Campos Vidal, 2003). Their helical arrangement, with a typical angular range in between 70 and 76 • with respect to the plane normal to the tendon axis evolution (Järvinen et al., 2004;Starborg et al., 2013), leads to a characteristic coupled axial and torsional behavior at the fascicle scale (Thorpe et al., 2013). The effective mechanical behavior at the tendon macroscale, arises as a rather complex function of the properties of its inner constituents (Reese et al., 2013). The latter have been commonly measured by elastic, uniaxial strain experiments, carried out independently at the different inner tendon scales (Figure 1). At the lower scale of fibers, elastic moduli values of E f = 0.57 ± 0.08 GPa and E f = 2.69 ± 0.42 GPa have been reported for wet and dry rat tail tendon fiber specimens, respectively (Kato et al., 1989). Independent experimental studies provided moduli values of E f = 1.17 ± 0.28 GPa for hydrated fiber tendon specimens (Gentleman et al., 2003), within the range reported by Kato et al. (1989). At the uppermost inner hierarchical scale of fascicles, experimental data suggest a substantially lower overall uncertainty. In particular, elastic moduli values of E fasc = 0.64 ± 0.03 GPa, E fasc = 0.48 ± 0.07 GPa, and E fasc = 0.55 ± 0.14 GPa have been reported (Derwin and Soslowsky, 1999;Lavagnino et al., 2005;Haraldsson et al., 2009;Svensson et al., 2010), exhibiting a maximum difference of 0.2 GPa with respect to their mean value.
The overall decrease observed in the elastic moduli of the upper scales is to be primarily attributed to the presence of the non-collagenous matrix substance (Figure 1). The latter has a quasi-negligible stiffness contribution with respect to the one provided by the fibrillar components. While no direct experimental measurement is available, analytical computations have estimated a matrix modulus E m below 1 MPa (Ault and Hoffman, 1992) that is more than two orders of magnitude lower than the one reported for any fibrillary component measurement (Thorpe et al., 2017); a value that has been typically employed in numerical studies (Reese et al., 2010). The content of fibers f r within fascicles (Figure 1) varies with the tendon's age and health state (Lavagnino et al., 2005;Svensson et al., 2010). Relevant studies provide fibrillar contents f r as low as 35% when a certain level of fiber swelling is accounted for (Svensson et al., 2012) and up to more than 60% (Hansen et al., 2010;Svensson et al., 2010).
The combination of high moduli differences amongst the embedding matrix and the fibrillar components (E f >> E m ), along with the helical arrangement of the tendon subunits have been shown to constitute the basic mechanisms responsible for the tendon's distinctive volumetric behavior (Reese et al., 2010;Swedberg et al., 2014;Karathanasopoulos et al., 2017). More specifically, experimental observations have reported Poisson's ratio values ν that well-exceed the ones observed in common engineering materials. In particular, mean Poisson's ratio values close to unity (Cheng and Screen, 2007) and up to 3 have been reported (Lynch, 2003), accompanied by substantial uncertainties.
The linear elastic response cannot however explicate the considerable stress and shock absorption capacities of tendons (Salathe and Arangio, 1990), which are directly associated to a viscous, time-dependent behavior. The tendon's viscoelastic nature allows for the attenuation of stress stimuli and is commonly characterized by a viscosity parameter η, which is directly related to a relaxation time τ (Lakes, 2009), as schematically depicted in Figure 2. Tendon inner scales have been shown to yield a viscoelastic, time-dependent response, which however has not been typically characterized with respect to its effective viscosity value η. In particular, at the tendon fascicle scale, relaxation curves with a time-dependent modulus evolution over a range of 200 s and up to more than 500 s have been provided in independent experimental studies (Screen, 2008;Davis and De Vita, 2012) (Figure 3D). Relaxation experiments suggest a loss of 40-70% of the initial elastic modulus at the tendon fascicle scale over this time interval (Screen, 2008).
The material properties of the tendon building blocks considerably affect its functionality and its load transfer efficiency (Rawson et al., 2015); parameters of primal importance for connective tissues (Galbusera et al., 2014). A thorough quantification of the material properties of the tendons' inner constituents is essential, not only for the understanding of its nature (Balint et al., 2014), but most importantly for any repair and restoration process to take place (Robinson et al., 2004;Laurent et al., 2014;Lee et al., 2017). Up to now, restoration processes have primarily used biological and synthetic scaffolds (Kuo et al., 2010;Kwon et al., 2014;Sandri et al., 2016). Biological treatments employed biodegradable silk-collagen scaffolds that were meant to provide advanced regeneration capabilities (Abdullah, 2015). Synthetic replacements were based on textile scaffolds selected out of a list of existing substitutes (Lomas et al., 2015). However, scaffolds of the kind have been reported to show limited mechanical biocompatibility (Ratcliffe et al., 2015;Tang et al., 2018). Statistical data on the success of surgical restorations reveal rather low success rates. More specifically, for tendons injuries with large and extensive damage, the mean post-surgical, operation success rate does not exceed a mere 50 and 40%, accordingly (Meimandi-Parizi et al., 2013).
A part of the low efficiency of current restoration practices needs to be attributed to the rather limited understanding of the tendon's inner material properties (Kuo et al., 2010), in particular with respect to its viscoelastic properties (Ganghoffer et al., 2016). The latter depend both on the physiology and on the loading type applied (e.g., quasi-static, high strain rate) (Oftadeh et al., 2018;Wu et al., 2018;Kuznetsov et al., 2019). The use of biocompatible materials, which mimic the physiological functionality of the native tissue constitutes a key aspect for any tendon restoration process (Ganghoffer et al., 2016). During the FIGURE 3 | Schematic representation of the tendon fascicle's helical geometry (A) that is composed of matrix-embedded fibers (B) (Goh et al., 2008), with a Poisson's ratio behavior (C) that is dependent on the fascicle's structural composition (B) and geometry (A). In subplot (D), the tendon fascicle's experimental relaxation curves as provided by Screen (2008) (Sc) and Davis and De Vita (2012) last years, substantial efforts have been devoted to the engineering of novel biomaterials with enhanced structural performance and improved biochemical compatibility (Zhang et al., 2017;Li et al., 2018a,b;Li S. et al., 2018;Liu et al., 2018). However, primal tendon inner material attributes, such as the effective viscosity at the tendon fiber scale remain unquantified.
In the current work, numerical models that describe the tendon's fascicle and fiber scale mechanics are combined with probabilistic inference models and experimental observations. By that means, a quantitative link between the tendon's viscoelastic mechanical response at the fascicle scale and experimental data is established, through a Bayesian inference framework (Papadopoulos et al., 2012;Taflanidis and Beck, 2013;Farrell et al., 2015). This allows for the quantification of uncertain tendon inner mechanical parameters that determine the experimentally observed fascicle-scale mechanics. The paper is structured as follows: in section Materials and Methods the parametrization of the fascicle's viscoelastic models is provided and the mathematical formalism of the Bayesian probabilistic framework is discussed (section Inference of the Fascicle's Viscoelastic Mechanical Properties), summarizing the experimental data used for the study (section Tendon Fascicle Experimental Data). In section Results, the results of the inference process are provided, namely the elastic and viscous fiber-scale material properties, which can reproduce the experimental observations reported at the fascicle scale. The work continuous with an overall discussion of the obtained results and a concluding summary in section Discussion and Conclusions.

Tendon Fascicle Geometry and Mechanical Models
We describe the fascicle's elastic and viscoelastic relaxation response in a two-step process. More specifically, we compute its elastic response with a dedicated composite helical fascicle finite element model, which is comprised of fibers immersed in a matrix substance, detailed in Karathanasopoulos et al. (2017). The fascicle geometry follows a helical angle θ , with respect to the plane perpendicular to the tendon evolution, as schematically depicted in Figure 3A. The fascicle contains fibers in different fiber contents f r defined as the ratio of the area A f covered by fibers (dark-gray) over the total fascicle cross sectional area A t ( Figure 3B). We allow for the composition of fascicles to vary so that different fiber contents are captured. More specifically, in order to take into consideration a large part of the experimentally observed fascicle compositions (Goh et al., 2008;Hansen et al., 2010;Svensson et al., 2012), we build fascicle models with fiber contents of 35% up to 65% with a spacing of 5%, accounting for a fiber swelling in the determination of the upper content limit (Hansen et al., 2010). For the linear finite element computations, the fascicles' domain covered by fibers is assigned a Young's modulus E f , while the matrix a modulus E m . The model computes the effective fascicle modulus E fasc and the fascicle's volumetric behavior ν fasc = R/R/ε z ( Figure 3C) for different fiber content values f r and angular arrangements θ (Figure 3A), which are considered to define distinct fascicle model classes M θ f r . Both parameters are complex functions (f ) of the tendon fascicle's inner material and geometric properties f f r , θ , E f , E m (Reese et al., 2010(Reese et al., , 2013Karathanasopoulos et al., 2017). We enumerate a total of 49 fascicle model classes The reader is referred to Karathanasopoulos and Hadjidoukas (2019) and Karathanasopoulos et al. (2017) for a detailed description of the numerical modeling specifications and for a quantification of the associated computational cost.
Following the linear computation, the relaxation response of the tendon fascicle is computed, using a Maxwell standard linear solid model (Figure 2). The timedependent relaxation curve is characterized by a relaxation experiment, using a single viscosity parameter η, as follows (Christensen, 1982;Lakes, 2009;Shen, 2010): where in Equation (1), E R is equal to the fiber's elastic modulus loss E R = E f (t = 0) − E ∞ during the relaxation process (Figure 2). For t = 0, the fascicles's elastic modulus is equal to its linear, non time-dependent value, so that the previous relation entails Noting that the matrix modulus E m has been shown to be more than two order of magnitudes lower than the one of its fibrillar components (E f >> E m ) (Ault and Hoffman, 1992;Reese et al., 2010;Karathanasopoulos et al., 2017), the primal contributors in the fascicle's relaxation curves ( Figure 3D) described by Equation (1) are its embedded fibers (E m (t) ≈ 0). As a result, the fascicle's time-dependent modulus loss is essentially characterized by the relaxation behavior of its inner, matrix-embedded fibers, so that the viscous parameter η entering Equation (1) describes the effective -homogenized-viscosity of its embedded fibers, here named as η f . We subsequently parametrize each fascicle model class M θ f r by a total of three parameters, namely by the elastic modulus of the fiber and of the matrix E f and E m and by the viscous modulus of its embedded fibers η f . We note that the fascicle's effective Poisson's ratio value ν fasc is a non-linear function of its fiber content f r , angle θ (thus M θ f r ) and fiber and matrix elastic moduli values E f and E m , as elaborated in Appendix section Fascicle Poisson Ratio. Accordingly, the fascicle's time-dependent response well differs, depending on the combination of the elastic and viscous properties entering Equation (1), as demonstrated in Appendix section Fascicle Relaxation Response. Using the previously defined parameters, we compute a total of three quantities of interest for each model class M θ f r , as follows: where in Equation (2),Ē fasc (t) stands for the time-evolution of the normalized fascicle modulus, the normalization carried out with respect to its initial elastic modulus E 0 fasc .

Inference of the Fascicle's Viscoelastic Mechanical Properties
The goal in parameter inference is to infer the parameters ϕ ∈ R N ϕ after observing the data D = {d i |i = 1, . . . , N}.
The computational model used to simulate the data D can be described through a function F(x; ϕ) ∈ R N , which depends on both the parameters ϕ and on input parameters x ∈ R N x that are not being inferred.
In Uncertainty Quantification we are interested in sampling the conditional posterior distribution p (ϕ |D ), using the Bayes' theorem, defined as follows: where in Equation (3)  The likelihood function is a measure of the probability that the data D are reproduced by the computational model employed, while the prior probability encodes all available information before observing any data. In the present work, we make the assumption that the data D i are independent and normally distributed around the observable of the model, with a proportional model error, so that: where the data D i are assigned a proportional error model. Using Equation (4), the likeli-hood function p (D|ϕ, M) takes the form: Where the model and the error parameters can be described by the vector ϑ = ϕ T , σ n T in a unified form. In most practical applications the posterior distribution (Equation 3) cannot be expressed analytically. Moreover, the model evidence p (D |M ) is typically not known, since it is given by integration of the nominator of Equation (3) over the potentially high dimensional domain of the parameters. In the current work, we rely on efficient sampling algorithms to identify the posterior Frontiers in Bioengineering and Biotechnology | www.frontiersin.org distribution. In particular, we use the Transitional Markov Chain Monte Carlo algorithm (TMCMC) (Ching and Chen, 2007;Farrell et al., 2015). The algorithm starts by sampling the prior distribution, which is usually trivial to be sampled, and through annealing steps it provides samples from the posterior distribution. One major advantage of this algorithm is its ability to sample multimodal distributions and provide low bias estimators of the model evidence (Ching and Chen, 2007). Here, we make use of uniform prior distributions U for each of the modeling parameters ϕ = E f , η f , E m which encompass and exceed the range reported in the corresponding experimental observations of section Introduction, indicating prior ignorance with respect to their actual value. In particular, we allow their values to vary in between 500 − 2500 MPa, 0.1 − 2500 GPa s and 0.01−5 MPa. A graphical representation of the probability model formulation is provided in Figure 4.

Tendon Fascicle Experimental Data
For the inference process we combine multiple available experimental data, pertaining to both the linear and the timedependent fascicle response. In particular, for the quantities of interest q 1 and q 2 of Equation (2), we made use of the elastic, non time-dependent experimental data reported in Derwin and Soslowsky (1999), Haraldsson et al. (2009) and Svensson et al. (2010);and Lynch (2003) and Cheng and Screen (2007), accordingly. What is more, for the time-dependent quantity q 3 of Equation (2), we digitalize the experimental relaxation curves provided in Screen (2008) and Davis and De Vita (2012) and summarized in Figure 3D, using ∼30 equally spaced timeintervals along the relaxation curve, as control points for our modeling predictions. We consider a total of five groups of experimental observations D 1 , D 2 , D 3 , D 4 , and D 5 . Each data set D i contains a set of normalized fascicle moduli ratiosĒ fasc (t), as schematically depicted in Figure 3D. In Table 1, we summarize the considered experimental measurements for each group of experimental observations.

RESULTS
We compute the posterior sample distributions using the TMCMC algorithm elaborated in Ching and Chen (2007) for all model classes M θ f r , for each of the data sets D i included in Table 1. In Figure 5 we provide the posterior frequency distribution for each of the inferred parameters ϕ = E f , η f , E m for the model M θ f r with a fiber content f r = 35% and an angle θ = 75 o . The results correspond to the data set D 1 of Table 1. Figure 5 depicts a posterior distribution of the modeling parameters ϕ with a clear clustering of values for each of the E f , η f , E m . In particular, for each parameter, the range of values with a non-zero posterior probability is a narrow subspace of the uniform prior used as initial modeling hypothesis and summarized in Figure 4. More specifically, the fiber modulus E f ranges in between 1,300 and 1, 500 MPa, while the matrix modulus E m within 0.35 and 0.5 MPa. Accordingly, the entire probability mass for the viscosity parameter η f ranges among 35 and 45 GPa s. The model proportional error (Equation 4) associated to the results of Figure 5 is restrained to a maximum of 8% for all model classes and data sets D i . Analogous posterior parameter frequency distributions with the one provided in Figure 5 are obtained for all model classes M θ f r and data sets D i of Table 1.
In Figures 6A,B we provide the values E f and E m that maximize the posterior PDF for each respective model class M θ f r for the data set D 1 . These values, named as the Maximum Figures 6C,D, we provide the values for E f and E m for the data set D 2 , pertaining to a fascicle Poisson's ratio value ν = 3 ( Table 1).
Figures 6A,C suggest a non-linear relation between the fiber modulus E f and the fiber content f r , so that the higher the fiber content value, the lower the most probable fiber modulus value obtained. In particular, a maximum and a minimum most probable fiber modulus of approximately 1550 MPa and 900 MPa is obtained for a 35% and a 65% fiber content, accordingly. The fiber angle θ affects the magnitude of the E f value: higher angle values θ correspond to lower fiber moduli E f , irrespective of the fascicle's fibrillary content f r . The associated difference lies however within a maximum range of 150 MPa for a given fiber content value (Figures 6A,C).
The different fascicle Poisson's ratio values (ν fasc in Table 1) used among the data sets D 1 and D 2 , appear to minorly affect the peaks of the posterior probability distribution for E f (Figures 6A,C), while they decisively affect the E m MAP values (Figures 6B,D). In particular, the higher fascicle Poisson's ratio value of data set D 2 results in considerably reduced E m MAP values with respect to the ones computed for the data set D 1 . However, the dependence of E m on f r and θ is analogous to the one obtained for the fiber modulus E f . The a-posteriori, most probable inferred fiber moduli values E f well compare to the ones suggested by the experimental study of Gentleman et al. (2003), while the matrix modulus values E m to the analytical modeling predictions of Ault and Hoffman (1992), being smaller than 0.7 MPa in the entire parametric space. For the data sets D 3 to D 5 of Table 1, fiber and matrix moduli values that differ by a maximum of 5% with respect to the ones provided in Figures 6A,B   1 | Fascicle experimental data sets D i encompassing the fascicle initial elastic modulus E fasc , Poisson's ratio value ν fasc , and the time-dependent response of its modulusĒ fasc (t).
E fasc (t = 0) (Derwin and Soslowsky, 1999;Haraldsson et al., 2009;Svensson et al., 2010) ν fasc (−)Ē fasc (t) (−) D 1 640| 480| 550 1 (Cheng and Screen, 2007) Figure 3D (Sc) (Screen, 2008) D 2 640| 480| 550 3 (Lynch, 2003) Figure 3D (Sc) (Screen, 2008) D 3 640| 480| 550 1 (Cheng and Screen, 2007) Figure 3D Table 1 relate to a range of η f values in between 40 and 80 GPa s. For all data sets, low fibrillar contents f r pair to high relaxation moduli values, while increased fiber angles θ to lower η f values for a given data set D i and content value f r . We note that the linear elastic, non time-dependent fascicle Poisson's ratio value ν fasc is independent of the relaxation moduli η f , so that no significant differences arise between the inferred η f values of data sets D 1 and D 2 . The inferred MAP values for the parameter E f , η f , E m are provided for all data sets D i and model classes M θ f r in the form of Supplementary Material. Each set of material properties   Table 1.
Frontiers in Bioengineering and Biotechnology | www.frontiersin.org provided in Figures 6, 7 can well-reproduce both the elastic and the experimentally observed time-dependent fascicle behavior summarized in Table 1.
In Figure 8, we provide for verification purposes a comparison of the fascicle relaxation response as computed by the previously inferred mechanical parameters and by the experimentally provided response. To that scope, we make use of the inferred MAP elastic and viscous mechanical properties part of them depicted in Figures 6C,D, 7B,D for the model class M 70 35 , which we compare with the experimental curves Dv 1 and Dv 3 , provided in Figure 3D. Figure 8 suggests a very good agreement between the experimental relaxation response and the relaxation behavior arising from the inferred mechanical properties. Analogous behavior is obtained for all model classes M θ f r and data sets D i of Table 1.We note that non-optimized sets of the parameters E f , E m , η f lead to utterly different fascicle relaxation behaviors, as elaborated and showcased in Appendix section Fascicle Relaxation Response.

DISCUSSION AND CONCLUSIONS
Knowledge of the tendon's inner material properties is a primal prerequisite for the application of any successful treatment or restoration process (Snedeker and Foolen, 2017;. However, experimental data on the material attributes of inner, lower tendon subunits are commonly highly uncertain, while fundamental parameters, such as the η f of tendon fibers remain unquantified (section Introduction). Combining upper (fascicle) and lower (fiber) tendon scale mechanics with experimental observations, provides a means to infer and quantify lower scale tendon mechanical properties, so that they are able to reproduce upper scale mechanics.
The inferred fiber scale mechanics indicate that the most probable fiber modulus values E f range between 900 and 1, 600 MPa (Figure 6), thus in a subspace of the literature reported range of moduli (Kato et al., 1989;Gentleman et al., 2003), when using experimental observations at the fascicle tendon scale. In particular, if information on the fibrillar content value f r is provided, the fiber modulus E f can be estimated, within a range of 200 MPa when the 95% of the mass of the posterior probability distribution is accounted for (Figure 5). Moreover, for the observed volumetric behaviors at the fascicle scale to be retrieved (Table 1), a matrix modulus E m below 1 MPa is required for all fibrillar content and fiber angle values model classes M θ f r . Increased fascicle lateral volumetric contractions are obtained by decreasing the matrix modulus E m (Figure 6), thus by increasing the relative contrast of elastic properties E f /E m at the fiber scale.
The viscous moduli η f values provided in Figure 7 constitute the first quantitative estimates of the embedded fiber viscosity that is based on experimental observations. For a given fascicle composition M θ f r and data set D i , the posterior distribution of the embedded fiber viscosity suggests η f values which at most differ by 10 GPa s (Figure 5), considering the entire mass of the posterior probability distribution. For fascicles with a low fibrillar content and fiber angle, most probable η f values of more than 40 GPa s are obtained, irrespective of the data-set D i considered (Figure 7). The values of Figure 7 combined with the most probable fiber elastic moduli values of Figure 6 can be viewed as a fundamental reference table in the design of artificial tendons. In particular, each pair of most probable E f , η f values can be used as a set of basic elastic and time-dependent material properties for the engineering of artificial fibers in scaffold-based tendon restoration processes (Kuo et al., 2010;Abdullah, 2015;Sandri et al., 2016).
Note that the differences observed among the most probable viscoelastic embedded fiber moduli values η f (Figure 7) reflect to a large extent the discrepancies among the experimentally reported relaxation curves of Figure 3, as well as the wide range of possible fibrillar content values f r . Figure 9 summarizes the range of values for the most probable elastic and viscous fiber parameters (Figures 6, 7), when all fiber content values f r , fiber angles θ , and data-sets D i are considered. In order to further delimit the range of probable elastic and viscous parameters, further relaxation experiments need to be conducted at the tendon fascicle scale. The latter need to provide mean and standard deviation values for the relaxation moduli at different time-frames throughout the relaxation process at different strain magnitudes and initial loading strain-rates. Data of the kind will allow for a series of secondary analysis to be carried out, which are as of now intractable. In particular, they will allow for the consideration of more complex relaxation models with multiple relaxation times or for the modeling of strain-rate and strain-magnitude effects (Oftadeh et al., 2018). Moreover, they would allow to account for the presence of geometric or material non-linearities, phenomena which typically play a role in the mechanical behavior of biological materials (Kuznetsov et al., 2019); objectives which are way beyond the current tendon mechanics-related data availability.

NK:
conception, design, results acquisition, and analysis. J-FG: modeling, analysis, and interpretation of results.

ACKNOWLEDGMENTS
NK acknowledges support of the Freenovation 2017 Grant, Switzerland and the support of the ETH CSE-Lab in the writing of the corresponding application.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.

2019.00085/full#supplementary-material
Data Sheet 1 | Summary of the MAP values for the datasets D 1 -D 5 of Table 1.

APPENDIX Fascicle Relaxation Response
The time-dependent relaxation response of tendon fascicles depends on both their elastic properties E ∞ , E R as indicated by Equation 1 on their viscosity η f . In Figure A1 we provide an insight in the role of the effective viscosity parameter η f on the form of the relaxation curve. To that scope, we compute the relaxation response of a tendon fascicle with an initial zerotime (t = 0) elastic modulus E fasc (0) = 600 MPa (with an E f = 1025 MPa E m = 0.4 MPa f r = 0.65 θ = 72 o ) and a relaxing elastic modulus E R = 720 MPa for viscosity values of η f = 10 GPa s η f = 50 GPa s and η f = 100 GPa s. Figure A1 suggests that increasing the viscosity parameter η results in higher time-dependent moduli values E fasc (t), with the maximum difference to amount to several hundreds of MPa. As an example, for η = 10 GPa s at t = 100 s the structure has relaxed to E fasc (100) = 240 MPa, while for η = 50 GPa s the corresponding value is E fasc (100) = 420 MPa.

Fascicle Poisson Ratio
The effective Poisson ratio of the fascicle ν fasc hinges on both the material properties of the fiber and the matrix (E f , E m ) and on its geometric arrangement (f r , θ ) [21, 25]. As a result, different volumetric fascicle behaviors are obtained depending on the set of parameters selected. In Figure A2, we provide the effective Poisson's ratio values for a case study with a fibrillar content f r = 0.6 over different helix angles θ and for different ratios of the FIGURE A1 | Relaxation response for an initial elastic modulus E fasc = 600 MPa with a relaxing elastic modulus E R = 720 MPa for different viscosity values η. fiber modulus to the matrix modulus E M ratio = E f /E m , computed for a matrix modulus value of E m = 1. Figure A2 suggests a highly non-linear relation between the moduli ratio value E M ratio and the fascicle's Poisson's ratio value, which depends on both the fibrillar content f r and on the fiber angle θ .