Journal of Biomechanics 47 (2014) 3598–3604

Contents lists available at ScienceDirect

Journal of Biomechanics journal homepage: www.elsevier.com/locate/jbiomech www.JBiomech.com

Subject-specific finite element analysis to characterize the influence of geometry and material properties in Achilles tendon rupture Vickie B. Shim a,b,n, Justin W. Fernandez a,c, Prasad B. Gamage a, Camille Regnery a,f, David W. Smith d, Bruce S. Gardiner d, David G. Lloyd b,e, Thor F. Besier a,c a

Auckland Bioengineering Institute, University of Auckland, Auckland, New Zealand Centre for Musculoskeletal Research, Griffith Health Institute, Griffith University, Gold Coast, QLD, Australia c Department of Engineering Science, University of Auckland, Auckland, New Zealand d School of Computer Science and Software Engineering, University of Western Australia, Nedlands, Australia e School of Sport Science, Exercise and Health, University of Western Australia, Nedlands, Australia f ENSEIRB, MATMECH Bordeaux, France b

art ic l e i nf o

a b s t r a c t

Article history: Accepted 6 October 2014

Achilles tendon injuries including rupture are one of the most frequent musculoskeletal injuries, but the mechanisms for these injuries are still not fully understood. Previous in vivo and experimental studies suggest that tendon rupture mainly occurs in the tendon mid-section and predominantly more in men than women due to reasons yet to be identified. Therefore we aimed to investigate possible mechanisms for tendon rupture using finite element (FE) analysis. Specifically, we have developed a framework for generating subject-specific FE models of human Achilles tendon. A total of ten 3D FE models of human Achilles tendon were generated. Subject-specific geometries were obtained using ultrasound images and a mesh morphing technique called Free Form Deformation. Tendon material properties were obtained by performing material optimization that compared and minimized difference in uniaxial tension experimental results with model predictions. Our results showed that both tendon geometry and material properties are highly subject-specific. This subject-specificity was also evident in our rupture predictions as the locations and loads of tendon ruptures were different in all specimens tested. A parametric study was performed to characterize the influence of geometries and material properties on tendon rupture. Our results showed that tendon rupture locations were dependent largely on geometry while rupture loads were more influenced by tendon material properties. Future work will investigate the role of microstructural properties of the tissue on tendon rupture and degeneration by using advanced material descriptions. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Achilles tendon Finite element analysis Tendon rupture Material property optimization Subject-specific finite element analysis

1. Introduction Achilles tendinopathy, including tendon rupture, occur at a rate of about 250,000 per year in the US alone (Jarvinen et al., 2005; Pennisi, 2002). The mechanism of tendinopathy and rupture is complex and thought to be influenced by tendon geometry, material-strength, sex, disease and genetics. Achilles tendon ruptures are typically reported to occur at 2–6 cm above the insertion to the calcaneus bone, in a region that is hypovascular (Theobald et al., 2005). It is not understood why this region receives poor blood supply and is prone to rupture. Tendon rupture predominantly occurs in males with the reported ratios between men and women ranging between 2:1 (Lemm et al., 1992) n Corresponding author at: Auckland Bioengineering Institute, 70 Symonds Street, Auckland, New Zealand. Tel.: þ 64 9 3737599x86932. E-mail address: [email protected] (V.B. Shim).

http://dx.doi.org/10.1016/j.jbiomech.2014.10.001 0021-9290/& 2014 Elsevier Ltd. All rights reserved.

and 12:1 (Bandak et al., 2001). Finite element (FE) analysis of the human Achilles tendon may provide insight into the possible mechanism of tendinopathy and tendon rupture as it can provide a standardized framework for systematic parametric and integrative analysis. Previous computational models of the human Achilles tendon include 1-dimensional (1D) line elements to describe the tendon as a part of a foot model (Bandak et al., 2001; Bayod et al., 2012). Although such models are useful in predicting joint kinematics and kinetics, they cannot predict spatially varying 3-dimensional (3D) stresses and strains. These variables are most likely important in soft tissues like tendon that display large deformation and may have heterogeneous stress distribution patterns leading to regions of localized stress concentrations. To assess the possible reasons for the location of tendinopathy and rupture, and the aforementioned gender differences, 3D FE analyses may need to include subjectspecific geometry and material properties.

V.B. Shim et al. / Journal of Biomechanics 47 (2014) 3598–3604

