3D Models Reveal the Influence of Achilles Subtendon Twist on Strain and Energy Storage

Introduction

The Achilles tendon (AT) plays an important but complex role in human movement. This unique tendon transmits forces from the three muscles of the triceps surae to the calcaneus during production of ankle plantarflexion torque. Due to muscle loading, the AT stores and returns energy as it stretches and recoils (Lichtwark and Wilson, 2005; Zelik and Franz, 2017). Each of the triceps surae muscles perform different functions in the generation of propulsion and vertical support during walking (Neptune et al., 2001; Anderson and Pandy, 2003; McGowan et al., 2008; Francis et al., 2013), and each muscle has a unique architecture (Ward et al., 2009; Handsfield et al., 2014; Bolsterlee et al., 2019). Therefore, the three triceps surae muscles apply different forces to the AT (Arndt et al., 1998), which means that the AT functions as three combined tendons, leading to a more complex relationship between structure and function than in a tendon attached to a single muscle.

The evidence of complex loading conditions can be visualized using ultrasound imaging during walking. For example, non-uniform displacements have been observed in the in vivo free tendon; the deep portion of the tendon displaces more than the superficial portion between toe-off and mid-stance, indicating greater elongation of this region (Franz et al., 2015). These kinematic results likely occur due to the combination of complex loading and the complicated structure of the AT. The internal structure of the Achilles free tendon comprises three subtendons (Handsfield et al., 2016), which are distinguishable groups of fascicles that originate from individual muscles: the gastrocnemius lateral head, medial head, and the soleus. These subtendons have been observed in cadavers to twist around each other before inserting into the calcaneus (Cummins and Anson, 1946; Szaro et al., 2009; Edama et al., 2015). The twisted structure occurs in all ATs and the direction of twist is consistent across individuals. However, the amount of subtendon twist varies between individuals, such that previous authors have used that variation to classify tendons into three groups (Cummins and Anson, 1946; Edama et al., 2015; Pękala et al., 2017). These detailed anatomical studies present many questions about the functional consequences of AT internal structure and how mechanical behavior may vary with differences in anatomy. For example, would strain experienced by the tendon during a given elongation vary with differences in the twisted structure of its subtendons? Would differences in subtendon twist change the amount of energy stored during that same stretch?

Computational models enable exploration into the relationships between tendon structure and function. A finite element (FE) model of Achilles subtendons (Handsfield et al., 2017b) demonstrated how sliding and differential loading are possible mechanisms underlying observed non-uniform displacements (Slane and Thelen, 2015). Shim et al. (2018) developed subject-specific models of the Achilles free tendon in which they incorporated variations in fascicle twist. Their study confirmed that varying tendon twist does impact mechanical behavior. However, tendon geometry and material properties also varied in this study and have been shown to be highly variable between individuals and contribute to differences in mechanical behavior (Shim et al., 2014). It remains unclear to what extent observed fascicle twist within the subtendons may influence tendon mechanics, independent of these other variations in the AT.

Variation in subtendon twisted morphology could also affect in vivo measurements of tendon mechanical behavior. For example, strains in the AT are generally estimated by tracking the distance between the distal and proximal endpoints (Kubo et al., 2002; Lichtwark and Wilson, 2005; Onambele et al., 2006), assuming that strain is a linear measure between these two endpoints. Similar methods are employed in combination with force estimation to determine tendon work (Lichtwark and Wilson, 2005; Zelik and Franz, 2017). It is possible that tendon fibers twisted along the length of the tendon may affect the relationship between the tendon tissue strain and energy storage and the longitudinal estimates of strain and energy storage.

The goal of this work was to investigate how differences in subtendon internal twisted structure influences AT fascicle morphology, strains, and energy storage. Our secondary goal was to explore how subtendon twisting may affect quantification of AT strain and energy storage made with in vivo measurements. In order to quantify the effects of twist geometry, independent of differences in free tendon shape and material properties, we built a FE model of the AT with three different internal structures. The different versions of the model represented the average of the three classifications of Achilles subtendon fascicle twisting observed in unloaded cadaver tendons (Pękala et al., 2017) and simulations were performed in which the tendon models were loaded to represent tendon displacement in walking. We analyzed the simulations to (1) determine how functional tendon behavior varied with differences in subtendon twisted morphology, and (2) assess how varied morphology contributed to errors in quantification of strain and energy storage in vivo.

Materials and Methods

Model Geometries

A three dimensional (3D) Achilles free tendon geometry was created in Autodesk Inventor (Autodesk Inc. San Rafael, CA, United States) based on measurements from Pękala et al. (2017). Two elliptical axial cross sections were defined a distance of 60 mm apart and a surface was lofted between them to create a 3D geometry with a volume of 5804.85 mm³ (Figure 1). The proximal cross section (major axis = 17.83 mm, minor axis = 5.76 mm, area = 80.66 mm²) represented the point on the tendon that is just distal to the soleus musculotendinous junction (MTJ). This location was chosen as the most proximal location that the free tendon can be assumed to be in series with all three triceps surae muscles (Epstein et al., 2006). The distal cross section (major axis = 22.04 mm, minor axis = 6.42 mm, area = 111.13 mm²) represented the point on the tendon that is just proximal to the superior portion of the insertion into the calcaneus.

FIGURE 1

