Regionalizing muscle activity causes changes to the magnitude and direction of the force from whole muscles—a modeling study

Skeletal muscle can contain neuromuscular compartments that are spatially distinct regions that can receive relatively independent levels of activation. This study tested how the magnitude and direction of the force developed by a whole muscle would change when the muscle activity was regionalized within the muscle. A 3D finite element model of a muscle with its bounding aponeurosis was developed for the lateral gastrocnemius, and isometric contractions were simulated for a series of conditions with either a uniform activation pattern, or regionally distinct activation patterns: in all cases the mean activation from all fibers within the muscle reached 10%. The models showed emergent features of the fiber geometry that matched physiological characteristics: with fibers shortening, rotating to greater pennation, adopting curved trajectories in 3D and changes in the thickness and width of the muscle belly. Simulations were repeated for muscle with compliant, normal and stiff aponeurosis and the aponeurosis stiffness affected the changes to the fiber geometry and the resultant muscle force. Changing the regionalization of the activity resulted to changes in the magnitude, direction and center of the force vector from the whole muscle. Regionalizing the muscle activity resulted in greater muscle force than the simulation with uniform activity across the muscle belly. The study shows how the force from a muscle depends on the complex interactions between the muscle fibers and connective tissues and the region of muscle that is active.

1. Introduction

Skeletal muscles can contain subunits called neuromuscular compartments that are spatially distinct regions that contain specific motor units and motor drive from the nervous system (English et al., 1993). In muscles with broad attachments, a relationship between anatomical compartments and function may appear logical, and this has shown to be the case for both the biceps femoris in the cat (English and Weeks, 1987; Chanaud et al., 1990) and the masseter muscle in the pig (Herring et al., 1989, 1991). However, functional regionalization in muscles with long tendons has also been reported (Carrasco et al., 1999), leading to the suggestion that activation of motor units in different compartments may result in differences to both the direction and the magnitude of force applied at the tendon (English et al., 1993). It is likely that asymmetry in the fascicle architecture combines with the location of the neuromuscular compartments to result in varied force vectors from a contracting muscle.

Unipennate muscle is asymmetrical in its architecture, and muscle fibers in different locations have different moment arms and may exert different torques about a joint. The way in which forces are transmitted from the contractile fibers to the tendon can involve myofascial pathways (Huijing, 2003), that in turn may modify the resultant force vector from the individual fibers. It has been shown that activity can differ between the neuromuscular compartments and spatial regions in the gastrocnemii in the cat during walking (English, 1984), and in man during both cycling and postural tasks (Wakeling, 2009; Hodson-Tole et al., 2013). Changes in activity between the neuromuscular compartments in the cat lateral gastrocnemius led to changes in the direction of the force vector along the tendon, altering the moments of yaw, pitch or roll about the ankle (Carrasco et al., 1999). However little has been reported about the mechanisms that link the varied forces developed by individual fibers to the net mechanical output of the whole muscle.

The different stiffnesses of connective tissues such as aponeurosis and tendon add to the complexity of muscle-tendon unit (MTU) behavior. In-vitro measurements of mechanical properties of tendon and aponeurosis (Wren et al., 2001; Azizi et al., 2009; Lake et al., 2010), and ultrasound-based in vivo measurements of these properties (Maganaris and Paul, 2000; Magnusson et al., 2003) suggest that tendon and aponeurosis may have different tensile elastic moduli which can be alterd with age (Onambele et al., 2006), training (Kubo et al., 2002) and injury (Arampatzis et al., 2007). The elastic properties of the aponeurosis can affect the extent to which muscle fibers rotate as they contract (Randhawa et al., 2013), and can thus affect the magnitude and direction of the forces developed by the whole muscle. However, it is beyond current experimental techniques to measure the effect of aponeurosis stiffness on the force outputs from muscle.

Some of the limitations in experimental studies can be addressed using in silica models of muscles. These models need to contain a realistic architecture and physiologically relevant connective tissue properties. Further, they should be able to support different activation levels in different regions. Implementing fundamental physiological concepts and material properties within sophisticated mathematical frameworks has moved muscle simulations from one-dimensional models (Zajac, 1989; Delp et al., 1990) to more architecturally and functionally detailed two—(Van Leeuwen and Spoor, 1992; Otten and Hulliger, 1994) and three-dimensional models (Oomens et al., 2003; Blemker et al., 2005; Bol and Reese, 2008). Despite the level of architectural details that current models include, such models rarely include other heterogeneities within the muscle such as material distribution (e.g., fiber-type or connective tissue properties) or differential patterns of activation. In this study we investigate how uneven patterns of activation across a unipennate muscle affect the magnitude and direction of the force developed by the whole muscle, using a conceptual in silica muscle mode. The in silica model has a simple geometry but is assymetric in architecture, and regionalized in activity. This model was also used to investigate how the aponeurosis properties affect the force development in the fibers. and its transmission to the external tendon.

2. Materials and Methods