Previous studies have investigated the relationships between the Achilles tendon geometry (e.g. length and cross-sectional area) and mechanical or structural properties of the tendon (Hess, 2010; Houghton et al., 2013; Jozsa et al., 1989; Kongsgaard et al., 2005). However, those studies were either in-vivo human or animal studies and could not comparatively analyze the influence of geometrical and/or material properties, and their inter-relationships, on tendon rupture. The aim of this study was to create subject-specific 3D FE models of the human Achilles tendon to evaluate how variation in geometry and material properties influence tendon mechanics. This included parameter optimization for a transversely isotropic hyperelastic material, which have been previously used to represent the material properties of different knee ligaments (Gardiner and Weiss, 2003; Quapp and Weiss, 1998; Weiss and Gardiner, 2001). Using these models, we simulated tendon rupture under uniaxial loading to investigate (i) the influence of the subject-specific geometry; and (ii) material properties on tendon 3D rupture patterns. We hypothesized that geometry and material properties of the tendon would influence the tendon's (i) rupture location, and (ii) rupture load.

2. Methods 2.1. Experimental data Data from a previous experimental study were used (Wren et al., 2003). Briefly, in that study fresh-frozen human Achilles tendons from donors were imaged with high-frequency ultrasound. Specifically, the cross-sectional area of each tendon was measured at 2, 4 and 6 cm proximal to the tendon's insertion, which corresponds to the mid-section of the tendon. The ends of each specimen were gripped so that the grip-to-grip length was approximately 10 cm. Two types of mechanical testing (cyclic and creep testing) performed using a MTS machine (MTS, Eden Parairie, MN). Each tendon specimen had six markers along the centerline of the tendon and a CCD camera recorded movements of these markers during mechanical testing to enable material properties to be determined. Ten tendons (eight female and two male, average age 68 yr) from 18 specimens that underwent cyclic testing were selected for subject-specific FE modeling. 2.2. Generation of subject-specific FE model Three cross-sectional ultrasound images of the ten selected tendons were segmented to describe subject-specific geometry. To create each FE tendon model, we used the computational framework developed as a part of the International Union of Physiological Society (IUPS) Physiome project (Fernandez et al., 2012; Hunter et al., 2002; Shim et al., 2011). A generic FE mesh of the human Achilles tendon was first developed using the Visible Human dataset. We used high order cubic Hermite elements, which preserve both the continuity of nodal values and their first derivatives. This type of element can describe smooth surfaces often found in biological structures such as the heart (Stevens et al., 2003), muscle (Oberhofer et al., 2010) and bone (Munro et al., 2013; Shim et al., 2012) with a lot fewer elements as each node has total 24 degrees of freedom. Therefore, our model has 32 elements and 72 nodes, which corresponds to 1728 degrees of freedom. This mesh density was chosen after convergence analysis. Free Form Deformation (FFD) (Fernandez et al., 2004) was then used to morph the generic tendon mesh to the segmented ultrasound images (Fig. 1). This involved morphing a host mesh to minimize the distance between identified landmark and target points. The embedded tendon mesh was then morphed to a subject-specific shape Fig. 2.

3599

The tendon was represented as an incompressible, transversely isotropic hyperelastic material. It had a strain energy density function (Gardiner and Weiss, 2003) that modeled the tendon as a composite of ground substance matrix with embedded collagen fibers, i.e. W ¼ F 1 ðI 1 Þ þ F 2 ðλÞ

ð1Þ

where I1 is the first invariant of the right Cauchy stretch tensor and λ is the stretch ratio along the local fibre direction. The function F1 describes the behavior of the ground substance matrix and F2 represents the behavior of collagen fibers. The ground substance was described with the neo-Hookian material model, i.e. F1 ¼

C1 ðI 1  3Þ 2

ð2Þ

The strain energy of the collagen fibers was represented as a piecewise function that characterized their non-linear stress/strain behavior using the following (Weiss et al., 1996): ∂F 2 ¼ 0; ∂λ

λ

for λ r 1;

h i ∂F 2 ¼ C 3 eC 4 ðλ  1Þ  1 ; ∂λ

λ

∂F 2 ¼ C 5 λ þ C 6; ∂λ

λ

for 1 r λ r λ

for λ Z λ

n

n

ð3Þ

where λ* is the stretch value where the collagen fibers become uncrimped, while coefficients C3, scales the exponential stress, C4 represents the rate of collagen fibre loading, and C5 modulus of the straightened collagen. To ensure a smooth transition between the 2nd and 3rd piecewise functions C6 was determined using the following:      n n C 6 ¼ C 3 exp C 4 λ  1  1  C 5 λ ð4Þ