Figure 1. Achilles free tendon geometry was developed with elliptical cross sections at the ends. Three unique models with differing amounts of subtendon twist in the undeformed configuration were created by subdividing the free tendon geometry using internal surfaces. These surfaces had the same proximal subtendon cross section divisions and differing distal cross section divisions. Flow guides were used to direct the amount of fascicle twist in each subtendon by creating surfaces that were parallel in the axial cross sections that divided the subtendon volume. Flow guide surfaces are shown in the axial plane for the proximal and distal ends of the lateral gastrocnemius (LG) subtendon in the low twist model. The angle between the flow guides at the ends defined the twist angle of the fascicles within the subtendon. Fascicle tracts in the undeformed models, reconstructed from the resulting fiber directions, are shown.

Using this 3D tendon geometry, three models were created with different internal structures characterized by the amount of subtendon twist. Surfaces divided each 3D model into three subtendons that corresponded with each of the triceps surae muscles: lateral gastrocnemius (LG), medial gastrocnemius (MG), and soleus (SOL). Subtendon divisions were defined such that the proximal cross sections were the same in each model while the distal cross sections varied, resulting in three unique undeformed models with differing sub-structure geometries: low twist, medium twist and high twist (Figure 1). The distribution of volume between the three subtendons was the same in all models (LG = 44%, MG = 27.5%, SOL = 28.5%). The three model geometries corresponded with the average of the classifications of AT torsion (Type I, Type II, and Type III) described by Pękala et al. (2017). These twisted model geometries were the undeformed or resting configuration of models where there was zero stress.

Model Fascicles

Each model was meshed automatically into 3D tetrahedral elements (AMPS Technologies, Pittsburgh, PA, United States). The low twist model contained 4,336 nodes and 16,595 elements, the medium twist model contained 4,147 nodes and 16,168 elements, and the high twist model contained 4,291 nodes and 16,355 elements. Each subtendon was meshed independently so that nodes on adjacent surfaces were not shared. An element convergence analysis was performed by repeating simulations of the low twist model with different mesh densities. The selected mesh was chosen such that the differences in the primary metrics of average along-fiber strain was less than 0.1% and in total strain energy was less than 3% when the number of elements increased by a factor of ten. Additionally, the maximum and minimum strain differed by 5 and 3%, respectively, and the first principle stress differed by less than 3%.

A local fiber direction (a₀) was defined for each element to represent the tendon fascicle structure, with a previously described method utilizing Laplacian simulations (Handsfield et al., 2017a). For each 3D subtendon, fibers were directed from the proximal cross section (inlet surface) to the distal cross section (outlet surface). When characterizing of three tendon twist types, Pękala et al. (2017) dissected the individual subtendons and quantified the average degree of twist of the fascicles in each unloaded subtendon. Based on these measurements, fascicle twist in each undeformed subtendon was enforced by implementing flow guides (Handsfield et al., 2017a) that were internal surfaces that divided the subtendons into four twisting portions. These surfaces were parallel in axial cross sections. The twist angle was defined as the angle in the axial plane between the proximal and distal edges of the flow guides (Figure 1). Unique flow guides resulted in different fascicle twist angles in each model (low twist model: LG = 107°, MG = 17°, SOL = 105°; medium twist model: LG = 157°, MG = 35°, SOL = 145°; high twist model: LG = 211°, MG = 68°, SOL = 200°).

To compute the lengths of the subtendon fascicles in the undeformed models, streamlines, generated at seed points on the proximal surface, were mapped through the field of local fiber direction vectors (a₀) in MATLAB (MathWorks Inc., Natick, MA, United States). These streamlines defined fascicle tracts, which were then truncated to not extend beyond the volume of the subtendon geometry or extrapolated to terminate on the distal surface using a method adapted from Bolsterlee et al. (2017). We defined fascicle lengths as the lengths of the adjusted fascicle tracts as they twisted from the proximal origin to distal insertion were calculated as previously described (Bolsterlee et al., 2017). At least 250 fascicle tracts were created for each model subtendon, representing the internal fascicle geometry in the undeformed condition with zero strain.

Constitutive Model

Subtendons were modeled as transversely isotropic, hyperelastic, quasi-incompressible material (Weiss et al., 1996; Criscione et al., 2001; Blemker et al., 2005). The constitutive model has been described in detail by Blemker et al. (2005) and has the strain energy density function defined in Eq. 1:

$$\text{Φ}(\mathbf{C},{\mathbf{\text{a}}}_0)=W_1+W_2+W_3+{\text{Φ}}^{vol}(1)$$

where a₀is the local fiber direction, C is the right Cauchy-Green deformation tensor. The dilatational portion of the strain energy (${\text{Φ}}^{vol}=\frac{K}{2}\text{ln}(J^2)$) relates to the volume change where $J=\sqrt{\text{det}(\mathbf{C})}$ and depends on a bulk modulus with a value set to K = 5e3 MPa. The strain energy associated with along-fiber shear ($W_1=G_1{(B_1({\overline{I}}_4,{\overline{I}}_5))}^2$) depends on a shear modulus set to G₁ = 3 MPa and the strain energy associated with cross-fiber shear ($W_2=G_2{(B_2({\overline{I}}_1,{\overline{I}}_4,{\overline{I}}_5))}^2$) depends on a shear modulus set to G₂ = 15 MPa (Fiorentino and Blemker, 2014). ${\overline{I}}_1,{\overline{I}}_4,{\overline{I}}_5$ are deviatoric invariants of C. The function for the strain energy associated with along-fiber stretch (W₃) characterizes the relationship between Cauchy stress in the tendon (σ) and the fiber stretch ($\mathrm{\lambda }=\sqrt{{\overline{I}}_4}$) and is defined to be consistent with a piece-wise tendon stress-strain relationship in Eq. 2:

