The application of standard nonlinear solid material models in modelling the tensile behaviour of the supraspinatus tendon

: Tendons transmit forces from muscles to bones through joints. Typically, tendons and muscles work together to innovate a particular type of motion. Therefore, in order for the tendons to find attachment to the bones, they are naturally adapted as much thinner strands than the muscles that they serve. Thus, they are often subjected to much higher stresses than the muscles that they actually serve in any given action. As a result, tendons are susceptible to injuries that may lead to a permanent dysfunction in joint mobility due to the fact that the scar tissue that forms after healing does not often have the same mechanical properties of the original tissue. It is, therefore, very important to understand the mechanical response of tendons. This paper examines the performances of two viscoelastic standard nonlinear models in modelling the elastic and plastic behaviour of the tendon in the light of a well-known hyperelastic Yeoh model. The use of the Yeoh model is more for validating the performances of the viscoelastic models within the elastic region rather than for com-parison purposes. Yeoh model’s selection was based on its superior performance in modelling the elastic phase of soft tissue as reported in previous studies combined with its simplicity. The results show that the two standard nonlinear solid models perform extremely well both in fitting accuracies and in correlating stress results. The most promising result is the fact that the two standard nonlinear models can model tendon behaviour in the nonlinear plastic region. It is also noted that the two standard nonlinear models are physically insightful since their optimisation parameters can be eas-ily interpreted in terms of tendon elasticity and viscoelastic parameters.