2.3. Estimating the material properties The values for the material coefficients were estimated using the cyclic experimental data. First, the value for C5, the modulus of the straightened collagen, was obtained by taking the gradient of the linear region in the experimental stress/strain curve of the whole tendon. λ* was kept constant at 1.055, the average value found from a previous study (Quapp and Weiss, 1998). The remaining material coefficients (C1, C3 and C4) were estimated by performing material property optimization using a FE analysis to simulate the experiments. Boundary conditions simulated the experimental set up where the bottom nodes were fixed to represent the tendon clamp, while a force was uniaxially and equally applied to the top nodes to simulate the movement of the upper clamp during the experiments. The model was solved using finite elasticity mechanics in our open-source bioengineering software CMISS (www.cmiss.org) from which the predicted positions of the six markers from the experiment were computed. The root mean square (RMS) error between predicted and actual experimental marker positions was then computed. The optimal material parameters were estimated by finding the minimum RMS marker position error via the non-linear leastsquares lsqnonlin algorithm in Matlab's Optimization Toolbox (Matlab Version 2009b, The MathWorks, USA). The lsqnonlin algorithm used has been previously detailed (Babarenda Gamage et al., 2011). For each specimen, 5 different randomly generated initial parameter guesses were tested to ensure global minimums were found, which was the case.

3600

V.B. Shim et al. / Journal of Biomechanics 47 (2014) 3598–3604

Fig. 1. Schematic of the generation of the subject-specific FE model of the Achilles tendon using free form deformation.

Fig. 2. Comparison of subject-specific meshes (brown) compared to the generic mesh (white). The green and yellow boxes are host meshes used in mesh morphing, where the yellow boxes are the initial host meshes and green boxes are the deformed host meshes. Subjects 1–4 are the first 4 listed in Table 1. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

2.4. Predicting stress patterns and tendon rupture using virtual experiments Once subject-specific material properties were obtained, stress patterns and rupture of the ten subject-specific Achilles tendons models was simulated. First, 5% uniaxial stretch was applied to the tendon models with a displacement boundary condition. The Von Mises stress distributions were recorded for each of the ten Achilles tendons simulations. The tendon rupture in the simulations was the performed based on previous in vitro tests that showed the failure stress of tendon is close to 100 MPa (Butler et al., 1984; Kongsgaard et al., 2005; Wren et al., 2001b). Therefore, we employed the same boundary conditions as our material property estimations, but instead increased the magnitude of the applied force in steps of 100 N until the maximum von Mises stress was greater than 100 MPa. Specifically, tendon rupture was defined to be in the region of the tendon where 15 or more consecutive Gauss points had von Mises stress greater than 100 MPa (Shim et al., 2010). The 15 Gauss points roughly corresponded to the length of 3 mm in

our model, which was considered large enough to lead to a rupture. Due to the quasi-static simulation, the influence of loading rate was not considered. The rupture load and location were recorded. The influence of the material properties and tendon geometry on the rupture load and location was then examined using three different test cases. In Case 1, rupture was estimated using the subject-specific tendon geometries and subject-specific properties obtained from the optimization. In Case 2 subject-specific tendon geometries with the average material properties from all ten subjects were used. Finally, in Case 3 the subject-specific material properties were used with the generic geometry from the Visible Human data. Rupture loads and locations were recorded during these virtual experiments. The Euclidean distances between the centroids of predicted rupture locations from the three cases were calculated; i.e. Case 1 versus Case 2; Case 1 versus Case 3, and Case 2 versus Case 3. The Euclidean distances and rupture loads were compared using one-way ANOVA to test for statistically significant differences (po 0.05) between the predicted rupture locations and loads from subject-specific and average models.

V.B. Shim et al. / Journal of Biomechanics 47 (2014) 3598–3604

3. Results Morphing of the generic mesh with cross-sectional ultrasound images resulted in subject-specific geometries with an average RMS error of 0.78 70.05 mm between the morphed mesh and the ultrasound image at the mid-section of the tendon (Fig. 1). Following the material property optimization on the ten specimens the average RMS error between the predicted and actual markers was less than 0.1 mm (0.092 70.008 mm). The final material properties coefficients and specimen details are shown in Table 1. Using the optimized material coefficients and uniaxial whole tendon stretch of 5%, higher stresses were observed in the midsection of all the tendons (Fig. 3, top row). However, the distribution of stress varied across subjects (Fig. 3 top row). The stress

Table 1 Optimized material parameters from the ten subjects. Subject Sex Age

Weight (kg)

Height (m)

C1

C3

C4

1 2 3 4 5 6 7 8 9 10

63.64 54.55 113.64 78.18 59.87 79.55 60.91 52.27 63.64 49.89