$$\mathrm{\lambda }\frac{\partial W_3}{\partial \mathrm{\lambda }}=\{\begin{array}{c}\hfill \mathrm{\sigma }\left(\mathrm{\lambda }\right)=P_1\left(e^{P_2\left(\mathrm{\lambda }-1\right)}-1\right)\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}1}<\mathrm{\lambda }<{\mathrm{\lambda }}^{*}\hfill \\ \hfill \mathrm{\sigma }\left(\mathrm{\lambda }\right)=P_3\mathrm{\lambda }+P_4\mathit{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} }\mathrm{\lambda }\ge {\mathrm{\lambda }}^{*}\mathit{\hspace{1em}\hspace{1em} }\hfill \end{array}(2)$$

where ${\mathrm{\lambda }}^{{}^{*}}$ represents the fiber stretch at which σ becomes linear and was set to ${\mathrm{\lambda }}^{{}^{*}}$ = 1.03. In the piece-wise equation, P₃ and P₄ were defined so σ is C⁰ and C¹ continuous at $\mathrm{\lambda }={\mathrm{\lambda }}^{{}^{*}}$ P₁ and P₂ were set to values of 1.75 MPa and 48.3, respectively, so that the slope in the linear region was 360 MPa (Onambele et al., 2006).

The constitutive model was implemented in the multi-physics FE analysis program, AMPSol (AMPS Technologies, Pittsburgh, PA, United States), by creating a user-defined hyperelastic material with explicit strain energy function specification.

Model Boundary Conditions

Frictionless sliding contact was assigned between the surfaces of adjacent subtendons in each model. The mechanics of the inter-subtendon matrix in the human AT are unknown and the interfascicular matrix has been shown to allow relative sliding of tendon fascicles in comparative studies (Thorpe et al., 2015), so this approach has been used previously to model Achilles subtendon interaction (Handsfield et al., 2017b). The distal end of each sub tendon was fixed, and the proximal end was constrained to move only in the proximal-distal direction.

To simulate uniaxial loading applied to the AT during walking, pressure boundaries were applied to the proximal surface of each subtendon. The applied pressure on each subtendon was tuned so that in each model the displacement of the proximal surfaces of the MG and LG subtendon were 7.6 mm and the displacement of the proximal SOL subtendon surface was 5.9 mm. These displacements were determined based on measurements of the maximum elongations measured in the superficial and deep portions of the AT during walking (Franz et al., 2015). These in vivo subtendon elongations were estimated from the change in distance between average nodal positions of tendon tissue measured with ultrasound speckle tracking and calcaneus marker positions from motion capture of healthy young adults walking at a speed of 1.25 m/s.

Calculating Strain and Energy Storage

We determined subtendon strain in the longitudinal direction, which we called longitudinal strain, by dividing the change in length by the original length of the subtendon measured in the proximal-distal direction, consistent with in vivo methods of measuring AT strain (Kubo et al., 2002; Lichtwark and Wilson, 2005; Onambele et al., 2006). Longitudinal strain was equivalent to the proximal surface displacement divided by the initial distance between the distal and proximal cross sections (60 mm). We determined strain the fiber direction at the tissue-level, which we called along-fiber strain, by calculating the average along-fiber stretch (λ) in each subtendon and subtracting 1.

We calculated the energy stored in the tendon during loading by integrating the average tendon work relationship for the full AT (force in all subtendons versus the average subtendon length change), which we called full tendon stored energy. This approach was chosen to be consistent with in vivo methods of calculating AT negative work (Lichtwark and Wilson, 2005; Zelik and Franz, 2017). We alternatively calculated the energy stored in the tendon during loading by integrating the average tendon work relationship for the individual subtendons (force in individual subtendon versus that subtendon’s length change) then summing the work done by each subtendon; we called this value summed-subtendon stored energy. We also determined the total strain energy in the AT by integrating the strain energy density (Φ) across all the elements in the subtendon; the strain energy density is directly calculated using the tendon constitutive model.

Kruskal–Wallis tests were used to determine whether differences in fascicle lengths and along-fiber strain occurred between subtendons of the three models. Post hoc Wilcoxon rank sum tests were used to test for differences in subtendon fascicle lengths and along-fiber strains between each pair of models. The Holm–Bonferroni method was used to correct for family wise error rate for tests repeated over the three subtendons. Significance was set at p = 0.05.

Results

Undeformed Subtendon Fascicle Lengths Are Increasingly Longer Than Free Tendon Length at Greater Twist Angles

All subtendon fascicle lengths were longer than the free tendon length of 60mm in the undeformed configuration (Figure 2). For the low twist model, subtendon fascicle lengths were: LG = 61.09 ± 0.19 mm, MG = 61.03 ± 0.12 mm, and SOL = 61.33 ± 0.30 mm. For the medium twist model, subtendon fascicle lengths were: LG = 61.53 ± 0.28 mm, MG = 62.02 ± 0.27 mm, and SOL = 61.85 ± 0.52 mm. For the high twist model, subtendon fascicle lengths were: LG = 62.39 ± 1.13 mm, MG = 62.24 ± 0.39 mm, and SOL = 63.18 ± 0.78 mm. The MG subtendon had the shortest fascicles and the SOL subtendon had the longest fascicles in low and high twist models. Fascicle lengths were more variable in subtendons with higher twist angles. For all subtendons, resting fascicle lengths were significantly shorter (p < 0.00001) in the low compared to the medium twist model, in the low compared to the high twist model, and in the medium compared to the high twist model. The percent differences in average fascicle lengths compared to the free tendon length were relatively small (low twist: LG = 1.80%, MG = 1.71%, and SOL = 2.19%; medium twist: LG = 2.52%, MG = 3.32%, and SOL = 3.04%; high twist: LG = 3.90%, MG = 3.66%, and SOL = 5.16%).