A unipennate muscle model was created in silica to test the effects of the activation being regionalized on the direction and magnitude of the force developed by the whole muscle. The model had a regularized shape to help constrain model variants and results to the conceptual questions which are the focus of this study, rather than allowing the model to respond to idiosyncrasies of individual geometries. The dimensions, as well as the structural and material properties of the model were styled to be consistent with those of the lateral gastrocnemius in humans. The soft tissues were treated as transversally isotropic hyperelastic materials. The muscle-aponeuroses complex was meshed by a grid of hexahedral elements. Displacements, stresses and forces were calculated using a three-field finite element formulation. The activation levels of the different muscle fibers within the model could be independently varied. The methods for each of these parts are described in more detail below.

2.1. Geometry, Mesh, Boundary, and Fiber Architecture

A simplistic geometry (Figure 1) was developed for human lateral gastrocnemius (LG) (Dimensions mostly from Randhawa et al., 2013). For the purpose of this study, the muscle belly was modeled without the attached external tendon. This allowed the belly to remain isometric and constant between conditions. The cross-section of the undeformed muscle belly was a parallelogram. The short side of this parallogram was 65 mm, the same as the length of a muscle fiber. The smaller angle within the parallelogram was 15°, which was also the the angle a fiber made with the aponeurosis at both ends of attachment. This could be considered as the initial pennation in the muscle. However, the angle between the fiber and the line of action of the muscle was smaller. This is because the line of action was calculated to be the diagonal of the parallelogram. The initial width of the muscle belly was set at 55 mm.

FIGURE 1

Figure 1. Simplified geometry of lateral gastrocnemius. Aponeurosis tissue in dark gray and muscle tissue in light gray. The origin of the muscle coordinate system was set to the bottom right corner of the deep aponeurosis, and the axes were aligned with x, width; y, thickness, z, belly length.

Aponeuroses covered the two long sides of the muscle tissue. Each aponeurosis was 210 mm in length and a constant thickness (depth) of 3 mm, with a width matching that of the muscle belly.

Integration points in the in silica model were assigned material properties appropriate to aponeurosis or muscle fibers, relative to their position in the belly. The geometry was fixed at one opposing end of each aponeurosis and the other boundary surfaces were traction free.

2.2. Material Properties of the in silica Muscle Model

For purposes of completeness and setting notation, we first describe the mathematical background for our model. The soft tissues (muscle and aponeurosis) in this project were mathematically described by a fiber-reinforced composite material (Spencer, 1984). Specifically, they are described as nearly incompressible (e.g., for muscle Baskin and Paolini, 1967), transversally isotropic hyperelastic materials. If we denote the position vector of a point in the deformed tissue by x, it can be expressed in terms of the position vector of the same point in the undeformed configuration, X, and the displacement vector u in the natural way: x = X + u. The response of a nearly-incompressible material to loading can be decomposed into a volume-changing (volumetric) and a volume-preserving (isochoric) part. Specifically, if $F:=\left[\frac{\partial x_i}{\partial X_j}\right]$ denotes the usual deformation gradient (a second-order tensor) and J : = det(F) denotes the dilation, then we can define the isochoric part of the deformation gradient as F:

$\begin{array}{ll}\overline{F}=\frac{1}{J^{1/3}}F.\hfill & (1)\hfill \end{array}$

The deformation gradient is used to calculate the left Cauchy-Green tensor (the Finger tensor) B and the Lagrangian finite strain tensor E as:

$\begin{array}{ll}B:=F{F}^T=\left[\frac{\partial x_i}{\partial X_k}\frac{\partial x_j}{\partial X_k}\right],E:=\frac{1}{2}(F^{T}F-I).\hfill & (2)\hfill \end{array}$

Using the decomposition in (1), we can obtain the isochoric part of the left Cauchy-Green tensor, which we denote as B:

$\begin{array}{ll}\overline{B}:=\overline{F}{\overline{F}}^T=\frac{1}{J^{2/3}}F{F}^T=\frac{1}{J^{2/3}}B.\hfill & (3)\hfill \end{array}$

As is standard for hyperelastic materials undergoing finite deformations, the mechanical response is described in terms of the strain energy density function W. As with the deformation gradients, this strain energy function can be decomposed into its volumetric and isochoric parts:

$\begin{array}{ll}W(B)=W_{vol}(J)+W_{iso}(\overline{B}).\hfill & (4)\hfill \end{array}$

The volumetric component W_vol of the strain energy describes the nearly incompressible behavior of the tissues,

$\begin{array}{ll}{W}_{vol}(J):=\frac{\kappa }{4}\left(J^2-1-2\log (J)\right).\hfill & (5)\hfill \end{array}$

The volumetric strain energies of muscle and the aponeurosis are distinguished by different choices of the (constant) bulk modulus, κ. In our experiments we used κ = 10⁶ Pa for muscle, and κ = 10⁸ Pa for the aponeurosis.

The isochoric contributions to the strain energy W_iso are described in terms of one (Yeoh, 1993) or more (Blemker et al., 2005) invariants (I₁…I₅) of the left Cauchy-Green tensor B, the initial orientation of fiber in the tissue (a₀) and the along-fiber stretch λ:

$\begin{array}{cc}\begin{array}{l}{I}_1:=tr(B),I_2:=\frac{1}{2}\left[{(tr(B))}^2-tr(B^2)\right],I_3:=det(B)=J^2,\\ I_4:=a_0·B·a_0={\lambda }^2,I_5:=a_0·B^2·a_0.\end{array}& (6)\end{array}$