1.80 1.62 1.83 1.83 1.80 1.77 1.58 1.65 1.92 1.71

46.52 43.13 157.98 164.18 38.07 115.99 40.20 35.73 25.38 33.12

15.29 2.21 65.31 31.73 11.57 44.62 4.98 6.22 7.90 4.12

26.80 928.15 97.37 1176.82 6.50 543.57 2.95 1416.44 30.65 1153.65 4.99 820.04 83.34 992.08 73.12 691.17 14.59 935.83 70.13 855.16

67.61 18.94

1.75 0.11

Mean SD

F F M M F F F F F F

66 82 69 71 60 50 81 91 62 81 71.3 12.44

70.03 19.40 41.04 54.19 21.13 36.20

C5

3601

pattern in the mid-section varied from being elongated and broad to short and narrow, for example, as shown in the tendons from subjects 3 and 4 in Fig. 3 (top row). Tendon rupture occurred in the mid-section, although the exact locations, as well as the rupture loads, were different from subject to subject (Fig. 3 bottom row). There were two major types of rupture patterns. Type 1 rupture appeared in the center of the midsection of the tendon while Type 2 predicted tendon rupture at one or both sides of the mid-section. Out of 10 models that we developed, 4 displayed Type 1 rupture (all female) while the other 6 showed Type 2 rupture (2 male and 4 female). Tendon rupture locations were more dependent on tendon geometry than tendon material properties (Fig. 4). When the same subject-specific geometries with different material properties were used, the predicted rupture locations did not change, as can be seen by comparing Case 1 and 2 (Fig. 4). But when different geometries were used with the same material properties, the rupture location changed to the center of the tendon mid-section (Fig. 4, Cases 1 and 3). These results were confirmed with the Euclidean distances between rupture locations (Fig. 5). The Euclidean distance between Case 1 and 2 locations was smaller than the

961.97 264.74

Fig. 3. Exemplar Von Mises stress distributions under 5% constant uniaxial tension (top row) and predicted rupture locations (bottom row). Subjects 1–4 are the first 4 listed in Table 1.

Fig. 4. Influence of geometry and material properties on tendon rupture locations. Subjects 1–4 are the first 4 listed in Table 1.

3602

V.B. Shim et al. / Journal of Biomechanics 47 (2014) 3598–3604

distances between Case 1 and 3 and Case 2 and 3 locations (p ¼0.013 and p¼0.018 respectively). This indicates that, for Cases 1 and 2 where only material properties were changed, the predicted rupture locations were very similar to each other. Rupture locations differed once geometry changed from subjectspecific to generic (i.e. Case 1 versus 3, and Case 2 versus 3). There was no statistically significant difference when geometry was not changed (p¼0.723). In contrast to predicted rupture location, predicted tendon rupture loads were influenced by material properties more than tendon geometry (Fig. 6). Changes in material properties and geometry both brought changes in predicted rupture loads. However, when material properties changed from subject-specific to average, the change in the rupture load was generally greater than when geometry changed from subject-specific to generic. The difference in rupture load was significantly different only when material properties changed (p ¼ 0.026). The relationship between geometry and material properties was further analyzed by comparing the predicted rupture loads in the two major rupture types, where Type 1 ruptured at the center of the mid-section while Type 2 ruptured at one or both sides of the mid-section. When rupture loads between the two types were compared, we found that Type 1 had a lower rupture load (2125 N)

Fig. 7. Representative ultrasound images of Type 1 and Type 2 tendons. Type 1 has more round and smaller cross-sectional areas than Type 2.

than Type 2 (3440 N) (p¼ 0.02). The main geometrical difference between Types 1 and 2 was the shape of the cross-sectional area in the mid-section of the tendon. Type 1 had a cross-sectional area that was smaller and rounder than Type 2, which had an elliptical shape (Fig. 7). Therefore our results indicate that wide band shape tendon (Type 2) withstand loads better than thin tubular shape tendon (Type 1). The stiffness between these two types of tendon was also compared. Type 1 had a higher stiffness value (998.7 MPa) (C5 in Eq. (3)) than Type 2 (827.3 MPa). However these differences were not statistically significant (p¼ 0.13).

4. Discussion

Fig. 5. Euclidean distances between predicted rupture locations from the three cases tested.

Fig. 6. Influence of material properties and geometry on tendon rupture loads. Case 1 is where the both the geometries and material properties were subjectspecific. Case 2 where geometries were subject-specific and material properties were subjects' average values. Case 3 where geometry was generic and material properties subject-specific.

