# A multiscale chemo-electro-mechanical skeletal muscle model to analyze muscle contraction and force generation for different muscle fiber arrangements

^{1}Continuum Biomechanics and Mechanobiology Research Group, Institute of Applied Mechanics (CE), University of Stuttgart, Stuttgart, Germany^{2}Stuttgart Research Center for Simulation Technology (SimTech), University of Stuttgart, Stuttgart, Germany

The presented chemo-electro-mechanical skeletal muscle model relies on a continuum-mechanical formulation describing the muscle's deformation and force generation on the macroscopic muscle level. Unlike other three-dimensional models, the description of the activation-induced behavior of the mechanical model is entirely based on chemo-electro-mechanical principles on the microscopic sarcomere level. Yet, the multiscale model reproduces key characteristics of skeletal muscles such as experimental force-length and force-velocity data on the macroscopic whole muscle level. The paper presents the methodological approaches required to obtain such a multiscale model, and demonstrates the feasibility of using such a model to analyze differences in the mechanical behavior of parallel-fibered muscles, in which the muscle fibers either span the entire length of the fascicles or terminate intrafascicularly. The presented results reveal that muscles, in which the fibers span the entire length of the fascicles, show lower peak forces, more dispersed twitches and fusion of twitches at lower stimulation frequencies. In detail, the model predicted twitch rise times of 38.2 and 17.2 ms for a 12 cm long muscle, in which the fibers span the entire length of the fascicles and with twelve fiber compartments in series, respectively. Further, the twelve-compartment model predicted peak twitch forces that were 19% higher than in the single-compartment model. The analysis of sarcomere lengths during fixed-end single twitch contractions at optimal length predicts rather small sarcomere length changes. The observed lengths range from 75 to 111% of the optimal sarcomere length, which corresponds to a region with maximum filament overlap. This result suggests that stability issues resulting from activation-induced stretches of non-activated sarcomeres are unlikely in muscles with passive forces appearing at short muscle length.

## 1. Introduction

The fascicles in parallel-fibered muscle are aligned with the muscle's line of action and run almost the entire length of the muscle (Loeb et al., 1987). The fascicles either consist of long fibers spanning the entire length of the fascicles (in the following termed “spanning-fibered muscle”), or they are composed of several shorter in-series arranged fiber compartments (in the following termed “series-fibered muscle”) (Richmond et al., 1985; Heron and Richmond, 1993; Young et al., 2000). The fiber compartments in series-fibered muscle can either be separated by tendinous inscriptions, as, for example, in cat and human semitendinosus muscle, or the muscle fibers are arranged in short overlapping arrays (Loeb et al., 1987; Paul, 2001; Woodley and Mercer, 2005).

The advantages and disadvantages of series-fibered and spanning-fibered muscle arrangements on the force generation have not yet been systematically analyzed. Experiments provide only limited information on which effects are due to the fiber arrangement, and which effects are due to other anatomical or physiological properties, e. g., the muscle geometry. Mathematical models instead can be used to investigate the influence of a specific property on the overall behavior. Previous modeling works focused on the influences of the muscle geometry and the fiber direction on the force generation (Zuurbier and Huijing, 1992; Sánchez et al., 2014). To investigate the effect of different fiber arrangements, one requires a model that unifies the following features: (i) The dynamics of the active force generation are determined at discrete locations (“sarcomeres”) along a muscle fiber. (ii) The model takes into account the subsequent activation of adjacent “sarcomeres” through the propagation of action potentials (APs) along the fibers. This is required since the AP propagation speed is rather slow, and hence, sarcomere activation is non-synchronized along a muscle fiber, and the asynchronism increases with increasing fiber length. (iii) The model accounts for the muscle tissue, in which the muscle fibers are embedded, and which shows resistance to applied loads. The tissue representation is required due to the fact that isolated muscle fibers or myofibrils do not behave like fibers within a muscle (Prado et al., 2005).

There exists no muscle model that can incorporate all of these requirements. Hill-type muscle models are typically used to describe whole muscle behavior (Zajac, 1989), although they have also been used to model single sarcomeres, and, by in-series arranging multiple Hill-type models, myofibrils and muscle fiber segments have been modeled (Morgan et al., 1982; Stoecker et al., 2009; Günther et al., 2012). While these approaches can describe local changes in sarcomere length, they cannot capture the behavior of a fiber within the three-dimensional (3D) muscle tissue. This is due to the fact that the passive forces in isolated myofibrils and single muscle fibers are mainly attributed to the titin filament (Horowits, 1992; Denoth et al., 2002), unlike in muscle tissue, where the extracellular matrix contributes additional passive forces (Prado et al., 2005). Further, in isolated myofibrils and single muscle fibers, force transmission can only take place along their length. In whole muscle, however, force transmission also occurs in lateral direction (Huijing, 1999). For an adequate representation of fibers within the muscle tissue and its mechanical implications on the behavior of the whole muscle, a 3D model based on continuum-mechanical principles is required.

Previous continuum-mechanical skeletal muscle models (Blemker et al., 2005; Röhrle and Pullan, 2007) include the active force-length (*F*-*ℓ*) and/or the active force-velocity (*F*-*v*) relations on the macroscale, which implies the assumption of an averaged sarcomere length and an averaged sarcomere shortening velocity. Therefore, these modeling approaches cannot represent local changes in sarcomere length and shortening velocity, which are required for above motivated cases. Furthermore, both the *F*-*ℓ* and the *F*-*v* relations should be modeled on the microscale, since they can be attributed to properties on the sarcomere level. The length dependence of the active force is due to changes in the overlap of the thick and thin filaments within the sarcomeres (Gordon et al., 1966), while the velocity dependence is attributed to (i) a lower tension of the cross bridges (XBs) that reattach in a shortened state, and (ii) an increased XB-detachment rate (Piazzesi et al., 2007; Telley and Denoth, 2007).

To overcome the limiting modeling assumption of homogenized sarcomere lengths and shortening velocities, and to analyze the effects of different fiber arrangements, in this contribution, the multiscale chemo-electro-mechanical skeletal muscle model of Heidlauf and Röhrle (2013) is extended to include the *F*-*ℓ* and *F*-*v* relations on the microscale.

## 2. Materials and Methods

To model the active *F*-*ℓ* and *F*-*v* relations on the microscale, the biophysical half-sarcomere model of Shorten et al. (2007) is extended to non-isometric conditions. The extended half-sarcomere model is coupled to (i) bioelectrical field equations describing the propagation of APs along muscle fibers, and (ii) a 3D continuum-mechanical description of the muscle tissue (Heidlauf and Röhrle, 2013).

### 2.1. Detailed Overview of the Multiscale Model

Depolarization of the membrane potential of the biophysical half-sarcomere model located at the innervation zone is induced by a current injection of short duration. The timing of the current injections is given by the stimulation frequency, which is prescribed in this study (e. g., 50 Hz or 100 Hz). The constant firing frequency can also be replaced by discrete motor unit discharge times resulting, for example, from the decomposition of an EMG signal (De Luca and Hostage, 2010) or from a phenomenological (Fuglevand et al., 1993) or biophysical (Heidlauf and Röhrle, 2013) model of the α motor neurons. Based on the respective stimulation, the biophysical half-sarcomere model provides, among many others, two quantities that are essential for the multiscale framework—the locally generated active stresses and the changes in membrane potential due to ionic and capacitive currents. To simulate the propagation of APs, the bioelectrical field equations are used to describe the diffusion of the membrane potential along the fibers. This results in a bi-directional coupling between the half-sarcomere model and the bioelectrical field equations through the membrane potential. The locally-generated, sarcomere-based active stresses are included in the formulation of the continuum-mechanical constitutive relation (relation between local deformation and resulting local stresses). The continuum-mechanical model predicts the deformation of the muscle geometry, the internal stress and strain distributions, and the forces that can be passed to adjacent structures such as tendon. The local strain is used to determine the new sarcomere length and the sarcomere shortening velocity, which are in turn inputs to the biophysical half-sarcomere model. Hence at a point in space, the half-sarcomere model and the continuum-mechanical model are bi-directionally coupled. Furthermore, since deformation changes geometrical properties of the fibres, the AP propagation along a muscle fiber is solved on the deformed geometry.