Using these invariants we can describe the along fiber isochoric strain energy (W_tissue) and the base material isochoric energy (W_base):

$\begin{array}{ll}{W}_{iso}=W_{tissue}+W_{base}.\hfill & (7)\hfill \end{array}$

The contribution of the base material to the strain energy, W_base, encapsulates the elastic properties of the connective tissue within both muscle and the aponeurosis (e.g., the extracellular connective tissue in muscle). These contributions are mathematically modeled as approximations of an exponential fit for the first invariant I₁ to real data. Both the muscle and the aponeurosis possess their own distinct base material properties. In this paper, the base material isochoric energy is described using a model originally due to Yeoh (1993), which is cubic in I₁:

$\begin{array}{ll}\begin{array}{l}{W}_{base,muscle}:={\displaystyle \sum _{i=1}^3{c}_i{(I_1-3)}^i},c_1=6.75\times {10}^{-3}\text{Pa},\\ c_2=2.78\times {10}^{-2}\text{Pa},c_3=-1.9745\times {10}^{-3}\text{Pa}.\end{array}\hfill & (8)\hfill \end{array}$

The base material isochoric energy for the aponeurosis is described using a Humphrey (exponential) model (Humphrey and Yin, 1987),

$\begin{array}{ll}\begin{array}{l}{W}_{base,apo:=}{c}_1\left(e^{c_2(I_1-3)}-1\right),c_1=4.3510\times {10}^4\text{Pa},\\ c_2=5.796\times {10}^2\text{Pa}.\end{array}\hfill & (9)\hfill \end{array}$

The isochoric strain energy contributions, W_tissue, arise from the stretching of fibers along their length. If we denote the Cauchy stress in the fiber caused by a stretch of λ as σ_tissue(λ), then the isochoric strain energy for both muscle and aponeurosis becomes

$\begin{array}{ll}\lambda \frac{\partial W_{tissue}(\lambda )}{\partial \lambda }={\sigma }_{tissue}(\lambda )\hfill & (10)\hfill \end{array}$

Muscle and aponeuroses tissues are distinguished by the specific form of their Cauchy stress σ_tissue(λ) used in (10). For this study, we constructed stress-stretch relationships from experimental data (see below), using the curve-fitting functions in (MATLAB, 2014).

For the aponeurosis, σ_tissue ≡ σ_Apo was obtained by a piecewise exponential fit of stress-stretch curves from Magnusson et al. (2003) where human gastrocnemius aponeurosis tensile properties were measured. See Figure 2A on page 5.