FIGURE 2

Figure 2. Subtendon fascicle lengths in the undeformed models are longer than resting free tendon length. Fascicle lengths are longer in subtendons of models with more twisted geometry. The average length of fascicle tracts in each subtendon are plotted with error bars indicating standard deviations. Significant differences (p < 0.00001) in lengths were found in subtendons compared to the low twist (^∗) and medium twist (+) models. The dotted line indicates the length of the free tendon, which would be consistent with fascicle lengths if no twisting among or within the subtendons occurred. The twist angle of the fascicles within each subtendon are illustrated in the axial plane for each model.

Subtendon Along-Fiber Strains Are Lower Than Longitudinal Strain During Elongation and Are Lower With Greater Subtendon Twist

The longitudinal strain was the same in all models with higher strain in the LG and MG subtendon (12.7%) than the SOL subtendon (9.8%) (Figures 3A,B), corresponding with differential displacement between the deep and superficial portions of the AT (Franz et al., 2015). Average along-fiber strains were lower than longitudinal strains in all models with values for each subtendon decreasing from the low twist model (LG = 9.9 ± 0.8%, MG = 10.5 ± 0.6%, and SOL = 7.9 ± 0.7%) to the medium twist model (LG = 9.2 ± 1.0%, MG = 9.0 ± 1.4%, and SOL = 7.1 ± 0.9%) to the high twist model (LG = 7.2 ± 1.2%, MG = 8.8 ± 0.9%, and SOL = 6.4 ± 1.2%), with significant differences (p < 0.002) in the LG of the low and medium twist models compared to high twist The MG subtendon had the highest average along-fiber strain in the low and high twist models and the SOL subtendon had the lowest average along-fiber strain in all three models (Figure 3A). Along-fiber strains were non-uniform throughout all subtendons in all of the models (Figure 3C). The percent differences in average along-fiber strains compared to the longitudinal strains ranged from 20-55% (low twist: LG = −25.7%, MG = −19.5%, and SOL = −22.5%; medium twist: LG = −33.0%, MG = −37.0%, and SOL = -32.3%; high twist: LG = −55.0%, MG = −38.6%, and SOL = −41.8%).

FIGURE 3

Figure 3. Subtendon along-fiber strains are lower than longitudinal strain and was lower in subtendons of models with greater twist. (A) The longitudinal strain was the same in all models and measured as the change in length in the proximal-distal direction divided by the initial length. The average along-fiber strain was determined in each full tendon model as well as in the individual subtendons, with error bars indicating standard deviations. Significant differences (p < 0.002) in along-fiber strain were found in high twist LG subtendon compared to the low twist (^∗) and medium twist (+) models. (B) All twist models experienced the same elongation in the proximal-distal direction, with greater displacement in the superficial (at the proximal end) subtendons (LG and MG) than the deep subtendon (SOL). (C) Posterior and anterior view of each model show the along-fiber stretch (λ) from each simulation, where a stretch value of one is the zero-strain or undeformed condition.

Energy Stored in the Tendon for Similar Elongations Is Lower With Greater Subtendon Twist

The average subtendon length change was the same in all models but the tendon force required to achieve the same displacement was lower in models with greater twist (low twist = 2.49 kN, medium twist = 2.18 kN, and high twist = 1.90 kN). Therefore, the full tendon stored energy was also lower in models with greater twist (low twist = 7.91 J, medium twist = 6.92 J, high twist = 6.04 J) (Figure 4). Similarly, the summed-subtendon stored energy was lower in models with greater twist (low twist = 8.10 J, medium twist = 7.09 J, and high twist = 6.22 J), though values were slightly higher than the full tendon stored energy. Total strain energy showed a similar trend as tendon stored energy and was lower in models with greater twist (low twist = 6.73 J, medium twist = 5.63 J, and high twist = 4.34 J), though strain energy in all models was less than the energy determined with both longitudinal methods. The strain energy density was non-uniform throughout all subtendons in all of the models, with areas of high energy corresponding with areas of high strain (Figure 3C). The percent differences in total strain energy from the full tendon stored energy increased with greater twist (low twist = −16.2%, medium twist = −20.5%, and high twist = −32.7%), while percent differences in total strain energy from the summed-subtendon stored energy displayed a similar trend and were slightly larger (low twist = −18.5%, medium twist = −23.0%, and high twist = −35.7%).

FIGURE 4

Figure 4. Energy stored in the tendon for similar elongations is lower in models with greater subtendon twist. (A) The total energy stored in each twist type was calculated longitudinally using two methods. Full tendon stored energy was calculated by integrating the force (F) from all subtendons over the average length change (ΔL) of the all subtendons. Summed-subtendon stored energy was calculated by integrating the force (F) from each subtendon over the average length change (ΔL) of that subtendon then summing the results. Total strain energy was also computed using the strain energy density (Φ) and volume (V) of each element. Colored bars show the proportion of energy stored in each subtendon. (B) The average tendon length change is plotted with the total force in black and the area under the curve was calculated to find full tendon stored energy. Subtendon length changes and forces are plotted with colored lines and area under each curve was summed to compute the summed subtendon stored energy. (C) Posterior and anterior view of each model show the strain energy density (Φ) from each simulation, which was multiplied with element volumes and summed to calculate total strain energy.

Discussion