Introduction
Mechanical properties of soft tissues are vital in the devevelopment of accurate and reliable computational models [1][2][3][4]. Tendons are normally understood to have strong tissues that are connected together that transfers muscle forces to the bones [5,6]. They have a complex hierarchical structure that spans several length scales from 300 nm long for tropocollagen to a few centimetres for the fascicles (tendon tissue). This structural organisation provides characteristics that enable proper functional properties. Tendons exhibit anisotropic, inhomogeneous and viscoelastic mechanical properties that are determined by their complicated hierarchical structure [7]. The natural structure of the tendons makes it possible for the tendons to generate large forces with small deformation. This natural phenomenon allows tendons to transmit required forces with less energy required for stretching. Because tendons are regarded as viscoelastic in nature, some energy that is transferred to cause deformation in the muscles is dissipated. However, they are most represented as elastic tissue that connect muscles to bones across joints. This makes them susceptible to quite a lot of tissue damage especially tendons that work around the glenohumeral (shoulder) joint which undertakes complex forms of locomotion in its day-to-day functions. Mechanically, tendons are like 'tie' structures, and they often act in tension, so the current study focuses on the tensile behaviour of a supraspinatus tendon.
The supraspinatus tendon forms a group of shoulder tendons that are most susceptible to injury called the rotator cuff tendons. These tendons insert into the head of the humerus, helping to achieve dynamic stability during abduction, adduction, and rotation of the arm (see Figure 1). However, the supraspinatus tendon itself ensures stability of the joint during arm abduction which involves it to act in tension. In order to understand its pathology and healing process, it is important that an accurate model of its tensile behaviour is formulated [6] . A recent study [8] shows that a correct modelling of tendon mechanical behaviour can significantly enhance the understanding of the mechanisms involved in its adaptability and mal-adaptability in response to mechanical loading, and to inform the design of patient-specific rehabilitation exercises addressing tendon pathologies. Generally, the problem of mechanical behaviour of soft tissues has been widely studied [9]. Researchers have approached the problem of modelling the mechanical behaviour of tendon tissue in two main ways: mathematical and computational approaches. Computational approaches are further classified as tissue-scale models or multiscale models. In most applications, these models are incorporated into finite element (FE) programs as solvers. Tissue-level computational models consider a tendon as an assemblage of elements with properties that are averaged over the tissue scale. [10] highlights the significance of these models by stating that contraction of fibroblasts at lower length scale is a sufficient mechanical impulse to building a planar wavy pattern and that the value of the crimp wavelength is determined by the mechanical properties of collagen fibrils and the interfibrillar matrix. Multiscale computational models disaggregate the tendon tissue into its different constituent length scales and solve its overall tissue response by translating responses over different length scales. This is obviously more involved than the previous approach and may require a lot of detailed material characterisation from nanoscale to macroscale. The rationale is that the mechanical response of a biological material to applied force reflects deformation mechanisms occurring within a hierarchical architecture extending over several distinct length scales [11][12][13][14][15][16]. Characterising, and in turn predicting the behaviour of such a material, requires an understanding of the mechanical properties of the substructures within the hierarchy, the interaction between the substructures, and the relative influence of each substructure on the overall behaviour. Computational models have been previously used successfully to predict various mechanisms of the diseases [17][18][19][20][21][22].
On the other hand, mathematical approaches are classified either as phenomenological or microstructural models [7]. Microstructural mathematical approaches account for material structure and attempt to model the tendon by generalising the mechanical behaviour of the different components in tendons. Research [23,24] shows that availability of such models may enhance the development of novel nano medicine and procedures for treatment of a broad class of collagen diseases, as well as lead to the design of novel biomaterials for regenerative medicine. Besides, these models may also lead to the development of potent therapies for tissue injury since a decline in mechanical response is a function of injury and age [25,26].
Phenomenological modelling is typically grouped into hyperelastic methods (those that do not display time-dependency yet exhibit a nonlinear response) and viscoelastic methods (those that incorporate time-dependence effects) [27]. Hyperelastic models are based on the minimisation of strain energy function with respect to strain invariants that are themselves expressed as functions of principal stretches and non-physical parameters. Their solutions involve determining the values of these non-physical parameters that minimise the sum of square errors between measured and model calculated stresses. Though they are typically simple numerical procedures, most of them require quite a lot of these non-physical parameters to achieve acceptable levels of correlations and fitting accuracies [28,29] . Some of the well-documented challenges in the development of these models have been the formulation of physically insightful models, incorporation of anisotropy, nonhomogeneity, geometrical multidimensionality and compressibility.
Viscoelastic models seem to possess better prospects for formulation of general theoretical frameworks, which could account for effects of nonlinearity, multiple-dimensionality, and finite viscoelasticity. These models provide a natural extension to modelling of anisotropy, complex deformations and more accurate prediction of mechanical behaviour under complex loading conditions [30][31][32] developed a general continuum model for the nonlinear viscoelastic behaviour of soft biological tissue, which was validated with test data from uniaxial extension of younger and older human patellar tendons and canine medial collateral ligaments. The model accurately predicted cyclic stresses for deformations too large to be adequately modelled by linear models. [33] demonstrated that nonlinear viscoelastic behaviour of the tendon and non-uniformity of the tendon crosssectional area should be taken into account for accurate modelling of tendons. [34] used a maximum dissipation function to formulate nonlinear viscoelastic evolution differential equations and successfully determined non-equilibrium transient responses (creep, stress relaxation and variable loading) in biological soft tissue. Leeman and Jones [31] particularly recommend the Maxwell type viscoelastic models to provide a good starting point for modelling of soft tissues. However, Shim and Mohr [30] find that viscoelastic constitutive models of the Zener type are not suitable for the modelling of rate dependent behaviour of polyurea. Their model predicts monotonic loading, material relaxation and viscous dissipation during loading and unloading cycle very well but it overestimates stiffness during unloading. Realising the power of a dissipation potential, [35] formulated a visco-hyperelastic constitutive model of short-term memory response for soft materials to capture both linear and nonlinear large deformation behaviour over a wide range of strain rates. The model displays great accuracy and avoids possible thermodynamic instabilities normally found in quasi-static hyperelastic models that corrupt dynamic stresses.
Thus, viscoelastic models seem to hold quite a lot of promise within the mathematical approaches for accurate modelling of soft tissue behaviour. It is, therefore, the main focus of the present study, and in order to bring out its advantages, in this paper, the performances of the two models in the elastic region have been shown in the light of a well-known Yeoh model. The choice of the hyperelastic Yeoh model was based on its performance and its simplicity in terms of the number of required optimisation parameters relative to the other models of comparable performance as reported in a study by Martins et al. [29].