The current study demonstrates the value of 3D subjectspecific Achilles tendon model in investigating tendon mechanics. Stress analysis performed using 10 subject-specific FE models showed that tendon geometry, material properties and the stress distribution patterns were all subject-specific. Moreover, this subject variability was also shown in our rupture predictions as both predicted tendon rupture load and location were different from specimen to specimen. Our results indicate that material properties mainly influenced rupture load while tendon geometry determined rupture locations. Our model predicted that tendon rupture almost always occurs in the mid-section of the tissue under uniaxial tension. This agrees with other in-vivo or epidemiological studies that also reported the same region as being most prone to rupture (Hess, 2010; Jarvinen et al., 2005). A number of possible reasons for this have been discussed in the literature such as poor blood supply or focal stress concentration (Kongsgaard et al., 2005). Since our results showed that geometry is the main predictor for rupture location, it may indicate that thinning of the tissue in the mid-section is the major factor that makes the mid-section particularly susceptible for rupture. However, there was variation in predicted rupture locations within the mid-section. In fact, the predicted rupture patterns were categorized into two major types according to their shape – Type 1 tendon had smaller and rounder cross-sectional area and ruptured at the center of the mid-section; Type 2 tendon had thicker and elliptical cross-sectional area and ruptured on one or both sides of the mid-section. When predicted rupture loads for these two types were compared, Type 2 has a significantly higher rupture load than Type 1. Since the main difference between these two types was the cross-sectional area, our study showed that thicker tendons that have a larger cross-sectional area provide a greater safety margin. A number of in vivo human and animal studies measured tendon cross-sectional areas and compared between two groups of different activity types and levels. In humans, cross-sectional areas of long distance runners were markedly greater than the control (Hansen et al., 2003;

V.B. Shim et al. / Journal of Biomechanics 47 (2014) 3598–3604

Magnusson and Kjaer, 2003). A similar observation was made in an animal study that compared trained and untrained horses (Kasashima et al., 2002). These studies indicate an increase in cross-sectional area would reduce the average stress of the tendon, thereby reducing the risk of tendon rupture. Our study supports these findings. An interesting point can also be drawn from the material optimization. First, predicted tendon rupture load was always higher when subject-specific material properties were used than using average material properties. This may indicate that tendon material property adapts differently depending on their geometry. This point is further supported by comparing material properties for Types 1 and 2 tendon. Type 1 tendon, whose geometry may be more susceptible to rupture, had higher stiffness than Type 2. This means that thin tendons like Type 1 will require a larger force to produce a given strain, effectively providing a means to mitigate their intrinsic rupture susceptibility. Numerous in vivo studies have also shown that mechanical properties of human tendons undergo substantial changes mainly with aging and disuse (Maganaris et al., 2006; Reeves et al., 2005). Our study showed that the intrinsic tendon geometry might also influence the way mechanical adaptation of tendon occurs. It was not possible to determine why there is a significant difference in the male female ratio of Achilles tendon rupture as we only had two male specimens (out of 10). However, the two male specimens had a Type 2 tendon and had higher predicted rupture load. This may indicate that the higher incidence of tendon rupture in males may be activity type and level related than mechanical or geometrical properties. However a further study will be required to fully clarify this. Our study also demonstrated the value of 3D FE models in investigating tendon pathophysiology. A number of researchers have developed computational models of tendon and ligaments using tissue material properties derived from experimental studies. However, the majority of these used simple linear or nonlinear springs to represent the tendon (Lichtwark and Wilson, 2006, 2007). Such models are useful to study whole joint mechanics, but they cannot predict stress patterns within the tissue as well as its interactions with surrounding structures of the joint. Previous 3D FE models of musculoskeletal soft tissues have been mainly focused on the muscle and ligaments. Blemker and Delp (2005) developed a novel framework for representing complex muscle architecture and geometries with 3D models, which also included tendons as well. They found that 3D representations were important in representing complex path motion and fiber moment arms. Pioletti et al. (1998) developed a 3D FE model of the human anterior cruciate ligament (ACL) using the isotropic and elastic material property. Gardiner and Weiss (2003) developed a subject-specific FE model of the human medial collateral ligament using a similar approach to the present study. They also used the transversely isotropic hyperelastic material description and performed mechanical testing to obtain optimal material parameters. Their model was able to adequately describe 3D strain distribution during valgus knee loading. This provided the inspiration to develop the current Achilles tendon model. Computational models of tendon by itself, however, are rather scarce. Wakabayashi et al. (2003) developed a 2D FE model of the supraspinatus tendon but used only a simple linear elastic material properties. Cheung et al. (2006) developed a 3D FE model of the foot to investigate loading from the Achilles tendon on plantar fascia tension but tendon loading was approximated with a force boundary condition and did not have a separate model for the tendon. A number of microstructural models of tendon exist such as the work by Reese et al. (2010) and Lavagnino et al. (2008). Recently, Fiorentino et al. (2014) developed a 3D model of the biceps femoris long head and found that local fibre strains were