The primary goal of this study was to use models of the Achilles free tendon to explore how mechanical behavior varies with morphology differences in subtendon and fascicle twist. Models predicted that, with more twisted resting geometry, the AT had longer fascicle lengths when undeformed and exhibited reduced along-fiber strain and lower energies stored during elongation, thus altering the free tendon’s response to loading.

Models predicted that increasing the amount of twist of the Achilles subtendon resting morphology effectively lowered stiffness in the free tendon even though the tissue level material properties were held constant. In the models, the tissue-level strain varied with twist while the elongation was the same. The average along-fiber strain was lower in models with greater twist (Figure 3A), leading to larger differences in the strain at the tissue-level compared to the tendon-level strain. The distribution of strains within the tissue differed as well (Figure 3C). Similar to our results, AT FE models with twisted geometry developed by Shim et al. (2018) predicted that greater geometric twist redistributed internal stresses during loading allowing larger loads to be applied in simulations before stress in the tissue reached a given rupture limit. Tendon twisted geometry could be a mechanism to reduce the tendon tissue strains experienced during a given muscle-tendon-unit (MTU) excursion. A more twisted morphology may therefore aid in avoiding rupture as failure strain of tendon has been shown to be highly conserved between different tendons in different species (LaCroix et al., 2013).

In our simulations where tendon elongation was the same, greater twist corresponded with lower stored energy at both the tendon level (longitudinal stored energy) and at the tissue level (strain energy) (Figure 4A). Longitudinal energy differences can be attributed to the varied forces applied in each model for the same elongation (Figure 4B), while the strain energy differences were also influenced by variation in tissue-level strain (Figure 3A). Twisting may actually improve the efficiency of energy return by the tendon, as tendon-level strain and energy storage requires less tissue deformation that could lead to conformational changes at the collagen level that would result in energy loss. In fact, equine energy storing tendons with helical substructures have been shown to experience less hysteresis loss than postural tendons with less fascicle rotation (Thorpe et al., 2013). Furthermore, the lower forces required in higher twist tendon models to elongate as much as the low twist model (Figure 4B), suggest that the amount of twist would impact series-elasticity of the musculotendon unit. A more compliant AT may enable greater force control by the triceps surae muscles (Alexander, 2002) and has the potential to alter the efficiency of these muscles (Lichtwark and Wilson, 2007; Uchida et al., 2016).

It is unknown if the variation in morphologic twist between individuals is an adaptation to mechanical stimuli or simply a product of anatomical variability. However, the amount of twist seems to provide tradeoffs in injury prevention and energetic efficiency. The twist of the subtendons and the twist angle of their fascicles may be an important consideration in the repair of AT ruptures. As the amount of twist does not vary much between the right and left sides (Pękala et al., 2017), twist in the contralateral tendon could serve as a reference in reconstructing the ruptured tendon. Alternatively, twist could be applied before suturing in an effort to preserve compliance as elasticity decreases (Karatekin et al., 2018) or to protect the repaired tendon from further injury. This approach to improving AT reconstruction is an exciting opportunity for further studies.

Each of the three subtendons had different twist angle in all undeformed models, resulting in different average fascicle lengths (Figure 2). The MG subtendon had lower twist angles and generally shorter fascicles compared to the LG subtendon, and therefore experienced higher along-fiber strains even though longitudinal strains were the same in both subtendons (Figure 3). These results suggest that the MG subtendon is more vulnerable to injury due to its lower twist angle, which may help explain how failure of a single subtendon occurs, resulting in a partial tear in the Achilles (Smigielski, 2008). Further studies could work to determine if tears occur more often in the less twisted MG subtendon. The SOL subtendon had the greatest twist angle in all models and generally had longer fascicles. Subtendon twist within the free tendon influenced fascicle lengths in addition to the fascicle twist angle, so in the case of the medium twist model, the MG subtendon had longer fascicles despite its lower twist angle. The difference in twist angles between these subtendons could be associated with the diversity in morphology of the muscles they attach to. Dean et al. (2007) demonstrated that a twisted tendon in the jaw of a spotted ratfish facilitated more uniform strains in the fibers of a muscle with a broad attachment, so that these fibers could operate at similar lengths on the force-length curve. The triceps surae muscles also have broad attachments though the attachment morphology of each muscle is different (Dalmau-Pastor et al., 2014), possibly necessitating different subtendon twist angles to achieve the same equalizing effect on fiber operating lengths.

Our secondary goal was to assess how morphological complexity of the AT impacts how in vivo measurements of tendon behavior should be interpreted. We hypothesized that varying lengths and geometry resulting from twisted morphology may lead to errors during in vivo calculations of strain and energy. We found that resting fascicle lengths were only slightly underestimated by 1.5-5.5% by using measurements of tendon length. However, strains and energy storage measured at the fascicle level were both overestimated by 20–55 and 15–35%, respectively when using longitudinal methods that are consistent with in vivo measurements. The amount of error increases with twist. Unfortunately, definitively calculating measurement errors due to subtendon twist is not possible as methods of determining the amount of twist in vivo have not been developed, to our knowledge. Hopefully the future development of imaging techniques will allow for correction of such measurement errors in the future. Additionally, the high twist group accounts for less than 6% of the population studied in Pękala et al. (2017) so most tendons will likely have a low (48% of population) or medium (46% of population) amount of twist. Therefore, the extent of error associated with the low (strain = 19–26%, energy = 16–19%) and medium (strain = 32–37%, energy = 20–23%) twist models will represent most of the population included in in vivo studies.