Materials and Methods
The research methodology in this study involved two main steps: formulation of the Yeoh model and the two standard nonlinear models and Sample preparation, measurement and preparation of the secondary test data. The following two sections present these steps.

Theoretical formulation of the models
In this section, the theoretical formulations of the Yeoh and standard nonlinear models are presented. The two viscoelastic models have different configurations derived from the Maxwell and Kelvin-Voigt architecture.

Yeoh material model
Hyperelasticity implies that the structure can elongate under load >0.5 of its original length: a behaviour which can only be correctly estimated by nonlinear methods that assume finite elongations rather than the infinitesimal elongations as assumed by linear mathematical methods. Hooke's law does not satisfy the mechanical behaviour of soft tissue as tendons. Therefore, Yeoh's material model [36] is chosen in this study because of its reported superior performance in a previous study [29], and because it requires only three optimisation parameters. This model assumes material isotropy and incompressibility.
For uniaxial loading, such as during uniaxial tensile testing, a strain energy function can be given as where I1 is a strain invariant in the direction of the load which can be further expressed as a function of the principal stretch in that direction, λ1 as Yeoh (1993) derived the strain energy function expression for hyperelastic behaviour as where ci are the parameters to be estimated through optimisation. Although these parameters may have units of stress, they cannot be constrained to fall within values of measured stress because as is shown in Eq (4), they are parameters within a rather implicit function that also contains higher orders of strain invariants I1 and principal stretches. The principal Cauchy stress was then determined as [29,36] where the strain invariant = − .