3603

non-uniform and the peak strain occurred at the region of the smallest cross-sectional area, which is in line with the findings from our study. However, to our best knowledge, a subject-specific 3D FE analysis of the whole human Achilles tendon has not been presented in the research literature. Having a 3D subject-specific models seems to very important in regard to tendon stress distribution patterns, and rupture load and location. There are several limitations in our study. First, our subjectspecific models were generated from ultrasound scans taken from three regions of the tendon mid-section. Therefore it is uncertain whether the other parts of the tendon, other than the mid-section, were customized accurately. However, we tested our method by customizing the generic model from the Visible Human Male scans to generate a model for the Visible Female only using scans from the mid-section selected from the full set of scans. We found that, although the model was customized using scans from the midsection, the overall geometry was close to the mesh generated from the full set of Visible Female scans (RMS error of 1.15 mm, Appendix 1). Therefore, we are confident that our FFD technique was capable of generating a suitable subject-specific FE model of the free Achilles tendon using these ultrasound data. The second limitation was that we only had experimental results from uniaxial tension experiments with free tendon. Therefore all simulations were performed under uniaxial loading conditions. Biaxial loading experiments that characterize tissue deformation transverse to the fiber direction would provide additional experimental data for material property estimation. Moreover we modeled the free tendon and did not have muscles or bones in the joint. Due to this limitation, we have not considered responses under torsion or other complicated loading cases. However, we already have a full foot model from our IUPS Physiome Anatomy database (www.physiomeproject.org). Therefore, future work includes extending the current uniaxial experimental set-up to biaxial one to capture the tissue deformation more accurately. It also needs to be mentioned that our rupture results are predictions from the FE model and are without absolute experimental validation. As such, they should be used and interpreted with caution. However, our FE stress predictions are based on experimental data and our predicted regions of high stress in the Achilles tendon are consistent with where tendinopathy and ruptures have been reported to occur (Kongsgaard et al., 2005; Theobald et al., 2005; Wren et al., 2003, 2001a), and this gives some confidence in our rupture predictions. Furthermore, we used a hyperelastic material model, which cannot take into account time and history dependent viscoelastic properties of tendon. But it has been reported that material behavior of tendon is relatively insensitive to strain rate and there is a minimum amount of hysteresis once the tissue is “preconditioned” (Weiss and Gardiner, 2001). Therefore it is not uncommon to concentrate on the non-linear elastic aspect of the tissue behavior, neglecting the time- and rate-dependent components of the tissue (Quapp and Weiss, 1998; Reese et al., 2010; Wakabayashi et al., 2003). Finally, parameters in our study were fit to a widely used constitutive law suitable at the macro continuum level and does not account for local collagen distribution or collagen fibre twists. The parameters are used for the whole tendon in our work and do not reflect spatial variation, which would require localized experimental data to fit. However, there is limited quantitative information about the microstructural details of human tendon. Moreover we used an optimization method to find material parameters with only six points on the tendons experiments to match, as it was considered prudent to not use too many material property parameters. In conclusion, we have developed a computational framework for subject-specific FE analysis of the human Achilles tendon and generated ten FE models from ultrasound images of in vitro

3604

V.B. Shim et al. / Journal of Biomechanics 47 (2014) 3598–3604