There are several limitations to this study that should be noted. Our goal was to characterize the mechanical consequences of varying subtendon twist, independent of variation in tendon geometry or material properties. To this end we developed a model with a generic geometry so that we could create variations in internal structure and examine the effects on tendon behavior when loaded uniaxially. Future work to develop more detailed models is needed to fully investigate how subtendon twist influences in vivo AT behavior. Methods enabling in vivo quantification of subtendon twist would improve subject-specific models. Twist could be incorporated with other model inputs, like geometry and material properties, that tendon behavior is sensitive to Shim et al. (2014). Furthermore, validation with in vivo data is required to better understand the predictions made in this study. We chose to use the measurements of subtendon structure reported by Pękala et al. (2017) because of the extent of quantitative anatomical data provided. We are aware that the subtendon proportions and cross section in that study (Pękala et al., 2017), deviate somewhat from those reported by previous authors (Szaro et al., 2009; Edama et al., 2015). Since measurements in these studies were performed ex vivo in dissected tendons, a major assumption of our model was that subtendon twist was a morphological characteristic that could exist independent of loading. Therefore, we implemented twist in the geometry of the undeformed models, as opposed to applying a torsional loading condition to achieve twist. AT rotation has been observed in response to in vivo loading (Obst et al., 2014), and model twist likely changes during simulated loading. All model twist angles are reported for the undeformed configuration.

All tendon measurements used to create model geometries were made in fresh frozen cadavers (Pękala et al., 2017), though simulation conditions were chosen to represent in vivo tendon loading (Franz et al., 2015). Displacements of subtendon proximal ends were controlled to match elongations estimated with speckle tracking in the free tendon. In the in vivo experiments, elongation was calculated from the proximal end to the calcaneal marker. The distal end of the models was located at the superior edge of the calcaneus, which resulted in shorter model tendons and therefore higher strains than in the in vivo study. Although AT strains of this magnitude have been reported during one-legged hopping (Lichtwark and Wilson, 2005), model predicted strains are likely larger than what occurs during walking. However, all models experienced the same longitudinal strains, allowing for comparison of tissue level strains due to differences in subtendon twist. Further, conclusions about strain in each subtendon are difficult to interpret from these results as the free tendon modeled here does not capture full external portion of proximal tendon associated with lateral or MG head. Simulations were performed quasi-statically and only included tendon loading. Further work is needed to incorporate more detailed viscoelastic behavior of AT in order to accurately simulate full tendon work loops. We assumed the tendon was unloaded prior to simulated elongation and, thus, we applied no initial stretch to the models. Estimates of in vivo AT loading show that the tendon is not stress-free at the beginning of the gait cycle (Keuler et al., 2019). As the start of our simulations correspond to this time point, it may be appropriate to apply a pre-stretch to the subtendons, which would affect when the tendon would transition from the toe region to the linear portion of the stress-strain curve, leading to higher stresses at the strains enforced in this study. While estimates of stress in the full tendon exist, it is unclear how loading may be distributed between the subtendons. Evidence of differences in slack angles of the triceps surae muscles (Hirata et al., 2015) suggests that the initial stretch at a given joint angle would vary between the subtendons, though current methods are unable to estimate what these loads should be. An exciting direction for future work would be to investigate how inhomogeneous subtendon loading in addition to variation in subtendon morphology influence the behavior of the AT.

In conclusion, the models developed in this paper of the AT with varied subtendon twisted geometry help us understand how this morphological characteristic can result in different amounts of tissue strain and energy storage within the tendon in response to similar loading. High twist angles in tendon fascicles can contribute to errors in quantifying these mechanical behaviors when methods that rely on 2D measurements at the endpoints are employed. The knowledge of this effect will aid in the interpretation of future studies of AT behavior and inspire future work to design methods that enable measurements of in vivo subtendon structure.

Data Availability Statement

The datasets generated for this study are available on request to the corresponding author.

Author Contributions

KK and SB: conceptualization, writing – review, and editing. KK: model development and analysis and writing – original draft preparation. SB: funding acquisition. Both authors contributed to the article and approved the submitted version.

Funding

Funding for the completion of this study was provided by the National Institutes of Health grant R01AG051748 and grant U01AR069393.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors would like to acknowledge Geoffrey Handsfield, Darryl Thelen, Jason Franz, and all other Tendonados for their assistance and insight related to this project.

References

Alexander, R. M. (2002). Tendon elasticity and muscle function. Comp. Biochem. Physiol. A Mol. Integr. Physiol. 133, 1001–1011. doi: 10.1016/S1095-6433(02)00143-5

CrossRef Full Text | Google Scholar

Anderson, F. C., and Pandy, M. G. (2003). Individual muscle contributions to support in normal walking. Gait Posture 17, 159–169. doi: 10.1016/S0966-6362(02)00073-5

CrossRef Full Text | Google Scholar

Arndt, A. N., Komi, P. V., Brüggemann, G. P., and Lukkariniemi, J. (1998). Individual muscle contributions to the in vivo achilles tendon force. Clin. Biomech. 13, 532–541. doi: 10.1016/S0268-0033(98)00032-1

CrossRef Full Text | Google Scholar

Blemker, S. S., Pinsky, P. M., and Delp, S. L. (2005). A 3D model of muscle reveals the causes of nonuniform strains in the Biceps Brachii. J. Biomech. 38, 657–665. doi: 10.1016/j.jbiomech.2004.04.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Bolsterlee, B., D’Souza, A., Gandevia, S. C., and Herbert, R. D. (2017). How does passive lengthening change the architecture of the human medial gastrocnemius muscle? J. Appl. Physiol. 122, 727–738. doi: 10.1152/japplphysiol.00976.2016

PubMed Abstract | CrossRef Full Text | Google Scholar