Standard nonlinear solid (SNS) material models
Although some studies have modelled tendons as hyperelastic materials, others have shown that the extent of potential damage to their tissue depends on both the rate of change of their length and on applied force [34]. This behaviour causes them to exhibit both elastic and viscous characteristics when undergoing deformation. This observation presents a motivation for investigation into viscoelastic modelling of tendons. The formulation employed in this study can be thought of as 'pseudo-viscoelastic' since it does not necessarily render itself to dynamic loading, rather it merely utilises the topology for a static application.
There are several viscoelastic model structures presented in the literature, however, this study employs two variants of the standard nonlinear solid (SNS) model owing to its simplicity and its ability to model both creep and stress relaxation effects considering nonlinearity effects. In this study, these two effects are exploited for modelling on plasticity in the initial phases of material yielding. These models are: (1) one spring parallel with a series connection of another spring and a viscous damper (a modified Maxwell's model); and (2) a parallel connection of a spring and damper in series with another spring (a modified Kelvin-Voigt model). The latter combination was used by [37] in predicting nonequilibrium responses in a lumbar intervertebral disk and was found to yield very good results especially under very low frequency cyclic loading. The study particularly focussed on approximation of stress relaxation up to a time period of 15,000 s.
In this study, both models are adopted to model tendon behaviour during linear and plastic failure at a uniform strain rate. Figure 2 shows the architectures of the SNS models under study.  (5) and (6), respectively [38].
where σ is the stress, ε is the engineering strain, Ei, i=1,2 denote the two elasticity parameters represented as moduli of elasticity, η is the damping characteristic, and the dot over the variables represents the time rate of those variables. Both the moduli of elasticities and the damping parameter are assumed to be nonlinear functions of strain. Eqs (5) and (6) show that under constant stress, the models will deform to some instantaneous strain but under elongated loading, the material will continue to deform and asymptotically approach some steady-state strain, which is the retarded elastic portion of strain. It is hypothesised that strain does not really remain constant during the yielding phase. Thus the material will only obey Hooke's law during the elastic phase.
Due to the complexity in structural hierarchy and randomness in helical organisation of tendon tissue, its stress-strain curve from the toe region to the point of the ultimate tensile stress is typically governed by an s-shaped process. This may imply that the material properties within the tendon tissue do not necessarily remain constant but rather that they continually vary during the testing process. For example, stiffness in the toe region is dominated by the straightening of fibril helical crimp: in the linear region it is dominated by stretching of the fibrils, and in the nonlinear plastic region, it is dominated by a combination of stretching and fracture of fibrils. Therefore, in order to account for these changes in material properties during testing, it has been assumed that both Eqs (5) and (6) are governed by the s-shaped curve for soft tissue stress-strain behaviour given by a function where , 1,2,3, b and k0 are the optimisation parameters. Then we use these optimisation parameters to derive expressions for the material properties. Upon basic mathematical calculations, the elasticity and damping parameters in Eqs (5) and (6) For both the Yeoh and SNS models, the parameters are solved using a nonlinear least squares curve fitting routine in MATLAB employing a default 'trustregion-reflective' algorithm. This algorithm is a subspace trust region algorithm for large-scale problems, which is based on the interior-reflective Newton method (Mathworks® Inc., 2019). Each iteration involves the approximate solution of a large linear system using the method of preconditioned conjugate gradients. However, this algorithm does not accept underdetermined systems, so in those cases, the Levenberg-Marquardt algorithm is recommended. It has the advantage of accepting bound constraints which are not permitted by the other algorithm. One can also use this algorithm in small-to medium-scale problems without computing the Jacobian or providing the Jacobian sparsity pattern, which makes it much faster than the routines that compute the Jacobian.

Sample preparation, measurement, and test data preparation
The three samples were harvested from frozen shoulders of cadavers of undisclosed age and features. The authors were merely granted consent for testing and reporting the results of their measurements. The specimens were of different cross-sectional areas: 90 mm2, 70 mm2 and 50 mm2 with lengths of 50 ±2 mm, which translates into three different aspect ratios. The specimens were temporarily preserved in saline solution before testing under room air temperature and relative humidity between 30 60 %. During testing the specimen were kept at 35-37oC to mimic body temperature.
They were tested one after another using an MTS EM Tensile Tester up to the point of total tissue rupture. The machine was operated in displacement control at constant strain rates on all the three tests. Strain rates were kept constant throughout the testing procedure on all tests at around 0.0173 per s with variations in the strain rates at different sample points being 3.5% for test 1, 3% for test 2, and 5.2% for test 3. The variations in strain rates from one test to another as well as within the same test were inevitable due to challenges in achieving perfect tissue clamping, measurement precision errors, and changes in material properties through the different stress-strain phases. The data from the toe region to the maximum stress value corresponding to ultimate tensile stress were selected for all tests. For the Yeh model, the test data range was only up to the upper limit of the linear elastic region. A simple data cleaning procedure by the use of a 4-point moving average filter was applied to all the test data.

Results and Discussions
In this section, a summary of the test results and a number of distinctive features are presented. The performances of the three models in terms of their correlation and fitting accuracies, as well as computational efficiencies are discussed. The section concludes with a discussion of the evolution of material properties over the strain process for the two SNS models.

Summary of test results
Three uniaxial tensile tests were conducted on tendons having 90 mm2, 70 mm2 and 50 mm2 in cross-sectional area. Table 1 summarises the test results. The following test parameters are extracted from the test results and presented in the table: moduli of elasticity in the linear region, ultimate tensile stress (UTS), and reserve tensile strength after UTS. Tendons in Tests 2 and 3 have larger moduli of elasticity and exit the toe region at much lower strain level than the tendon in Test 1. The tendon in Test 2 expectedly reaches its UTS at much lower value of 2.621 N/mm2 however, the tendon in Test 3 reaches its UTS at almost equal UTS as in Test 1. But after reaching the UTS, both Tests 2 and 3 elongate by more than 50% of their original lengths, while Test 1 only elongates by about 28% of its original length before rupture. Apart from the large differences in UTS, Tests 2 and 3 bear some marked similarities in almost every aspect. The similarities between Tests 2 and 3 are also evident in the graphs of material property trends plotted in Figures 7 and 8 to be discussed later. Both tendons exited the toe region at strains of around 0.1 and at equal stresses; they had the same moduli of elasticity in the linear region; and they also had similar amounts of reserve strength after UTS. Furthermore, they displayed similar stress-strain behaviour after UTS. After UTS, they both exhibited two consecutive stages of 'stress relaxation' type of behaviour after which they transited into a strain softening phase before they finally ruptured. The tendon in Test 1 exited the toe region at a much higher strain and stress of 0.3 and 1.25 N/mm2, respectively, exhibited the lowest value of modulus of elasticity in the linear region, and reached its UTS at 4.393 N/mm2. Its behaviour after UTS showed two rather sharp consecutive tissue 'stress relaxation' stages before final rupture. In contrast to Tests 2 and 3, the tendon in Test 1 did not exhibit a distinct strain softening type behaviour before failure. Besides, its 'stress relaxation' stages were much steeper than in Tests 2 and 3. These differences can be attributed to inherent tissue material property differences. From the calculated differences in the moduli of elasticity, and behaviour after UTS, it is clear that the tendon in Test 1 was less elastic and lacked enough toughness to sustain it for long after reaching its UTS. This might have been caused by differences in tissue treatment or operating conditions prior to or during testing.

Model performances
Figures 3 -5 show correlation plots (in frames (a)) and error plots (in frames (b)) of models' outputs for all three tests. Included on frames (a) on all these figures, are the calculated correlation and fitting error percentages. The two versions of the SNS model perform similarly on each of the three tests. They range from correlation percentages of 99.8% to 100% having fitting errors below 5% on all tests. The two models perform comparably to the Yeoh model in the linear elastic range on all tests. The other advantage is that they are able to model tissue behaviour up to the nonlinear plastic region up to the UTS.
The error plots in frames (b) show that on average the SNS models perform better than the Yeoh model in the toe region. The toe region errors are clearly the largest contributor to the fitting errors for the Yeoh model. Martins et al., 2006 have attributed these errors to insufficient capacity for modelling material non-homogeneity. We also hypothesize that measurement errors in this region play a role. The error plots on all three tests consistently display spikes in errors at region boundaries, first between toe region and linear region, second between linear region and nonlinear plastic region, and finally at the UTS. Both the SNS models and the Yeoh model fail to model these discontinuities and have further contributed to the fitting errors. For the SNS models, the errors at the UTS are clearly the largest contributor to their fitting errors. At the UTS the material apparently behaves in an unstable manner and the model is incapable of fitting that behaviour. As already noted, the material transits into a stress relaxation type of behaviour after this point so that the test data for the SNS models were limited to the UTS on all the three tests. Another important observation is that the performances of these models were influenced by the ranges of the selected data. This is normal characteristic of regression type models, in which their accuracies are influenced by how well the test data represents the underlying theory. For example, the Yeoh model cannot yield good correlations in the plastic region.  Table 2 summarizes the optimisation results in terms of the optimisation parameters, number of iterations and function evaluations for the SNS models and Yeoh model for all three tests. For both models the optimisation parameters are non-physical with the only exception that those for the SNS models are parameters of tissue properties such as moduli of elasticity and tissue damping. For the optimisation process, Yeoh model clearly performs much better than the SNS models with fewer parameters, number of iterations and function evaluations. The Yeoh model also converges much faster at 0.07 s as compared to 1 s on average for the SNS models using an Intel(R) Core(TM) i7-3720QM CPU @ 2.60GHz with 4 cores and 8 GB (RAM) running on 64-bit Microsoft Windows 10 Enterprise.

Material property behaviour
The optimisation parameters for the SNS models are applied to Eqs (8), (9) and (10) to calculate the moduli of elasticity E1 and E2, and the tissue damping factor η for the Maxwell version and Eqs (11), (12) and (13)    (a) A closer look at the evolution of the material properties E1, E2 and η plotted in Figures 6 -8 reveals more detailed information about the tendon behaviour in each test. Despite a number of aesthetic similarities especially between Tests 2 and 3, it is important to note that the differences in the actual values of the properties imply that these tendons were not exactly the same. Mechanically, the properties serve different physiological functions such as material compliance, force isolation, and energy dissipation. The two SNS models distribute these functions differently between its two elastic components (E1 and E2) depending on their different architectures. In the following sections, we attempt to explain these functions in the light of the property trends as plotted in Figures 6 -8.

Material compliance
Material compliance is represented by E1 and E2 in Maxwell and Kelvin-Voigt models respectively. In the Kelvin-Voigt model, the trend of E2 on each strain level is the same for each test as shown in Figures 6(e), 7(e) and 8(e). The tissues have low elasticities in the toe region, then increase toward a maximum at the mid-point of the linear region and finally decrease towards the UTS. This behaviour is rather expected from a practical point of view since the tissue is expected to be more compliant in the toe region and towards its UTS but exhibit higher stiffness in the middle of its linear region. For the Maxwell model, material compliance is represented by E1 and it is shown in Figures 6(a), 7(a) and 8(a). However, the previously discussed tissue behaviour in Kelvin-Voigt model is not represented by the Maxwell model in these figures except the characteristic decrease in elasticity towards the UTS and the maximum at the mid-point of the linear region in Test 1. From these results, it may be concluded that the Kelvin-Voigt model is more stable in modelling the tendon elasticity than the Maxwell model.

Force isolation
This property measures how the tendon transmits external forces to the adjoining parts. It is dominated by negative elasticity in the models, represented by E1 for Kelvin-Voigt model and E2 for Maxwell model. In the Maxwell model, Figures 6(b), 7(b) and 8(b) typically show similar tissue behaviour in all tests with the only difference being in their exact values of elasticity. However, the Kelvin-Voigt model show similarities only in Tests 2 and 3 (Figures 7(d) and 8(d)). Figure 6(d) shows effective force isolation only in the toe region, which is actually tendon's physiological range. Therefore, it may be concluded that the Maxwell model is more stable in modelling force isolation than the Kelvin-Voigt model.

Damping or energy dissipation
Typically, this property is represented by η in Figures 6(f), 7(f) and 8(f) and gives an indication of tendon's energy dissipation capacity. A desirable damping property in the tendon tissue should therefore be either very low or negative values especially in the physiological state though this may cause material behaviour instabilities under certain loading conditions. For both Kelvin-Voigt and Maxwell models the tissues exhibit similar behaviour for Tests 2 and 3. Test 1 exhibits higher levels of damping in both models in the linear and plastic regions. This further confirms the similarities between Tests 2 and 3.
However, it must be mentioned that there are quite numerous areas around tendon material property evolution during tensile testing that require further study. Some of these issues are the exact mechanisms within the tendon that cause the complicated nature of its variations over the strain process and the exact evolution of the material properties in the toe region.

Discussion
Authors should discuss the results and how they can be interpreted from the perspective of previous studies and of the working hypotheses. The findings and their implications should be discussed in the broadest context possible. Future research directions may also be highlighted.

Conclusions
In this paper two models whose structures are derived from viscoelastic topology are evaluated in the light of a well-known Yeoh model. The rationale behind this is to exploit viscoelasticity to model unstable yield behaviour only up to the ultimate tensile stress, and also to elucidate the material property evolution during the strain process. It has been observed that the two viscoelastic models are able to accurately fit the tendon stress-strain behaviour in the linear region quite comparable to the Yeoh model. It has further been shown that the two viscoelastic models are capable of predicting the stress-strain beyond the linear elastic limit up to the ultimate tensile stress point.
On understanding the material property behaviour, the two viscoelastic models promise to be physically insightful. It has been shown that the Kelvin-Voigt architecture is more stable in modelling tendon stiffness while the Maxwell architecture is more stable in modelling the tendon's force isolation property. Although it is hard at this stage to cascade the obtained information into microscale tendon behaviour, it can be satisfactorily employed to understand several influences of the ground substance that influence tendon's dissipative characteristics on the macroscale.
There are still few more areas that require further study. For example, the variation of material with tissue treatment and different testing conditions; more stable ways of modelling damping in tendon tissue; and designing of reliable testing procedures of tendons in the toe region. Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.
Data Availability Statement: In this section, please provide details regarding where data supporting reported results can be found, including links to publicly archived datasets analyzed or generated during the study. Please refer to suggested Data Availability Statements in section "MDPI Research Data Policies" at https://www.mdpi.com/ethics. You might choose to exclude this statement if the study did not report any data.