Due to the complexity of the model, Table 1 lists the model's variables including their dependencies, while Table 2 summarizes the parameters of the model.

### 2.2. The Continuum-Mechanical Muscle Model

Since the physiological working range of many muscles involves changes in length of 50% and more (Burkholder and Lieber, 2001), a continuum-mechanical analysis must be based on the finite elasticity theory (Holzapfel, 2000; Bonet and Wood, 2008). In continuum mechanics, the placement function χ assigns a material point with position **X** in the reference (undeformed) configuration at time *t*_{0} to a position in the actual (deformed) configuration **x** at time *t*, i. e., **x** = **χ**(**X**, *t*). The material deformation gradient tensor **F** is defined as the derivative of the placement function with respect to the material coordinates, i. e., ${F}{=}\frac{{\partial}{\chi}}{{\partial}{X}}{=}\frac{{\partial}{x}}{{\partial}{X}}$. Local deformations and strains are conveniently described by the right Cauchy-Green deformation tensor **C** = **F**^{T} **F** and the Green's strain tensor ${E}{=}\frac{{1}}{{2}}{(}{C}{-}{I}{)}$, respectively, where **I** is the second-order identity tensor.

Considering the stress equilibrium in the actual configuration and neglecting inertia and body forces, the momentum balance equation reduces to div **T** = **0**, where **T** denotes the Cauchy stress tensor. To characterize the material behavior, a constitutive equation is required that relates the local deformations or strains to the resulting local stresses. This is conveniently done in the reference configuration. The Cauchy stress tensor of the actual configuration is related to the second Piola-Kirchhoff stress tensor, **S**, of the reference configuration via a scaled covariant push forward operation: **T** = (det **F**)^{−1} **F S F**^{T}. Muscle tissue can actively generate tension and in the passive state, it exhibits transversal isotropic material behavior. This is reflected in **S**, which consists of an isotropic part based on the Mooney-Rivlin material, **S**^{iso}, a term appealing to stretches in the fiber direction, **S**^{ani} (cf. Markert et al., 2005), which together with **S**^{iso} characterizes the transversal isotropic passive behavior of muscle tissue, and a term representing the muscle's ability to actively generate tension, **S**^{act}. The form of **S** is derived in Heidlauf and Röhrle (2013), and is given by

where *p* is the hydrostatic pressure, *I*_{1} = tr **C** is the first principal invariant of **C**, and **a**_{0} is a unit vector in fiber direction defined in the reference configuration. Further, ${{\lambda}}_{{f}}{=}\sqrt{{{I}}_{{4}}}$ denotes the fiber stretch with *I*_{4} = **a**_{0} · **Ca**_{0} being the fourth (mixed) invariant of **C**.

While the fiber stretch is a (spatially varying) macroscopic quantity, it can be related to the corresponding quantity on the sarcomere level, i. e., the sarcomere length, ℓ_{S}, by λ_{f} = ℓ_{S}/ℓ^{0}_{S} with ℓ^{0}_{S} = 2.0 μm being the sarcomere resting length. Finally, *P ^{act}* represents a scalar-valued active nominal stress, which is the product of the maximum active stress at optimal fiber length and under isometric conditions,

*P*, and the normalized active stress

^{max}*γ*:

Therein, *γ* depends on the stimulation frequency *f _{s}*, the fiber length (represented through the fiber stretch λ

_{f}), and the contraction velocity, ${\dot{{\lambda}}}_{{f}}$. Note that previous models (Johansson et al., 2000; Röhrle et al., 2008; Heidlauf and Röhrle, 2013) employ the

*F*-

*ℓ*and/or

*F*-

*v*relations on the macroscopic continuum level in the form

*P*=

^{act}*P*

^{max}

*f*

_{1}(α)

*f*

_{2}(λ

_{f})

*f*

_{3}(${\dot{{\lambda}}}_{{f}}$) with α ϵ [0, 1] being an internal activation parameter. In contrast to these models, the present work provides novel contributions to the field of multiscale skeletal muscle modeling by determining

*γ*as part of the biophysical model on the microscale (see next section).

The macroscopic material parameters *c*_{10} and *c*_{01} in Equation (2) have been fitted in a least-squares sense to the uniaxial compression experiments of Zheng et al. (1999). Further, *b*_{1} and *d*_{1} have been determined similarly from the passive experimental data of Hawkins and Bey (1994), from which also the value of *P ^{max}* in Equation (2) is adopted. The parameters are summarized in Table 2.

### 2.3. The Biophysical Half-Sarcomere Model