Bolsterlee, B., D’Souza, A., and Herbert, R. D. (2019). Reliability and robustness of muscle architecture measurements obtained using diffusion tensor imaging with anatomically constrained tractography. J. Biomech. 86, 71–78. doi: 10.1016/J.JBIOMECH.2019.01.043

PubMed Abstract | CrossRef Full Text | Google Scholar

Criscione, J. C., Douglas, A. S., and Hunter, W. C. (2001). Physically based strain invariant set for materials exhibiting transversely isotropic behavior. J. Mech. Phys. Solids 49, 871–897. doi: 10.1016/S0022-5096(00)00047-8

CrossRef Full Text | Google Scholar

Cummins, E. J., and Anson, B. J. (1946). The structure of the calcaneal tendon (of Achilles) in relation to orthopedic surgery, with additional observations on the plantaris muscle. Surg. Gynecol. Obstet. 83, 107–116.

Google Scholar

Dalmau-Pastor, M., Fargues-Polo, B., Casanova-Martínez, D., Vega, J., and Golanó, P. (2014). Anatomy of the triceps surae. Foot Ankle Clin. 19, 603–635. doi: 10.1016/j.fcl.2014.08.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Dean, M. N., Azizi, E., and Summers, A. P. (2007). Uniform strain in broad muscles: active and passive effects of the twisted tendon of the spotted ratfish hydrolagus colliei. J. Exp. Biol. 210, 3395–3406. doi: 10.1242/jeb.007062

PubMed Abstract | CrossRef Full Text | Google Scholar

Edama, M., Kubo, M., Onishi, H., Takabayashi, T., Inai, T., Yokoyama, E., et al. (2015). The twisted structure of the human achilles tendon. Scand. J. Med. Sci. Sports 25, e497–e503. doi: 10.1111/sms.12342

PubMed Abstract | CrossRef Full Text | Google Scholar

Epstein, M., Wong, M., and Herzog, W. (2006). Should tendon and aponeurosis be considered in series? J. Biomech. 39, 2020–2025. doi: 10.1016/j.jbiomech.2005.06.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Fiorentino, N. M., and Blemker, S. S. (2014). Musculotendon variability influences tissue strains experienced by the biceps femoris long head muscle during high-speed running. J. Biomech. 47, 3325–3333. doi: 10.1016/j.jbiomech.2014.08.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Francis, C. A., Lenz, A. L., Lenhart, R. L., and Thelen, D. G. (2013). The modulation of forward propulsion, vertical support, and center of pressure by the plantarflexors during human walking. Gait Posture 38, 993–997. doi: 10.1016/j.gaitpost.2013.05.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Franz, J. R., Slane, L. C., Rasske, K., and Thelen, D. G. (2015). Non-uniform in vivo deformations of the human achilles tendon during walking. Gait Posture 41, 192–197. doi: 10.1016/j.gaitpost.2014.10.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Handsfield, G. G., Bolsterlee, B., Inouye, J. M., Herbert, R. D., Besier, T. F., and Fernandez, J. W. (2017a). Determining skeletal muscle architecture with laplacian simulations: a comparison with diffusion tensor imaging. Biomech. Model. Mechanobiol. 16, 1845–1855. doi: 10.1007/s10237-017-0923-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Handsfield, G. G., Inouye, J. M., Slane, L. C., Thelen, D. G., Miller, G. W., and Blemker, S. S. (2017b). A 3D model of the achilles tendon to determine the mechanisms underlying nonuniform tendon displacements. J. Biomech. 51, 17–25. doi: 10.1016/j.jbiomech.2016.11.062

PubMed Abstract | CrossRef Full Text | Google Scholar

Handsfield, G. G., Meyer, C. H., Hart, J. M., Abel, M. F., and Blemker, S. S. (2014). Relationships of 35 lower limb muscles to height and body mass quantified using MRI. J. Biomech. 47, 631–638. doi: 10.1016/j.jbiomech.2013.12.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Handsfield, G. G., Slane, L. C., and Screen, H. R. C. (2016). Nomenclature of the tendon hierarchy: an overview of inconsistent terminology and a proposed size-based naming scheme with terminology for multi-muscle tendons. J. Biomech. 49, 3122–3124. doi: 10.1016/j.jbiomech.2016.06.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Hirata, K., Kanehisa, H., Miyamoto-Mikami, E., and Miyamoto, N. (2015). Evidence for intermuscle difference in slack angle in human triceps surae. J. Biomech. 48, 1210–1213. doi: 10.1016/j.jbiomech.2015.01.039

PubMed Abstract | CrossRef Full Text | Google Scholar

Karatekin, Y. S., Karaismailoglu, B., Kaynak, G., Ogut, T., Dikici, A. S., Esmerer, E. U., et al. (2018). Does elasticity of achilles tendon change after suture applications? Evaluation of repair area by acoustic radiation force impulse elastography. J. Orthop. Surg. Res. 13:45. doi: 10.1186/s13018-018-0751-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Keuler, E. M., Loegering, I. F., Martin, J. A., Roth, J. D., and Thelen, D. G. (2019). Shear wave predictions of achilles tendon loading during human walking. Sci. Rep. 9:13419. doi: 10.1038/s41598-019-49063-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Kubo, K., Kawakami, Y., Kanehisa, H., and Fukunaga, T. (2002). Measurement of viscoelastic properties of tendon structures in vivo. Scand. J. Med. Sci. Sports 12, 3–8. doi: 10.1034/j.1600-0838.2002.120102.x

PubMed Abstract | CrossRef Full Text | Google Scholar