specimens. Our model predictions of tendon rupture showed large variation across specimens, which was consistent with the variation in tendon size and macro scale tissue response. Predicted rupture location was more dependent on geometry while rupture load was more dependent on material properties. Future work includes implementing microstructural features to our model to investigate how tissue microstructure influences tendon rupture and degeneration. Conflict of interest There is no conflict of interest with any financial organizations regarding the material discussed in the manuscript. Acknowledgments This work was funded by the Australian Research Council Linkage Grant (LP110100581). The authors would like to thank Prof. Wren for proving the experimental dataset and for the guidance of Professor Ming Hao Zheng, Professor Jiake Xu, Professor Brett Kirk, Professor Peter Hunter, Dr. Paul Anderson, Dr. Allan Wang, and Assistant Prof. Jonas Rubenson in this research. Appendix A. Supporting information Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.jbiomech.2014.10.001. References Babarenda Gamage, T.P., Rajagopal, V., Ehrgott, M., Nash, M.P., Nielsen, P.M.F., 2011. Identification of mechanical properties of heterogeneous soft bodies using gravity loading. Int. J. Numer. Methods Biomed. Eng. 27, 391–407. Bandak, F.A., Tannous, R.E., Toridis, T., 2001. On the development of an osseoligamentous finite element model of the human ankle joint. Int. J. Solids Struct. 38, 1681–1697. Bayod, J., Becerro-de-Bengoa-Vallejo, R., Losa-Iglesias, M.E., Doblare, M., 2012. Mechanical stress redistribution in the calcaneus after autologous bone harvesting. J. Biomech. 45, 1219–1226. Blemker, S.S., Delp, S.L., 2005. Three-dimensional representation of complex muscle architectures and geometries. Ann. Biomed. Eng. 33, 661–673. Butler, D.L., Grood, E.S., Noyes, F.R., Zernicke, R.F., Brackett, K., 1984. Effects of structure and strain measurement technique on the material properties of young human tendons and fascia. J. Biomech. 17, 579–596. Cheung, J.T., Zhang, M., An, K.N., 2006. Effect of Achilles tendon loading on plantar fascia tension in the standing foot. Clin. Biomech. 21, 194–203. Fernandez, J., Hunter, P., Shim, V., Mithraratne, K., 2012. A subject-specific framework to inform musculoskeletal modeling: outcomes from the IUPS physiome project. In: Calvo Lopez, B., Peña, E. (Eds.), Patient-Specific Computational Modeling. Springer, Netherlands, pp. 39–60. Fernandez, J.W., Mithraratne, P., Thrupp, S.F., Tawhai, M.H., Hunter, P.J., 2004. Anatomically based geometric modelling of the musculo-skeletal system and other organs. Biomech. Model. Mechanobiol. 2, 139–155. Fiorentino, N.M., Rehorn, M.R., Chumanov, E.S., Thelen, D.G., Blemker, S.S., 2014. Computational models predict larger muscle tissue strains at faster sprinting speeds. Med. Sci. Sports. Exerc. 46, 776–786. Gardiner, J.C., Weiss, J.A., 2003. Subject-specific finite element analysis of the human medial collateral ligament during valgus knee loading. J. Orthop. Res.: Off. Publ. Orthop. Res. Soc. 21, 1098–1106. Hansen, P., Aagaard, P., Kjaer, M., Larsson, B., Magnusson, S.P., 2003. Effect of habitual running on human Achilles tendon load-deformation properties and cross-sectional area. J. Appl. Physiol. (1985) 95, 2375–2380. Hess, G.W., 2010. Achilles tendon rupture: a review of etiology, population, anatomy, risk factors, and injury prevention. Foot Ankl. Spec. 3, 29–32. Houghton, L., Dawson, B., Rubenson, J., 2013. Achilles tendon mechanical properties after both prolonged continuous running and prolonged intermittent shuttle running in cricket batting. J. Appl. Biomech. 29, 453–462.