The basis for modeling the subcellular level in this contribution is the model of Shorten et al. (2007), which describes the complex, nonlinear, biophysical processes leading from electrical excitation to contraction and force generation. To model the excitation-contraction coupling, Shorten et al. (2007) combined several component models describing (a) membrane electrophysiology, (b) calcium release from the sarcoplasmic reticulum and (c) calcium dynamics, (d) cross-bridge (XB) dynamics, and (e) fatigue. The model of Shorten et al. (2007) can be freely accessed and downloaded from the CellML website (http://www.cellml.org/).

The modifications of this contribution to the model of Shorten et al. (2007) are restricted to the eight-state XB-dynamics component model (d), which is based on the four-state XB-dynamics model of Razumova et al. (1999, 2000) and Campbell et al. (2001a,b). A schematic representation of the eight-state model is shown in Figure 1. In six of the eight states, the XBs are in a detached state with zero, one or two Ca^{2+} ions bound to troponin (denoted by indices 0, 1, and 2, respectively) and with tropomyosin in either a blocking (*B*) or non-blocking (*D*) position. Only in the case when two Ca^{2+} ions are bound to troponin and the tropomyosin block is in a non-blocking position (the *D*_{2} state), the detached XB can move to a state where the myosin head is attached. Two attached states are distinguished—the pre-power stroke state *A*_{1} and the post-power stroke state *A*_{2}. The transition from the *A*_{1} to the *A*_{2} state represents the power stroke, i. e., the force producing step, for which the forward and backward reaction rates, *h*_{0} and *h*′, respectively, apply. The forward and backward reaction rates changing the *D*_{2} to the *A*_{1} state and vice versa (XB attachment and detachment) are denoted by *f*_{0} and *f*′, respectively. Finally, the detachment of XBs from state *A*_{2} to state *D*_{2} is described by reaction rate *g*_{0}, see Figure 1. Shorten et al. (2007) provide a slow-twitch (type-I) and a fast-twitch (type-II) parametrization of the model to simulate isometric contractions of mouse soleus and EDL muscle, respectively. These parameter sets are adopted in the present work, i. e., no reparametrization is required.

In the present contribution, the model of Shorten et al. (2007) is extended to non-isometric conditions. This is done, first, by incorporating changes in the myofilament overlap, and further, by adding a distortion dependence and cooperative effects to the XB dynamics component model. These extensions are based on the works of Razumova et al. (1999) and Campbell et al. (2001b).

The force that can be exerted by a sarcomere depends on the number of XB connections between the actin and myosin filaments (Huxley, 1957). The number of possible XB connections depends on the filament overlap, and hence, on the sarcomere length (Gordon et al., 1966). Based on analytical considerations of the filament overlap, Campbell et al. (2001b) proposed a piecewise linear relation between the sarcomere half-length and the number of possible XB connections. The relation is depicted in Figure 2 (green dashed line) assuming a direct relation between the number of possible XB connections and the isometric active force at full activation. Experiments on single sarcomeres, however, suggest a steeper decline of the force on the ascending limb of the active force-length curve at sarcomere lengths below 1.7 μm, and no active force production at lengths below 1.27 μm (Gordon et al., 1966). This is attributed to an interaction of the myosin filament with the Z-disks at low sarcomere lengths. The red solid line in Figure 2 shows the experimentally determined relation between the sarcomere length and the isometric active stress at full activation. In the present work, a fourth-order polynomial is used, cf. Figure 2 (dot-dashed blue line):

where Ξ is the normalized isometric active force at full activation, and ℓ_{S} denotes the sarcomere length. The polynomial in (3) is symmetric with respect to the optimal sarcomere length ℓ^{opt}_{S} = 2.4 μm (Burkholder and Lieber, 2001), and can be seen as an approximation to the experimentally determined force-sarcomere length relation, where the largest deviations occur at very long sarcomere lengths. In this contribution, the behavior at very long sarcomere lengths, however, is dominated by the passive stiffness of the muscle tissue, and therefore, the implications of the deviations will be limited. Note that the fourth-order polynomial in (3) is a generic description of a muscle's *F*-*ℓ* behavior (cf. Zuurbier et al., 1995). This approximation can be easily replaced by a different *F*-*ℓ* curve that was fitted to experimental data of a specific muscle. Furthermore, the optimal sarcomere length, which is invariant for all presented simulations, can be changed to simulate a specific muscle.

**Figure 2. The relation between the normalized maximum isometric active stress and the sarcomere length**. Plotted is the experimentally determined force-sarcomere length relation for cat skeletal muscle (Rassier et al., 1999) (red solid line), the piecewise linear relation of Campbell et al. (2001b) (dashed green line), and the fourth-order polynomial of this work (dot-dashed blue line).

To account for length changes during contraction, average distortions (or elastic deformations) of XBs in a sarcomere are introduced into the XB-dynamics component model according to Campbell et al. (2001b). The average distortion induced by the power stroke during an isometric contraction of a half-sarcomere is denoted by *x*_{0}. The average elastic deformations among XBs in the pre-power stroke state *A*_{1} and post-power stroke state *A*_{2} induced through filament sliding during non-isometric contractions are denoted by *x*_{1} and *x*_{2}, respectively. Note that the term average refers in this context to the spatial average over all XBs of that sarcomere in the respective state. Figure 3 illustrates the different distortions. While *x*_{0} is assumed to be constant, *x*_{1} and *x*_{2} account for distortions entering and leaving due to XB cycling and for distortions imposed by shearing between thick and thin filaments. From the distortional balances, Campbell et al. (2001b) derived the following ODEs, which are included in the present model:

**Figure 3. (A)** The average distortion *x*_{0} induced through the power stroke in an isometric contraction. In the pre-power stroke state *A*_{1} the cross-bridge is attached to the myosin binding site (small filled circle) and does not experience an elastic distortion. The power stroke converts the *A*_{1} to the *A*_{2} state by transducing chemically stored energy into mechanical energy, which is stored in the elastically distorted cross-bridges. **(B)** Average distortion *x*_{1} induced through filament sliding during non-isometric contractions on the cross-bridges in the *A*_{1} state. **(C)** Average distortion *x*_{2} induced through filament sliding during non-isometric contractions on the cross-bridges in the *A*_{2} state.

Therein, ${\dot{{\ell}}}_{{S}}$ denotes the sarcomere contraction velocity. Further, quantities in square brackets denote concentrations of XBs in the respective state. The differential equations describing the concentrations of XBs in the different states are part of the biophysical half-sarcomere model of Shorten et al. (2007).

The force exerted by a half-sarcomere is proportional to the product of the stiffness of all parallel XBs and their average distortions (Razumova et al., 1999), i. e.,

The normalized sarcomere-based active stress γ is defined to be the product of the force-length relation and the ratio between *B* and the value of *B* at maximum stimulation *f ^{max}_{s}* and under isometric conditions ${\dot{{\ell}}}_{{S}}$, corrected for the value of

*B*at zero activation:

Before, however, including γ in the macroscopic continuum-mechanical constitutive equation (2), the normalized sarcomere-based active stress values are homogenized, cf. Section 2.5 and Heidlauf and Röhrle (2013) for details.

To reproduce the hyperbolic *F*-*v* relation (Hill, 1938), Razumova et al. (1999) proposed two modifications to their four-state XB-dynamics model: (i) The forward rate of XB attachment *f*_{0} contains now nearest-neighbor cooperative effects, i. e., increased XB-attachment probabilities due to neighboring XBs in the force-bearing state. (ii) A distortion dependence is now incorporated in the XB-detachment rate *g*_{0} accounting for an increasing probability of XB detachment with an increasing XB distortion:

where *g* is the XB-detachment rate of an isometric contraction, ϑ controls the distortion dependence, *f* denotes the forward rate of XB attachment if no neighbor is in the force-bearing state, and ν controls the influence of the cooperative effects. Further, *T _{tot}* is the total number of possible XB connections at optimal filament overlap. Equations (7) and (8) are added to the XB-dynamics component model within the model of Shorten et al. (2007). Note that in the method of Campbell et al. (2001b),

*T*depends on the sarcomere length. This approach, however, bears some problems, for example, in the case when sarcomeres are stretched to beyond myofilament overlap (

_{tot}*T*→ 0). Therefore, the present approach includes the

_{tot}*F*-

*ℓ*relation in Equation (6) in a form that is inspired by Hill-type models (Siebert et al., 2008). However, in contrast to Hill-type models that include the

*F*-

*ℓ*relation at the macroscopic whole muscle level, the present approach contains this relation at the microscopic sarcomere level.

### 2.4. Action Potential Propagation

Previous electro-mechanical muscle models (Fernandez et al., 2005; Böl et al., 2011) describe the AP propagation as a continuous 3D wave front moving through the entire muscle domain. However, the macroscopic electrical conductivity of skeletal muscle tissue perpendicular to the fiber direction is up to one magnitude lower than the conductivity along the fiber direction (Epstein and Foster, 1983; Gielen et al., 1984), and electrical stimulation from one fiber to adjacent ones is not observed. Therefore, the propagation of an AP along a skeletal muscle fiber is modeled as a one-dimensional (1D) problem (cf. Röhrle et al., 2008; Heidlauf and Röhrle, 2013). The AP propagation can be described by the monodomain equation, which is in 1D identical to the cable equation, see e. g., Hodgkin and Huxley (1952); Pullan et al. (2005):

Therein, *s* is the spatial variable describing the position along the path of the fiber, σ denotes the conductivity, *V _{m}* represents the membrane voltage,

*A*reflects the surface-area-to-volume ratio of the cell, and

_{m}*C*is the capacitance of the cell membrane per unit area. The monodomain equation links through the ionic currents crossing the cell membrane,

_{m}*I*, to the half-sarcomere model described in the previous section, i. e.,

_{ion}*I*=

_{ion}*I*(

_{ion}*t*,

*V*,

_{m}**y**) with

**y**denoting the state variables of the half-sarcomere model. The term on the left-hand side of Equation (9) describes the diffusion of membrane potential along a muscle fiber. For details, the reader is referred to Heidlauf and Röhrle (2013).

### 2.5. Computational Framework

Due to interactions between the half-sarcomere model, the AP propagation model, and the continuum-mechanical model, a fully coupled system needs to be solved in an integrated fashion. In this contribution, a staggered solution scheme is employed (Heidlauf and Röhrle, 2013), which allows usage of different solution methods and different time step sizes for the solution of the individual subsystems. Moreover, due to computational efficiency, an approach that uses different finite element discretizations for the 1D bioelectrical and the 3D continuum-mechanical subsystems is adopted. Since different meshes are used for different subsystems, transfer operations for sharing variables between different meshes are required. For example, the normalized active stress γ determined in the half-sarcomere models, cf. Equation (6), needs to be homogenized to the coarser 3D continuum-mechanical mesh (Γ: γ → *γ*) to be included in the evaluation of the stress tensor in Equations (1) and (2). Without loss of generality, a geometrically based homogenization is used in this contribution. The convergence behavior of this approach is investigated in Röhrle et al. (2008), and shows good results. Further details on the computational framework can be found in Heidlauf and Röhrle (2013) and Bradley et al. (2011).

## 3. Results

Before comparing series-fibered and spanning-fibered muscles, the behavior of the extended half-sarcomere model and the new fully coupled chemo-electro-mechanical skeletal muscle model is investigated.

### 3.1. Half-Sarcomere Model

To show that the extended half-sarcomere model (Shorten et al., 2007) exhibits a *F*-*v* relation as muscle fibers do, the sensitivity of the model to the newly introduced parameters ν and ϑ is analyzed first. To do so, experiments are carried out using a single extended half-sarcomere model at a stimulation frequency of *f _{s}* = 100 Hz. For different prescribed constant velocities, the corresponding normalized active stresses γ are computed at optimal sarcomere length.

The model predicts a linear *F*-*v* relation for constant rate coefficients *f*_{0} and *g*_{0} (ϑ = 0, ν = 1), cf. Figure 4. When considering nearest-neighbor cooperative effects in *f*_{0} (ϑ = 0, ν = 3.4), the model is able to predict a hyperbolic relation for shortening contractions, but unreasonable high forces occur for lengthening contractions. The distortion dependence in *g*_{0} (ν = 3.4, ϑ = 1000, 2000) mainly influences lengthening contractions. Note that Figure 4 only depicts results for the type-II parametrization (Shorten et al., 2007). Similar results are obtained for type-I fibers as demonstrated in Section 3.2.

**Figure 4. F-v relation of an isolated half-sarcomere model**. Shown is the relation for constant rate coefficients (red dash-dotted line), for variability in

*f*

_{0}only (blue dotted line), and for additional variability in

*g*

_{0}(blue dashed line and turquoise solid line). To depict the Hill relation in its typical form, the x-axis is inverted to show shortening contractions on the right.

Further, three shortening contractions are simulated to demonstrate the influence of the *F*-*ℓ* and the *F*-*v* relations on the active stress profiles. To this end, a single half-sarcomere model is stimulated at a frequency of 100 Hz. After 500 ms of isometric contraction at optimal sarcomere length, the sarcomere shortens at a constant prescribed velocity. Three different velocities are considered: 5, 10, and 15% of the maximum shortening velocity *v _{max}*.

Figure 5 shows the evolution of the normalized sarcomere-based active stresses (top) and the sarcomere length (ℓ_{S}, bottom). The profiles show, first, an increase in the active stress due to the stimulation, which is identical for all three traces. After 500 ms, when the active stress approximately saturates and the shortening starts, the model shows an instantaneous stress drop which is due to the shortening velocity. As expected, the magnitude of the stress drop increases with the shortening velocity, cf. Figure 4. The model further predicts a decrease in the stress, which is due to the *F*-*ℓ* relation, i. e., as the sarcomere shortens, it moves along the ascending limb of the *F*-*ℓ* relation (Figure 2) from the optimal length toward smaller sarcomere lengths.

**Figure 5. Evolution of the normalized sarcomere-based active stress of an isolated half-sarcomere model for three different shortening velocities at a stimulation frequency of 100 Hz (top)**. The shortening contraction is preceded by an isometric contraction at optimal length of 500 ms duration. Additionally, the actual sarcomere length (ℓ_{S}) is shown for each of the stress profiles (bottom).

### 3.2. The Chemo-Electro-Mechanical Model Compared to Experimental Data

The chemo-electro-mechanical model is first compared to experimental *F*-*ℓ* data to demonstrate that the multiscale muscle model that includes the entire active behavior on the microscale can reproduce typical mechanical behavior of whole muscle on the macroscale. For the comparison, the experimental *F*-*ℓ* data of Hawkins and Bey (1994) are used.

Hawkins and Bey (1994) analyzed the rat tibialis anterior (TA) muscle, which consists of about 97.5% type-II fibers (Staron et al., 1999). Therefore, in the model all fibers are assumed to be of type II. The numerical specimen used for the comparison is chosen as a rectangular cuboid with dimensions 4 cm × 2 cm × 2 cm. The fibers are aligned with the long edge of the cuboid. Starting from the stress-free reference configuration, the muscle specimen is passively stretched along the fiber direction to the desired muscle length. After passively stretching, displacement in the direction of the fibers is constrained at both ends of the specimen in order to simulate fixed-end contractions. Moreover, displacement at two further non-parallel faces of the specimen is constrained in the direction perpendicular to the respective face (symmetry boundary conditions). Note that the lengths of the individual half-sarcomeres are not constrained but only the total length of the muscle. A stimulation frequency of *f _{s}* = 100 Hz is applied to the central half-sarcomere model of each muscle fiber model. The simulation output is the nominal stress, which is defined as the ratio of the resulting reaction forces in fiber direction and the initial cross-sectional area of the specimen. The peak nominal stress of the chemo-electro-mechanical model induced through the passive stretch and the applied stimulation provides the value of the total model. The determined passive and total nominal stresses at different muscle stretches are shown in Figure 6, together with the experimental stress-stretch data of Hawkins and Bey (1994). Note that Hawkins and Bey (1994) used an unrealistic high stimulation frequency of 250 Hz. The biophysical half-sarcomere model can not account for such high frequencies. However, in the model force saturation occurs at a stimulation frequency of about 100 Hz. For this stimulation frequency complete fusion of twitches occurs, which was also reported for the experiment.

**Figure 6. Muscle stress-stretch relation**. Shown are the passive and total stresses computed using the coupled chemo-electro-mechanical model, and the experimental data of rat TA muscle (Hawkins and Bey, 1994). Simulations are carried out at stretches varying from 0.8 to 1.4 in steps of size 0.1, and at λ_{f} = 0.75, 0.76, 1.35, and 1.375.

After establishing realistic mechanical behavior under isometric conditions, the coupled chemo-electro-mechanical model is now tested for its capacity to reproduce experimental *F*-*v* data of whole muscle. The hyperbolic *F*-*v* relation of Hill (1938) can be expressed by

where κ is a dimensionless parameter, *F _{iso}* denotes the maximum isometric force, and

*v*is the maximum shortening velocity at

_{max}*F*= 0.

In the literature, κ ranges from 0.15 to 0.25 (McMahon, 1984). For example, Ranatunga (1984) reports a mean value of κ = 0.24 for rat soleus muscle. Rat soleus muscle consists mainly of type-I fibers (Soukup et al., 2002), and hence, all half-sarcomere models in the chemo-electro-mechanical model use now the type-I parametrization of Shorten et al. (2007). The parameters quantifying the cooperative effects and the distortion dependence are set to ν = 3.0 and ϑ = 1700, respectively.

Within the numerical experiments, the model specimen is first passively stretched to the optimal length. Then, the length of the specimen is kept fixed, and all fibers are fully activated (*f _{s}* = 100 Hz). For a prescribed velocity, the corresponding reaction force is computed. The resulting

*F*-

*v*data are depicted in Figure 7, where the force values have been normalized to the value at isometric conditions and the velocity has been normalized to the maximum shortening velocity.

**Figure 7. F-v data computed using the multiscale chemo-electro-mechanical model**. Shown are the

*F*-

*v*data of the model (black crosses), the corresponding fit of Hill's hyperbolic relation (κ = 0.241, blue line), and the region of typical muscle

*F*-

*v*curves (0.15≤ κ ≤ 0.25, light-blue shaded area).

Fitting the parameter κ in Equation (10) to the simulation results obtained for shortening contractions in a least-squares sense yields κ = 0.241, cf. Figure 7. For lengthening contractions, the chemo-electro-mechanical model predicts a maximum force of 1.77 times the isometric force, cf. Figure 7.

### 3.3. Compartmentalization

After verifying that the multiscale model is capable of predicting experimental *F*-*ℓ* and *F*-*v* data of whole muscle, the chemo-electro-mechanical skeletal muscle model is used to compare series-fibered and spanning-fibered muscles. The aim of this comparison is to reveal differences in the mechanical behavior of the different muscle fiber arrangements.

In all of the following numerical experiments, a rectangular cuboid with dimensions 12 cm × 2 cm × 2 cm is considered. The fascicle direction is assumed to be aligned with the cuboid's long edge. To mimic series-fibered skeletal muscle arrangements, the long side of the muscle specimen is subdivided into compartments of equal length. The fibers in adjacent compartments are aligned end-to-end, and do not interdigitate with each other. As in real muscle, electrical activation from one fiber to adjacent ones does not occur, neither between adjacent compartments, nor in lateral direction within a compartment. The neuromuscular junction of each fiber is assumed to be located in the middle of the respective fiber. All half-sarcomeres are assumed to be of type II. The mechanical behavior of the chemo-electro-mechanical muscle model is investigated for simultaneously stimulating all fibers. Before stimulating the muscle specimen, it is passively stretched to the optimal length (λ^{opt}_{f} = 1.2, ℓ^{opt}_{S} = 2.4 μm).

First, fixed-end contractions and shortening contractions at 10% of the maximum shortening velocity at *f _{s}* = 50 and 100 Hz are considered. A muscle model with fibers that span the entire length of the fascicles (referred to as

*SPA*) and a model consisting of four fiber compartments in series (referred to as

*SER*·

*4*) are compared to each other. The resulting nominal stresses are depicted in Figure 8. Fixed-end contractions predict differences of almost up to 80% between the different muscle fiber arrangements. The largest differences occur at the beginning of the contraction, i. e., during the first twitch, but decline rapidly to approximately 10% and less. Moreover, the results show that the initial differences are less pronounced in shortening contractions independent of the stimulation frequency. At

*f*= 50 Hz, twitches tend to be more fused for model

_{s}*SPA*than for model

*SER*·

*4*. This applies to both fixed-end and shortening contractions. Completely fused twitches are observed for both models for

*f*= 100 Hz.

_{s}**Figure 8. Comparison of a spanning-fibered muscle model (SPA) and a series-fibered muscle model consisting of four in-series arranged fiber compartments (SER·4) stimulated with f_{s} = 100 Hz (top row) and f_{s} = 50 Hz (bottom row) in fixed-end (left column) and shortening contractions at v = 0.1 v_{max} (right column), and their differences in percent**.

Independent of the stimulation frequency, model *SER*·*4* shows higher peak forces than model *SPA* in fixed-end and shortening contractions. At *f _{s}* = 100 Hz, the maximum force of model

*SER*·

*4*is 3.29% and 6.61% higher than the maximum force of model

*SPA*, in fixed-end and shortening contractions, respectively. The observed decrease after reaching the maximal value in all simulations with

*f*= 100 Hz is due to fatigue, which is contained in the half-sarcomere model of Shorten et al. (2007).

_{s}The results reveal that the largest differences between spanning-fibered and series-fibered muscle models occur during the first twitch in fixed-end contractions. Hence, fixed-end single twitch experiments are further investigated in the following. The aim is to reveal a potential relation between the twitch shape and the fiber length.

In addition to the model with spanning fibers (termed *SPA*), muscle specimens consisting of two, four, six, and twelve fiber compartments of equal length are considered. The series-fibered models are termed *SER*·*2*, *SER*·*4*, *SER*·*6*, and *SER*·*12* indicating the respective number of compartments. Furthermore, two different scenarios are considered. In the first scenario, all fibers in all compartments receive a stimulus at the same time to simulate a coordinated single twitch contraction. The second scenario appeals to the model with six in-series arranged compartments, in which only the fibers within the first compartment are stimulated. (Note that the choice which of the compartments is stimulated does not influence the resulting reaction forces.) This model is referred to as *SER*·*6a*.

Figure 9 shows the distribution of the membrane potential and the contraction-induced resulting deformation of the muscle in the different models of the first scenario.

**Figure 9. Distribution of membrane potential, V_{m} in [mV], and contraction-induced deformation during single twitch contractions of models SPA (t = 22 ms after stimulation), SER·2 (t = 10 ms), SER·6 (t = 5 ms), and SER·12 (t = 2 ms) (from top to bottom)**.

Figure 10 demonstrates that the twitch rise time of a muscle depends on the length of its fibers, i. e., the twitch rise time increases with increasing fiber length. Thus, model *SER*·*12* has the lowest twitch rise time of 17.2 ms, while the maximum twitch rise time occurs in model *SPA*, where the peak stress occurs 38.2 ms after stimulation. The computed AP propagation speed of the models is 2.186 m/s. In model *SPA*, where the AP propagates 6 cm from the motor end-plates to each end of the fibers, this propagation speed yields an AP propagation time of 27.45 ms. In comparison, a half-sarcomere model considered in isolation shows a twitch rise time of 16.1 ms. Hence, the AP propagation time in model *SPA* exceeds the twitch rise time of a single half-sarcomere. In other words, the sarcomeres located at the motor end-plates reach their peak twitch force before the sarcomeres located at the ends of the fibers are activated.

**Figure 10. Comparison of single twitch contractions in a spanning-fibered model and in series-fibered models with different fiber lengths and number of compartments**. The reader is referred to the text for model definitions.

While the twitch rise time increases, the peak twitch stress of the muscle model decreases with increasing fiber length. In detail, the peak twitch stresses are 0.82 and 0.98 N/cm^{2} in models *SPA* and *SER*·*12*, respectively, which corresponds to an increase of 19.4%. Integrating the area below the stress curve over 200 ms, i. e., to a point where the active stress has declined and only passive stress components remain, yields 84.95 N·ms/cm^{2} for model *SPA*, and 83.25 N·ms/cm^{2} for model *SER*·*12*.

Deducting from the total stresses the respective passive stresses, which are due to the initial stretch to optimal length, the peak twitch force obtained in model *SER*·*6a* is 6.5 times smaller than the peak twitch force of model *SER*·*6*.

Further, changes in local sarcomere length during fixed-end single twitch contractions are analyzed. The aim is to investigate if activation-induced stretches of passive sarcomeres to beyond myofilament overlap occur. The resulting maximum and minimum sarcomere lengths are reported in Table 3.

**Table 3. Minimum and maximum sarcomere lengths in fixed-end single twitch contractions absolute and in percent of their length prior to stimulation, ℓ ^{opt}_{S} = 2.4 μm**.

Considering the first scenario, the shortest and largest sarcomere lengths of 1.81 and 2.66 μm, respectively, occur for model *SPA*. Changes in sarcomere length decrease with an increasing number of in-series fiber compartments. In the second scenario, a minimum sarcomere length of 1.74 μm is observed for model *SER*·*6a*.

## 4. Discussion

A multiscale skeletal muscle model was presented that includes the description of the active behavior entirely on the microscopic sarcomere level. Yet, the model proved to be able to reproduce experimentally determined data of whole muscle on the macroscale. The multiscale model was used to investigate differences in the muscle contraction and force generation caused by different muscle fiber arrangements.

### 4.1. From Isometric Half-Sarcomere Model to Non-Isometric Whole Muscle Simulations

The *F*-*ℓ* and *F*-*v* relationships of skeletal muscle originate from properties on the microscopic filament level (Gordon et al., 1966; Piazzesi et al., 2007). For example, Winters et al. (2011) point out that the active *F*-*ℓ* relation of a whole muscle is very similar to the *F*-*ℓ* relation of a single sarcomere. Likewise, the *F*-*v* relation shows very similar characteristics on the cell level (Edman, 1988) and on the whole muscle level *in situ* (Devrome and MacIntosh, 2007). Based on these findings, the proposed model has the advantage to contain the active *F*-*ℓ* and *F*-*v* relations on the microscopic half-sarcomere level.

Extending the half-sarcomere model of Shorten et al. (2007) to non-isometric contractions introduces two more parameters to the model. The additional uncertainty due to the introduction of these parameters is minor, since both of them can easily be determined by comparing computational results to experimental *F*-*v* data. The extended half-sarcomere model can reproduce the hyperbolic *F*-*v* relation of shortening contractions and the bounded force increase in lengthening contractions known from experiments (Hill, 1938; Zajac, 1989). Similar results are reported by Razumova et al. (1999) using a different approach. Razumova et al. (1999) assumed quasi-static conditions and rearranged their XB-dynamics model such that they could compute the corresponding velocity for a prescribed force.

It is noteworthy that, in contrast to previous macroscopic models (Zajac, 1989), the hyperbolic *F*-*v* relation is not explicitly prescribed in the model but results from the XB-dynamics component model formulation. Thus, the model can be used to reveal the underlying mechanisms leading to the characteristic *F*-*v* behavior (Hernández-Gascón et al., 2013).

The active behavior on the macroscopic whole muscle level is modeled to be entirely determined by the extended half-sarcomere model. The presented results demonstrate that the multiscale model is capable of reproducing microscopic properties of the sarcomere level on the macroscopic whole muscle level. This applies likewise to the *F*-*ℓ* and the *F*-*v* relationships.

In the literature, different behaviors are reported for lengthening contractions of skeletal muscles (cf. Morgan, 1990). Zajac (1989) report a bounded increase up to 1.8 times the isometric force, which is adopted in this contribution. Since the model behavior for lengthening contractions proved to be sensitive to a single parameter, the presented model can easily be adapted to a different shape. However, the fact that experimental *F*-*v* relations show a non-continuously differentiable behavior at the transition from shortening to lengthening contractions (Katz, 1939) is not predicted by the model. Once the origin of this unique feature is completely understood, it could potentially be included in the XB-dynamics component model.

### 4.2. Compartmentalization

First the computational results obtained for the different muscle fiber arrangements are discussed, before using this data to analyze its implications on stability.

The presented model predicts the largest differences between series-fibered and spanning-fibered muscles in the rise time, shape and peak force of single twitches. During sustained contractions, twitches tended to fuse at lower stimulation frequencies in spanning-fibered muscles, while series-fibered muscles showed higher peak forces. Since the basic descriptions of passive and active material behavior are identical in the different models, the observed differences in the force responses must result from the differences in the muscle fiber arrangement. Although the same half-sarcomere model is used in all simulations, single twitches are more dispersed in muscle models with longer fibers, which can be explained by longer AP propagation times. Experimentally observed differences in the twitch shape in different fibers of the same twitch type might therefore be largely governed by the fiber length. This might explain the different twitch shapes observed in different species. For example, the twitch rise time in mouse soleus muscle consisting purely of type-I fibers is approximately 35 ms (Shorten et al., 2007), while 90 ms are observed in human type-I motor units (Fuglevand et al., 1993). Further, the simulations demonstrated that a fascicle consisting of end-to-end terminating fibers does functionally not perform like a single muscle fiber of equivalent length, as hypothesized by Lieber and Fridén (2000).

According to Harris et al. (2005), long fibers are less efficient than short fibers since sarcomere shortening cannot be well synchronized along the length of a fiber. Harris et al. (2005) speculate that a twitch in a long fiber will produce much less force than a more synchronous contraction of the sarcomeres. The presented results confirm that the peak twitch force in spanning-fibered muscle is lower than in series-fibered muscle of the same length, however, it is also more dispersed, such that the stress induced through a single twitch integrated over time is similar in series-fibered and spanning-fibered muscles. This can be attributed to the fact that the number of sarcomeres contributing to the active force is identical in both models. The non-activated parts of the fibers behave as series elastic elements, i. e., they store contractile energy. It is believed that the minor differences observed in the integrated stress values stem from local changes in sarcomere length due to the *F*-*ℓ* relation and from different sarcomere contraction velocities due to the *F*-*v* relation. At this point, however, one has to bear in mind that the modeling assumption of hyperelastic passive material behavior neglects viscous effects, which exist in passive muscle (Hoyt et al., 2008; Van Loocke et al., 2008).

The model further predicts that the peak force exerted by a synchronous activation of all in-series arranged compartments exceeds the product of the number of in-series arranged compartments and the peak force produced when stimulating only the fibers in one compartment. This might be explained by the fact that an additional series compliance is introduced through inactive compartments against which the activated fibers contract (Botterman et al., 1983). It is hypothesized that the effect will be more pronounced at shorter muscle lengths than at the optimal length (at which the numerical experiments are carried out) (cf. Mutungi and Ranatunga, 2000), or in muscles with passive forces appearing only at long muscle length (see further below).

Changes in sarcomere length due to the contraction of activated parts of the fibers against non-activated parts are reported for spanning-fibered and series-fibered muscle models. Fixed-end single twitch contractions, in which the fibers of all compartments are simultaneously activated, show that changes in sarcomere length increase with increasing fiber length. Shorter sarcomere lengths are only observed if one out of six compartments is activated (model *SER*·*6a*). This is not surprising as the five non-activated compartments act as series elastic elements. Comparing the extreme values of the sarcomere length with Figure 2 reveals that the range of sarcomere lengths of the numerical experiments is limited to a rather narrow region with considerable filament overlap. Mutungi and Ranatunga (2000) report experimental sarcomere length changes in fixed-end single twitch contractions that are considerably smaller than those found in the present numerical investigations. The difference can be explained based on the fact that Mutungi and Ranatunga (2000) simultaneously stimulated the entire fiber bundle using plate electrodes, and hence, almost all sarcomeres shortened concurrently against a small region at the fiber ends.

The fact that the model predicts rather small changes in sarcomere length during fixed-end single twitch contractions might be explained by the following considerations. A resting sarcomere length of 2.0 μm (Edman, 1979) is assigned to the model's stress-free reference configuration (λ_{f} = 1). Thus, the longest observed sarcomere length of 2.66 μm corresponds to a local fiber stretch of λ_{f} = 1.33. Comparing this value with the *F*-*ℓ* relation in Figure 6, one observes that considerable passive forces start to appear at this fiber stretch. This can be explained by the fact that at every instant in time, the contractile forces in the activated parts of the muscle need to be matched by the stretch-induced passive forces in the non-activated parts, since they are in-series arranged. Sarcomere length changes will therefore be more pronounced in muscles with passive forces appearing at long whole muscle length.

A description of tendon was not included in the model. Since tendinous tissue is much stiffer than passive muscle tissue (Hawkins and Bey, 1997), the series compliance added to the system by including tendon is small. Therefore, the effect of neglecting tendon in this study is expected to have a minor effect on the force generation and the sarcomere length changes. It should be noted, however, that this only applies to parallel-fibered muscles. In general, tendons and aponeuroses are crucial to muscle-joint dynamics. Therefore, future models should incorporate tendon and aponeurosis compliance to better link sarcomere dynamics to joint dynamics during movement.

The study of compartmentalization is particularly interesting with regard to stability issues. The model results demonstrate that activated parts of a muscle can contract against non-activated parts. It has been hypothesized that in long spanning-fibered muscle, in which the AP propagation time exceeds the twitch rise time, activation-induced stresses might stretch non-activated sarcomeres to beyond myofilament overlap potentially leading to instabilities (Loeb et al., 1987). Loeb et al. (1987) therefore speculate that the twitch rise time might impose a limit on the length of the fibers. The presented results, however, demonstrate that a muscle model, in which the AP propagation time exceeds the twitch rise time of a single sarcomere, does not necessarily show any instabilities. In series-fibered muscle, a similar stability problem is believed to exist when activation of series-arranged compartments is unbalanced or asynchronous, i. e., if fibers in an activated compartment shorten against fibers in non-activated compartments (Richmond et al., 1985; Loeb et al., 1987). This instability was not observed either in the numerical experiments (model *SER*·*6a*) using the presented model settings.

The fact that instabilities are observed neither in the spanning-fibered model nor in the series-fibered model might be due to the fact that in the present model passive forces appear already at short muscle length. According to Hawkins and Bey (1994), this corresponds to the behavior of rat TA muscle, which shows even at full activation a monotonically increasing isometric *F*-*ℓ* relation, cf. Figure 6. The passive stiffness of the muscle tissue might therefore prevent an overextension of non-activated sarcomeres. However, in muscles with passive forces appearing at long muscle length, sarcomere extensions to beyond myofilament overlap might be possible, and this might lead to stability problems and damage (Loeb et al., 1987).

In the future, the proposed multiscale model can be used, for example, to study sarcomere length changes in muscles, in which passive forces appear at long muscle length and the associated potential instabilities. Furthermore, the presented framework can be used to study the implications of the task-specific activation of sub-volumes of a muscle on the muscle contraction and force generation.

### 4.3. Summary

A chemo-electro-mechanical skeletal muscle model has been developed to reveal differences between parallel-fibered muscles, in which the muscle fibers either span the entire length of the fascicles or terminate intrafascicularly. The multiscale model proved to be able to reveal differences in the muscle contraction and force generation that result from the muscle fiber arrangement. The largest differences in the mechanical behaviors due to the different arrangements have been found during fixed-end single twitch contractions. Spanning-fibered muscles showed lower but more dispersed twitch forces than series-fibered muscles of the same length. Similarly, series-fibered muscles showed significantly higher peak forces during sustained contractions. Further, sarcomere length changes during fixed-end single twitch contractions of the multiscale muscle model at optimal sarcomere length have been analyzed. It was found that the sarcomere length changes were limited to a rather narrow region with considerable filament overlap. Stability issues resulting from activation-induced stretches of non-activated sarcomeres to beyond myofilament overlap were not observed. It is concluded that in muscles with passive forces appearing at short muscle length these stability problems do not exist.

## 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.

## Acknowledgments

This work was supported by the German Research Foundation (DFG) within the funding programme Open Access Publishing. Further, the authors would like to thank the German Research Foundation (DFG) for financial support of the project within the Cluster of Excellence in Simulation Technology (EXC 310/1) at the University of Stuttgart. Moreover, the research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement n. 306757 (LEAD).

## References

Böl, M., Weikert, R., and Weichert, C. (2011). A coupled electromechanical model for the excitation-dependent contraction of skeletal muscle. *J. Mech. Behav. Biomed. Mater*. 4, 1299–1310. doi: 10.1016/j.jmbbm.2011.04.017

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

Bonet, J., and Wood, R. D. (2008). *Nonlinear Continuum Mechanics for Finite Element Analysis*. Cambridge: Cambridge University Press.

Botterman, B. R., Hamm, T. M., Reinking, R. M., and Stuart, D. G. (1983). Distribution of monosynaptic Ia excitatory post-synaptic potentials in the motor nucleus of the cat semitendinosus muscle. *J. Physiol*. 338, 379–393.

Bradley, C. P., Bowery, A., Britten, R., Budelmann, V., Camara, O., Christie, R., et al. (2011). OpenCMISS: a multi-physics & multi-scale computational infrastructure for the VPH/Physiome project. *Prog. Biophys. Mol. Biol*. 107, 32–47. doi: 10.1016/j.pbiomolbio.2011.06.015

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Burkholder, T. J., and Lieber, R. L. (2001). Sarcomere length operating range of vertebrate muscles during movement. *J. Exp. Biol*. 204, 1529–1536.

Campbell, K. B., Razumova, M. V., Kirkpatrick, R. D., and Slinker, B. K. (2001a). Myofilament kinetics in isometric twitch dynamics. *Ann. Biomed. Eng*. 29, 384–405. doi: 10.1114/1.1366669

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Campbell, K. B., Razumova, M. V., Kirkpatrick, R. D., and Slinker, B. K. (2001b). Nonlinear myofilament regulatory processes affect frequency-dependent muscle fiber stiffness. *Biophys. J*. 81, 2278–2296. doi: 10.1016/S0006-3495(01)75875-4

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

De Luca, C. J., and Hostage, E. C. (2010). Relationship between firing rate and recruitment threshold of motoneurons in voluntary isometric contractions. *J. Neurophysiol*. 104, 1034–1046. doi: 10.1152/jn.01018.2009

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Denoth, J., Stüssi, E., Csucs, G., and Danuser, G. (2002). Single muscle fiber contraction is dictated by inter-sarcomere dynamics. *J. Theor. Biol*. 216, 101–122. doi: 10.1006/jtbi.2001.2519

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Devrome, A., and MacIntosh, B. (2007). The biphasic force-velocity relationship in whole rat skeletal muscle *in situ*. *J. Appl. Physiol*. 102, 2294–2300. doi: 10.1152/japplphysiol.00276.2006

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Edman, K. (1979). The velocity of unloaded shortening and its relation to sarcomere length and isometric force om vertebrate muscle fibres. *J. Physiol*. 291, 143–159.

Edman, K. (1988). Double-hyperbolic force-velocity relation in frog muscle fibres. *J. Physiol*. 404, 301–321.

Epstein, B. R., and Foster, K. R. (1983). Anisotropy in the dielectric properties of skeletal muscle. *Med. Biol. Eng. Comput*. 21, 51–55. doi: 10.1007/BF02446406

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Fernandez, J. W., Buist, M. L., Nickerson, D. P., and Hunter, P. J. (2005). Modelling the passive and nerve activated response of the rectus femoris muscle to a flexion loading: a finite element framework. *Med. Eng. Phys*. 27, 862–870. doi: 10.1016/j.medengphy.2005.03.009

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Fuglevand, A. J., Winter, D. A., and Patla, A. E. (1993). Models of recruitment and rate coding organization in motor unit pools. *J. Neurophysiol*. 70, 2470–2488.

Günther, M., Röhrle, O., Häufle, D. F., and Schmitt, S. (2012). Spreading out muscle mass within a hill-type model: a computer simulation study. *Comput. Math. Methods Med*. 2012, 1–13. doi: 10.1155/2012/848630

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Gielen, F. L. H., Wallinga-de Jonge, W., and Boon, K. L. (1984). Electrical conductivity of skeletal muscle tissue: experimental results from different muscles *in vivo*. *Med. Biol. Eng. Comput*. 22, 569–577. doi: 10.1007/BF02443872

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

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

Harris, A. J., Duxson, M. J., Butler, J. E., Hodges, P. W., Taylor, J. L., and Gandevia, S. C. (2005). Muscle fiber and motor unit behavior in the longest human skeletal muscle. *J. Neurosci*. 25, 8528–8533. doi: 10.1523/JNEUROSCI.0923-05.2005

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Hawkins, D., and Bey, M. (1994). A comprehensive approach for studying muscle-tendon mechanics. *J. Biomech. Eng*. 116, 51–55. doi: 10.1115/1.2895704

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Hawkins, D., and Bey, M. (1997). Muscle and tendon force-length properties and their interactions *in vivo*. *J. Biomech*. 30, 63–70. doi: 10.1016/S0021-9290(96)00094-2

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Heidlauf, T., and Röhrle, O. (2013). Modeling the chemoelectromechanical behavior of skeletal muscle using the parallel open-source software library OpenCMISS. *Comput. Math. Methods Med*. 2013, 1–14. doi: 10.1155/2013/517287

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Hernández-Gascón, B., Grasa, J., Calvo, B., and Rodríguez, J. (2013). A 3D electro-mechanical continuum model for simulating skeletal muscle contraction. *J. Theor. Biol*. 335, 108–118. doi: 10.1016/j.jtbi.2013.06.029

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Heron, M. I., and Richmond, F. J. (1993). In-series fiber architecture in long human muscles. *J. Morphol*. 216, 35–45. doi: 10.1002/jmor.1052160106

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Hill, A. V. (1938). The heat of shortening and the dynamic constants of muscle. *Proc. R. Soc. Lond. Ser. B Biol. Sci*. 126, 136–195. doi: 10.1098/rspb.1938.0050

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Hodgkin, A. L., and Huxley, A. F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. *J. Physiol*. 117, 500–544.

Horowits, R. (1992). Passive force generation and titin isoforms in mammalian skeletal muscle. *Biophys. J*. 61, 392–398. doi: 10.1016/S0006-3495(92)81845-3

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Hoyt, K. H., Kneezel, T., Castaneda, B., and Parker, K. J. (2008). Quantitative sonoelastography for the *in vivo* assessment of skeletal muscle viscoelasticity. *Phys. Med. Biol*. 53, 4063–4080. doi: 10.1088/0031-9155/53/15/004

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Huijing, P. A. (1999). Muscle as a collagen fiber reinforced composite: a review of force transmission in muscle and whole limb. *J. Biomech*. 32, 329–345. doi: 10.1016/S0021-9290(98)00186-9

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Huxley, A. F. (1957). Muscle structure and theories of contraction. *Prog. Biophys. Biophys. Chem*. 7, 255–318.

Johansson, T., Meier, P., and Blickhan, R. (2000). A finite-element model for the mechanical analysis of skeletal muscles. *J. Theor. Biol*. 206, 131–149. doi: 10.1006/jtbi.2000.2109

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Katz, B. (1939). The relation between force and speed in muscular contraction. *J. Physiol*. 96, 45–64.

Lieber, R., and Fridén, J. (2000). Functional and clinical significance of skeletal muscle architecture. *Muscle Nerve* 23, 1647–1666. doi: 10.1002/1097-4598(200011)23:11<1647::AID-MUS1>3.0.CO;2-M

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Loeb, G., Pratt, C., Chanaud, C., and Richmond, F. (1987). Distribution and innervation of short, interdigitated muscle fibers in parallel-fibered muscles of the cat hindlimb. *J. Morphol*. 191, 1–15. doi: 10.1002/jmor.1051910102

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Markert, B., Ehlers, W., and Karajan, N. (2005). A general polyconvex strain-energy function for fiber-reinforced materials. *Proc. Appl. Math. Mech*. 5, 245–246. doi: 10.1002/pamm.200510099

Morgan, D., Mochon, S., and Julian, F. (1982). A quantitative model of intersarcomere dynamics during fixed-end contractions of single frog muscle fibers. *Biophys. J*. 39, 189–196. doi: 10.1016/S0006-3495(82)84507-4

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Morgan, D. (1990). New insights into the behavior of muscle during active lengthening. *Biophys. J*. 57, 209–221. doi: 10.1016/S0006-3495(90)82524-8

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Mutungi, G., and Ranatunga, K. (2000). Sarcomere length changes during end-held (isometric) contractions in intact mammalian (rat) fast and slow muscle fibres. *J. Muscle Res. Cell Motil*. 21, 565–575. doi: 10.1023/A:1026588408907

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Paul, A. C. (2001). Muscle length affects the architecture and pattern of innervation differently in leg muscles of mouse, guinea pig, and rabbit compared to those of human and monkey muscles. *Anat. Rec*. 262, 301–309. doi: 10.1002/1097-0185(20010301)262:3<301::AID-AR1045>3.0.CO;2-H

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Piazzesi, G., Reconditi, M., Linari, M., Lucii, L., Bianco, P., Brunello, E., et al. (2007). Force transmission in skeletal muscle: from actomyosin to external tendons. *Cell* 131, 784–795. doi: 10.1016/j.cell.2007.09.045

Prado, L. G., Makarenko, I., Andresen, C., Krüger, M., Opitz, C. A., and Linke, W. A. (2005). Isoform diversity of giant proteins in relation to passive and active contractile properties of rabbit skeletal muscles. *J. Gen. Physiol*. 126, 461–480. doi: 10.1085/jgp.200509364

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Pullan, A. J., Buist, M. L., and Cheng, L. K. (2005). *Mathematically Modelling the Electrical Activity of the Heart: From Cell to Body Surface and Back Again*. Singapore: World Scientific Publishing Company.

Röhrle, O., and Pullan, A. J. (2007). Three-dimensional finite element modelling of muscle forces during mastication. *J. Biomech*. 40, 3363–3372. doi: 10.1016/j.jbiomech.2007.05.011

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Röhrle, O., Davidson, J. B., and Pullan, A. J. (2008). Bridging scales: a three-dimensional electromechanical finite element model of skeletal muscle. *SIAM J. Sci. Comput*. 30, 2882–2904. doi: 10.1137/070691504

Röhrle, O., Davidson, J. B., and Pullan, A. J. (2012). A physiologically based, multi-scale model of skeletal muscle structure and function. *Front. Physiol*. 3:358. doi: 10.3389/fphys.2012.00358

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Ranatunga, K. W. (1984). The force-velocity relation of rat fast- and slow-twitch muscles examined at different temperatures. *J. Physiol*. 351, 517–529.

Rassier, D., MacIntosh, B., and Herzog, W. (1999). Length dependence of active force production in skeletal muscle. *J. Appl. Physiol*. 86, 1445–1457.

Razumova, M. V., Bukatina, A. E., and Campbell, K. B. (1999). Stiffness-distortion sarcomere model for muscle simulation. *J. Appl. Physiol*. 87, 1861–1876.

Razumova, M. V., Bukatina, A. E., and Campbell, K. B. (2000). Different myofilament nearest neighbor interactions have distinctive effects on contractile bahavior. *Biophys. J*. 78, 3120–3137. doi: 10.1016/S0006-3495(00)76849-4

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Richmond, F. J., MacGillis, D. R., and Scott, D. A. (1985). Muscle-fiber compartmentalization in cat splenius muscles. *J. Neurophysiol*. 53, 868–885.

Sánchez, C. A., Lloyd, J. E., Fels, S., and Abolmaesumi, P. (2014). Embedding digitized fibre fields in finite element models of muscles. *Comput. Methods Biomech. Biomed. Eng. Imaging Vis*. 2, 223–236. doi: 10.1080/21681163.2013.862861

Shorten, P. R., O'Callaghan, P., Davidson, J. B., and Soboleva, T. K. (2007). A mathematical model of fatigue in skeletal muscle force contraction. *J. Muscle Res. Cell Motil*. 28, 293–313. doi: 10.1007/s10974-007-9125-6

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Siebert, T., Rode, C., Herzog, W., Till, O., and Blickhan, R. (2008). Nonlinearities make a difference: comparison of two common Hill-type models with real muscle. *Biol. Cybern*. 98, 133–143. doi: 10.1007/s00422-007-0197-6

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Soukup, T., Zachařová, G., and Smerdu, V. (2002). Fibre type composition of soleus and extensor digitorum longus muscles in normal female inbred Lewis rats. *Acta Histochem*. 104, 399–405. doi: 10.1078/0065-1281-00660

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Staron, R. S., Kraemer, W. J., Hikida, R. S., Fry, A. C., Murray, J. D., and Campos, G. E. (1999). Fiber type composition of four hindlimb muscles of adult fisher 344 rats. *Histochem. Cell Biol*. 111, 117–123. doi: 10.1007/s004180050341

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Stoecker, U., Telley, I. A., Stüssi, E., and Denoth, J. (2009). A multisegmental cross-bridge kinetics model of the myofibril. *J. Theor. Biol*. 259, 714–726. doi: 10.1016/j.jtbi.2009.03.032

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Telley, I. A., and Denoth, J. (2007). Sarcomere dynamics during muscular contraction and their implications to muscle function. *J. Muscle Res. Cell Motil*. 28, 89–104. doi: 10.1007/s10974-007-9107-8

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Van Loocke, M., Lyons, C., and Simms, C. (2008). Viscoelastic properties of passive skeletal muscle in compression: stress-relaxation behaviour and constitutive modelling. *J. Biomech*. 41, 1555–1566. doi: 10.1016/j.jbiomech.2008.02.007

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Winters, T. M., Takahashi, M., Lieber, R. L., and Ward, S. R. (2011). Whole muscle length-tension relationships are accurately modeled as scaled sarcomeres in rabbit hindlimb muscles. *J. Biomech*. 44, 109–115. doi: 10.1016/j.jbiomech.2010.08.033

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Woodley, S., and Mercer, S. (2005). Hamstring muscles: architecture and innervation. *Cells Tissues Organs* 179, 125–141. doi: 10.1159/000085004

Young, M., Paul, A., Rodda, J., Duxson, M., and Sheard, P. (2000). Examination of intrafascicular muscle fiber terminations: implications for tension delivery in series-fibered muscles. *J. Morphol*. 245, 130–145. doi: 10.1002/1097-4687(200008)245:2<130::AID-JMOR4>3.0.CO;2-R

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

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

Zheng, Y., Mak, A. F., and Lue, B. (1999). Objective assessment of limb tissue elasticity: development of a manual indentation procedure. *J. Rehabil. Res. Dev*. 36, 71–85.

Zuurbier, C. J., and Huijing, P. A. (1992). Influence of muscle geometry on shortening speed of fibre, aponeurosis and muscle. *J. Biomech*. 25, 1017–1026. doi: 10.1016/0021-9290(92)90037-2

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Zuurbier, C. J., Heslinga, J. W., Lee-de Groot, M. B., and Van der Laarse, W. J. (1995). Mean sarcomere length-force relationship of rat muscle fibre bundles. *J. Biomech*. 28, 83–87. doi: 10.1016/0021-9290(95)80009-3

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Keywords: spanning-fibered, series-fibered, sarcomere stretch, sarcomere instability, biophysical cell model, non-isometric

Citation: Heidlauf T and Röhrle O (2014) A multiscale chemo-electro-mechanical skeletal muscle model to analyze muscle contraction and force generation for different muscle fiber arrangements. *Front. Physiol*. **5**:498. doi: 10.3389/fphys.2014.00498

Received: 02 October 2014; Accepted: 02 December 2014;

Published online: 23 December 2014.

Edited by:

Taian Mello Martins Vieira, Universidade Federal do Rio de Janeiro, BrazilReviewed by:

Chris Richards, Royal Veterinary College, UKNiccolo Fiorentino, University of Utah, USA

Copyright © 2014 Heidlauf and Röhrle. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Oliver Röhrle, Institute of Applied Mechanics (CE), University of Stuttgart, Pfaffenwaldring 7, 70569 Stuttgart, Germany e-mail: roehrle@simtech.uni-stuttgart.de

## COMMENTARY