LaCroix, A. S., Duenwald-Kuehl, S. E., Lakes, R. S., and Vanderby, R. (2013). Relationship between tendon stiffness and failure: a metaanalysis. J. Appl. Physiol. 115, 43–51. doi: 10.1152/japplphysiol.01449.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Lichtwark, G. A., and Wilson, A. M. (2005). In vivo mechanical properties of the human achilles tendon during one-legged hopping. J. Exp. Biol. 208, 4715–4725. doi: 10.1242/jeb.01950

PubMed Abstract | CrossRef Full Text | Google Scholar

Lichtwark, G. A., and Wilson, A. M. (2007). Is achilles tendon compliance optimised for maximum muscle efficiency during locomotion? J. Biomech. 40, 1768–1775. doi: 10.1016/j.jbiomech.2006.07.025

PubMed Abstract | CrossRef Full Text | Google Scholar

McGowan, C. P., Neptune, R. R., and Kram, R. (2008). Independent effects of weight and mass on plantar flexor activity during walking: implications for their contributions to body support and forward propulsion. J. Appl. Physiol. 105, 486–494. doi: 10.1152/japplphysiol.90448.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Neptune, R. R., Kautz, S. A., and Zajac, F. E. (2001). Contributions of the individual ankle plantar flexors to support, forward progression and swing initiation during walking. J. Biomech. 34, 1387–1398. doi: 10.1016/S0021-9290(01)00105-1

CrossRef Full Text | Google Scholar

Obst, S. J., Renault, J.-B., Newsham-West, R., and Barrett, R. S. (2014). Three-dimensional deformation and transverse rotation of the human free achilles tendon in vivo during isometric plantarflexion contraction. J. Appl. Physiol. 116, 376–384. doi: 10.1152/japplphysiol.01249.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Onambele, G. L., Narici, M. V., and Maganaris, C. N. (2006). Calf muscle-tendon properties and postural balance in old age. J. Appl. Physiol. 100, 2048–2056. doi: 10.1152/japplphysiol.01442.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Pękala, P. A., Henry, B. M., Ochała, A., Kopacz, P., Tatoń, G., Młyniec, A., et al. (2017). The twisted structure of the achilles tendon unraveled: a detailed quantitative and qualitative anatomical investigation. Scand. J. Med. Sci. Sports 27, 1705–1715. doi: 10.1111/sms.12835

PubMed Abstract | CrossRef Full Text | Google Scholar

Shim, V. B., Fernandez, J. W., Gamage, P. B., Regnery, C., Smith, D. W., Gardiner, B. S., et al. (2014). Subject-specific finite element analysis to characterize the influence of geometry and material properties in achilles tendon rupture. J. Biomech. 47, 3598–3604. doi: 10.1016/j.jbiomech.2014.10.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Shim, V. B., Handsfield, G. G., Fernandez, J. W., Lloyd, D. G., and Besier, T. F. (2018). Combining in silico and in vitro experiments to characterize the role of fascicle twist in the achilles tendon. Sci. Rep. 8:13856. doi: 10.1038/s41598-018-31587-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Slane, L. C., and Thelen, D. G. (2015). Achilles tendon displacement patterns during passive stretch and eccentric loading are altered in middle-aged adults. Med. Eng. Phys. 37, 712–716. doi: 10.1016/j.medengphy.2015.04.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Smigielski, R. (2008). Management of partial tears of the gastro-soleus complex. Clin. Sports Med. 27, 219–229. doi: 10.1016/j.csm.2007.10.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Szaro, P., Witkowski, G., Śmigielski, R., Krajewski, P., and Ciszek, B. (2009). Fascicles of the adult human achilles tendon – an anatomical study. Ann. Anat. 191, 586–593. doi: 10.1016/j.aanat.2009.07.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Thorpe, C. T., Godinho, M. S. C., Riley, G. P., Birch, H. L., Clegg, P. D., and Screen, H. R. C. (2015). The interfascicular matrix enables fascicle sliding and recovery in tendon, and behaves more elastically in energy storing tendons. J. Mech. Behav. Biomed. Mater. 52, 85–94. doi: 10.1016/j.jmbbm.2015.04.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Thorpe, C. T., Klemt, C., Riley, G. P., Birch, H. L., Clegg, P. D., and Screen, H. R. C. (2013). Helical sub-structures in energy-storing tendons provide a possible mechanism for efficient energy storage and return. Acta Biomater. 9, 7948–7956. doi: 10.1016/j.actbio.2013.05.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Uchida, T. K., Hicks, J. L., Dembia, C. L., and Delp, S. L. (2016). Stretching your energetic budget: how tendon compliance affects the metabolic cost of running. PLoS One 11:e0150378. doi: 10.1371/journal.pone.0150378

PubMed Abstract | CrossRef Full Text | Google Scholar

Ward, S. R., Eng, C. M., Smallwood, L. H., and Lieber, R. L. (2009). Are current measurements of lower extremity muscle architecture accurate? Clin. Orthop. Relat. Res. 467, 1074–1082. doi: 10.1007/s11999-008-0594-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Weiss, J. A., Maker, B. N., and Govindjee, S. (1996). Finite element implementation of incompressible, transversely isotropic hyperelasticity. Comput. Methods Appl. Mech. Eng. 135, 107–128. doi: 10.1016/0045-7825(96)01035-3

CrossRef Full Text | Google Scholar

Zelik, K. E., and Franz, J. R. (2017). It’s positive to be negative: achilles tendon work loops during human locomotion. PLoS One 12:e0179976. doi: 10.1371/journal.pone.0179976

PubMed Abstract | CrossRef Full Text | Google Scholar