Hunter, P., Robbins, P., Noble, D., 2002. The IUPS human Physiome Project. Pflug. Archiv.: Eur. J. Physiol. 445, 1–9. Jarvinen, T.A., Kannus, P., Maffulli, N., Khan, K.M., 2005. Achilles tendon disorders: etiology and epidemiology. Foot Ankl. Clin. 10, 255–266. Jozsa, L., Lehto, M., Kvist, M., Balint, J.B., Reffy, A., 1989. Alterations in dry mass content of collagen fibers in degenerative tendinopathy and tendon-rupture. Matrix 9, 140–146. Kasashima, Y., Smith, R.K., Birch, H.L., Takahashi, T., Kusano, K., Goodship, A.E., 2002. Exercise-induced tendon hypertrophy: cross-sectional area changes during growth are influenced by exercise. Equine Vet. J. Supplement, S264–S268. Kongsgaard, M., Aagaard, P., Kjaer, M., Magnusson, S.P., 2005. Structural Achilles tendon properties in athletes subjected to different exercise modes and in Achilles tendon rupture patients. J. Appl. Physiol. 99, 1965–1971. Lavagnino, M., Arnoczky, S.P., Kepich, E., Caballero, O., Haut, R.C., 2008. A finite element model predicts the mechanotransduction response of tendon cells to cyclic tensile loading. Biomech. Model. Mechanobiol. 7, 405–416. Lemm, M., Blake, R.L., Colson, J.P., Ferguson, H., 1992. Achilles peritendinitis. A literature review with case report. J. Am. Podiatr. Med. Assoc. 82, 482–490. Lichtwark, G.A., Wilson, A.M., 2006. Interactions between the human gastrocnemius muscle and the Achilles tendon during incline, level and decline locomotion. J. Exp. Biol. 209, 4379–4388. Lichtwark, G.A., Wilson, A.M., 2007. Is Achilles tendon compliance optimised for maximum muscle efficiency during locomotion? J. Biomech. 40, 1768–1775. Maganaris, C.N., Reeves, N.D., Rittweger, J., Sargeant, A.J., Jones, D.A., Gerrits, K., De Haan, A., 2006. Adaptive response of human tendon to paralysis. Muscle Nerve 33, 85–92. Magnusson, S.P., Kjaer, M., 2003. Region-specific differences in Achilles tendon cross-sectional area in runners and non-runners. Eur. J. Appl. Physiol. 90, 549–553. Munro, J.T., Anderson, I.A., Walker, C.G., Shim, V.B., 2013. Finite element analysis of retroacetabular osteolytic defects following total hip replacement. J. Biomech. 46, 2529–2533. Oberhofer, K., Stott, N.S., Mithraratne, K., Anderson, I.A., 2010. Subject-specific modelling of lower limb muscles in children with cerebral palsy. Clin. Biomech. 25, 88–94. Pennisi, E., 2002. Tending tender tendons. Science 295, 1011. Pioletti, D.P., Rakotomanana, L., Benvenuti, J., Leyvraz, P., 1998. Finite element model of the human anterior cruciate ligament. Comput. Methods Biomech. Biomed. Eng. 21, 221–226. Quapp, K.M., Weiss, J.A., 1998. Material characterization of human medial collateral ligament. J. Biomech. Eng. 120, 757–763. Reese, S.P., Maas, S.A., Weiss, J.A., 2010. Micromechanical models of helical superstructures in ligament and tendon fibers predict large Poisson's ratios. J. Biomech. 43, 1394–1400. Reeves, N.D., Maganaris, C.N., Narici, M.V., 2005. Plasticity of dynamic muscle performance with strength training in elderly humans. Muscle Nerve 31, 355–364. Shim, V., Bohme, J., Vaitl, P., Klima, S., Josten, C., Anderson, I., 2010. Finite element analysis of acetabular fractures – development and validation with a synthetic pelvis. J. Biomech. 43, 1635–1639. Shim, V.B., Hunter, P.J., Pivonka, P., Fernandez, J.W., 2011. A multiscale framework based on the physiome markup languages for exploring the initiation of osteoarthritis at the bone–cartilage interface. IEEE Trans. Biomed. Eng. 58, 3532–3536. Shim, V.B., Pitto, R.P., Anderson, I.A., 2012. Quantitative CT with finite element analysis: towards a predictive tool for bone remodelling around an uncemented tapered stem. Int. Orthop. 36, 1363–1369. Stevens, C., Remme, E., LeGrice, I., Hunter, P., 2003. Ventricular mechanics in diastole: material parameter sensitivity. J. Biomech. 36, 737–748. Theobald, P., Benjamin, M., Nokes, L., Pugh, N., 2005. Review of the vascularisation of the human Achilles tendon. Injury 36, 1267–1272. Wakabayashi, I., Itoi, E., Sano, H., Shibuya, Y., Sashi, R., Minagawa, H., Kobayashi, M., 2003. Mechanical environment of the supraspinatus tendon: a twodimensional finite element model analysis. J. Shoulder Elb. Surg./Am. Shoulder Elb. Surg. 12, 612–617. Weiss, J.A., Gardiner, J.C., 2001. Computational modeling of ligament mechanics. Crit. Rev. Biomed. Eng. 29, 303–371. Weiss, J.A., Maker, B.N., Govindjee, S., 1996. Finite element implementation of incompressible, transversely isotropic hyperelasticity. Comput. Method Appl. M 135, 107–128. Wren, T.A., Lindsey, D.P., Beaupre, G.S., Carter, D.R., 2003. Effects of creep and cyclic loading on the mechanical properties and failure of human Achilles tendons. Ann. Biomed. Eng. 31, 710–717. Wren, T.A., Yerby, S.A., Beaupre, G.S., Carter, D.R., 2001a. Influence of bone mineral density, age, and strain rate on the failure mode of human Achilles tendons. Clin. Biomech. 16, 529–534. Wren, T.A., Yerby, S.A., Beaupre, G.S., Carter, D.R., 2001b. Mechanical properties of the human achilles tendon. Clin. Biomech. 16, 245–251.