$\begin{array}{ll}{\sigma }_{Apo}(\lambda )=\{\begin{array}{l}{10}^6\times 3.053({\lambda }^{124.6}-1)\text{Pa},1\le \lambda \le 1.025,\hfill \\ \begin{array}{l}{10}^6\times 3.053[17374.99\hfill \\ (\lambda -1.025)+20.687]\text{Pa},\hfill \end{array}1.025<\lambda .\hfill \end{array}\hfill & (11)\hfill \end{array}$

FIGURE 2

Figure 2. Along-fiber stress-stretch curves used by the model. (A) Aponeurosis (Equation 11). (B) Passive muscle (Equation 13). (C) Active muscle (Equation 16).

Muscle fibers can be passively or actively stretched, and we can therefore decompose σ_Muscle(λ) = σ_active(λ) + σ_passive(λ). We can also express the active and passive stresses normalized by the maximum isometric stress σ₀ = 200 KPa at the optimal length of the fiber. Experimental data is commonly given in terms of normalized stresses,

$\begin{array}{ll}{\sigma }_{muscle}(\lambda )={\sigma }_0\left({\widehat{\sigma }}_{Active}(\lambda )+{\widehat{\sigma }}_{Passive}(\lambda )\right)\hfill & (12\hfill \end{array})$

The normalized passive stress-stretch relation, $\widehat{\sigma }$_Passive(λ), was based on a piecewise exponential fit (Figure 2B) of the experimental data from Zajac (1989).

$\begin{array}{ll}{\widehat{\sigma }}_{Passive}(\lambda )=\{\begin{array}{l}0\lambda \le 1.0,\hfill \\ {10}^{-5}(38.495e^{5.339\lambda }-7945)1.0<\lambda .\hfill \end{array}\hfill & (13)\hfill \end{array}$

To obtain the (normalized) stress-stretch relation $\widehat{\sigma }$_Active(λ) for active muscle, we used a Hill-type model (Hill, 1938) of muscle fiber contraction. If we denote the normalized activation level by α(t) which ranges from 0 to 1, and the normalized stress-stretch rate by ${\widehat{\sigma }}_{\dot{\lambda }}$, then the (normalized) stress-stretch relation for the active fiber is given by

$\begin{array}{ll}{\widehat{\sigma }}_{Active}(\lambda )=\alpha (t){\widehat{\sigma }}_{\lambda }{\widehat{\sigma }}_{\dot{\lambda }}\hfill & (14)\hfill \end{array}$

In fact, in order to define different activation levels in different regions of the muscle, we modify (14) to make the dependence on location explicit:

$\begin{array}{ll}{\widehat{\sigma }}_{Active}(X,\lambda )=\alpha (X,t){\widehat{\sigma }}_{\lambda }{\widehat{\sigma }}_{\dot{\lambda }}\hfill & (15)\hfill \end{array}$

In the simplest instance, we can set α(X, t) = 0 for inactive regions, and linearly increasing over the simulation time from zero to a selected maximum level of activity α_max in the remaining. However, this may result in abrupt transitions in activity within the muscle, which may not be physiological. Instead, we used combinations of arctan (X) functions to vary activity smoothly between regions which were active and those which were not.

In this study, we allowed α(X, t) to vary linearly from 0 to α_max of 0.2 during 0.4 s of simulation. In order to describe $\widehat{\sigma }$_λ, we used experimental data from Gordon et al. (1966). In this experiment, the stresses were measured for fully active and isometric muscle fibers with isometric contractions where ${\widehat{\sigma }}_{\dot{\lambda }}$ = 1 at steady state. We fit the experimental data with a trigonometric polynomial (Figure 2C) in order to capture the complexity of the data while allowing for smooth derivatives.

$\begin{array}{ll}\begin{array}{l}{\widehat{\sigma }}_{\lambda }=0.534+0.229\cos (\omega \lambda )-0.095\cos (2\omega \lambda )\\ +0.024\cos (3\omega \lambda )-0.021\cos (4\omega \lambda )+0.013\cos (5\omega \lambda )\\ -0.421\sin (\omega \lambda )+0.079\sin (2\omega \lambda )-0.029\sin (3\omega \lambda )\\ +0.013\sin (4\omega \lambda )+0.002\sin (5\omega \lambda ),\end{array}\hfill & (16)\hfill \end{array}$

with ω = 4.957.

2.3. Numerical Methods

A three-field finite element formulation was used to solve for the displacement u, internal pressure $\tilde{p}:=\frac{\delta W_{vol}(\tilde{J})}{\tilde{J}}$ and the dilation constraint $\tilde{J}$ = J(u). We obtained these quantities by using the principle of stationary potential energy, i.e., minima of

$\begin{array}{ll}\begin{array}{l}U(u,\tilde{J},\tilde{p})={\displaystyle {\int }_{\Omega }{W}_{vol}}+\tilde{p}(J(u)-\tilde{J})dv+{\displaystyle {\int }_{\Omega }{W}_{iso}}\overline{B}(u)dv\\ -{\displaystyle {\int }_{\Omega }{\text{f}}_b}·udv-{\displaystyle {\int }_{\partial \Omega }{\text{f}}_t}·uda\end{array}\hfill & (17)\hfill \end{array}$

Here Ω, ∂Ω, v and a are the system's domain, boundary, volume and boundary area respectively and f_b, f_t are body and traction forces acting on the domain and boundary of the system respectively. In this study the applied body force f_b ≡ 0.

We used a discontinuous Galerkin method for u, $\tilde{p}$ and $\tilde{J}$. The resultant non-linear system was solved using Newton-Raphson iterations, and the linear solves within each Newton step were performed using a preconditioned conjugate gradient method. The discretization and solution was performed within the deal.II finite element library (Bangerth et al., 2013), and is a modification of a code (http://www.dealii.org/developer/doxygen/deal.II/step_44.html) by Pelteret and McBride to compute the material response of a neo-Hookean material.

2.4. Numerical Simulations

A set of simulations were designed to investigate the effect of regionalized activation, as well as different aponeurosis stiffnesses, on the magnitude and the direction of force developed by an isometrically activated muscle. The average activation of the whole muscle tissue was set to be 10% but the distribution of activation was changed between the simulations. For initial undeformed (relaxing) state, the muscle fibers were considered to be at optimal length and along-fiber strain in aponeurosis was set to zero.

The different distributions of activation are shown in Table 1. All activation distribution scenarios were repeated for three different elasticity moduli for aponeurosis. The aponeurosis was considered with maximum strains of 2, 5 and 10 % when the muscle was developing maximum isometric force, where the 5% case is given in Equation 11.

TABLE 1

Table 1. Level and regionalization of activation in different simulations.

The model was run on an eight-core machine with multi-threading over the cores (16 threads). The average CPU time for each simulation was approximately 10 min. This included the time needed to initialize the mesh, assemble matrices and iteratively solve the system.

3. Results

During the isometric contractions simulated in this study, the aponeuroses stretched, allowing the muscle fibers to shorten and rotate to greater pennation angles (Figure 3) than the initial value of 15°. A common end-point for the contractions was defined as the time where there was a mean 10% activation across the muscle tissue. The end-point of a contraction can be seen in Figure 4A for the condition with a uniform activation and compliant aponeurosis; this figure shows the total strain for each element in the tissues. However, note that the maximum along-fiber strain in the muscle was 26%. For this condition the aponeurosis stretched up to 1.5%, and the greatest shortening of the muscle fibers occurred in the center of the muscle belly. The muscle belly bulged in its width (x-direction) by approximately 12%, and decreased in the thickness between the aponeuroses. The muscle fibers curved during contraction, with the greatest curvatures occurring close to the aponeusoses, and the fibers following S-shaped paths. The initial condition had the fibers arranged in plane (parallel to the yz plane), and these curved outwards as the muscle belly width increased during contraction.

FIGURE 3

Figure 3. Deformed (active) and undeformed (relaxed) geometries for (A) the Uniform activation pattern and (B) the proximal-distal activation pattern. These geometries are shown with pale areas and blue lines for the undeformed states, and darker areas and gray lines for the deformed states. Note that in the deformed states the pennation angle for the proximal-distal activation pattern (18.37°) is larger than for the Uniform activation pattern. Transverse sections through the muscles are shown for the (C) Midline activation pattern, and (D) Medial-lateral activation pattern. In these panels the undeformed shape is shown by the rectangular and dark red area. The colored elements show the magnitude of the strain in the model tissues in their deformed state, ranging from low strains (blue) in the aponeurosis to greatest strains (red) in the muscle belly. Note how the muscle belly thickness between the aponeuroses is least over the active region of fibers, and the width of the muscles has increased beyond the undeformed state. Also note that in the Medial-lateral activation pattern the maximum strains have moved laterally (to the left) within the muscle.

FIGURE 4

Figure 4. Simulation results for the uniform activation condition with compliant aponeurosis and 10% activation. (A) Magnitude of strain. (B) “xz” and “yz” shear, and “zz” tensile stress contours on the plane connecting the aponeurosis and tendon (z = 0). (C) The direction cosines (dark gray) and force magnitude for the resultant force (light gray) acting on the z = 0 plane.

The results comparing all 12 simulations can be seen in Figure 5. The force vectors for the whole muscle were calculated from the shear and tensile stresses developed across the transverse plane bounding the deep aponeurosis (z = 0, for contour details see Figure 4B). The force vectors were described by the x-, y- and z- direction cosines of the force vector (δx, δy, and δz, respectively), and the resultant force magnitude (for details see Figure 4C).

FIGURE 5

Figure 5. Stress contours and force magnitudes and directions for the 12 test conditions. The scales are shown in Figures 4B,C.

In general, an increase in aponeurosis stiffness caused an increase in the magnitude of force and a change to its direction (see δy in Figure 5). The stretch in the aponeurosis was reduced for increased aponeurosis stiffness, and this led to a reduction in the shortening of the muscle fibers and a reduction in their rotation to higher pennation angles. Additionally, as the aponeurosis stiffness increased, the changes in muscle belly thickness and width became smaller. For the example of the uniformly activated muscle with a stiff aponeurosis, the width increased by only 7% during contraction.

The conditions with uniform activation had each muscle fiber activated to 10%. For the other 4 conditions with regionalization of the activity, the mean activity level across the muscle was kept at 10%, but this was concentrated in half the fibers each being activated to 20%. The magnitude of the muscle force was similar for the simulations with uniform, midline and medial-lateral distributions of activity, however the proximal-distal activation pattern resulted in greater muscle force.

The active fibers in the conditions with heterogeneous activation patterns (proximal-distal, midline and medial-lateral) contracted to a shorter length and rotated to a greater pennation angle than the uniform pattern. Additionally, in the conditions with proximal-distal activation patterns the thickness of the belly changed non-uniformly along the length of the muscle with a greater reduction in thickness in the active region.

For the compliant conditions, the y-component of the center of force moved from a position midway down the aponeurosis to a level closer to the deep surface for the stiff aponeurosis condition (Table 2). Conditions with the uniform, proximal-distal and midline activation patterns are all symmetrical about the midplane of the muscle (x = 27.5 mm). For these conditions there was a negligible x-component to the whole-muscle force, with δx < −10⁻⁶. The medial-lateral condition has the region of activity displaced to one side (Table 2), with the activity centered about the plane x = 32.1 mm. The x-component for the whole-muscle force was increased for this condition (δx ≅ −3 × 10⁻³): this value was still small, however, there was a more substantial increase in the x-component of the center of force (Table 2) acting at the end of the aponeurosis (to x = 31.3 mm).

TABLE 2

Table 2. x and y components of the center of Force (COF) on the z = 0 plane.

4. Discussion

The in silica isometric contractions of the lateral gastrocnemius in this study show decreases in muscle fiber length and increases in pennation (Figures 4A,B) in a similar manner to in vivo reports (Maganaris et al., 1998). However, the in silica muscle belly thickness decreased during contraction in a manner more representative of in vivo measurements from the medial gastrocnemius (Maganaris et al., 1998). The simulated muscle belly thickness decreased because a component of the contractile force of the muscle fibers acts to compress the muscle between the aponeuroses, and this was balanced by increases in the width of the muscle belly (Figures 4C,D) to maintain the incompressibility of the muscle tissue. Compressing the muscle belly would cause increases in intramuscular pressure, and this in turn drives increases in the curvature of the muscle fibers (Van Leeuwen and Spoor, 1992). Indeed, in our simulations the initial configurations of the fibers were straight, and this changed to curved S-shapes during contraction. These emergent features of the muscle fibers are consistent with recent reports of S-shaped fascicle trajectories that develop in vivo during isometric contractions of the medial gastrocnemius (Namburete et al., 2011). The width of the simulated muscle belly increased during contractions more than the width of the aponeuroses (Figures 4C,D), and so the belly bulged outwards to the sides. A consequence of this bulging was that the planes across which the fibers were initially aligned (parallel to the yz plane) transformed to curved sheets during contraction. Sejerested and co-workers (Sejersted et al., 1984) observed curved fascicle sheets in the vastus lateralis of human cadavers and have suggested the arrangement of fascicle sheets in an onion-like arrangement, and recent 3D reconstructions of fiber curvatures have also demonstrated that the fiber trajectories in both the medial and lateral gastrocnemius curve around concentric sheets during contraction (Rana, 2012). Despite the simplistic initial geometry (Figure 1) of the muscle in these simulations with straight fibers and rectangular aponeuroses, an emergent feature from the simulations was for the fibers to develop curved trajectories in 3D, and this has been suggested to be an important property to maintain mechanical stability within the muscle (Van Leeuwen and Spoor, 1992).

All the in silica conditions in this study had an equivalent level of muscle activity with an average of 10% activity across the muscle fibers (see Table 1) in the contracting state. However, both the magnitude and direction of the force (Figure 5) varied with the different regionalizations of the muscle activity. All conditions with regionalized activity resulted in greater muscle forces than for the condition with uniform activity across all fibers, and they also showed localized decreases in thickness, fascicle length and larger pennations around their active regions. These results are consistent with findings by Otten and Hulliger (1994) who compared the level of force in a uniformly fully active muscle with muscles where half the fibers were fully active: they activated the fibers in one end of the muscle and noted that the muscle force was 14% greater than expected from the average activity level. There was a pronounced increase in the muscle force for the non-uniform pattern of activation in the proximal-distal model (Figure 5); this could be due to specific changes in the shape of the muscle belly compared to other activity patterns, where the active region is decreasing in thickness more than the inactive region, and the muscle is bending around the region where the activity levels transitioned. We additionally show that changing the regionalization of the activity results to changes in the direction of the force vector from the whole muscle (Table 2), and this would result in different directions of the muscle torque in the joints that it spans, explaining experimental observations for the lateral gastrocnemius in the cat (Carrasco et al., 1999).

The increases in aponeurosis stiffness in these in silica contractions resulted in decreased aponeurosis stretch and muscle fiber strain, and smaller increases in pennation of the fibers. These changes resulted to increases in the magnitude of the force (Figure 5) developed by the whole muscle. These results are consistent with previous finite element models of the biceps femoris in which aponeurosis dimensions were changed with conditions that had stiffer aponeuroses, resulting in reduced fiber strains (Rehorn and Blemker, 2010). In another study, a finite element model of the medial gastrocnemius from the cat demonstrated that increased aponeurosis stiffness resulted in greater muscle forces for isometric contractions (Lemos et al., 2008), again consistent with our results. The combined findings from these three modeling studies show that the interactions between the muscle fibers and the connective tissue are important for shaping the mechanical output from the whole muscle. Our simulations additionally demonstrate how the direction of the muscle force is affected by the stiffness of the aponeurosis (Figure 5), presumably due to the changes in the shape of the muscle belly and force transmission that are caused by the different aponeurosis properties.

In our simulations the aponeurosis was included in an unrealistically thick (3 mm) state for computational simplicity, however the material properties of the aponeurosis was scaled so that its overall stiffness matched that expected for the in vivo condition. However, the thickness in the aponeurosis in our simulations allowed us to observe gradients of stress in the thickness direction (y-direction) that change with the aponeurosis stiffness. The simulations showed a reduced muscle fiber strain near the ends of the muscle where the fibers are structurally close to the fixed-end boundary condition and pull against the stiff aponeurosis. It is possible that removing the fixed boundary constraints, for instance by including the external tendon, will reduce this effect in future studies. Our in silica results show that the stress in the aponeurosis increased for the stiffest conditions, with the highest stress being more concentrated toward the outer layers of the aponeurosis (Figure 5), and it is possible that this indicates increases in the risk of injury initiation at areas closer to the outer surface of the aponeurotic sheets.

The results from this study highlight that the mechanical output of a whole muscle should not simply be considered to be a scaled-up muscle fiber that matches the size of the whole muscle (Wakeling and Lee, 2011), or a simple sum of actions of all the individual muscle fibers (Hill, 1970), but instead depends on the complex interactions between the muscle fibers and connective tissues that is brought about by the 3D structure of the muscle. In particular we show how the effect of regionalizing the muscle activity to a particular volume of muscle fibers causes changes to both the magnitude and direction of the whole muscle force, even when the mean level of muscle activity remains unchanged.

In summary, our simulations indicate that muscles with stiffer aponeuroses would result in smaller aponeurosis stretches and muscle fiber shortening. This effect would place the muscle fibers at a longer length on the ascending limb of their force-length curves, allowing them to develop greater stress and force. Additionally, as the stretch in the aponeurosis is reduced, the muscle fibers did not increase in their pennation angle as much during contraction. The simulations of regionalized and non-uniform activation patterns caused local differences in the shape of the muscle belly, strains and orientations of the muscle fibers. These factors affect both the magnitude and direction of the resultant muscle force.

Funding

We gratefully acknowledge funding from Natural Sciences and Engineering Research Council of Canada (Nilima Nigam and James M. Wakeling) and the Canada Research Chairs Program (Nilima Nigam).

Conflict of Interest Statement

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.

References

Arampatzis, A., Karamanidis, K., Morey-Klapsing, G., De Monte, G., and Stafilidis, S. (2007). Mechanical properties of the triceps surae tendon and aponeurosis in relation to intensity of sport activity. J. Biomech. 40, 1946–1952. doi: 10.1016/j.jbiomech.2006.09.005

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Azizi, E., Halenda, G. M., and Roberts, T. J. (2009). Mechanical properties of the gastrocnemius aponeurosis in wild turkeys. Integr. Comp. Biol. 49, 51–58. doi: 10.1093/icb/icp006

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bangerth, W., Heister, T., Heltai, L., Kanschat, G., Kronbichler, M., Maier, M., et al. (2013). The Deal.ii Library, Version 8.1. arXiv preprint: http://arxiv.org/abs/1312.2266v4.

Baskin, R. J., and Paolini, P. J. (1967). Volume change and pressure development in muscle during contraction. Am. J. Physiol. 213, 1025–1030.

Pubmed Abstract | Pubmed Full Text

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 | Pubmed Full Text | CrossRef Full Text

Bol, M., and Reese, S. (2008). Micromechanical modelling of skeletal muscles based on the finite element method. Comp. Methods Biomech. Biomed. Eng. 11, 489–504. doi: 10.1080/10255840701771750

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Carrasco, D. I., Lawrence, J., and English, A. W. (1999). Neuromuscular compartments of cat lateral gastrocnemius produce different torques about the ankle joint. Motor control 3, 436–446.

Pubmed Abstract | Pubmed Full Text

Chanaud, C., Pratt, C., and Loeb, G. (1990). Functionally complex muscles of the cat hindlimb. ii. mechanical and architectural heterogenity within the biceps femoris. Exp. Brain Res. 85, 257–270. doi: 10.1007/BF00229405

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Delp, S. L., Loan, J. P., Hoy, M. G., Zajac, F. E., Topp, E. L., and Rosen, J. M. (1990). An interactive graphics-based model of the lower extremity to study orthopaedic surgical procedures. IEEE Trans. Biomed. Eng. 37, 757–767. doi: 10.1109/10.102791

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

English, A. W. (1984). An electromyographic analysis of compartments in cat lateral gastrocnemius muscle during unrestrained locomotion. J. Neurophysiol. 52, 114–125.

Pubmed Abstract | Pubmed Full Text

English, A. W., and Weeks, O. I. (1987). An anatomical and functional analysis of cat biceps femoris and semitendinosus muscles. J. Morphol. 191, 161–175. doi: 10.1002/jmor.1051910207

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

English, A. W., Wolf, S. L., and Segal, R. L. (1993). Compartmentalization of muscles and their motor nuclei: the partitioning hypothesis. Phys. Ther. 73, 857–867.

Pubmed Abstract | Pubmed Full Text

Gordon, A., Huxley, A. F., and Julian, F. (1966). The variation in isometric tension with sarcomere length in vertebrate muscle fibres. J. Physiol. 184, 170–192.

Pubmed Abstract | Pubmed Full Text

Herring, S., Anapol, F., and Wineski, L. (1991). Motor-unit territories in the masseter muscle of infant pigs. Arch. Oral Biol. 36, 867–873. doi: 10.1016/0003-9969(91)90116-C

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Herring, S. W., Anapol, F. C., and Wineski, L. E. (1989). Neural organization of the masseter muscle in the pig. J. Comp. Neurol. 280, 563–576. doi: 10.1002/cne.902800407

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hill, A. V. (1938). The heat of shortening and the dynamic constants of muscle. Proc. R. Soc. Lond. B Biol. Sci. 126, 136–195. doi: 10.1152/advan.00072.2005

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hill, A. V. (1970). First and Last Experiments in Muscle Mechanics. Cambridge: Cambridge University Press.

Hodson-Tole, E. F., Loram, I. D., and Vieira, T. M. (2013). Myoelectric activity along human gastrocnemius medialis: Different spatial distributions of postural and electrically elicited surface potentials. J. Electromyogr. Kinesiol. 23, 43–50. doi: 10.1016/j.jelekin.2012.08.003

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Huijing, P. A. (2003). Muscular force transmission necessitates a multilevel integrative approach to the analysis of function of skeletal muscle. Exerc. Sport Sci. Rev. 31, 167–175. doi: 10.1097/00003677-200310000-00003

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Humphrey, J., and Yin, F. (1987). On constitutive relations and finite deformations of passive cardiac tissue: I. a pseudostrain-energy function. J. Biomech. Eng. 109, 298–304. doi: 10.1115/1.3138684

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kubo, K., Kanehisa, H., and Fukunaga, T. (2002). Effects of resistance and stretching training programmes on the viscoelastic properties of human tendon structures in vivo. J. Physiol. 538, 219–226. doi: 10.1113/jphysiol.2001.012703

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lake, S. P., Miller, K. S., Elliott, D. M., and Soslowsky, L. J. (2010). Tensile properties and fiber alignment of human supraspinatus tendon in the transverse direction demonstrate inhomogeneity, nonlinearity, and regional isotropy. J. Biomecha. 43, 727–732. doi: 10.1016/j.jbiomech.2009.10.017

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lemos, R. R., Epstein, M., and Herzog, W. (2008). Modeling of skeletal muscle: the influence of tendon and aponeuroses compliance on the force–length relationship. Med. Biol. Eng. Comput. 46, 23–32. doi: 10.1007/s11517-007-0259-x

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Maganaris, C. N., Baltzopoulos, V., and Sargeant, A. J. (1998). in vivo measurements of the triceps surae complex architecture in man: implications for muscle function. J. Physiol. 512(Pt 2), 603–614. doi: 10.1111/j.1469-7793.1998.603be.x

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Maganaris, C. N., and Paul, J. P. (2000). in vivo human tendinous tissue stretch upon maximum muscle force generation. J. Biomech. 33, 1453–1459. doi: 10.1016/S0021-9290(00)00099-3

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Magnusson, S. P., Hansen, P., Aagaard, P., Brnd, J., Dyhre-Poulsen, P., Bojsen-Moller, J., et al. (2003). Differential strain patterns of the human gastrocnemius aponeurosis and free tendon, in vivo. Acta Physiol. Scand. 177, 185–195. doi: 10.1046/j.1365-201X.2003.01048.x

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

MATLAB. (2014). Version 8.3.0 (R2014a). Natick, MA: The MathWorks Inc.

Namburete, A. I., Rana, M., and Wakeling, J. M. (2011). Computational methods for quantifying in vivo muscle fascicle curvature from ultrasound images. J. Biomech. 44, 2538–2543. doi: 10.1016/j.jbiomech.2011.07.017

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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 | Pubmed Full Text | CrossRef Full Text

Oomens, C. W., Maenhout, M., van Oijen, C. H., Drost, M. R., and Baaijens, F. P. (2003). Finite element modelling of contracting skeletal muscle. Philos. Trans. R. Soc. Lond. B Biol. Sci. 358, 1453–1460. doi: 10.1098/rstb.2003.1345

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Otten, E., and Hulliger, M. (1994). A finite-elements approach to the study of functional architecture in skeletal-muscle. Zool. Anal. Comp. Syst. 98, 233–242.

Rana, M. (2012). 3D Muscle Architecture in the Triceps Surae Muscle: 3D Ultrasound Methods and Maps of Fascicle Orientation and Curvature. PhD thesis, Science: Department of Biomedical Physiology and Kinesiology. Burnaby, BC, Simon Fraser University.

Randhawa, A., Jackman, M. E., and Wakeling, J. M. (2013). Muscle gearing during isotonic and isokinetic movements in the ankle plantarflexors. Eur. J. Appl. Physiol. 113, 437–447. doi: 10.1007/s00421-012-2448-z

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Rehorn, M. R., and Blemker, S. S. (2010). The effects of aponeurosis geometry on strain injury susceptibility explored with a 3d muscle model. J. Biomech. 43, 2574–2581. doi: 10.1016/j.jbiomech.2010.05.011

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Sejersted, O., Hargens, A. R., Kardel, K. R., Blom, P., Jensen, O., and Hermansen, L. (1984). Intramuscular fluid pressure during isometric contraction of human skeletal muscle. J. Appl. Physiol. 56, 287–295.

Pubmed Abstract | Pubmed Full Text

Spencer, A. J. M. (1984). Continuum Theory of the Mechanics of Fibre-Reinforced Composites. New York, NY: Springer-Verlag. doi: 10.1007/978-3-7091-4336-0

CrossRef Full Text

Van Leeuwen, J. L., and Spoor, C. W. (1992). Modelling mechanically stable muscle architectures. Philos. Trans. R. Soc. Lond. B Biol. Sci. 336, 275–292. doi: 10.1098/rstb.1992.0061

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wakeling, J. M. (2009). Patterns of motor recruitment can be determined using surface emg. J. Electromyogr. Kinesiol. 19, 199–207. doi: 10.1016/j.jelekin.2007.09.006

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wakeling, J. M., and Lee, S. S. (2011). Modelling muscle forces: from scaled fibres to physiological task-groups. Proc. IUTAM 2, 317–326. doi: 10.1016/j.piutam.2011.04.028

CrossRef Full Text

Wren, T. A., Yerby, S. A., Beaupre, G. S., and Carter, D. R. (2001). Mechanical properties of the human achilles tendon. Clin. Biomech. 16, 245–251. doi: 10.1016/S0268-0033(00)00089-9

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Yeoh, O. (1993). Some forms of the strain energy function for rubber. Rubb. Chem. Technol. 66, 754–771. doi: 10.5254/1.3538343

CrossRef Full Text

Zajac, F. E. (1989). Muscle and tendon - properties, models, scaling, and application to biomechanics and motor control. Crit. Rev. Biomed. Eng. 17, 359–411.

Pubmed Abstract | Pubmed Full Text