Size, History-Dependent, Activation and Three-Dimensional Effects on the Work and Power Produced During Cyclic Muscle Contractions

doi:10.1093/icb/icy021

Size, History-Dependent, Activation and Three-Dimensional Effects on the Work and Power Produced During Cyclic Muscle Contractions

2018-08-01 00:00:00 Abstract Muscles undergo cycles of length change and force development during locomotion, and these contribute to their work and power production to drive body motion. Muscle fibers are typically considered to be linear actuators whose stress depends on their length, velocity, and activation state, and whose properties can be scaled up to explain the function of whole muscles. However, experimental and modeling studies have shown that a muscle’s stress additionally depends on inactive and passive tissues within the muscle, the muscle’s size, and its previous contraction history. These effects have not been tested under common sets of contraction conditions, especially the cyclic contractions that are typical of locomotion. Here we evaluate the relative effects of size, history-dependent, activation and three-dimensional effects on the work and power produced during cyclic contractions of muscle models. Simulations of muscle contraction were optimized to generate high power outputs: this resulted in the muscle models being largely active during shortening, and inactive during lengthening. As such, the history-dependent effects were dominated by force depression during simulated active shortening rather than force enhancement during active stretch. Internal work must be done to deform the muscle tissue, and to accelerate the internal muscle mass, resulting in reduced power and work that can be done on an external load. The effect of the muscle mass affects the scaling of muscle properties, with the inertial costs of contraction being relatively greater at larger sizes and lower activation levels. Introduction Muscles typically undergo cycles of activation, length change, and force development during locomotion (Dickinson et al. 2000). An understanding of how muscles operate during these cyclical regimes is thus necessary to determine the muscle design criteria for achieving locomotor tasks in a manner that is metabolically efficient or mechanically advantageous in terms of work and power. Pioneering studies on insect flight muscle in the 1950s characterized muscle forces during high frequency bursts of excitation coupled with cyclic length changes by expressing muscle force as a function of muscle length in a manner now known as a “work-loop” (Boettinger 1957; Machin and Pringle 1960). Boettinger (1957) noted that the slope of the major axis of the loop represented the compliance of the excited muscle, the area enclosed in the loop represented the work done against damping forces in the mechanical system, and the position of the loop in the force–length area represented the magnitude of contraction velocity. These early studies on asynchronous muscle, where the excitation is not synchronized with the contraction cycles, were later extended in the 1980s to synchronous insect flight muscle (Josephson 1985) and crustacean muscle (Josephson and Stokes 1989), where the timing of the excitation or action potentials is related to phase of the contraction cycle. This work-loop approach has since been generalized to a range of vertebrate muscles. Understanding the determinants of muscle performance during vertebrate locomotion requires measures of both muscle force and length changes during in vivo activities. Prior to the 1980s when sonomicrometry was introduced to measure in vivo lengths and tendon buckles to measure in-series tendon forces (Loeb et al. 1985), work-loop studies were primarily done in situ or in vitro with imposed artificial length changes and excitation patterns. These new techniques to measure in vivo function fueled a proliferation of studies applying the work-loop technique in fish muscle (Altringham and Johnston 1990) to mimic the deformations observed in vivo (Franklin and Johnston 1997), and subsequently, in vivo work-loops have been determined for a range of species including pigeons (Biewener et al. 1998a), turkeys (Roberts et al. 1997), guinea fowl (Daley and Biewener 2003), mallard ducks (Williamson et al. 2001), wallabies (Biewener et al. 1998b), and goats (McGuigan et al. 2009). Nonetheless, we still know very little about the in vivo behavior of contracting muscle for a wide range of locomotor activities in the majority of muscles in the majority of species, as well as the mechanisms that underlie this behavior. Hill-type muscle models are widely used in comparative biomechanics to estimate in vivo muscle forces where direct measures are not feasible (Van Leeuwen 1992; Wakeling and Johnston 1999; Hodson-Tole and Wakeling 2009; Lee et al. 2013; Richards and Clemente 2013), and because of our limited understanding of in vivo whole muscle mechanics, the intrinsic properties of these models are often derived from experiments on single fibers or fiber bundles that are maximally activated and constrained to contract in a steady manner. As a consequence, the few studies that have validated Hill-type models against in vivo or in situ measures have found that these models perform best during maximal, steady-state contractions and suffer under more functionally-relevant conditions where activations and muscle lengths and velocities are not constant. Perreault et al. (2003) found that model errors (percent root-mean-square) can exceed 50% when large changes in muscle length are modeled in situ. More recently, we have achieved modeling errors of 32% for in vivo locomotion in goats (Lee et al. 2013), and 13.5% in humans (Dick et al. 2017), however, correlation between measured and predicted forces only reached a maximum of r2 ≈ 0.5. Therefore, assuming that whole muscles behave as single fibers during controlled, steady contractions fails to account for the complex mechanisms that underlie whole muscle behavior in vivo. Hill-type muscle models typically assume that muscles are composed of a single fiber-type and contain no internal mass or inertial effects, and so can be scaled to any size without any consequences on the kinematic and dynamic behavior of the model (Ross et al. 2018). However, studies have shown that the presence of different fiber-types is important for muscle force production (Wakeling et al. 2012; Lee et al. 2013), and that muscle mass provides inertial damping within the muscle (Günther et al. 2012; Ross and Wakeling 2016) that is dependent on fiber-type and activation level (Holt et al. 2014). The mass within muscle tissue takes time to move when the muscle forces cause accelerations, and this has been shown to slow contraction times (Günther et al. 2012) and maximum effective contraction speeds (Ross and Wakeling 2016) in isotonic shortening contractions. Although the implications of these inertial effects due to muscle tissue mass on muscle performance during cyclic contractions are not yet known, we would expect the inclusion of internal muscle mass to reduce the shortening velocity of the muscle, and thus the power generated. Most Hill-type muscle models calculate muscle force as a function of purely the instantaneous activation state, muscle length, and velocity, and do not consider previous contraction history. However, history-dependent effects, which are distinct from time-dependent effects due to excitation-activation dynamics, may be important for muscle contractile behavior. In particular, when compared with isometric forces at a given length, muscle forces may be elevated or enhanced when the muscle reaches that length with prior active lengthening (Abbott and Aubert 1952; Cavagna and Citterio 1974; Edman et al. 1982; Herzog and Leonard 2002; Hisey et al. 2009), and may be reduced or depressed when the muscle reaches that length with prior active shortening (Abbott and Aubert 1952; Maréchal and Plaghki 1979; Meijer et al. 1998; Herzog et al. 2000). The few experimental studies examining history-dependent effects during cyclic contractions have been limited to isokinetic stretch-shortening and shortening-stretch regimes with constant, maximal activation in situ in cats (Herzog and Leonard 2000; Lee et al. 2001) or in vivo in humans (Seiberl et al. 2015; Fortuna et al. 2017, 2018). Furthermore, while a number of muscle models have incorporated history-dependent effects (Forcinito et al. 1998; Meijer et al. 1998; Wu and Herzog 1999; Ettema and Meijer 2000; Ettema 2002; Tamura et al. 2005; McGowan et al. 2010, 2013; Williams 2010; Nishikawa et al. 2012), none have examined the contribution of these effects during cyclic contractions with unconstrained length changes and time-varying activation. Muscle is nearly isovolumetric (Baskin and Paolini 1967) and hence when fibers shorten they must expand in girth to maintain their volume. The extent to which such transverse deformations occur is resisted by stiffness both within and external to the muscle fibers (Huijing 1999). Transverse stiffness in the muscle resists the transverse expansion, and the volumetric constraints result in this being seen as a resistance of the fibers to shortening (Azizi et al. 2017). These three-dimensional (3D) phenomena are also not considered by typical Hill-type models (Zajac 1989) that treat muscles as linear actuators with only 1D properties. However, we would expect that the inclusion of base material stiffness within the muscle, with properties in all three dimensions, would decrease the power available to work against an external load. Despite the wealth of experimental and modeling evidence that these additional effects are important for in vivo muscle behavior, they have only yet been evaluated in separate studies with specific contraction regimes to elucidate their specific effects, so their individual and relative contributions are still not known. Therefore, the purpose of this study was to evaluate the contributions of internal mass, fiber-type, history-dependent, and 3D effects to the work and power output of muscle during cyclic contractions. To achieve this, we used a modeling framework that drives cyclic contractions of a Hill-type muscle model in series with a damped harmonic oscillator to create muscle work-loops (Ross et al. 2018). These additional features were implemented within the muscle model, and so could be tested under common contractile regimes with the same intrinsic muscle properties to determine their individual and relative contributions to cyclic muscle performance. Methods System and muscle properties The modeling framework we used to test the effects of different contractile and mechanical properties on muscle behavior was composed of a harmonic oscillator, that represented an external load placed on the muscle, in series with the muscle model (Fig. 1A). The whole system (muscle and oscillator) was held at constant length by fixing both ends. The system dynamics were driven by a cyclical activation signal a that caused the muscle to generate active force at the start of each activation cycle when the magnitude of activation increased. The muscle force Fm in turn accelerated an inertial point mass m, within the harmonic oscillator, that was considered to be external to the muscle, and this acceleration was resisted by a linear spring with stiffness k and a viscous damper with damping coefficient b that both acted only along the length (x-axis) of the system. When the muscle activation decreased at the end of each activation cycle, the harmonic oscillator force dominated the muscle force and caused the muscle to increase in length. Thus, the system allowed the muscle to cyclically contract with unconstrained muscle lengths and forces to mimic the contraction cycles seen during in vivo locomotion. Fig. 1 View largeDownload slide The different muscle models used for the simulated work-loops. (A) The spring, mass, and damper of the harmonic oscillator were in series connected to a muscle model (cylinder), with the ends of the whole system being held at fixed positions. The muscle model was evaluated using a series of 1D models: (B) a Hill-type model with a contractile element (CE) and parallel elastic element (PEE); (C) a Hill-type model with an additional history-dependent element acting in parallel; (D) a mass enhanced model where the muscle mass was distributed into 16 point masses (p1 – p16), separated by Hill-type actuator. A 3D model was also tested within the oscillating system (E) that was evaluated using a fiber-reinforced composite biomaterial and solved using the finite element method. Fig. 1 View largeDownload slide The different muscle models used for the simulated work-loops. (A) The spring, mass, and damper of the harmonic oscillator were in series connected to a muscle model (cylinder), with the ends of the whole system being held at fixed positions. The muscle model was evaluated using a series of 1D models: (B) a Hill-type model with a contractile element (CE) and parallel elastic element (PEE); (C) a Hill-type model with an additional history-dependent element acting in parallel; (D) a mass enhanced model where the muscle mass was distributed into 16 point masses (p1 – p16), separated by Hill-type actuator. A 3D model was also tested within the oscillating system (E) that was evaluated using a fiber-reinforced composite biomaterial and solved using the finite element method. The activation signal a that drove the active muscle forces was derived from a pulsed square wave excitation signal u with prescribed maximum amplitude umax, frequency f, and duty cycle D. Excitation u was converted to a via a first-order bilinear differential equation (Zajac 1989), the characteristics of which have been described in detail elsewhere (Ross et al. 2018). The activity of muscle in vivo can be increased by recruiting more muscle fibers, and also by increasing the firing rate and thus the activity of the individual fibers. However, in our simplified model, typical of the ubiquitous Hill-type muscle models, the muscle contains only one fiber that can be scaled to any size. Thus, the activity is merely a scalar value that is applied to the muscle force. The muscle simulated in this study was selected to contain fibers that were parallel to the line-of-action of the muscle, and has no aponeurosis or internal tendon. As such, the model ignores the structural effects that come from pennation and muscle fiber gearing. Muscle fibers are slender with an aspect ratio (length:diameter) that can approach 1000, however, whole muscles are much less slender with a lower aspect ratio. In order to create a shape for the model that was a compromise between a single fiber and a whole muscle, we selected an intermediate geometry with an aspect ratio of 100. This may represent the geometry of a bundle of muscle fibers, or fascicle, but is still more slender than typical whole muscles. The optimal muscle length l0 was set to 0.02 m and the muscle cross-sectional area A was set to 0.0314 mm2. The maximum isometric stress σ0 was held constant at 225 kPa. A typical Hill-type muscle model without internal mass can be scaled to any size without any consequences on the kinematic and dynamic behavior of the model (Ross et al. 2018). However, the inclusion of internal muscle mass results in scale-dependent distortions where the load due to the muscle mass mm is proportional to the muscle volume and the muscle force Fm is proportional to A (Ross and Wakeling 2016). We tested these mass effects by comparing muscles of different sizes with identical geometric proportions and aspect ratios. For example, a muscle with a length scale factor λl0 of 10 would have 10 times greater l0, 100 times greater A, and 1000 times greater volume than a muscle with λl0 of 1. To control for the dynamics of the harmonic oscillator and achieve kinematic and dynamic similarity across different values of λl0, the parameters k, b, and m for the damped harmonic oscillator were scaled according to Ross et al. (2018). The average mass-specific power per cycle P* was calculated as the product of the net mass-specific work per cycle W* and f. The harmonic oscillator properties were chosen to allow the muscle to maximize P*. As such, the values for k, b, and m, as well as D for the excitation signal, were determined using nonlinear global optimization with the objective of maximizing P* for the control condition with umax of 1.0, f of 1 Hz, and fast fiber-type properties. These parameters were then kept constant, except for D, which was optimized for each f and fiber-type. Description of models Root model The root model has been previously described elsewhere (Ross et al. 2018), and is shown in Fig. 1B. In brief, the muscle in the root model consisted of a contractile element (CE) and a parallel elastic element (PEE). The CE developed forces when activated, and the force was proportional to the activation state, the force–length and force–velocity properties within the muscle (Fig. 2), as well as the instantaneous length and velocity of the muscle. Although the force–velocity properties of the CE result in implicit damping effects, we did not incorporate any additional damping element within the muscle. The force–length and force–velocity properties were defined by their Bezier curve descriptions (Ross et al. 2018). The root model did not include internal mass, 3D, or history-dependent effects. Fig. 2 View largeDownload slide Force–velocity (A) and force–length (B) properties for the Hill-type models of the muscle fibers. The curves were parameterized using Bezier curves (Ross et al. 2018). Fig. 2 View largeDownload slide Force–velocity (A) and force–length (B) properties for the Hill-type models of the muscle fibers. The curves were parameterized using Bezier curves (Ross et al. 2018). Mass model We tested this by subdividing the muscle into a series of smaller CEs and PEEs, divided by point masses that accounted for the muscle mass, following the scheme presented by Ross and Wakeling (2016). The total muscle mass was divided into 16 CEs that were separated by 16 point masses (Fig. 1D). The intrinsic properties of the CEs and PEEs remained the same as for the root model. The muscle mass mm was calculated as the product of the muscle volume and the density ρ, and as ρ was assumed to be constant at 1060 kg m−3, mm scaled in proportion to the muscle volume. As mm was distributed evenly along the length of the muscle at rest, the mass of each point mass was equal to mm divided by the number of point masses. History-dependent model To test for the effect of history-dependent phenomena on the resulting cyclic work-loops in this study, we added a history-dependent element in parallel to the CE of the root model (Fig. 1C). The force of this history-dependent element was determined by the sum of the stretch-induced force enhancement FFE and the shortening-induced force depression FFD. The method to include FFE was styled from a recent theory (Rode et al. 2009; Herzog et al. 2012; Herzog 2014). FFE was calculated using an empirical description (in a similar manner to the force–length and force–velocity parameters of the Hill-type model) and represents the additional titin forces that are believed to develop after active stretch. When the activation of the model increased beyond a threshold (50% maximum activation for each trial), the instantaneous muscle length at the time the threshold was reached was retained as a state constant le. If the instantaneous muscle length lm was stretched beyond le, then FFE was calculated from the normalized elastic modulus and the normalized muscle stretch: FFE=F0lm-leleE^lel0, (1) where F0 is the maximum isometric force for the muscle, and Ê is the normalized elastic modulus that varies as a function of le normalized to l0 in a non-linear fashion (Fig. 3A) that was determined from experimental data for the ascending (Hisey et al. 2009) and descending (Edman et al. 1982) limbs of the force–length relationship. This empirical description of force enhancement after stretch results in predicted enhanced forces that match experimental values for steady-state contractions (Fig. 3B). When the activation state returned below the threshold, the FFE was reduced to half its value to estimate the residual force enhancement, consistent with experimental data (Hisey et al. 2009). The force enhancement FFE reduced to zero, and the muscle was returned to its non-enhanced state, when lm shortened to less than le. This process would then be repeated for the subsequent contraction cycles. For simulations where the excitation was increased to a sub-maximal level with umax of 0.1 or 0.5, the maximum activation in each burst was <1. We assumed that only a portion of the fibers in the muscle are recruited during these sub-maximal simulations, and therefore, only the active portion of fibers would have enhanced titin forces and contribute to FFE. Thus, FFE was scaled by the maximum activation achieved in each contraction cycle. Fig. 3 View largeDownload slide Method for calculating force enhancement that occurs during periods of active stretch. (A) The modulus of the enhanced force (relative to the length that the muscle was activated) was determined for experimental data from both the ascending (Hisey et al. 2009) and descending (Edman et al. 1982) limbs of the force–length relation. (B) The predicted enhanced forces are shown in the context of steady forces after stretch in comparison to experimental values from Edman et al. (1982) and Hisey et al. (2009) (circles). Fig. 3 View largeDownload slide Method for calculating force enhancement that occurs during periods of active stretch. (A) The modulus of the enhanced force (relative to the length that the muscle was activated) was determined for experimental data from both the ascending (Hisey et al. 2009) and descending (Edman et al. 1982) limbs of the force–length relation. (B) The predicted enhanced forces are shown in the context of steady forces after stretch in comparison to experimental values from Edman et al. (1982) and Hisey et al. (2009) (circles). The magnitude of force depression has been described as a linear function of the net work done during each active contraction (Herzog et al. 2000). The force depression FFD was initially modeled following McGowan et al. (2010): FFD=-CFDWc, (2) where Wc is the work done during the active contraction and CFD is the coefficient for force depression that was set to a value of 9.43 N J−1 (McGowan et al. 2010 from Herzog et al. 2000). However, the coefficient CFD had dimensions and needed to be non-dimensionalized so that it could be applied across scales. Assuming that the magnitude of force depression per sarcomere is the same across all muscle sizes, the non-dimensional force depression coefficient was calculated as: F^FD=0.3387F^mW^c, (3) where the non-dimensional work W^c and non-dimensional force F^m were given by: W^c=WcWmax (4) and F^m=FmF0. (5) Therefore, C^FD=0.3387FmWmaxF0Wc, (6) where Wmax is the product of l0 and F0, calculated using data from Herzog et al. (2000). Thus, FFD was calculated as: FFD=-F0C^FDW^c. (7) If the muscle started to lengthen while still active, Wc would decrease and so too would FFD, however, FFD was set to zero if Wc became negative. In a similar manner to the force enhancement, the contribution of the force depression was assumed to only occur in active fibers. Thus, the FFD was also scaled by the maximum activation achieved in each contraction cycle. 3D finite element model To study the effects of the 3D structure of muscle on its behavior during cyclic contractions, we extended the isometric 3D finite element model (FEM) of muscle from Rahemi et al. (2014) to allow for dynamic behavior (Fig. 1E). This FEM is described in the following section, with additional details in Appendix 1. In this scheme, the parallel muscle fibers were represented by contractile tissue surrounded by base material. The intrinsic force–velocity and force–length properties of the fibers, as well as the excitation u and activation a profiles, were the same as those of the Hill-based root model (Fig. 2). The 3D FEM muscle had the same optimal length l0 as the root model, but a smaller aspect ratio of 10 with a square A. Following the works of Weiss et al. (1996) and Blemker and Delp (2005), we modeled the muscle as a fiber-reinforced, nonlinearly-elastic transversely-isotropic material. The constitutive laws for the materials have been described in detail elsewhere (Rahemi et al. 2014), and include the effects of the muscle being nearly incompressible to prevent large changes in volume, as well as the isotropic nonlinearly-elastic behavior associated with tissue base material properties in all directions. As in Rahemi et al. (2014), a Yeoh-type model (Yeoh 1993) was used to define the base material properties, which was scaled by a positive factor sbase to control the stiffness of the base material. The force–velocity and force–length properties of the muscle only operated along the direction of the fibers, and were incorporated into a description of the strain energy. The effect of a uniform muscle density ρ throughout the material was incorporated into a kinematic energy contribution. The FEM system was defined by: dharmonicoscillatorE+muscleEdt=-bdxdt2, (8) where E denotes energy, x denotes the displacement of the point mass m of the harmonic oscillator, and b denotes the damping coefficient in the harmonic oscillator, both as for the root model. As in the root model, we have not explicitly incorporated damping effects within the muscle, although they do occur implicitly as a consequence of the force–velocity properties of the CE. The energy associated with the point mass in the harmonic oscillator m is equal to the sum of its kinetic and potential energies, while the energy associated with the muscle is equal to the sum of the muscle kinetic energy, the volumetric strain energy, the base material strain energy, and the contractile strain energy. Note that only the contractile strain energy, and not the other three terms that make up the muscle energy, is accounted for in the standard 1D Hill-type model. The volumetric strain energy captures the material’s resistance to volume changes (nearly incompressible), while the contractile strain energy describes the energetic contributions due to the passive (PEE) and active (CE) elements of the fiber. The system dynamics are described by the changes in how the total energy is distributed among these different forms of energy over time. Thus, if the total energy is fixed, the relative contributions of the harmonic oscillator energy will be higher in the 1D root model than they are in the 3D model. The overall action integral of the muscle system describes how the total energy evolves over time as a function of the local displacement, dilation, and pressure. The principle of stationarity of the action integral was used to derive the equations in these unknowns, following the work of Simo et al. (1985) and Simo and Taylor (1991). This method allows the muscle fiber stretch to vary along its length, depending on the local forces it experiences, and moreover, the fibers are permitted to curve, unlike in the 1D Hill-based models. The dynamics of the muscle were coupled through a balance of forces to the external damped harmonic oscillator system described earlier, and the fibers in the 3D model were activated using the same u as was used for the root, mass, and history-dependent models. This resulted in a mathematical system that was dynamic and highly nonlinear. The FEM muscle was subdivided into 128 smaller elements (parallelepipeds), and each element had local activation, base material properties, Hill-type contractile properties of the fibers, and specified initial fiber orientation relative to the geometry. Although these properties can be varied from voxel to voxel, for the purposes of this study, we assumed they were uniform throughout the tissue. We used semi-implicit time stepping to discretize in time, and the Newton method to solve the nonlinear equations describing the displacement, dilation, and pressure at each time step. Simulations We tested the root, history-dependent, and mass models under a range of maximum excitations umax, frequencies f, length scale factors λl0, and fiber-type properties to replicate the contraction regimes seen during in vivo locomotion. The muscle behaved as if it was composed of either entirely fast or entirely slow fibers, and this was implemented into the models by setting the maximum shortening strain-rate v0 to either 10 s−1 or 5 s−1 and the activation rate constant τact to either 0.025 s or 0.045 s, respectively. These values are selected to represent typical properties for fast and slow fibers (Dick et al. 2017). umax was set to either 0.1, 0.5, or 1.0, and f was set to either 1, 2, 4, or 8 Hz. The duty cycle D was re-optimized for the different values of f, and decreased with increased f. For the root model with v0 of 10 s−1, D ranged from 0.52 at f of 0.5 Hz to 0.31 at f of 8 Hz; for the slower-fibered model with v0 of 5 s−1, D ranged from 0.55 at f of 0.5 Hz to 0.29 at f of 8 Hz. The mass of the harmonic oscillator m was optimized for the control condition and scaled with λl0 across the range of muscle sizes tested (Ross et al. 2018). We also tested a range of muscle masses from 0.67 mg at λl0 of 1, which is typical of a muscle fascicle, to 670 g at λl0 of 100, which is less than the plantar flexor muscle mass in humans (Handsfield et al. 2014), for the mass model. When muscle is sub-maximally activated, not all the muscle fibers need to be recruited, and the intrinsic properties (fiber-type) of the activated fibers affect the mechanical output of the whole muscle. If the output of the whole muscle depended only on the combined output of the active fibers, for instance following Hill (1970), then recruiting slower fibers would result in lower power output than recruiting faster fibers. However, when the internal mass of the muscle is considered, the presence of inactive fibers provides inertial resistance, even when these inactive fibers are not developing active force (Holt et al. 2014; Ross and Wakeling 2016), and this may occlude differences due purely to the fiber-type. Thus, the root model would be expected to develop different power outputs depending on its fiber-type, with these differences seen across activation levels, however, the fiber-type contribution would be smaller for the model with internal mass, particularly at larger muscle sizes. This was tested by evaluating the muscle model with slower intrinsic properties with v0 of 5 s−1 and τact of 0.045 s, and at the lower levels of excitation with umax of 0.1 and 0.5. 3D simulations were achieved at umax of 0.1 and 1.0, and for increased base material stiffness (133% sbase), decreased density (50% ρ), and slower fibers than for a control simulation. The control simulation used fast intrinsic properties and sbase of 250. The numerical solution of the 3D model involves choosing several algorithmic parameters for the linear and nonlinear iterations. These in turn affect the maximal level to which the muscle can be activated in a numerically stable fashion, and these parameter choices depend on the stiffness of the base material. To simulate the 3D model at umax of 1.0, it was necessary to use a base material stiffness sbase that was higher than used by Rahemi et al. (2014). However, in that work the maximal activation achieved was considerably lower. Further mathematical investigation into the relations between the algorithmic parameters, base material stiffness, and activation levels is needed in addition to more extensive experimental measurements of this stiffness. Nonetheless, simulations for the 3D FEM were repeated for excitation frequencies of 1, 2, and 4 Hz. Results The muscle excitation u for each simulation was represented by a repeating square wave that took a value of either 0 or umax. The resulting activation a increased rapidly at the beginning of each cycle of u, plateaued across mid-burst, and finally decayed at a slower rate than the rate of activation when u returned to zero. The muscle generated active forces when a was greater than 0, and this caused the muscle to shorten and accelerate the point mass m within the harmonic oscillator (Fig. 1A). When the muscle was inactive, the restoring force due to the spring in the harmonic oscillator exceeded the muscle force, and caused m to accelerate to lengthen the muscle. The muscle length excursions were thus in-phase with the timing of u. However, there was a short time lag between when the muscle was excited and when it started to shorten due to both the time required for the muscle force to accelerate and change the direction of the velocity of m, and the time required to reach peak a after u changes to umax. Similarly, the muscle started to lengthen before it was fully deactivated due to the restoring force exceeding the muscle force toward the end of the activation period. As a consequence of the non-zero passive muscle and harmonic oscillator forces when the muscle was at l0, the values of the initial position and velocity of m required to solve the system were not immediately evident. However, because only the transient solution of a damped system, and not the steady-state solution, depends on the choice of initial conditions, we were able to avoid finding the ideal initial conditions by examining only the steady-state solution. Thus, herein we report only steady-state results where the behavior of the system repeats each cycle. The greatest average mass-specific muscle power output P* across all models and simulations was 42 W kg−1, and was achieved for the control simulation in the root model where f was 1 Hz, umax was 1, and the muscle had fast fiber-type properties. Rather than being of physiological importance, this was due to selecting the harmonic oscillator parameters and D to maximize P* for only these control conditions in the root model. For these conditions, the maximum shortening velocity reached was only 11% of v0, and the maximum stress was 66% of σ0. However, the theoretical maximum P* for the control simulation in the root model would have been 92 W kg−1, given these force–velocity and force–length properties, σ0, and v0 (following Weis-Fogh and Alexander 1977), so the control condition does not reflect the optimal conditions for maximum steady-state power production. The main effects of f, umax, and fiber-type on the Hill-based root, mass, and history-dependent models are shown in Fig. 4. The net mass-specific work per cycle W* was lower for higher cycle frequencies, and at the highest f of 8 Hz there was insufficient time for the muscle to fully activate or deactivate. This resulted in muscle force oscillations that were too fast to drive substantial displacements of the harmonic oscillator, and thus, there was little length excursion in the muscle for the high frequencies (Fig. 5A). The value of P* was greatest at the intermediate f of 1 Hz where the harmonic oscillator was calibrated to allow the muscle to generate maximal power. As umax decreased from 1.0 to 0.5 or 0.1, the activation and muscle force decreased. This led to the muscle operating at longer lengths with smaller excursions, as it was not able to generate large enough forces to contract to shorter lengths. Consequently, both W* and P* were smaller in magnitude for lower levels of umax (Figs. 4 and 5B). The fiber-type properties of the muscle were determined by v0 and τact, where the fast muscles had v0 of 10 s−1 and τact of 0.025 s, and the slow muscles had v0 of 5 s−1 and τact of 0.045 s. Muscles composed of the slower fiber-type developed lower active force during shortening leading to smaller length excursions and lower W* and P* (Figs. 4 and 5C). Fig. 4 View largeDownload slide Main effects of frequency, fiber-type, and excitation level on the mass-specific work W* and power P* from the 1D model during cycling contractions. The main effects were calculated from a set of 166 simulations across this parameter range, and included data from the mass-enhanced, history-dependent, and root Hill-type models. Fig. 4 View largeDownload slide Main effects of frequency, fiber-type, and excitation level on the mass-specific work W* and power P* from the 1D model during cycling contractions. The main effects were calculated from a set of 166 simulations across this parameter range, and included data from the mass-enhanced, history-dependent, and root Hill-type models. Fig. 5 View largeDownload slide Work-loops for muscle contractions demonstrating all tested parameters. (A–E) Work-loops are compared with the control simulation at 1 Hz frequency, maximal excitation, fast-fiber type, and small muscle size. Comparisons are shown that vary: (A) the frequency from 1 to 8 Hz; (B) the maximal excitation from 1.0 to 0.1; (C) the fiber type from fast to slow; (D) the muscle size from 0.67 mg to 670 g using the mass-enhanced model; and (E) history-dependent effects. (F) Results from the finite element simulations of the 3D model. Again, curves are compared with their control simulation at 1 Hz frequency, maximal excitation, and fast-fiber type with comparisons made to a model with greater base material stiffness, reduced muscle density, or slower fiber-type. The temporal direction of all work-loops is counterclockwise. Fig. 5 View largeDownload slide Work-loops for muscle contractions demonstrating all tested parameters. (A–E) Work-loops are compared with the control simulation at 1 Hz frequency, maximal excitation, fast-fiber type, and small muscle size. Comparisons are shown that vary: (A) the frequency from 1 to 8 Hz; (B) the maximal excitation from 1.0 to 0.1; (C) the fiber type from fast to slow; (D) the muscle size from 0.67 mg to 670 g using the mass-enhanced model; and (E) history-dependent effects. (F) Results from the finite element simulations of the 3D model. Again, curves are compared with their control simulation at 1 Hz frequency, maximal excitation, and fast-fiber type with comparisons made to a model with greater base material stiffness, reduced muscle density, or slower fiber-type. The temporal direction of all work-loops is counterclockwise. The non-dimensional dynamics of the muscle were identical across muscle length scales λl0 for the root Hill-type model that had no mass (i.e., stress, strain, W*, and P*). By contrast, inclusion of internal mass within the muscle model resulted in reduced net active muscle force during contraction, as well as increased passive force during lengthening due to the inertial costs of this mass. These changes in force were subtle, but were accentuated for simulations with greater λl0, with the larger mass-enhanced models generating lower W* and P*. Larger muscles and lower umax resulted in smaller differences in P* between fast and slow fiber-type conditions in the mass model compared with the massless root model. There was minimal active stretch in the muscle after the start of each cycle of u, so the contribution of force enhancement was minimal across all conditions. However, the muscle was largely active during shortening, resulting in significant force depression effects. The lower forces in the history-dependent model due to force depression resulted in the muscle not contracting to such short lengths, so both W* and P* were lower for the history-dependent model compared with the root model (Fig. 5E). The main effects of the different contractile and muscle parameters were mirrored in the 3D model. The 3D FEM was evaluated at umax of 0.1 and 1.0, and at f of 1, 2, and 4 Hz. There were large increases in W* and P* when umax was increased from 0.1 to 1.0, and decreases in W* and P* when the frequency increased from 1 Hz to 4 Hz. The effects of the base material stiffness, muscle density, and fiber-type were similar for both the low and high activation levels, with increased base material stiffness resulting in reduced W* and P*, decreased muscle density ρ (and thus mass) resulting in greater W* and P*, and the slower-fibered model resulting in lower W* and P* (Figs. 5F and 7). Overall, the greatest effects of the 1D mass-enhanced and history-dependent models can be observed for the contractile conditions with f of 1 Hz and umax of 1.0 (Fig. 5D, E). The mass-enhanced model at the largest scale tested (λl0 of 100; equivalent to a 670 g muscle) generated 12% less work and power with fast fiber-type properties and 4% less with slow properties compared with the root model with the same conditions. The history-dependent model was dominated by the force depression effect for these contractions and resulted in a 7% reduction in work and power compared with the root model. For the 3D model, there were additional substantial effects of the base material stiffness on W* and P*, and the mass effect was pronounced and important even at its 67 mg muscle mass. Discussion Muscle mass and scaling of contractile properties The root model was styled from a typical Hill-type muscle model, and as such, there were no mechanisms within the model to distort the intrinsic properties of the muscle as it scaled to different sizes. This implicitly means that simulations of muscles at any scale would result in the same non-dimensional muscle performance provided the external load from the harmonic oscillator was scaled appropriately. Indeed, this is how Hill-type muscle models are typically used: they are defined using intrinsic properties measured during experiments on single fibers or small muscles, which are in turn used to predict the function of whole muscles that may be of larger size (Wakeling and Lee 2011). The results from the root model in this study confirm that the function predicted by typical Hill-type models is conserved across scales. However, these models ignore the inertial effects of internal muscle mass on the dynamics of contraction. Previous studies have shown that internal mass acts to slow the rate of force development (Günther et al. 2012), and the maximum contraction velocity of larger muscles (Ross and Wakeling 2016). This study confirms these findings and also demonstrates that the internal mass results in a reduction in the work and power produced by whole muscles during cyclic contractions (Figs. 5 and 6). Fig. 6 View largeDownload slide Relative changes in performance between the 1D root model with the history-dependent and mass-enhanced 1D models. (A) Muscle mass-specific work W* and (B) muscle mass-specific power P* are shown relative to W* and P* for the control simulation (at 1 Hz frequency and maximal excitation). Simulations with different fiber-types were compared with their fiber-type control. Fig. 6 View largeDownload slide Relative changes in performance between the 1D root model with the history-dependent and mass-enhanced 1D models. (A) Muscle mass-specific work W* and (B) muscle mass-specific power P* are shown relative to W* and P* for the control simulation (at 1 Hz frequency and maximal excitation). Simulations with different fiber-types were compared with their fiber-type control. The mass of the muscle tissue alters the properties and function of muscle depending on the scale. Muscle mass varies in proportion to its volume, whereas muscle force varies with its cross-sectional area. This means that increasing muscle size, and decreasing the area to volume ratio, would result in the internal mass increasing faster than the available forces to accelerate that mass. Thus, the mass effects will predominate at larger muscle sizes. Indeed, the 1D models predicted a negligible effect of the internal mass for a scale of 1 where the muscle was 0.67 mg, whereas the mass effect was large and important for scale of 100 where the muscle was 670 g. However, the 3D model resulted in the internal mass being accelerated in all three dimensions, and here the mass effects were apparent even at its muscle size of 67 mg. The effect of the internal mass of the muscle tissue was important in this study due to the dynamic nature of the contractions. It should be noted that models of isometric contractions where the accelerations are negligible would show little difference if mass were included. However, internal work is required to accelerate the mass of the muscle tissue when the muscle changes length. Additionally, the 3D finite element representation of the muscle allows transverse deformations of the tissue and therefore transverse accelerations of the muscle mass, so the effect of the internal mass becomes even more important for the fully 3D model. Differences in muscle fiber-types result in different forces for given muscle lengths and velocities, and so we could expect the mass effect to depend on the intrinsic speed of the contracting fibers. Indeed, model results show that the maximum (normalized) shortening velocity of a muscle is reduced for larger sizes, and this effect is exacerbated for muscles with faster intrinsic speeds of their CEs (Ross and Wakeling 2016). Here we additionally show greater reductions in the mass-specific work and power output during cyclic contractions for muscles with faster intrinsic speeds (Fig. 6), which are caused by the dynamics of the muscle mass. Additional mass effects are apparent in muscles that are not fully active. Inactive muscle fibers have been proposed to contribute drag to a contracting whole muscle, resulting in a mechanical cost that would diminish the mechanical output of the muscle (Josephson and Edman 1988; Holt et al. 2014). The effects of the internal muscle mass would predominate at lower activations in a similar manner to which they predominate at larger muscle sizes, since at low activations there would be reduced contractile force available to accelerate the internal mass of the muscle tissue. Indeed, the results of this study support this notion. Muscle function and history-dependent properties When a muscle generates high power, it must contract with high force while shortening, and generate negligible force when it is lengthening. This would require the muscle to be activated close to its maximum length, and to be deactivated close to its minimum length, and this pattern emerged in our contraction cycles. However, activating close to maximum length provides minimal history-dependent force enhancement during active stretch. Indeed, the magnitude of force enhancement was negligible for the stretch-shortening contraction cycles in this study (Fig. 5E). By contrast, a substantial effect of force depression was observed. Thus, the majority of the history-dependent effects from these power-producing cyclic contractions occurred due to the force depression during shortening. In this study, the history-dependent effects resulted in a reduction in the net work and average power output per cycle of the muscle of about 7% for the 1 Hz frequency and maximal excitation condition. Although few studies have examined history-dependent effects during cyclic contractions, similar magnitudes of force depression have been previously reported. McGowan et al. (2010) found a 6% reduction in isometric torque during simulated knee extension and a 28% reduction in average crank power during simulated cycling when force depression was incorporated into a lower-body musculoskeletal model. This magnitude of joint torque reduction during knee extension is comparable to a 6.5% reduction in torque reported experimentally with the same conditions (Lee et al. 1999). Despite these findings, these musculoskeletal simulations involved more than one joint and multiple muscle actuators, and the Hill-type formulation of these actuators accounted for the effects of series elasticity. Thus, the contribution of force depression to the behavior of individual muscles in these simulations is difficult to resolve, so the results of our simulations are not directly comparable to the results of these previous studies. However, as our control simulation was designed to maximize muscle work and power output per cycle, the magnitude of force depression reported here is likely as large as would be achieved by a single muscle for contraction cycles during locomotion. Intrinsic speeds and fiber-type effects The maximum power output from a muscle depends on the intrinsic speed (maximum strain-rate) of its muscle fibers (Weis-Fogh and Alexander 1977). Not only can faster muscle fibers contract at faster speeds than slow fibers, but they generate greater stress at any given speed. This study used forward dynamic simulations of contractile cycles, whereby the range of motion, shortening velocity, and muscle force could all vary as a function of the excitation parameters used to drive the simulations. The simulations showed that the muscle with faster intrinsic speed (v0 of 10 versus 5 s−1) generated greater forces, allowing for greater length changes and thus higher shortening velocities than the slower muscle. The resultant work-loops showed that the faster muscle generated greater work per cycle and power output for an equivalent excitation and frequency (Fig. 5C). Mammalian muscle is commonly comprised of a heterogeneous population of fiber-types. If a heterogeneous muscle is sub-maximally activated then the intrinsic properties of the fibers that are activated may govern the contractile properties of the whole muscle. However, experimental studies have shown that there may be less difference in the shortening velocity of whole muscle than would be predicted by fiber-type alone (Holt et al. 2014). This effect has been attributed to the “drag” of the inactive tissue (Josephson and Edman 1988; Holt et al. 2014), and may be caused by the inertial properties of the internal muscle mass (Ross and Wakeling 2016), as discussed above. The pattern of recruitment between different types of motor units is therefore important in determining the mechanical output of whole muscle. The usual recruitment pattern follows the “size principle” whereby slower motor units are recruited first and for lower activations (Henneman et al. 1965a, 1965b). However, alternative recruitment patterns have been described that may favor recruitment of greater proportions of faster motor units at sub-maximal activations (see review: Hodson-Tole and Wakeling 2009) at fast contraction speeds and high contraction frequencies (Blake and Wakeling 2014). The simulations from this study show that, even at low excitation (10%), the type of fibers recruited alters the work and power output of the whole muscle. Additionally, our findings confirm that inertial effects from the muscle mass cause slower muscle shortening speeds and result in lower W* and P* than in smaller muscle or muscle without mass (Fig. 6). 3D structural, volumetric, and base material effects The contractile forces generated in the muscle are developed within the myofilaments in each muscle fiber from the cross-bridge forces between the actins and myosins, and from stretch in the PEEs such as the titins. The myofilaments are aligned with the longitudinal axis of the muscle fibers, and thus give rise to the magnitude and direction of the fiber force in this study. These forces are represented by the CE and the PEE in these models (Fig. 1). The fiber forces act to shorten the muscle fibers in their longitudinal direction. Muscle fibers are typically assumed to be isovolumetric (Baskin and Paolini 1967) and so must increase in girth to maintain their volume as they shorten. In the 3D model in this study, the volumetric constraints redistribute the longitudinal fiber forces to act in all three dimensions. This isovolumetric property is governed by incompressibility of the fluids that make up the intra and extracellular muscle tissue (muscle contains about 80% water; Van Loocke et al. 2008). Typical Hill-type muscle models (Zajac 1989) are 1D and contain no volume; however, in our 3D model the volumetric constraints act to maintain a nearly constant volume (changing by no >0.03% in these simulations). It should be noted that small changes in muscle volume have been seen at the fascicle level attributed to compressibility of the extracellular matrix (ECM) (Smith et al. 2011), and at the whole muscle level as the intramuscular pressure restricts the blood flow and blood volume (Ashton 1975), and so relaxation of the volumetric constraints may be warranted for future investigations. The 3D model represents the muscle as a fiber-reinforced composite material, where the fiber properties are governed by the myofilament contractile forces, and the composite properties represent the base material of the tissue. This base material contains properties across all muscle scales, including intracellular stiffness and properties from the ECM. Similar approaches have been employed in other 3D models of muscle where a variety of strain energy functions have been used to characterize the base material properties (e.g., Johansson et al. 2000; Blemker and Delp 2005). The volumetric constraints in the 3D model convert longitudinal shortening into a tendency for transverse expansion of the muscle. However, this expansion is resisted by the base material stiffness and so internal work must be done to achieve transverse expansion to allow for longitudinal shortening. Results from this study show that as the base material stiffness is increased, a greater magnitude of internal work is required for tissue deformation, resulting in reduced muscle external work and power (Figs. 5F and 7). Fig. 7 View largeDownload slide Changes in muscle mass-specific power P* between different 3D FEM simulations. The simulations were evaluated for 1 Hz cycle frequencies and maximal excitation. Models varied by increasing the scaling of their base material stiffness sbase, decreasing their tissue density ρ and thus mass, and by changing their fiber-type properties which were fast-fibered in the control model compared with a slow-fibered simulation. Note that being at a 1 Hz cycle frequency, the percentage changes in muscle mass specific work W* were the same as for these changes in P* for these simulations. Fig. 7 View largeDownload slide Changes in muscle mass-specific power P* between different 3D FEM simulations. The simulations were evaluated for 1 Hz cycle frequencies and maximal excitation. Models varied by increasing the scaling of their base material stiffness sbase, decreasing their tissue density ρ and thus mass, and by changing their fiber-type properties which were fast-fibered in the control model compared with a slow-fibered simulation. Note that being at a 1 Hz cycle frequency, the percentage changes in muscle mass specific work W* were the same as for these changes in P* for these simulations. The base material properties were represented in the model independently from the volumetric constraints, and indeed these were both independent from the density of the muscle tissue. However, the resultant properties of the whole muscle contraction emerged from the nonlinear interactions of these effects. Thus, the model provides an ideal framework for testing the specific effects of base material stiffness (or indeed ECM stiffness), volume changes, and mass properties. Such models inevitably benefit from thorough validation, however it is experimentally very challenging to independently vary any of these individual effects. Nonetheless, the mechanisms relating muscle shortening to its transverse expansion are consistent with a range of recent experimental studies: (A) transversely loading rat muscle results in reduced external force during isometric contractions (Siebert et al. 2014b), smaller transverse expansion, and greater internal work being done for transverse tissue deformation (Siebert et al. 2014a); (B) increasing the volume and intramuscular pressure within bullfrog muscle fibers by bathing them in hypotonic solution results in greater passive forces in the longitudinal direction (Sleboda and Roberts 2017). In this case the increased intramuscular pressure would result in greater transverse strain in the base material that would be redistributed by the volumetric constraints of the intracellular fluid, to result in greater resistance to longitudinal stretch due to the base material stiffness; and (C) resisting transverse expansion of leopard frog muscle by external physical constraint (to augment the ECM stiffness) causes the contracting muscle to shorten less in the longitudinal direction, and thus do less external work (Azizi et al. 2017). This effect is presumably due to the additional internal work that must be done to achieve the transverse deformation against the augmented base material stiffness. The 3D model in this study allowed us to test the effect of tissue density, and again, this is a difficult parameter to specifically vary in experiments. For instance, muscle density does vary with the fat content of the tissue, however, the adiposity also changes volumetric constraints and base material stiffness. Directly altering the density and thus tissue mass of the 3D model in this study showed that increased internal mass, and therefore additional internal work required for accelerating this mass, results in reduced external work and power from the muscle (Figs. 5F and 7). These mass effects were the same as those achieved in the mass-enhanced 1D Hill-type model and support the conclusion that the relative reductions in muscle work and power due to internal mass are more important at larger sizes and lower activation levels. The models in this study simulated contractions where the fibers were parallel to the line-of-action (longitudinal direction) of the muscle (Fig. 1). This allowed our study to focus on the contractile mechanisms in the fibers themselves rather than adding the geometric complexities that come with pennate architectures. However, the mechanisms and processes identified in the 3D model in this study will exist in any muscle fiber, regardless of the architecture of the muscle that contains it. Adjacent fibers interact between their shared boundaries, so transverse forces and deformations will propagate across the tissue. In pennate muscle, transverse fiber expansion must be balanced by increasing pennation angle in order for the muscles to maintain their spatial packing, with the aponeuroses that envelope the muscle influencing overall tissue deformation and changes in pennation. Indeed, we have previously used this 3D finite element framework (Rahemi et al. 2014, 2015) to show that: (A) stiffer aponeuroses lead to smaller increases in pennation, consistent with the reduction in changes in rat muscle gearing as it stiffens with age (Holt et al. 2016); (B) the curvatures that develop in contracting fibers are greatest near the aponeuroses, consistent with imaging studies in man (Namburete and Wakeling 2012); and (C) asymmetries in stress distribution across contracting muscles also lead to asymmetries in the transverse deformations of the muscle fibers and muscle belly (Rahemi et al. 2014). These findings are consistent with other experimental studies on muscle belly deformations (Azizi et al. 2008; Siebert et al. 2014b; Wakeling and Randhawa 2014), and provide a framework for understanding how the fundamental contractile properties of the fibers are related to tissue deformations, gearing, and the dynamics of whole muscle contraction. Assumptions, limitations, and future work The purpose of this study was to determine the effect of contractile features that are not accounted for in most Hill-type muscle models on the force and power of muscle during cyclic contractions. As such, we chose to orient the line-of-action of the fibers in these simulations along the longitudinal x-axis of the system (Fig. 1A). This enabled us to directly interpret the effects of the contractile properties on muscle performance because the entirety of the force vector from the fibers acted to accelerate the point mass m of the harmonic oscillator. However, not all muscles have their fibers arranged in parallel, and many have their fibers arranged at an oblique, pennate angle from the line-of-action of the muscle (Wickiewicz et al. 1983). As muscle fibers shorten they must increase in girth to maintain their volume. Such increases in girth allow pennate fibers to shorten and increase their pennation angle so that they can maintain their packing within the muscle belly (Lieber and Ward 2011). When the fibers rotate, their shortening velocity becomes uncoupled from that of the muscle belly in a process termed belly gearing (Wakeling et al. 2011), which lowers the fiber velocity relative to the belly velocity and alters muscle force, work, and power. Thus, architectural features of muscle, including pennation, add a layer of complexity to our understanding of muscle function in vivo. This understanding is further challenged when we consider that the exact nature of how muscles bulge in different directions (Maganaris et al. 1998; Azizi et al. 2008; Randhawa et al. 2013), and how variations in such bulging alters gearing (Azizi et al. 2008) is not fully known. Hence, focusing on parallel-oriented fibers allows us to probe fiber-level questions without being confounded by the additional complexity due to muscle architecture. We also chose to exclude tendinous, or series elastic, elements from the model for similar reasons. Tendons stretch when loaded, and if a tendon is long or compliant, this stretch can uncouple the velocity of the muscle belly from that of the muscle-tendon unit in a process known as tendon gearing (Wakeling et al. 2011). To make more direct interpretations of the effect of muscle (fiber or belly) shortening on the motion of the point mass within the harmonic oscillator, we assumed that tendon stretch was negligible and that the muscle was connected to m via a rigid link. However, it should be acknowledged that where significant tendinous tissue occurs in vivo, the tendon can further modify the behavior of the muscle-tendon-unit by amplifying muscle power output or attenuating muscle power input (Roberts and Azizi 2011). As the primary goal of this study was to examine the relative effect of different contractile phenomena on the work and power production of a muscle during cyclic contractions, the muscle model was designed to generate its maximum power for one of the intermediate contractile regimes with frequency f of 1 Hz, maximum excitation umax of 1, and fast fiber-type properties. Generating maximum power is an example of the muscle acting as a motor, and indeed our test condition operated at 50% of the maximum theoretical average mass-specific power per cycle (calculated for maximal activation and constant velocity) (Weis-Fogh and Alexander 1977). However, muscles can also act as springs, struts, or brakes during locomotion depending on the requirements of the task (Dickinson et al. 2000), and muscles may have specific adaptations to their properties depending on their primary function (Biewener 1998; Mörl et al. 2016). The muscle model in this study was coupled to a spring-mass-damper system that acted as a damped harmonic oscillator to provide an oscillating external load to allow for cyclic contractions of the muscle. Damped harmonic oscillators, by their nature, are only able to store and return or dissipate energy, and as a result, cannot be used to make the muscle behave like a brake. Therefore, the results of this study should be considered in the context of the power producing type of contractions that they were intended to simulate. As a consequence of only examining power-producing contraction cycles where the muscle is primarily active during shortening, force depression had a far greater effect than force enhancement in this study. However, these history-dependent effects are particularly challenging to represent mathematically due to a lack of consensus on the features of, and mechanisms behind, these effects, particularly during cyclic contractions with non-constant activations and length changes. Although force enhancement following active lengthening has been reported at sub-maximal levels of activation during voluntary (Oskouei and Herzog 2005, 2006a, 2006b; Pinniger and Cresswell 2007; Altenburg et al. 2008; Seiberl et al. 2012, 2013) and electrically stimulated (Abbott and Aubert 1952; Meijer 2002) contractions, the relationship between level of activation and magnitude of force enhancement is still poorly characterized. Furthermore, Meijer (2002) found that the magnitude of force enhancement was lower for lower levels of activation, and Oskouei and Herzog (2006b) only observed this phenomenon at sub-maximal activation in a portion of their study sample. Therefore, we chose to use 50% of maximum activation for the given simulation as a threshold value to indicate the start of force enhancement. However, as the excitation is prescribed as a square wave, the activation rises and decays rapidly so this threshold value would have little consequence on the magnitude of force enhancement for a whole cycle. We also chose to model force enhancement using a normalized elastic modulus that varied as a function of the instantaneous length at which this activation threshold was reached (Fig. 3A). To determine this relationship, we assumed that the enhanced muscle force as a function of muscle length (Fig. 3B) was linear and that the modulus depended only on the length where the activation threshold was reached and not on the current muscle length. However, further experimental studies examining force enhancement and depression during sub-maximal cyclic contractions are needed to determine if these assumptions are reasonable. In this study, we chose to model force depression as a function of muscle work, consistent with experimental evidence (De Ruiter et al. 1998; Herzog et al. 2000; Kosterina et al. 2008, 2009) and previous model formulations (McGowan et al. 2010, 2013). However, force depression has been shown to depend on contraction speed (Maréchal and Plaghki 1979; Meijer et al. 1997), and may decay during the period of time when the muscle is shortening (Till et al. 2014). More recently, Fortuna et al. (2018) found the magnitude of isometric, steady-state force depression following a stretch-shortening cycle to be independent of work. Thus, there is still debate as to whether or not force depression is accurately expressed as a function of only mechanical work. To include force depression across muscle model scales it was necessary to non-dimensionalize the force depression effect (Equation 7). Assuming force depression is work-dependent, the force would scale with muscle area, but work would scale with muscle volume in a similar manner to the mass scaling arguments discussed previously. Thus, a large muscle would generate a higher ratio of work to force, and at some limit, the contractile force would be depressed more than the force the muscle could generate. Hence, we assumed that the force depression per sarcomere would be independent of muscle size so that the force depression would result in an equivalent non-dimensional effect across muscle sizes. However, future experimental studies should be conducted to test the validity of this assumption. Hill-type muscle models typically do not use muscle-specific parameters but rather use parameter values taken from different studies on single fibers that vary widely depending on the muscle, species, temperature, preparation method (skinned versus intact), and stimulation protocol. The geometric dimensions and intrinsic properties of the muscle models in this study were chosen to reflect typical or average values reported in the literature so that our results could be more widely applicable than if we were to select values to represent a particular muscle. While we expect that our conclusions would hold if these parameter values were changed, this has not yet been evaluated and so should be considered when interpreting the results of this study. Using computational muscle models in this study allowed us to change features of the muscle that are difficult or impossible to manipulate experimentally. Although doing so provided a means to explore theoretical questions about muscle design and performance across muscles of different sizes under a range of conditions, it also limited our ability to directly validate our models against empirical data. Although we would expect that including more parameters in the muscle model formulation would provide more accurate predictions of in vivo muscle behavior, we have not explicitly tested this assumption. However, changing different model parameters, such as the fiber-type properties, led to changes in the muscle work and power output that were in a direction consistent with experimental evidence. Furthermore, the primary goal of this study was not to improve muscle model predictions, but rather to improve our understanding of the physiological mechanisms that may underlie muscle mechanical behavior during cyclic contractions. Conclusions This study examined a range of contractile properties that may contribute to the work and power output of muscle during cycling contractions. In our simulations, substantial history-dependent effects occurred for the power-producing contractions in these simulations, and were dominated by the force depression effect that occurs during active shortening. Our results show that the mechanical behavior of whole muscle during sub-maximal contractions depends on the fiber-type of the active fibers, even when the inertial effect of inactive muscle tissue is accounted for. There is an inertial cost to moving the tissue mass that results in reduced mass-specific work and power, and this is most pronounced for large muscles and low activations. 3D constraints also affect the mechanics of the muscle across all size scales. As the muscle changes length, the volumetric constraints cause it to undergo transverse deformations, and internal work against the base material stiffness must be done to deform the muscle tissue during contraction. Therefore, internal work must be done both to achieve tissue deformation and allow length changes, and to accelerate the muscle bulk. This study shows that the 3D effects, base material stiffness, and muscle mass all contribute to the dynamics of the muscle tissue, and also that these effects are size-dependent. These properties should thus be considered for a more complete understanding of muscle function during locomotor activities, where muscle activations are often sub-maximal and time-varying, muscles undergo unsteady cycles of length change, and whole muscles may be much larger than the single fibers from which their intrinsic properties have been determined. Acknowledgments The authors thank Dr. Emma Hodson-Tole for insightful discussions. Funding This work was supported by the National Institutes of Health [grant R01AR055648] and National Sciences and Engineering Research Council of Canada (NSERC) Discovery Grants [to both N.N. and J.M.W]. S.D. thanks the support of CONICYT-Chile through Becas Chile. References Abbott BC , Aubert XM. 1952 . The force exerted by active striated muscle during and after change of length . J Physiol 117 : 77 – 86 . Google Scholar CrossRef Search ADS PubMed Altenburg TM , de Ruiter CJ , Verdijk PW , Van Mechelen W , De Haan A. 2008 . Vastus lateralis surface and single motor unit EMG following submaximal shortening and lengthening contractions . Appl Physiol Nutr Metab 33 : 1086 – 95 . Google Scholar CrossRef Search ADS PubMed Altringham JD , Johnston IA. 1990 . Scaling effects on muscle function: power output of isolated fish muscle fibres performing oscillatory work . J Exp Biol 151 : 453 – 67 . Ashton H. 1975 . The effect of increased tissue pressure on blood flow . Clin Orthop Relat Res 113 : 15 – 26 . Google Scholar CrossRef Search ADS Azizi E , Brainerd EL , Roberts TJ. 2008 . Variable gearing in pennate muscles . Proc Natl Acad Sci U S A 105 : 1745 – 50 . Google Scholar CrossRef Search ADS PubMed Azizi E , Deslauriers AR , Holt NC , Eaton CE. 2017 . Resistance to radial expansion limits muscle strain and work . Biomech Model Mechanobiol 16 : 1633 – 43 . Google Scholar CrossRef Search ADS PubMed Baskin RJ , Paolini PJ. 1967 . Volume change and pressure development in muscle during contraction . Am J Physiol 213 : 1025 – 30 . Google Scholar PubMed Biewener AA. 1998 . Muscle function in vivo: a comparison of muscles used for elastic energy savings versus muscles used to generate mechanical power . Am Zool 38 : 703 – 17 . Google Scholar CrossRef Search ADS Biewener AA , Corning WR , Tobalske BW. 1998a . In vivo pectoralis muscle force–length behavior during level flight in pigeons (Columba livia) . J Exp Biol 201 : 3293 – 307 . Biewener AA , Konieczynski DD , Baudinette RV. 1998b . In vivo muscle force–length behavior during steady-speed hopping in tammar wallabies . J Exp Biol 201 : 1681 – 94 . Blake OM , Wakeling JM. 2014 . Early deactivation of slower muscle fibres at high movement frequencies . J Exp Biol 217 : 3528 – 34 . Google Scholar CrossRef Search ADS PubMed Blemker SS , Delp SL. 2005 . Three-dimensional representation of complex muscle architectures and geometries . Ann Biomed Eng 33 : 661 – 73 . Google Scholar CrossRef Search ADS PubMed Boettinger EG. 1957 . The machinery of insect flight. In: Scheer BT , editor. Recent advances in invertebrate physiology . Eugene : University of Oregon Publications . p. 117 – 42 . Cavagna GA , Citterio G. 1974 . Effect of stretching on the elastic characteristics and the contractile component of frog striated muscle . J Exp Biol 239 : 1 – 43 . Daley MA , Biewener AA. 2003 . Muscle force–length dynamics during level versus incline locomotion: a comparison of in vivo performance of two guinea fowl ankle extensors . J Exp Biol 206 : 2941 – 58 . Google Scholar CrossRef Search ADS PubMed Davydov D , Arndt D , Bangerth W , Heister T , Heltai L , Kronbichler M , Maier M , Pelteret JP , Turcksin B , Wells D. 2017 . The deal. II library, version 8.5 . J Numer Math 25 : 137 – 45 . De Ruiter CJ , De Haan AD , Jones DA , Sargeant AJ. 1998 . Shortening‐induced force depression in human adductor pollicis muscle . J Physiol 507 : 583 – 91 . Google Scholar CrossRef Search ADS PubMed Dick TJ , Biewener AA , Wakeling JM. 2017 . Comparison of human gastrocnemius forces predicted by Hill-type muscle models and estimated from ultrasound images . J Exp Biol 220 : 1643 – 53 . Google Scholar CrossRef Search ADS PubMed Dickinson MH , Farley CT , Full RJ , Koehl MAR , Kram R , Lehman S. 2000 . How animals move: an integrative view . Science 288 : 100 – 6 . Google Scholar CrossRef Search ADS PubMed Edman KAP , Elizinga G , Noble MIM. 1982 . Residual force enhancement after stretch of contracting frog single muscle fibers . J Gen Physiol 80 : 769 – 84 . Google Scholar CrossRef Search ADS PubMed Ettema GJ. 2002 . Effects of contraction history on control and stability in explosive actions . J Electromyogr Kinesiol 12 : 455 – 61 . Google Scholar CrossRef Search ADS PubMed Ettema GJ , Meijer K. 2000 . Muscle contraction history: modified Hill versus an exponential decay model . Biol Cybern 83 : 491 – 500 . Google Scholar CrossRef Search ADS PubMed Forcinito M , Epstein M , Herzog W. 1998 . Can a rheological muscle model predict force depression/enhancement? J Biomech 31 : 1093 – 9 . Google Scholar CrossRef Search ADS PubMed Fortuna R , Groeber M , Seiberl W , Power GA , Herzog W. 2017 . Shortening‐induced force depression is modulated in a time‐ and speed‐dependent manner following a stretch-shortening cycle . Physiol Rep published online (doi:10.14814/phy2.13279). Fortuna R , Kirchhuebel H , Seiberl W , Power GA , Herzog W. 2018 . Force depression following a stretch-shortening cycle is independent of stretch peak force and work performed during shortening . Sci Rep published online 8 (doi:10.1038/s41598-018-19657-8). Franklin CE , Johnston IA. 1997 . Muscle power output during escape responses in an Antarctic fish . J Exp Biol 200 : 703 – 12 . Google Scholar PubMed Günther M , Röhrle O , Haeufle DFB , Schmitt S. 2012 . Spreading out muscle mass within a Hill-type model: a computer simulation study. Comput Math Methods Med 2012 : 1 – 13 . Google Scholar CrossRef Search ADS Handsfield GG , Meyer CH , Hart JM , Abel MF , Blemker SS. 2014 . Relationships of 35 lower limb muscles to height and body mass quantified using MRI . J Biomech 47 : 631 – 8 . Google Scholar CrossRef Search ADS PubMed Henneman E , Somjen G , Carpenter DO. 1965 . Excitability and inhibitability of motoneurons of different sizes . J Neurophysiol 28 : 599 – 620 . Google Scholar CrossRef Search ADS PubMed Henneman E , Somjen G , Carpenter DO. 1965 . Functional significance of cell size in spinal motoneurons . J Neurophysiol 28 : 560 – 80 . Google Scholar CrossRef Search ADS PubMed Herzog W. 2014 . The role of titin in eccentric muscle contraction . J Exp Biol 217 : 2825 – 33 . Google Scholar CrossRef Search ADS PubMed Herzog W , Leonard TR. 2000 . The history dependence of force production in mammalian skeletal muscle following stretch-shortening and shortening-stretch cycles . J Biomech 33 : 531 – 42 . Google Scholar CrossRef Search ADS PubMed Herzog W , Leonard TR. 2002 . Force enhancement following stretching of skeletal muscle a new mechanism . J Exp Biol 205 : 1275 – 83 . Google Scholar PubMed Herzog W , Leonard T , Joumaa V , DuVall M , Panchangam A. 2012 . The three filament model of skeletal muscle stability and force production . Mol Cell Biomech 9 : 175 – 91 . Google Scholar PubMed Herzog W , Leonard TR , Wu JZ. 2000 . The relationship between force depression following shortening and mechanical work in skeletal muscle . J Biomech 33 : 659 – 68 . Google Scholar CrossRef Search ADS PubMed Hill AV. First and last experiments in muscle mechanics . UK : Cambridge University Press . p. 52 – 5 . Hisey B , Leonard TR , Herzog W. 2009 . Does residual force enhancement increase with increasing stretch magnitudes? J Biomech 42 : 1488 – 92 . Google Scholar CrossRef Search ADS PubMed Hodson-Tole EF , Wakeling JM. 2009 . Motor unit recruitment for dynamic tasks: current understanding and future directions . J Comp Physiol B 179 : 57 – 66 . Google Scholar CrossRef Search ADS PubMed Holt NC , Danos N , Roberts TJ , Azizi E. 2016 . Stuck in gear: age-related loss of variable gearing in skeletal muscle . J Exp Biol 219 : 998 – 1003 . Google Scholar CrossRef Search ADS PubMed Holt NC , Wakeling JM , Biewener AA. 2014 . The effect of fast and slow motor unit activation on whole-muscle mechanical performance: the size principle may not pose a mechanical paradox . Proc R Soc B 281 : 20140002. Google Scholar CrossRef Search ADS PubMed Huijing PA. 1999 . Muscle as a collagen fiber reinforced composite: a review of force transmission in muscle and whole limb . J Biomech 32 : 329 – 45 . Google Scholar CrossRef Search ADS PubMed Johansson T , Meier P , Blickhan R. 2000 . A finite-element model for the mechanical analysis of skeletal muscles . J Theor Biol 206 : 131 – 49 . Google Scholar CrossRef Search ADS PubMed Josephson RK. 1985 . Mechanical power output from striated muscle during cyclic contraction . J Exp Biol 114 : 493 – 512 . Josephson RK , Edman KAP. 1988 . The consequences of fibre heterogeneity on the force–velocity relation of skeletal muscle . Acta Physiol Scand 132 : 341 – 52 . Google Scholar CrossRef Search ADS PubMed Josephson RK , Stokes DR. 1989 . Strain, muscle length and work output in a crab muscle . J Exp Biol 145 : 45 – 61 . Kosterina N , Westerblad H , Eriksson A. 2009 . Mechanical work as predictor of force enhancement and force depression . J Biomech 42 : 1628 – 34 . Google Scholar CrossRef Search ADS PubMed Kosterina N , Westerblad H , Lännergren J , Eriksson A. 2008 . Muscular force production after concentric contraction . J Biomech 41 : 2422 – 9 . Google Scholar CrossRef Search ADS PubMed Lee SS , Arnold AS , de Boef Miara M , Biewener AA , Wakeling JM. 2013 . Accuracy of gastrocnemius muscles forces in walking and running goats predicted by one-element and two-element Hill-type models . J Biomech 46 : 2288 – 95 . Google Scholar CrossRef Search ADS PubMed Lee HD , Herzog W , Leonard T. 2001 . Effects of cyclic changes in muscle length on force production in in-situ cat soleus . J Biomech 34 : 979 – 87 . Google Scholar CrossRef Search ADS PubMed Lee HD , Suter E , Herzog W. 1999 . Force depression in human quadriceps femoris following voluntary shortening contractions . J Appl Physiol 87 : 1651 – 5 . Google Scholar CrossRef Search ADS PubMed Lieber RL , Ward SR. 2011 . Skeletal muscle design to meet functional demands . Philos Trans R Soc Lond B Biol Sci 366 : 1466 – 76 . Google Scholar CrossRef Search ADS PubMed Loeb GE , Hoffer JA , Pratt CA. 1985 . Activity of spindle afferents from cat anterior thigh muscles. I. Identification and patterns during normal locomotion . J Neurophysiol 54 : 549 – 64 . Google Scholar CrossRef Search ADS PubMed Machin KE , Pringle JW. 1960 . The physiology of insect fibrillar muscle. III. The effect of sinusoidal changes of length on a beetle flight muscle . Proc R Soc B 152 : 311 – 30 . Google Scholar CrossRef Search ADS Maganaris CN , Baltzopoulos V , Sargeant AJ. 1998 . In vivo measurements of the triceps surae complex architecture in man: implications for muscle function . J Physiol 512 : 603 – 14 . Google Scholar CrossRef Search ADS PubMed Maréchal G , Plaghki L. 1979 . The deficit of the isometric tetanic tension redeveloped after a release of frog muscle at a constant velocity . J Gen Physiol 73 : 453 – 67 . Google Scholar CrossRef Search ADS PubMed McGowan CP , Neptune RR , Herzog W. 2010 . A phenomenological model and validation of shortening-induced force depression during muscle contractions . J Biomech 43 : 449 – 54 . Google Scholar CrossRef Search ADS PubMed McGowan CP , Neptune RR , Herzog W. 2013 . A phenomenological muscle model to assess history dependent effects in human movement . J Biomech 46 : 151 – 7 . Google Scholar CrossRef Search ADS PubMed McGuigan MP , Yoo E , Lee DV , Biewener AA. 2009 . Dynamics of goat distal hind limb muscle-tendon function in response to locomotor grade . J Exp Biol 212 : 2092 – 104 . Google Scholar CrossRef Search ADS PubMed Meijer K. 2002 . History dependence of force production in submaximal stimulated rat medial gastrocnemius muscle . J Electromyogr Kinesiol 12 : 463 – 70 . Google Scholar CrossRef Search ADS PubMed Meijer K , Grootenboer HJ , Koopman BF , Huijing PA. 1997 . Fully isometric length–force curves of rat muscle differ from those during and after concentric contractions . J Appl Biomech 13 : 164 – 81 . Google Scholar CrossRef Search ADS Meijer K , Grootenboer HJ , Koopman HF , Van Der Linden BJ , Huijing PA. 1998 . A Hill type model of rat medial gastrocnemius muscle that accounts for shortening history effects . J Biomech 31 : 555 – 63 . Google Scholar CrossRef Search ADS PubMed Mörl F , Siebert T , Häufle D. 2016 . Contraction dynamics and function of the muscle–tendon complex depend on the muscle fibre–tendon length ratio: a simulation study . Biomech Model Mechanobiol 15 : 245 – 58 . Google Scholar CrossRef Search ADS PubMed Namburete AI , Wakeling JM. 2012 . Regional variations in fascicle curvatures within a muscle belly change during contraction . J Biomech 45 : 2835 – 40 . Google Scholar CrossRef Search ADS PubMed Nishikawa KC , Monroy JA , Uyeno TE , Yeo SH , Pai DK , Lindstedt SL. 2012 . Is titin a ‘winding filament’? A new twist on muscle contraction . Proc R Soc B 279 : 981 – 90 . Google Scholar CrossRef Search ADS PubMed Oskouei AE , Herzog W. 2005 . Observations on force enhancement in submaximal voluntary contractions of human adductor pollicis muscle . J Appl Physiol 98 : 2087 – 95 . Google Scholar CrossRef Search ADS PubMed Oskouei AE , Herzog W. 2006 . Force enhancement at different levels of voluntary contraction in human adductor pollicis . Eur J Appl Physiol 97 : 280 – 7 . Google Scholar CrossRef Search ADS PubMed Oskouei AE , Herzog W. 2006 . The dependence of force enhancement on activation in human adductor pollicis . Eur J Appl Physiol 98 : 22 – 9 . Google Scholar CrossRef Search ADS PubMed Perreault EJ , Heckman CJ , Sandercock TG. 2003 . Hill muscle model errors during movement are greatest within the physiologically relevant range of motor unit firing rates . J Biomech 36 : 211 – 8 . Google Scholar CrossRef Search ADS PubMed Pinniger GJ , Cresswell AG. 2007 . Residual force enhancement after lengthening is present during submaximal plantar flexion and dorsiflexion actions in humans . J Appl Physiol 102 : 18 – 25 . Google Scholar CrossRef Search ADS PubMed Rahemi H , Nigam N , Wakeling JM. 2014 . Regionalizing muscle activity causes changes to the magnitude and direction of the force from whole muscles . Front Physiol 5 : 298. Google Scholar CrossRef Search ADS PubMed Rahemi H , Nigam N , Wakeling JM. 2015 . The effect of intramuscular fat on skeletal muscle mechanics: implications for the elderly and obese . J R Soc Interface 12 : 20150365. Google Scholar CrossRef Search ADS PubMed Randhawa A , Jackman M , Wakeling JM. 2013 . Muscle gearing during isotonic and isokinetic movements in the ankle plantarflexors . Eur J Appl Physiol 113 : 437 – 47 . Google Scholar CrossRef Search ADS PubMed Richards CT , Clemente CJ. 2013 . Built for rowing: frog muscle is tuned to limb morphology to power swimming . J R Soc Interface 10 : 20130236. Google Scholar CrossRef Search ADS PubMed Roberts TJ , Azizi E. 2011 . Flexible mechanisms: the diverse roles of biological springs in vertebrate movement . J Exp Biol 214 : 353 – 61 . Google Scholar CrossRef Search ADS PubMed Roberts TJ , Marsh RL , Weyand PG , Taylor CR. 1997 . Muscular force in running turkeys: the economy of minimizing work . Science 275 : 1113 – 5 . Google Scholar CrossRef Search ADS PubMed Rode C , Siebert T , Blickhan R. 2009 . Titin-induced force enhancement and force depression: a ‘sticky-spring’ mechanism in muscle contractions? . J Theor Biol 259 : 350 – 60 . Google Scholar CrossRef Search ADS PubMed Ross SA , Nigam N , Wakeling JM. 2018 . A modelling approach for exploring muscle dynamics during cyclic contractions . PLoS Comput Biol 14 : e1006123. Google Scholar CrossRef Search ADS PubMed Ross SA , Wakeling JM. 2016 . Muscle shortening velocity depends on tissue inertia and level of activation during submaximal contractions . Biol Lett 12 : 20151041. Google Scholar CrossRef Search ADS PubMed Seiberl W , Hahn D , Herzog W , Schwirtz A. 2012 . Feedback controlled force enhancement and activation reduction of voluntarily activated quadriceps femoris during sub-maximal muscle action . J Electromyogr Kinesiol 22 : 117 – 23 . Google Scholar CrossRef Search ADS PubMed Seiberl W , Paternoster F , Achatz F , Schwirtz A , Hahn D. 2013 . On the relevance of residual force enhancement for everyday human movement . J Biomech 46 : 1996 – 2001 . Google Scholar CrossRef Search ADS PubMed Seiberl W , Power GA , Herzog W , Hahn D. 2015 . The stretch‐shortening cycle (SSC) revisited: residual force enhancement contributes to increased performance during fast SSCs of human m. adductor pollicis . Physiol Rep published online (doi:10.14814/phy2.12401). Siebert T , Till O , Blickhan R. 2014a . Work partitioning of transversally loaded muscle: experimentation and simulation . Comp Methods Biomech Biomed Eng 17 : 217 – 29 . Google Scholar CrossRef Search ADS Siebert T , Till O , Stutzig N , Günther M , Blickhan R. 2014b . Muscle force depends on the amount of transversal muscle loading . J Biomech 47 : 1822 – 8 . Google Scholar CrossRef Search ADS Simo JC , Taylor RL , Pister KS. 1985 . Variational and projection methods for the volume constraints in finite deformation elasto-plasticity . Comput Methods Appl Mech Eng 51 : 177 – 208 . Google Scholar CrossRef Search ADS Simo JC , Taylor RL. 1991 . Quasi-incompressible finite elasticity in principal stretches . Comput Methods Appl Mech Eng 85 : 273 – 310 . Google Scholar CrossRef Search ADS Sleboda DA , Roberts TJ. 2017 . Incompressible fluid plays a mechanical role in the development of passive muscle tension . Biol Lett 13 : 20160630. Google Scholar CrossRef Search ADS PubMed Smith LR , Gerace-Fowler L , Lieber RL. 2011 . Muscle extracellular matrix applies a transverse stress on fibers with axial strain . J Biomech 44 : 1618 – 20 . Google Scholar CrossRef Search ADS PubMed Tamura Y , Saito M , Nagato R. 2005 . A new motor model representing the stretch-induced force enhancement and shortening-induced force depression in skeletal muscle . J Biomech 38 : 877 – 84 . Google Scholar CrossRef Search ADS PubMed Till O , Siebert T , Blickhan R. 2014 . Force depression decays during shortening in the medial gastrocnemius of the rat . J Biomech 47 : 1099 – 103 . Google Scholar CrossRef Search ADS PubMed Van Leeuwen JL. 1992 . Muscle function in locomotion. In: Alexander RMcN , editors. Mechanics of animal locomotion . Berlin : Springer . p. 191 – 250 . Van Loocke M , Lyons CG , Simms CK. 2008 . Viscoelastic properties of passive skeletal muscle in compression: stress-relaxation behaviour and constitutive modelling . J Biomech 41 : 1555 – 66 . Google Scholar CrossRef Search ADS PubMed Wakeling JM , Blake OM , Wong I , Rana M , Lee SS. 2011 . Movement mechanics as a determinate of muscle structure, recruitment and coordination . Philos Trans R Soc B 366 : 1554 – 64 . Google Scholar CrossRef Search ADS Wakeling JM , Johnston IA. 1999 . Predicting muscle force generation during fast-starts for the common carp Cyprinus carpio . J Comp Physiol B 169 : 391 – 401 . Google Scholar CrossRef Search ADS Wakeling JM , Lee SS. 2011 . Modelling muscle forces: from scaled fibres to physiological task-groups . Procedia IUTAM 2 : 317 – 26 . Google Scholar CrossRef Search ADS Wakeling JM , Lee SSM , de Boef Miara M , Arnold AS , Biewener AA. 2012 . A muscle’s force depends on the recruitment patterns of its fibers . Ann Biomed Eng 40 : 1708 – 20 . Google Scholar CrossRef Search ADS PubMed Wakeling JM , Randhawa A. 2014 . Transverse strains in muscle fascicles during voluntary contraction: a 2D frequency decomposition of B-mode ultrasound images . Int J Biomed Imaging 2014 : 352910. Google Scholar CrossRef Search ADS PubMed Weis-Fogh T , Alexander RM. 1977 . The sustained power output from striated muscle. In: Pedley TJ , editor. Scale effects in animal locomotion . London : Academic Press . p. 511 – 25 . Weiss JA , Maker BN , Govindjee S. 1996 . Finite element implementation of incompressible, transversely isotropic hyperelasticity . Comput Methods Appl Mech Eng 135 : 107 – 28 . Google Scholar CrossRef Search ADS Wickiewicz TL , Roy RR , Powell PL , Edgerton VR. 1983 . Muscle architecture of the human lower limb . Clin Orthopaed Rel Res 179 : 275 – 83 . Google Scholar CrossRef Search ADS Williams TL. 2010 . A new model for force generation by skeletal muscle, incorporating work-dependent deactivation . J Exp Biol 213 : 643 – 50 . Google Scholar CrossRef Search ADS PubMed Williamson MR , Dial KP , Biewener AA. 2001 . Pectoralis muscle performance during ascending and slow level flight in mallards (Anas platyrhynchos) . J Exp Biol 204 : 495 – 507 . Google Scholar PubMed Wu JZ , Herzog W. 1999 . Modelling concentric contraction of muscle using an improved cross-bridge model . J Biomech 32 : 837 – 48 . Google Scholar CrossRef Search ADS PubMed Yeoh OH. 1993 . Some forms of the strain energy function for rubber . Rubber Chem Technol 66 : 754 – 71 . Google Scholar CrossRef Search ADS Zajac FE. 1989 . Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control . Crit Rev Biomed Eng 17 : 359 – 411 . Google Scholar PubMed Appendix 1 A mixed (P2–P1–P1) Finite Element Method was used to approximate the displacement, dilation, and pressure of the muscle in three dimensions. Piecewise continuous quadratic polynomials (P2) were used to approximate the displacement, whereas piecewise linear polynomials (P1) were used to approximate the pressure and dilation. The deformations of the model were governed by the balance of forces in the system, derived from the system’s total energy. The dilation was added as an unknown using a Lagrange multiplier to represent the internal pressure within the muscle. A three-field formulation was implemented in the library deal.II (Davydov et al. 2017), where the tutorial step-44 was taken as an initial starting point to solve the three resultant dynamic, nonlinear equations. The boundary conditions of the system were set as follows: we imposed zero displacement on some parts on the boundary, and set a traction (superficial stress in the normal direction) free condition on the rest of the boundary. However, a prescribed traction in the x-direction was used on part of this boundary to couple the 3D muscle system with the forces produced by the 1D oscillator. The initial displacement and pressure of the system were set to zero and the initial dilation was set to 1. Prior 3D muscle models using the Finite Element Method have assumed that the inertial effects within muscle are negligible by using quasi-static analysis (e.g., Rahemi et al. 2014). However, to examine the effects of inertia due to tissue mass, we accounted for the muscle’s density and velocity through its kinetic energy. The governing equations were solved using semi-implicit time stepping coupled with a finite element discretization of the space variables. The nonlinear nature of the equations implies that a nonlinear iteration must be applied at every time step to find the subsequent state. As the nonlinear iteration requires an initial guess to start, the converged solution found in the previous time step by the nonlinear solver was used to start the nonlinear iteration in the current state. For the first time step, the initial conditions for the displacement, pressure, and dilation were used as the initial guess. © The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Integrative and Comparative Biology. All rights reserved. For permissions please email: [email protected]. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Integrative and Comparative Biology Oxford University Press http://www.deepdyve.com/lp/oxford-university-press/size-history-dependent-activation-and-three-dimensional-effects-on-the-VbayDRqlkp

Loading next page...

References (114)

R. Alexander (2011)
Mechanics of Animal Locomotion
A. Namburete, J. Wakeling (2012)
Regional variations in fascicle curvatures within a muscle belly change during contraction.
Journal of biomechanics, 45 16
C. Richards, C. Clemente (2013)
Built for rowing: frog muscle is tuned to limb morphology to power swimming
Journal of The Royal Society Interface, 10
J. Wakeling, Avleen Randhawa (2014)
Transverse Strains in Muscle Fascicles during Voluntary Contraction: A 2D Frequency Decomposition of B-Mode Ultrasound Images
International Journal of Biomedical Imaging, 2014
Lucas Smith, L. Fowler-Gerace, R. Lieber (2011)
Muscle extracellular matrix applies a transverse stress on fibers with axial strain.
Journal of biomechanics, 44 8
J. Leeuwen (1992)
Muscle Function in Locomotion
(2000)
Muscle contraction history
H. Ashton (1975)
The effect of increased tissue pressure on blood flow.
Clinical orthopaedics and related research, 113
G. Ettema (2002)
Effects of contraction history on control and stability in explosive actions.
Journal of electromyography and kinesiology : official journal of the International Society of Electrophysiological Kinesiology, 12 6
E. Henneman, G. Somjen, D. Carpenter (1965)
FUNCTIONAL SIGNIFICANCE OF CELL SIZE IN SPINAL MOTONEURONS.
Journal of neurophysiology, 28
John Wu, Walter Herzog (1999)
Modelling concentric contraction of muscle using an improved cross-bridge model.
Journal of biomechanics, 32 8
(1977)
The sustained power output from striated muscle. In: Pedley TJ, editor. Scale effects in animal locomotion
Cavagna (1974)
Effect of stretching on the elastic characteristics and the contractile component of frog striated muscle
J Exp Biol, 239
G. Pinniger, A. Cresswell (2007)
Residual force enhancement after lengthening is present during submaximal plantar flexion and dorsiflexion actions in humans.
Journal of applied physiology, 102 1
(1998)
Can a rheological
S. Blemker, S. Delp (2005)
Three-Dimensional Representation of Complex Muscle Architectures and Geometries
Annals of Biomedical Engineering, 33
Youjiro Tamura, M. Saito, Rie Nagato (2005)
A new motor model representing the stretch-induced force enhancement and shortening-induced force depression in skeletal muscle.
Journal of biomechanics, 38 4
G. Handsfield, C. Meyer, J. Hart, M. Abel, S. Blemker (2014)
Relationships of 35 lower limb muscles to height and body mass quantified using MRI.
Journal of biomechanics, 47 3
G. Maréchal, L. Plaghki (1979)
The deficit of the isometric tetanic tension redeveloped after a release of frog muscle at a constant velocity
The Journal of General Physiology, 73
H. lee, W. Herzog, T. Leonard (2001)
Effects of cyclic changes in muscle length on force production in in-situ cat soleus.
Journal of biomechanics, 34 8
W. Seiberl, G. Power, W. Herzog, D. Hahn (2015)
The stretch-shortening cycle (SSC) revisited: residual force enhancement contributes to increased performance during fast SSCs of human m. adductor pollicis
Physiological Reports, 3
M. Günther, O. Röhrle, D. Haeufle, S. Schmitt (2012)
Spreading out Muscle Mass within a Hill-Type Model: A Computer Simulation Study
Computational and Mathematical Methods in Medicine, 2012
M. Dickinson, C. Farley, R. Full, M. Koehl, R. Kram, S. Lehman (2000)
How animals move: an integrative view.
Science, 288 5463
K. Machin, J. Pringle (1960)
The physiology of insect fibrillar muscle III. The effect of sinusoidal changes of length on a beetle flight muscle
Proceedings of the Royal Society of London. Series B. Biological Sciences, 152
O. Till, T. Siebert, R. Blickhan (2014)
Force depression decays during shortening in the medial gastrocnemius of the rat.
Journal of biomechanics, 47 5
J. Wakeling, Sabrina Lee, A. Arnold, Maria Miara, A. Biewener (2012)
A Muscle’s Force Depends on the Recruitment Patterns of Its Fibers
Annals of Biomedical Engineering, 40
R. Lieber, S. Ward (2011)
Skeletal muscle design to meet functional demands
Philosophical Transactions of the Royal Society B: Biological Sciences, 366
Abbott (1952)
The force exerted by active striated muscle during and after change of length
J Physiol, 117
T. Wickiewicz, R. Roy, P. Powell, V. Edgerton (1983)
Muscle architecture of the human lower limb.
Clinical orthopaedics and related research, 179
G. Ettema, K. Meijer (2000)
Muscle contraction history: modified Hill versus an exponential decay model
Biological Cybernetics, 83
C. Maganaris, V. Baltzopoulos, A. Sargeant (1998)
In vivo measurements of the triceps surae complex architecture in man: implications for muscle function
The Journal of Physiology, 512
C. Rode, T. Siebert, R. Blickhan (2009)
Titin-induced force enhancement and force depression: a 'sticky-spring' mechanism in muscle contractions?
Journal of theoretical biology, 259 2
A. Biewener, WR Corning, B. Tobalske (1998)
In vivo pectoralis muscle force-length behavior during level flight in pigeons (Columba livia)
The Journal of experimental biology, 201 (Pt 24)
F. Mörl, T. Siebert, Daniel Häufle (2015)
Contraction dynamics and function of the muscle-tendon complex depend on the muscle fibre-tendon length ratio: a simulation study
Biomechanics and Modeling in Mechanobiology, 15
M. Williamson, K. Dial, A. Biewener (2001)
Pectoralis muscle performance during ascending and slow level flight in mallards (Anas platyrhynchos).
The Journal of experimental biology, 204 Pt 3
M. Mcguigan, E. Yoo, David Lee, A. Biewener (2009)
Dynamics of goat distal hind limb muscle–tendon function in response to locomotor grade
Journal of Experimental Biology, 212
J. Altringham, I. Johnston (1990)
Scaling Effects on Muscle Function: Power Output of Isolated Fish Muscle Fibres Performing Oscillatory Work
The Journal of Experimental Biology, 151
I. Johnston, J. Altringham (1988)
Muscle function in locomotion
Nature, 335
G. Loeb, J. Hoffer, C. Pratt (1985)
Activity of spindle afferents from cat anterior thigh muscles. I. Identification and patterns during normal locomotion.
Journal of neurophysiology, 54 3
J. Simo, R. Taylor (1991)
Quasi-incompressible finite elasticity in principal stretches. Continuum basis and numerical algorithms
Applied Mechanics and Engineering, 85
R. Fortuna, Hannah Kirchhuebel, W. Seiberl, G. Power, W. Herzog (2018)
Force depression following a stretch-shortening cycle is independent of stretch peak force and work performed during shortening
Scientific Reports, 8
E. Henneman, G. Somjen, D. Carpenter (1965)
Excitability and inhibitability of motoneurons of different sizes.
Journal of neurophysiology, 28 3
E. Azizi, E. Brainerd, T. Roberts (2008)
Variable gearing in pennate muscles
Proceedings of the National Academy of Sciences, 105
N. Kosterina, H. Westerblad, J. Lännergren, A. Eriksson (2008)
Muscular force production after concentric contraction.
Journal of biomechanics, 41 11
W. Herzog, T. Leonard (2000)
The history dependence of force production in mammalian skeletal muscle following stretch-shortening and shortening-stretch cycles.
Journal of biomechanics, 33 5
Hadi Rahemi, N. Nigam, J. Wakeling (2015)
The effect of intramuscular fat on skeletal muscle mechanics: implications for the elderly and obese
Journal of The Royal Society Interface, 12
J. Murray, T. Pedley (1978)
Scale Effects in Animal Locomotion.
Biometrics, 34
E. Hodson-Tole, J. Wakeling (2008)
Motor unit recruitment for dynamic tasks: current understanding and future directions
Journal of Comparative Physiology B, 179
W. Herzog (2014)
The role of titin in eccentric muscle contraction
Journal of Experimental Biology, 217
S. Ross, J. Wakeling (2016)
Muscle shortening velocity depends on tissue inertia and level of activation during submaximal contractions
Biology Letters, 12
R. Josephson, Stokes Darrell (1989)
Strain, Muscle Length and Work Output in a Crab Muscle
The Journal of Experimental Biology, 145
J. Simo, R. Taylor, K. Pister (1985)
Variational and projection methods for the volume constraint in finite deformation elasto-plasticity
Computer Methods in Applied Mechanics and Engineering, 51
M. Loocke, C. Lyons, C. Simms (2008)
Viscoelastic properties of passive skeletal muscle in compression: stress-relaxation behaviour and constitutive modelling.
Journal of biomechanics, 41 7
W. Herzog, T. Leonard, V. Joumaa, M. Duvall, A. Panchangam (2012)
The three filament model of skeletal muscle stability and force production.
Molecular & cellular biomechanics : MCB, 9 3
(1957)
The machinery of insect flight
J. Wakeling, O. Blake, Iris Wong, Manku Rana, Sabrina Lee (2011)
Movement mechanics as a determinate of muscle structure, recruitment and coordination
Philosophical Transactions of the Royal Society B: Biological Sciences, 366
E. Perreault, C. Heckman, T. Sandercock (2003)
Hill muscle model errors during movement are greatest within the physiologically relevant range of motor unit firing rates.
Journal of biomechanics, 36 2
Brandon Hisey, T. Leonard, W. Herzog (2009)
Does residual force enhancement increase with increasing stretch magnitudes?
Journal of biomechanics, 42 10
Avleen Randhawa, Meghan Jackman, J. Wakeling (2013)
Muscle gearing during isotonic and isokinetic movements in the ankle plantarflexors
European Journal of Applied Physiology, 113
A. Biewener (1998)
Muscle Function in vivo: A Comparison of Muscles used for Elastic Energy Savings versus Muscles Used to Generate Mechanical Power1
Integrative and Comparative Biology, 38
Craig McGowan, R. Neptune, Walter Herzog (2013)
A phenomenological muscle model to assess history dependent effects in human movement.
Journal of biomechanics, 46 1
T. Roberts, R. Marsh, P. Weyand, C. Taylor (1997)
Muscular Force in Running Turkeys: The Economy of Minimizing Work
Science, 275
R. Baskin, P. Paolini (1967)
Volume change and pressure development in muscle during contraction.
The American journal of physiology, 213 4
G. Cavagna, G. Citterio (1974)
Effect of stretching on the elastic characteristics and the contractile component of frog striated muscle
The Journal of Physiology, 239
T. Williams (2010)
A new model for force generation by skeletal muscle, incorporating work-dependent deactivation
Journal of Experimental Biology, 213
T. Siebert, O. Till, N. Stutzig, M. Günther, R. Blickhan (2014)
Muscle force depends on the amount of transversal muscle loading.
Journal of biomechanics, 47 8
Ali Oskouei, W. Herzog (2006)
Force enhancement at different levels of voluntary contraction in human adductor pollicis
European Journal of Applied Physiology, 97
Franklin, Johnston (1997)
Muscle power output during escape responses in an Antarctic fish
The Journal of experimental biology, 200 Pt 4
W. Seiberl, D. Hahn, W. Herzog, A. Schwirtz (2012)
Feedback controlled force enhancement and activation reduction of voluntarily activated quadriceps femoris during sub-maximal muscle action.
Journal of electromyography and kinesiology : official journal of the International Society of Electrophysiological Kinesiology, 22 1
P. Huijing (1999)
Muscle as a collagen fiber reinforced composite: a review of force transmission in muscle and whole limb.
Journal of biomechanics, 32 4
T. Altenburg, C. Ruiter, P. Verdijk, W. Mechelen, A. Haan (2008)
Vastus lateralis surface and single motor unit EMG following submaximal shortening and lengthening contractions.
Applied physiology, nutrition, and metabolism = Physiologie appliquee, nutrition et metabolisme, 33 6
Ali Oskouei, W. Herzog (2005)
Observations on force enhancement in submaximal voluntary contractions of human adductor pollicis muscle.
Journal of applied physiology, 98 6
E. Azizi, A. Deslauriers, N. Holt, C. Eaton (2017)
Resistance to radial expansion limits muscle strain and work
Biomechanics and Modeling in Mechanobiology, 16
R. Josephson (1985)
Mechanical Power output from Striated Muscle during Cyclic Contraction
The Journal of Experimental Biology, 114
David Sleboda, T. Roberts (2017)
Incompressible fluid plays a mechanical role in the development of passive muscle tension
Biology Letters, 13
(1977)
The sustained power output from striated muscle
W. Bangerth, T. Heister, L. Heltai, G. Kanschat, M. Kronbichler, Matthias Maier, Bruno Turcksin, Toby Young (2013)
The deal.II library, version 8.5
Journal of Numerical Mathematics, 25
K. Meijer (2002)
History dependence of force production in submaximal stimulated rat medial gastrocnemius muscle.
Journal of electromyography and kinesiology : official journal of the International Society of Electrophysiological Kinesiology, 12 6
Kenneth Meijer, H. Grootenboer, H. Koopman, B. Linden, P. Huijing (1998)
A Hill type model of rat medial gastrocnemius muscle that accounts for shortening history effects.
Journal of biomechanics, 31 6
S. Ross, N. Nigam, J. Wakeling (2018)
A modelling approach for exploring muscle dynamics during cyclic contractions
PLoS Computational Biology, 14
M. Daley, A. Biewener (2003)
Muscle force-length dynamics during level versus incline locomotion: a comparison of in vivo performance of two guinea fowl ankle extensors
Journal of Experimental Biology, 206
Hae-Dong Lee, E. Suter, Walter Herzog (1999)
Force depression in human quadriceps femoris following voluntary shortening contractions.
Journal of applied physiology, 87 5
O. Blake, J. Wakeling (2014)
Early deactivation of slower muscle fibres at high movement frequencies
Journal of Experimental Biology, 217
C. McGowan, R. Neptune, W. Herzog (2010)
A phenomenological model and validation of shortening-induced force depression during muscle contractions.
Journal of biomechanics, 43 3
C. Ruiter, A. Haan, D. Jones, A. Sargeant (1998)
Shortening‐induced force depression in human adductor pollicis muscle
The Journal of Physiology, 507
T. Dick, A. Biewener, J. Wakeling (2017)
Comparison of human gastrocnemius forces predicted by Hill-type muscle models and estimated from ultrasound images
Journal of Experimental Biology, 220
J. Wakeling, I. Johnston (1999)
Predicting muscle force generation during fast-starts for the common carp Cyprinus carpio
Journal of Comparative Physiology B, 169
Mechanical effects on muscle work-loops 249
J. Weiss, B. Maker, S. Govindjee (1996)
Finite element implementation of incompressible, transversely isotropic hyperelasticity
Computer Methods in Applied Mechanics and Engineering, 135
A. Hill (1970)
First and Last Experiments in Muscle Mechanics
K. Nishikawa, J. Monroy, Theodore Uyeno, S. Yeo, D. Pai, S. Lindstedt (2012)
Is titin a ‘winding filament’? A new twist on muscle contraction
Proceedings of the Royal Society B: Biological Sciences, 279
Mario Forcinito, Marcelo Epstein, Walter Herzog (1998)
Can a rheological muscle model predict force depression/enhancement?
Journal of biomechanics, 31 12
(1974)
Effect of stretching on
F. Zajac (1989)
Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control.
Critical reviews in biomedical engineering, 17 4
Walter Herzog, Tim Leonard (2002)
Force enhancement following stretching of skeletal muscle: a new mechanism.
The Journal of experimental biology, 205 Pt 9
Sabrina Lee, A. Arnold, Maria Miara, A. Biewener, J. Wakeling (2013)
Accuracy of gastrocnemius muscles forces in walking and running goats predicted by one-element and two-element Hill-type models.
Journal of biomechanics, 46 13
Walter Herzog, Tim Leonard, John Wu (2000)
The relationship between force depression following shortening and mechanical work in skeletal muscle.
Journal of biomechanics, 33 6
T. Roberts, E. Azizi (2011)
Flexible mechanisms: the diverse roles of biological springs in vertebrate movement
Journal of Experimental Biology, 214
R. Josephson, K. Edman (1988)
The consequences of fibre heterogeneity on the force-velocity relation of skeletal muscle.
Acta physiologica Scandinavica, 132 3
N. Holt, J. Wakeling, A. Biewener (2014)
The effect of fast and slow motor unit activation on whole-muscle mechanical performance: the size principle may not pose a mechanical paradox
Proceedings of the Royal Society B: Biological Sciences, 281
K. Meijer, H. Grootenboer, Bart Koopman, P. Huijing (1997)
Fully Isometric Length-Force Curves of Rat Muscle Differ from those during and after Concentric Contractions
Journal of Applied Biomechanics, 13
T. Johansson, P. Meier, Reinhard Blickhan (2000)
A finite-element model for the mechanical analysis of skeletal muscles.
Journal of theoretical biology, 206 1
Ali Oskouei, W. Herzog (2006)
The dependence of force enhancement on activation in human adductor pollicis
European Journal of Applied Physiology, 98
N. Kosterina, H. Westerblad, A. Eriksson (2009)
Mechanical work as predictor of force enhancement and force depression.
Journal of biomechanics, 42 11
Hadi Rahemi, N. Nigam, J. Wakeling (2014)
Regionalizing muscle activity causes changes to the magnitude and direction of the force from whole muscles—a modeling study
Frontiers in Physiology, 5
R. Fortuna, Martin Groeber, W. Seiberl, G. Power, W. Herzog (2017)
Shortening‐induced force depression is modulated in a time‐ and speed‐dependent manner following a stretch–shortening cycle
Physiological Reports, 5
W. Seiberl, F. Paternoster, F. Achatz, A. Schwirtz, D. Hahn (2013)
On the relevance of residual force enhancement for everyday human movement.
Journal of biomechanics, 46 12
K. Edman, G. Elzinga, M. Noble (1982)
Residual force enhancement after stretch of contracting frog single muscle fibers
The Journal of General Physiology, 80
O. Yeoh (1993)
Some Forms of the Strain Energy Function for Rubber
Rubber Chemistry and Technology, 66
T. Siebert, O. Till, R. Blickhan (2014)
Work partitioning of transversally loaded muscle: experimentation and simulation
Computer Methods in Biomechanics and Biomedical Engineering, 17
N. Holt, N. Danos, T. Roberts, E. Azizi (2016)
Stuck in gear: age-related loss of variable gearing in skeletal muscle
Journal of Experimental Biology, 219
Andrew Biewener, D. Konieczynski, R. Baudinette (1998)
In vivo muscle force-length behavior during steady-speed hopping in tammar wallabies.
The Journal of experimental biology, 201 Pt 11
J. Wakeling, Sabrina Lee (2011)
Modelling muscle forces: from scaled fibres to physiological task-groups
Procedia IUTAM, 2
By Alexandrowicz (1953)
Recording work done at the Plymouth Laboratory
Journal of the Marine Biological Association of the United Kingdom, 31

Publisher: Oxford University Press
Copyright: © The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Integrative and Comparative Biology. All rights reserved. For permissions please email: [email protected].
ISSN: 1540-7063
eISSN: 1557-7023
DOI: 10.1093/icb/icy021
Publisher site: See Article on Publisher Site

Abstract

Abstract Muscles undergo cycles of length change and force development during locomotion, and these contribute to their work and power production to drive body motion. Muscle fibers are typically considered to be linear actuators whose stress depends on their length, velocity, and activation state, and whose properties can be scaled up to explain the function of whole muscles. However, experimental and modeling studies have shown that a muscle’s stress additionally depends on inactive and passive tissues within the muscle, the muscle’s size, and its previous contraction history. These effects have not been tested under common sets of contraction conditions, especially the cyclic contractions that are typical of locomotion. Here we evaluate the relative effects of size, history-dependent, activation and three-dimensional effects on the work and power produced during cyclic contractions of muscle models. Simulations of muscle contraction were optimized to generate high power outputs: this resulted in the muscle models being largely active during shortening, and inactive during lengthening. As such, the history-dependent effects were dominated by force depression during simulated active shortening rather than force enhancement during active stretch. Internal work must be done to deform the muscle tissue, and to accelerate the internal muscle mass, resulting in reduced power and work that can be done on an external load. The effect of the muscle mass affects the scaling of muscle properties, with the inertial costs of contraction being relatively greater at larger sizes and lower activation levels. Introduction Muscles typically undergo cycles of activation, length change, and force development during locomotion (Dickinson et al. 2000). An understanding of how muscles operate during these cyclical regimes is thus necessary to determine the muscle design criteria for achieving locomotor tasks in a manner that is metabolically efficient or mechanically advantageous in terms of work and power. Pioneering studies on insect flight muscle in the 1950s characterized muscle forces during high frequency bursts of excitation coupled with cyclic length changes by expressing muscle force as a function of muscle length in a manner now known as a “work-loop” (Boettinger 1957; Machin and Pringle 1960). Boettinger (1957) noted that the slope of the major axis of the loop represented the compliance of the excited muscle, the area enclosed in the loop represented the work done against damping forces in the mechanical system, and the position of the loop in the force–length area represented the magnitude of contraction velocity. These early studies on asynchronous muscle, where the excitation is not synchronized with the contraction cycles, were later extended in the 1980s to synchronous insect flight muscle (Josephson 1985) and crustacean muscle (Josephson and Stokes 1989), where the timing of the excitation or action potentials is related to phase of the contraction cycle. This work-loop approach has since been generalized to a range of vertebrate muscles. Understanding the determinants of muscle performance during vertebrate locomotion requires measures of both muscle force and length changes during in vivo activities. Prior to the 1980s when sonomicrometry was introduced to measure in vivo lengths and tendon buckles to measure in-series tendon forces (Loeb et al. 1985), work-loop studies were primarily done in situ or in vitro with imposed artificial length changes and excitation patterns. These new techniques to measure in vivo function fueled a proliferation of studies applying the work-loop technique in fish muscle (Altringham and Johnston 1990) to mimic the deformations observed in vivo (Franklin and Johnston 1997), and subsequently, in vivo work-loops have been determined for a range of species including pigeons (Biewener et al. 1998a), turkeys (Roberts et al. 1997), guinea fowl (Daley and Biewener 2003), mallard ducks (Williamson et al. 2001), wallabies (Biewener et al. 1998b), and goats (McGuigan et al. 2009). Nonetheless, we still know very little about the in vivo behavior of contracting muscle for a wide range of locomotor activities in the majority of muscles in the majority of species, as well as the mechanisms that underlie this behavior. Hill-type muscle models are widely used in comparative biomechanics to estimate in vivo muscle forces where direct measures are not feasible (Van Leeuwen 1992; Wakeling and Johnston 1999; Hodson-Tole and Wakeling 2009; Lee et al. 2013; Richards and Clemente 2013), and because of our limited understanding of in vivo whole muscle mechanics, the intrinsic properties of these models are often derived from experiments on single fibers or fiber bundles that are maximally activated and constrained to contract in a steady manner. As a consequence, the few studies that have validated Hill-type models against in vivo or in situ measures have found that these models perform best during maximal, steady-state contractions and suffer under more functionally-relevant conditions where activations and muscle lengths and velocities are not constant. Perreault et al. (2003) found that model errors (percent root-mean-square) can exceed 50% when large changes in muscle length are modeled in situ. More recently, we have achieved modeling errors of 32% for in vivo locomotion in goats (Lee et al. 2013), and 13.5% in humans (Dick et al. 2017), however, correlation between measured and predicted forces only reached a maximum of r2 ≈ 0.5. Therefore, assuming that whole muscles behave as single fibers during controlled, steady contractions fails to account for the complex mechanisms that underlie whole muscle behavior in vivo. Hill-type muscle models typically assume that muscles are composed of a single fiber-type and contain no internal mass or inertial effects, and so can be scaled to any size without any consequences on the kinematic and dynamic behavior of the model (Ross et al. 2018). However, studies have shown that the presence of different fiber-types is important for muscle force production (Wakeling et al. 2012; Lee et al. 2013), and that muscle mass provides inertial damping within the muscle (Günther et al. 2012; Ross and Wakeling 2016) that is dependent on fiber-type and activation level (Holt et al. 2014). The mass within muscle tissue takes time to move when the muscle forces cause accelerations, and this has been shown to slow contraction times (Günther et al. 2012) and maximum effective contraction speeds (Ross and Wakeling 2016) in isotonic shortening contractions. Although the implications of these inertial effects due to muscle tissue mass on muscle performance during cyclic contractions are not yet known, we would expect the inclusion of internal muscle mass to reduce the shortening velocity of the muscle, and thus the power generated. Most Hill-type muscle models calculate muscle force as a function of purely the instantaneous activation state, muscle length, and velocity, and do not consider previous contraction history. However, history-dependent effects, which are distinct from time-dependent effects due to excitation-activation dynamics, may be important for muscle contractile behavior. In particular, when compared with isometric forces at a given length, muscle forces may be elevated or enhanced when the muscle reaches that length with prior active lengthening (Abbott and Aubert 1952; Cavagna and Citterio 1974; Edman et al. 1982; Herzog and Leonard 2002; Hisey et al. 2009), and may be reduced or depressed when the muscle reaches that length with prior active shortening (Abbott and Aubert 1952; Maréchal and Plaghki 1979; Meijer et al. 1998; Herzog et al. 2000). The few experimental studies examining history-dependent effects during cyclic contractions have been limited to isokinetic stretch-shortening and shortening-stretch regimes with constant, maximal activation in situ in cats (Herzog and Leonard 2000; Lee et al. 2001) or in vivo in humans (Seiberl et al. 2015; Fortuna et al. 2017, 2018). Furthermore, while a number of muscle models have incorporated history-dependent effects (Forcinito et al. 1998; Meijer et al. 1998; Wu and Herzog 1999; Ettema and Meijer 2000; Ettema 2002; Tamura et al. 2005; McGowan et al. 2010, 2013; Williams 2010; Nishikawa et al. 2012), none have examined the contribution of these effects during cyclic contractions with unconstrained length changes and time-varying activation. Muscle is nearly isovolumetric (Baskin and Paolini 1967) and hence when fibers shorten they must expand in girth to maintain their volume. The extent to which such transverse deformations occur is resisted by stiffness both within and external to the muscle fibers (Huijing 1999). Transverse stiffness in the muscle resists the transverse expansion, and the volumetric constraints result in this being seen as a resistance of the fibers to shortening (Azizi et al. 2017). These three-dimensional (3D) phenomena are also not considered by typical Hill-type models (Zajac 1989) that treat muscles as linear actuators with only 1D properties. However, we would expect that the inclusion of base material stiffness within the muscle, with properties in all three dimensions, would decrease the power available to work against an external load. Despite the wealth of experimental and modeling evidence that these additional effects are important for in vivo muscle behavior, they have only yet been evaluated in separate studies with specific contraction regimes to elucidate their specific effects, so their individual and relative contributions are still not known. Therefore, the purpose of this study was to evaluate the contributions of internal mass, fiber-type, history-dependent, and 3D effects to the work and power output of muscle during cyclic contractions. To achieve this, we used a modeling framework that drives cyclic contractions of a Hill-type muscle model in series with a damped harmonic oscillator to create muscle work-loops (Ross et al. 2018). These additional features were implemented within the muscle model, and so could be tested under common contractile regimes with the same intrinsic muscle properties to determine their individual and relative contributions to cyclic muscle performance. Methods System and muscle properties The modeling framework we used to test the effects of different contractile and mechanical properties on muscle behavior was composed of a harmonic oscillator, that represented an external load placed on the muscle, in series with the muscle model (Fig. 1A). The whole system (muscle and oscillator) was held at constant length by fixing both ends. The system dynamics were driven by a cyclical activation signal a that caused the muscle to generate active force at the start of each activation cycle when the magnitude of activation increased. The muscle force Fm in turn accelerated an inertial point mass m, within the harmonic oscillator, that was considered to be external to the muscle, and this acceleration was resisted by a linear spring with stiffness k and a viscous damper with damping coefficient b that both acted only along the length (x-axis) of the system. When the muscle activation decreased at the end of each activation cycle, the harmonic oscillator force dominated the muscle force and caused the muscle to increase in length. Thus, the system allowed the muscle to cyclically contract with unconstrained muscle lengths and forces to mimic the contraction cycles seen during in vivo locomotion. Fig. 1 View largeDownload slide The different muscle models used for the simulated work-loops. (A) The spring, mass, and damper of the harmonic oscillator were in series connected to a muscle model (cylinder), with the ends of the whole system being held at fixed positions. The muscle model was evaluated using a series of 1D models: (B) a Hill-type model with a contractile element (CE) and parallel elastic element (PEE); (C) a Hill-type model with an additional history-dependent element acting in parallel; (D) a mass enhanced model where the muscle mass was distributed into 16 point masses (p1 – p16), separated by Hill-type actuator. A 3D model was also tested within the oscillating system (E) that was evaluated using a fiber-reinforced composite biomaterial and solved using the finite element method. Fig. 1 View largeDownload slide The different muscle models used for the simulated work-loops. (A) The spring, mass, and damper of the harmonic oscillator were in series connected to a muscle model (cylinder), with the ends of the whole system being held at fixed positions. The muscle model was evaluated using a series of 1D models: (B) a Hill-type model with a contractile element (CE) and parallel elastic element (PEE); (C) a Hill-type model with an additional history-dependent element acting in parallel; (D) a mass enhanced model where the muscle mass was distributed into 16 point masses (p1 – p16), separated by Hill-type actuator. A 3D model was also tested within the oscillating system (E) that was evaluated using a fiber-reinforced composite biomaterial and solved using the finite element method. The activation signal a that drove the active muscle forces was derived from a pulsed square wave excitation signal u with prescribed maximum amplitude umax, frequency f, and duty cycle D. Excitation u was converted to a via a first-order bilinear differential equation (Zajac 1989), the characteristics of which have been described in detail elsewhere (Ross et al. 2018). The activity of muscle in vivo can be increased by recruiting more muscle fibers, and also by increasing the firing rate and thus the activity of the individual fibers. However, in our simplified model, typical of the ubiquitous Hill-type muscle models, the muscle contains only one fiber that can be scaled to any size. Thus, the activity is merely a scalar value that is applied to the muscle force. The muscle simulated in this study was selected to contain fibers that were parallel to the line-of-action of the muscle, and has no aponeurosis or internal tendon. As such, the model ignores the structural effects that come from pennation and muscle fiber gearing. Muscle fibers are slender with an aspect ratio (length:diameter) that can approach 1000, however, whole muscles are much less slender with a lower aspect ratio. In order to create a shape for the model that was a compromise between a single fiber and a whole muscle, we selected an intermediate geometry with an aspect ratio of 100. This may represent the geometry of a bundle of muscle fibers, or fascicle, but is still more slender than typical whole muscles. The optimal muscle length l0 was set to 0.02 m and the muscle cross-sectional area A was set to 0.0314 mm2. The maximum isometric stress σ0 was held constant at 225 kPa. A typical Hill-type muscle model without internal mass can be scaled to any size without any consequences on the kinematic and dynamic behavior of the model (Ross et al. 2018). However, the inclusion of internal muscle mass results in scale-dependent distortions where the load due to the muscle mass mm is proportional to the muscle volume and the muscle force Fm is proportional to A (Ross and Wakeling 2016). We tested these mass effects by comparing muscles of different sizes with identical geometric proportions and aspect ratios. For example, a muscle with a length scale factor λl0 of 10 would have 10 times greater l0, 100 times greater A, and 1000 times greater volume than a muscle with λl0 of 1. To control for the dynamics of the harmonic oscillator and achieve kinematic and dynamic similarity across different values of λl0, the parameters k, b, and m for the damped harmonic oscillator were scaled according to Ross et al. (2018). The average mass-specific power per cycle P* was calculated as the product of the net mass-specific work per cycle W* and f. The harmonic oscillator properties were chosen to allow the muscle to maximize P*. As such, the values for k, b, and m, as well as D for the excitation signal, were determined using nonlinear global optimization with the objective of maximizing P* for the control condition with umax of 1.0, f of 1 Hz, and fast fiber-type properties. These parameters were then kept constant, except for D, which was optimized for each f and fiber-type. Description of models Root model The root model has been previously described elsewhere (Ross et al. 2018), and is shown in Fig. 1B. In brief, the muscle in the root model consisted of a contractile element (CE) and a parallel elastic element (PEE). The CE developed forces when activated, and the force was proportional to the activation state, the force–length and force–velocity properties within the muscle (Fig. 2), as well as the instantaneous length and velocity of the muscle. Although the force–velocity properties of the CE result in implicit damping effects, we did not incorporate any additional damping element within the muscle. The force–length and force–velocity properties were defined by their Bezier curve descriptions (Ross et al. 2018). The root model did not include internal mass, 3D, or history-dependent effects. Fig. 2 View largeDownload slide Force–velocity (A) and force–length (B) properties for the Hill-type models of the muscle fibers. The curves were parameterized using Bezier curves (Ross et al. 2018). Fig. 2 View largeDownload slide Force–velocity (A) and force–length (B) properties for the Hill-type models of the muscle fibers. The curves were parameterized using Bezier curves (Ross et al. 2018). Mass model We tested this by subdividing the muscle into a series of smaller CEs and PEEs, divided by point masses that accounted for the muscle mass, following the scheme presented by Ross and Wakeling (2016). The total muscle mass was divided into 16 CEs that were separated by 16 point masses (Fig. 1D). The intrinsic properties of the CEs and PEEs remained the same as for the root model. The muscle mass mm was calculated as the product of the muscle volume and the density ρ, and as ρ was assumed to be constant at 1060 kg m−3, mm scaled in proportion to the muscle volume. As mm was distributed evenly along the length of the muscle at rest, the mass of each point mass was equal to mm divided by the number of point masses. History-dependent model To test for the effect of history-dependent phenomena on the resulting cyclic work-loops in this study, we added a history-dependent element in parallel to the CE of the root model (Fig. 1C). The force of this history-dependent element was determined by the sum of the stretch-induced force enhancement FFE and the shortening-induced force depression FFD. The method to include FFE was styled from a recent theory (Rode et al. 2009; Herzog et al. 2012; Herzog 2014). FFE was calculated using an empirical description (in a similar manner to the force–length and force–velocity parameters of the Hill-type model) and represents the additional titin forces that are believed to develop after active stretch. When the activation of the model increased beyond a threshold (50% maximum activation for each trial), the instantaneous muscle length at the time the threshold was reached was retained as a state constant le. If the instantaneous muscle length lm was stretched beyond le, then FFE was calculated from the normalized elastic modulus and the normalized muscle stretch: FFE=F0lm-leleE^lel0, (1) where F0 is the maximum isometric force for the muscle, and Ê is the normalized elastic modulus that varies as a function of le normalized to l0 in a non-linear fashion (Fig. 3A) that was determined from experimental data for the ascending (Hisey et al. 2009) and descending (Edman et al. 1982) limbs of the force–length relationship. This empirical description of force enhancement after stretch results in predicted enhanced forces that match experimental values for steady-state contractions (Fig. 3B). When the activation state returned below the threshold, the FFE was reduced to half its value to estimate the residual force enhancement, consistent with experimental data (Hisey et al. 2009). The force enhancement FFE reduced to zero, and the muscle was returned to its non-enhanced state, when lm shortened to less than le. This process would then be repeated for the subsequent contraction cycles. For simulations where the excitation was increased to a sub-maximal level with umax of 0.1 or 0.5, the maximum activation in each burst was <1. We assumed that only a portion of the fibers in the muscle are recruited during these sub-maximal simulations, and therefore, only the active portion of fibers would have enhanced titin forces and contribute to FFE. Thus, FFE was scaled by the maximum activation achieved in each contraction cycle. Fig. 3 View largeDownload slide Method for calculating force enhancement that occurs during periods of active stretch. (A) The modulus of the enhanced force (relative to the length that the muscle was activated) was determined for experimental data from both the ascending (Hisey et al. 2009) and descending (Edman et al. 1982) limbs of the force–length relation. (B) The predicted enhanced forces are shown in the context of steady forces after stretch in comparison to experimental values from Edman et al. (1982) and Hisey et al. (2009) (circles). Fig. 3 View largeDownload slide Method for calculating force enhancement that occurs during periods of active stretch. (A) The modulus of the enhanced force (relative to the length that the muscle was activated) was determined for experimental data from both the ascending (Hisey et al. 2009) and descending (Edman et al. 1982) limbs of the force–length relation. (B) The predicted enhanced forces are shown in the context of steady forces after stretch in comparison to experimental values from Edman et al. (1982) and Hisey et al. (2009) (circles). The magnitude of force depression has been described as a linear function of the net work done during each active contraction (Herzog et al. 2000). The force depression FFD was initially modeled following McGowan et al. (2010): FFD=-CFDWc, (2) where Wc is the work done during the active contraction and CFD is the coefficient for force depression that was set to a value of 9.43 N J−1 (McGowan et al. 2010 from Herzog et al. 2000). However, the coefficient CFD had dimensions and needed to be non-dimensionalized so that it could be applied across scales. Assuming that the magnitude of force depression per sarcomere is the same across all muscle sizes, the non-dimensional force depression coefficient was calculated as: F^FD=0.3387F^mW^c, (3) where the non-dimensional work W^c and non-dimensional force F^m were given by: W^c=WcWmax (4) and F^m=FmF0. (5) Therefore, C^FD=0.3387FmWmaxF0Wc, (6) where Wmax is the product of l0 and F0, calculated using data from Herzog et al. (2000). Thus, FFD was calculated as: FFD=-F0C^FDW^c. (7) If the muscle started to lengthen while still active, Wc would decrease and so too would FFD, however, FFD was set to zero if Wc became negative. In a similar manner to the force enhancement, the contribution of the force depression was assumed to only occur in active fibers. Thus, the FFD was also scaled by the maximum activation achieved in each contraction cycle. 3D finite element model To study the effects of the 3D structure of muscle on its behavior during cyclic contractions, we extended the isometric 3D finite element model (FEM) of muscle from Rahemi et al. (2014) to allow for dynamic behavior (Fig. 1E). This FEM is described in the following section, with additional details in Appendix 1. In this scheme, the parallel muscle fibers were represented by contractile tissue surrounded by base material. The intrinsic force–velocity and force–length properties of the fibers, as well as the excitation u and activation a profiles, were the same as those of the Hill-based root model (Fig. 2). The 3D FEM muscle had the same optimal length l0 as the root model, but a smaller aspect ratio of 10 with a square A. Following the works of Weiss et al. (1996) and Blemker and Delp (2005), we modeled the muscle as a fiber-reinforced, nonlinearly-elastic transversely-isotropic material. The constitutive laws for the materials have been described in detail elsewhere (Rahemi et al. 2014), and include the effects of the muscle being nearly incompressible to prevent large changes in volume, as well as the isotropic nonlinearly-elastic behavior associated with tissue base material properties in all directions. As in Rahemi et al. (2014), a Yeoh-type model (Yeoh 1993) was used to define the base material properties, which was scaled by a positive factor sbase to control the stiffness of the base material. The force–velocity and force–length properties of the muscle only operated along the direction of the fibers, and were incorporated into a description of the strain energy. The effect of a uniform muscle density ρ throughout the material was incorporated into a kinematic energy contribution. The FEM system was defined by: dharmonicoscillatorE+muscleEdt=-bdxdt2, (8) where E denotes energy, x denotes the displacement of the point mass m of the harmonic oscillator, and b denotes the damping coefficient in the harmonic oscillator, both as for the root model. As in the root model, we have not explicitly incorporated damping effects within the muscle, although they do occur implicitly as a consequence of the force–velocity properties of the CE. The energy associated with the point mass in the harmonic oscillator m is equal to the sum of its kinetic and potential energies, while the energy associated with the muscle is equal to the sum of the muscle kinetic energy, the volumetric strain energy, the base material strain energy, and the contractile strain energy. Note that only the contractile strain energy, and not the other three terms that make up the muscle energy, is accounted for in the standard 1D Hill-type model. The volumetric strain energy captures the material’s resistance to volume changes (nearly incompressible), while the contractile strain energy describes the energetic contributions due to the passive (PEE) and active (CE) elements of the fiber. The system dynamics are described by the changes in how the total energy is distributed among these different forms of energy over time. Thus, if the total energy is fixed, the relative contributions of the harmonic oscillator energy will be higher in the 1D root model than they are in the 3D model. The overall action integral of the muscle system describes how the total energy evolves over time as a function of the local displacement, dilation, and pressure. The principle of stationarity of the action integral was used to derive the equations in these unknowns, following the work of Simo et al. (1985) and Simo and Taylor (1991). This method allows the muscle fiber stretch to vary along its length, depending on the local forces it experiences, and moreover, the fibers are permitted to curve, unlike in the 1D Hill-based models. The dynamics of the muscle were coupled through a balance of forces to the external damped harmonic oscillator system described earlier, and the fibers in the 3D model were activated using the same u as was used for the root, mass, and history-dependent models. This resulted in a mathematical system that was dynamic and highly nonlinear. The FEM muscle was subdivided into 128 smaller elements (parallelepipeds), and each element had local activation, base material properties, Hill-type contractile properties of the fibers, and specified initial fiber orientation relative to the geometry. Although these properties can be varied from voxel to voxel, for the purposes of this study, we assumed they were uniform throughout the tissue. We used semi-implicit time stepping to discretize in time, and the Newton method to solve the nonlinear equations describing the displacement, dilation, and pressure at each time step. Simulations We tested the root, history-dependent, and mass models under a range of maximum excitations umax, frequencies f, length scale factors λl0, and fiber-type properties to replicate the contraction regimes seen during in vivo locomotion. The muscle behaved as if it was composed of either entirely fast or entirely slow fibers, and this was implemented into the models by setting the maximum shortening strain-rate v0 to either 10 s−1 or 5 s−1 and the activation rate constant τact to either 0.025 s or 0.045 s, respectively. These values are selected to represent typical properties for fast and slow fibers (Dick et al. 2017). umax was set to either 0.1, 0.5, or 1.0, and f was set to either 1, 2, 4, or 8 Hz. The duty cycle D was re-optimized for the different values of f, and decreased with increased f. For the root model with v0 of 10 s−1, D ranged from 0.52 at f of 0.5 Hz to 0.31 at f of 8 Hz; for the slower-fibered model with v0 of 5 s−1, D ranged from 0.55 at f of 0.5 Hz to 0.29 at f of 8 Hz. The mass of the harmonic oscillator m was optimized for the control condition and scaled with λl0 across the range of muscle sizes tested (Ross et al. 2018). We also tested a range of muscle masses from 0.67 mg at λl0 of 1, which is typical of a muscle fascicle, to 670 g at λl0 of 100, which is less than the plantar flexor muscle mass in humans (Handsfield et al. 2014), for the mass model. When muscle is sub-maximally activated, not all the muscle fibers need to be recruited, and the intrinsic properties (fiber-type) of the activated fibers affect the mechanical output of the whole muscle. If the output of the whole muscle depended only on the combined output of the active fibers, for instance following Hill (1970), then recruiting slower fibers would result in lower power output than recruiting faster fibers. However, when the internal mass of the muscle is considered, the presence of inactive fibers provides inertial resistance, even when these inactive fibers are not developing active force (Holt et al. 2014; Ross and Wakeling 2016), and this may occlude differences due purely to the fiber-type. Thus, the root model would be expected to develop different power outputs depending on its fiber-type, with these differences seen across activation levels, however, the fiber-type contribution would be smaller for the model with internal mass, particularly at larger muscle sizes. This was tested by evaluating the muscle model with slower intrinsic properties with v0 of 5 s−1 and τact of 0.045 s, and at the lower levels of excitation with umax of 0.1 and 0.5. 3D simulations were achieved at umax of 0.1 and 1.0, and for increased base material stiffness (133% sbase), decreased density (50% ρ), and slower fibers than for a control simulation. The control simulation used fast intrinsic properties and sbase of 250. The numerical solution of the 3D model involves choosing several algorithmic parameters for the linear and nonlinear iterations. These in turn affect the maximal level to which the muscle can be activated in a numerically stable fashion, and these parameter choices depend on the stiffness of the base material. To simulate the 3D model at umax of 1.0, it was necessary to use a base material stiffness sbase that was higher than used by Rahemi et al. (2014). However, in that work the maximal activation achieved was considerably lower. Further mathematical investigation into the relations between the algorithmic parameters, base material stiffness, and activation levels is needed in addition to more extensive experimental measurements of this stiffness. Nonetheless, simulations for the 3D FEM were repeated for excitation frequencies of 1, 2, and 4 Hz. Results The muscle excitation u for each simulation was represented by a repeating square wave that took a value of either 0 or umax. The resulting activation a increased rapidly at the beginning of each cycle of u, plateaued across mid-burst, and finally decayed at a slower rate than the rate of activation when u returned to zero. The muscle generated active forces when a was greater than 0, and this caused the muscle to shorten and accelerate the point mass m within the harmonic oscillator (Fig. 1A). When the muscle was inactive, the restoring force due to the spring in the harmonic oscillator exceeded the muscle force, and caused m to accelerate to lengthen the muscle. The muscle length excursions were thus in-phase with the timing of u. However, there was a short time lag between when the muscle was excited and when it started to shorten due to both the time required for the muscle force to accelerate and change the direction of the velocity of m, and the time required to reach peak a after u changes to umax. Similarly, the muscle started to lengthen before it was fully deactivated due to the restoring force exceeding the muscle force toward the end of the activation period. As a consequence of the non-zero passive muscle and harmonic oscillator forces when the muscle was at l0, the values of the initial position and velocity of m required to solve the system were not immediately evident. However, because only the transient solution of a damped system, and not the steady-state solution, depends on the choice of initial conditions, we were able to avoid finding the ideal initial conditions by examining only the steady-state solution. Thus, herein we report only steady-state results where the behavior of the system repeats each cycle. The greatest average mass-specific muscle power output P* across all models and simulations was 42 W kg−1, and was achieved for the control simulation in the root model where f was 1 Hz, umax was 1, and the muscle had fast fiber-type properties. Rather than being of physiological importance, this was due to selecting the harmonic oscillator parameters and D to maximize P* for only these control conditions in the root model. For these conditions, the maximum shortening velocity reached was only 11% of v0, and the maximum stress was 66% of σ0. However, the theoretical maximum P* for the control simulation in the root model would have been 92 W kg−1, given these force–velocity and force–length properties, σ0, and v0 (following Weis-Fogh and Alexander 1977), so the control condition does not reflect the optimal conditions for maximum steady-state power production. The main effects of f, umax, and fiber-type on the Hill-based root, mass, and history-dependent models are shown in Fig. 4. The net mass-specific work per cycle W* was lower for higher cycle frequencies, and at the highest f of 8 Hz there was insufficient time for the muscle to fully activate or deactivate. This resulted in muscle force oscillations that were too fast to drive substantial displacements of the harmonic oscillator, and thus, there was little length excursion in the muscle for the high frequencies (Fig. 5A). The value of P* was greatest at the intermediate f of 1 Hz where the harmonic oscillator was calibrated to allow the muscle to generate maximal power. As umax decreased from 1.0 to 0.5 or 0.1, the activation and muscle force decreased. This led to the muscle operating at longer lengths with smaller excursions, as it was not able to generate large enough forces to contract to shorter lengths. Consequently, both W* and P* were smaller in magnitude for lower levels of umax (Figs. 4 and 5B). The fiber-type properties of the muscle were determined by v0 and τact, where the fast muscles had v0 of 10 s−1 and τact of 0.025 s, and the slow muscles had v0 of 5 s−1 and τact of 0.045 s. Muscles composed of the slower fiber-type developed lower active force during shortening leading to smaller length excursions and lower W* and P* (Figs. 4 and 5C). Fig. 4 View largeDownload slide Main effects of frequency, fiber-type, and excitation level on the mass-specific work W* and power P* from the 1D model during cycling contractions. The main effects were calculated from a set of 166 simulations across this parameter range, and included data from the mass-enhanced, history-dependent, and root Hill-type models. Fig. 4 View largeDownload slide Main effects of frequency, fiber-type, and excitation level on the mass-specific work W* and power P* from the 1D model during cycling contractions. The main effects were calculated from a set of 166 simulations across this parameter range, and included data from the mass-enhanced, history-dependent, and root Hill-type models. Fig. 5 View largeDownload slide Work-loops for muscle contractions demonstrating all tested parameters. (A–E) Work-loops are compared with the control simulation at 1 Hz frequency, maximal excitation, fast-fiber type, and small muscle size. Comparisons are shown that vary: (A) the frequency from 1 to 8 Hz; (B) the maximal excitation from 1.0 to 0.1; (C) the fiber type from fast to slow; (D) the muscle size from 0.67 mg to 670 g using the mass-enhanced model; and (E) history-dependent effects. (F) Results from the finite element simulations of the 3D model. Again, curves are compared with their control simulation at 1 Hz frequency, maximal excitation, and fast-fiber type with comparisons made to a model with greater base material stiffness, reduced muscle density, or slower fiber-type. The temporal direction of all work-loops is counterclockwise. Fig. 5 View largeDownload slide Work-loops for muscle contractions demonstrating all tested parameters. (A–E) Work-loops are compared with the control simulation at 1 Hz frequency, maximal excitation, fast-fiber type, and small muscle size. Comparisons are shown that vary: (A) the frequency from 1 to 8 Hz; (B) the maximal excitation from 1.0 to 0.1; (C) the fiber type from fast to slow; (D) the muscle size from 0.67 mg to 670 g using the mass-enhanced model; and (E) history-dependent effects. (F) Results from the finite element simulations of the 3D model. Again, curves are compared with their control simulation at 1 Hz frequency, maximal excitation, and fast-fiber type with comparisons made to a model with greater base material stiffness, reduced muscle density, or slower fiber-type. The temporal direction of all work-loops is counterclockwise. The non-dimensional dynamics of the muscle were identical across muscle length scales λl0 for the root Hill-type model that had no mass (i.e., stress, strain, W*, and P*). By contrast, inclusion of internal mass within the muscle model resulted in reduced net active muscle force during contraction, as well as increased passive force during lengthening due to the inertial costs of this mass. These changes in force were subtle, but were accentuated for simulations with greater λl0, with the larger mass-enhanced models generating lower W* and P*. Larger muscles and lower umax resulted in smaller differences in P* between fast and slow fiber-type conditions in the mass model compared with the massless root model. There was minimal active stretch in the muscle after the start of each cycle of u, so the contribution of force enhancement was minimal across all conditions. However, the muscle was largely active during shortening, resulting in significant force depression effects. The lower forces in the history-dependent model due to force depression resulted in the muscle not contracting to such short lengths, so both W* and P* were lower for the history-dependent model compared with the root model (Fig. 5E). The main effects of the different contractile and muscle parameters were mirrored in the 3D model. The 3D FEM was evaluated at umax of 0.1 and 1.0, and at f of 1, 2, and 4 Hz. There were large increases in W* and P* when umax was increased from 0.1 to 1.0, and decreases in W* and P* when the frequency increased from 1 Hz to 4 Hz. The effects of the base material stiffness, muscle density, and fiber-type were similar for both the low and high activation levels, with increased base material stiffness resulting in reduced W* and P*, decreased muscle density ρ (and thus mass) resulting in greater W* and P*, and the slower-fibered model resulting in lower W* and P* (Figs. 5F and 7). Overall, the greatest effects of the 1D mass-enhanced and history-dependent models can be observed for the contractile conditions with f of 1 Hz and umax of 1.0 (Fig. 5D, E). The mass-enhanced model at the largest scale tested (λl0 of 100; equivalent to a 670 g muscle) generated 12% less work and power with fast fiber-type properties and 4% less with slow properties compared with the root model with the same conditions. The history-dependent model was dominated by the force depression effect for these contractions and resulted in a 7% reduction in work and power compared with the root model. For the 3D model, there were additional substantial effects of the base material stiffness on W* and P*, and the mass effect was pronounced and important even at its 67 mg muscle mass. Discussion Muscle mass and scaling of contractile properties The root model was styled from a typical Hill-type muscle model, and as such, there were no mechanisms within the model to distort the intrinsic properties of the muscle as it scaled to different sizes. This implicitly means that simulations of muscles at any scale would result in the same non-dimensional muscle performance provided the external load from the harmonic oscillator was scaled appropriately. Indeed, this is how Hill-type muscle models are typically used: they are defined using intrinsic properties measured during experiments on single fibers or small muscles, which are in turn used to predict the function of whole muscles that may be of larger size (Wakeling and Lee 2011). The results from the root model in this study confirm that the function predicted by typical Hill-type models is conserved across scales. However, these models ignore the inertial effects of internal muscle mass on the dynamics of contraction. Previous studies have shown that internal mass acts to slow the rate of force development (Günther et al. 2012), and the maximum contraction velocity of larger muscles (Ross and Wakeling 2016). This study confirms these findings and also demonstrates that the internal mass results in a reduction in the work and power produced by whole muscles during cyclic contractions (Figs. 5 and 6). Fig. 6 View largeDownload slide Relative changes in performance between the 1D root model with the history-dependent and mass-enhanced 1D models. (A) Muscle mass-specific work W* and (B) muscle mass-specific power P* are shown relative to W* and P* for the control simulation (at 1 Hz frequency and maximal excitation). Simulations with different fiber-types were compared with their fiber-type control. Fig. 6 View largeDownload slide Relative changes in performance between the 1D root model with the history-dependent and mass-enhanced 1D models. (A) Muscle mass-specific work W* and (B) muscle mass-specific power P* are shown relative to W* and P* for the control simulation (at 1 Hz frequency and maximal excitation). Simulations with different fiber-types were compared with their fiber-type control. The mass of the muscle tissue alters the properties and function of muscle depending on the scale. Muscle mass varies in proportion to its volume, whereas muscle force varies with its cross-sectional area. This means that increasing muscle size, and decreasing the area to volume ratio, would result in the internal mass increasing faster than the available forces to accelerate that mass. Thus, the mass effects will predominate at larger muscle sizes. Indeed, the 1D models predicted a negligible effect of the internal mass for a scale of 1 where the muscle was 0.67 mg, whereas the mass effect was large and important for scale of 100 where the muscle was 670 g. However, the 3D model resulted in the internal mass being accelerated in all three dimensions, and here the mass effects were apparent even at its muscle size of 67 mg. The effect of the internal mass of the muscle tissue was important in this study due to the dynamic nature of the contractions. It should be noted that models of isometric contractions where the accelerations are negligible would show little difference if mass were included. However, internal work is required to accelerate the mass of the muscle tissue when the muscle changes length. Additionally, the 3D finite element representation of the muscle allows transverse deformations of the tissue and therefore transverse accelerations of the muscle mass, so the effect of the internal mass becomes even more important for the fully 3D model. Differences in muscle fiber-types result in different forces for given muscle lengths and velocities, and so we could expect the mass effect to depend on the intrinsic speed of the contracting fibers. Indeed, model results show that the maximum (normalized) shortening velocity of a muscle is reduced for larger sizes, and this effect is exacerbated for muscles with faster intrinsic speeds of their CEs (Ross and Wakeling 2016). Here we additionally show greater reductions in the mass-specific work and power output during cyclic contractions for muscles with faster intrinsic speeds (Fig. 6), which are caused by the dynamics of the muscle mass. Additional mass effects are apparent in muscles that are not fully active. Inactive muscle fibers have been proposed to contribute drag to a contracting whole muscle, resulting in a mechanical cost that would diminish the mechanical output of the muscle (Josephson and Edman 1988; Holt et al. 2014). The effects of the internal muscle mass would predominate at lower activations in a similar manner to which they predominate at larger muscle sizes, since at low activations there would be reduced contractile force available to accelerate the internal mass of the muscle tissue. Indeed, the results of this study support this notion. Muscle function and history-dependent properties When a muscle generates high power, it must contract with high force while shortening, and generate negligible force when it is lengthening. This would require the muscle to be activated close to its maximum length, and to be deactivated close to its minimum length, and this pattern emerged in our contraction cycles. However, activating close to maximum length provides minimal history-dependent force enhancement during active stretch. Indeed, the magnitude of force enhancement was negligible for the stretch-shortening contraction cycles in this study (Fig. 5E). By contrast, a substantial effect of force depression was observed. Thus, the majority of the history-dependent effects from these power-producing cyclic contractions occurred due to the force depression during shortening. In this study, the history-dependent effects resulted in a reduction in the net work and average power output per cycle of the muscle of about 7% for the 1 Hz frequency and maximal excitation condition. Although few studies have examined history-dependent effects during cyclic contractions, similar magnitudes of force depression have been previously reported. McGowan et al. (2010) found a 6% reduction in isometric torque during simulated knee extension and a 28% reduction in average crank power during simulated cycling when force depression was incorporated into a lower-body musculoskeletal model. This magnitude of joint torque reduction during knee extension is comparable to a 6.5% reduction in torque reported experimentally with the same conditions (Lee et al. 1999). Despite these findings, these musculoskeletal simulations involved more than one joint and multiple muscle actuators, and the Hill-type formulation of these actuators accounted for the effects of series elasticity. Thus, the contribution of force depression to the behavior of individual muscles in these simulations is difficult to resolve, so the results of our simulations are not directly comparable to the results of these previous studies. However, as our control simulation was designed to maximize muscle work and power output per cycle, the magnitude of force depression reported here is likely as large as would be achieved by a single muscle for contraction cycles during locomotion. Intrinsic speeds and fiber-type effects The maximum power output from a muscle depends on the intrinsic speed (maximum strain-rate) of its muscle fibers (Weis-Fogh and Alexander 1977). Not only can faster muscle fibers contract at faster speeds than slow fibers, but they generate greater stress at any given speed. This study used forward dynamic simulations of contractile cycles, whereby the range of motion, shortening velocity, and muscle force could all vary as a function of the excitation parameters used to drive the simulations. The simulations showed that the muscle with faster intrinsic speed (v0 of 10 versus 5 s−1) generated greater forces, allowing for greater length changes and thus higher shortening velocities than the slower muscle. The resultant work-loops showed that the faster muscle generated greater work per cycle and power output for an equivalent excitation and frequency (Fig. 5C). Mammalian muscle is commonly comprised of a heterogeneous population of fiber-types. If a heterogeneous muscle is sub-maximally activated then the intrinsic properties of the fibers that are activated may govern the contractile properties of the whole muscle. However, experimental studies have shown that there may be less difference in the shortening velocity of whole muscle than would be predicted by fiber-type alone (Holt et al. 2014). This effect has been attributed to the “drag” of the inactive tissue (Josephson and Edman 1988; Holt et al. 2014), and may be caused by the inertial properties of the internal muscle mass (Ross and Wakeling 2016), as discussed above. The pattern of recruitment between different types of motor units is therefore important in determining the mechanical output of whole muscle. The usual recruitment pattern follows the “size principle” whereby slower motor units are recruited first and for lower activations (Henneman et al. 1965a, 1965b). However, alternative recruitment patterns have been described that may favor recruitment of greater proportions of faster motor units at sub-maximal activations (see review: Hodson-Tole and Wakeling 2009) at fast contraction speeds and high contraction frequencies (Blake and Wakeling 2014). The simulations from this study show that, even at low excitation (10%), the type of fibers recruited alters the work and power output of the whole muscle. Additionally, our findings confirm that inertial effects from the muscle mass cause slower muscle shortening speeds and result in lower W* and P* than in smaller muscle or muscle without mass (Fig. 6). 3D structural, volumetric, and base material effects The contractile forces generated in the muscle are developed within the myofilaments in each muscle fiber from the cross-bridge forces between the actins and myosins, and from stretch in the PEEs such as the titins. The myofilaments are aligned with the longitudinal axis of the muscle fibers, and thus give rise to the magnitude and direction of the fiber force in this study. These forces are represented by the CE and the PEE in these models (Fig. 1). The fiber forces act to shorten the muscle fibers in their longitudinal direction. Muscle fibers are typically assumed to be isovolumetric (Baskin and Paolini 1967) and so must increase in girth to maintain their volume as they shorten. In the 3D model in this study, the volumetric constraints redistribute the longitudinal fiber forces to act in all three dimensions. This isovolumetric property is governed by incompressibility of the fluids that make up the intra and extracellular muscle tissue (muscle contains about 80% water; Van Loocke et al. 2008). Typical Hill-type muscle models (Zajac 1989) are 1D and contain no volume; however, in our 3D model the volumetric constraints act to maintain a nearly constant volume (changing by no >0.03% in these simulations). It should be noted that small changes in muscle volume have been seen at the fascicle level attributed to compressibility of the extracellular matrix (ECM) (Smith et al. 2011), and at the whole muscle level as the intramuscular pressure restricts the blood flow and blood volume (Ashton 1975), and so relaxation of the volumetric constraints may be warranted for future investigations. The 3D model represents the muscle as a fiber-reinforced composite material, where the fiber properties are governed by the myofilament contractile forces, and the composite properties represent the base material of the tissue. This base material contains properties across all muscle scales, including intracellular stiffness and properties from the ECM. Similar approaches have been employed in other 3D models of muscle where a variety of strain energy functions have been used to characterize the base material properties (e.g., Johansson et al. 2000; Blemker and Delp 2005). The volumetric constraints in the 3D model convert longitudinal shortening into a tendency for transverse expansion of the muscle. However, this expansion is resisted by the base material stiffness and so internal work must be done to achieve transverse expansion to allow for longitudinal shortening. Results from this study show that as the base material stiffness is increased, a greater magnitude of internal work is required for tissue deformation, resulting in reduced muscle external work and power (Figs. 5F and 7). Fig. 7 View largeDownload slide Changes in muscle mass-specific power P* between different 3D FEM simulations. The simulations were evaluated for 1 Hz cycle frequencies and maximal excitation. Models varied by increasing the scaling of their base material stiffness sbase, decreasing their tissue density ρ and thus mass, and by changing their fiber-type properties which were fast-fibered in the control model compared with a slow-fibered simulation. Note that being at a 1 Hz cycle frequency, the percentage changes in muscle mass specific work W* were the same as for these changes in P* for these simulations. Fig. 7 View largeDownload slide Changes in muscle mass-specific power P* between different 3D FEM simulations. The simulations were evaluated for 1 Hz cycle frequencies and maximal excitation. Models varied by increasing the scaling of their base material stiffness sbase, decreasing their tissue density ρ and thus mass, and by changing their fiber-type properties which were fast-fibered in the control model compared with a slow-fibered simulation. Note that being at a 1 Hz cycle frequency, the percentage changes in muscle mass specific work W* were the same as for these changes in P* for these simulations. The base material properties were represented in the model independently from the volumetric constraints, and indeed these were both independent from the density of the muscle tissue. However, the resultant properties of the whole muscle contraction emerged from the nonlinear interactions of these effects. Thus, the model provides an ideal framework for testing the specific effects of base material stiffness (or indeed ECM stiffness), volume changes, and mass properties. Such models inevitably benefit from thorough validation, however it is experimentally very challenging to independently vary any of these individual effects. Nonetheless, the mechanisms relating muscle shortening to its transverse expansion are consistent with a range of recent experimental studies: (A) transversely loading rat muscle results in reduced external force during isometric contractions (Siebert et al. 2014b), smaller transverse expansion, and greater internal work being done for transverse tissue deformation (Siebert et al. 2014a); (B) increasing the volume and intramuscular pressure within bullfrog muscle fibers by bathing them in hypotonic solution results in greater passive forces in the longitudinal direction (Sleboda and Roberts 2017). In this case the increased intramuscular pressure would result in greater transverse strain in the base material that would be redistributed by the volumetric constraints of the intracellular fluid, to result in greater resistance to longitudinal stretch due to the base material stiffness; and (C) resisting transverse expansion of leopard frog muscle by external physical constraint (to augment the ECM stiffness) causes the contracting muscle to shorten less in the longitudinal direction, and thus do less external work (Azizi et al. 2017). This effect is presumably due to the additional internal work that must be done to achieve the transverse deformation against the augmented base material stiffness. The 3D model in this study allowed us to test the effect of tissue density, and again, this is a difficult parameter to specifically vary in experiments. For instance, muscle density does vary with the fat content of the tissue, however, the adiposity also changes volumetric constraints and base material stiffness. Directly altering the density and thus tissue mass of the 3D model in this study showed that increased internal mass, and therefore additional internal work required for accelerating this mass, results in reduced external work and power from the muscle (Figs. 5F and 7). These mass effects were the same as those achieved in the mass-enhanced 1D Hill-type model and support the conclusion that the relative reductions in muscle work and power due to internal mass are more important at larger sizes and lower activation levels. The models in this study simulated contractions where the fibers were parallel to the line-of-action (longitudinal direction) of the muscle (Fig. 1). This allowed our study to focus on the contractile mechanisms in the fibers themselves rather than adding the geometric complexities that come with pennate architectures. However, the mechanisms and processes identified in the 3D model in this study will exist in any muscle fiber, regardless of the architecture of the muscle that contains it. Adjacent fibers interact between their shared boundaries, so transverse forces and deformations will propagate across the tissue. In pennate muscle, transverse fiber expansion must be balanced by increasing pennation angle in order for the muscles to maintain their spatial packing, with the aponeuroses that envelope the muscle influencing overall tissue deformation and changes in pennation. Indeed, we have previously used this 3D finite element framework (Rahemi et al. 2014, 2015) to show that: (A) stiffer aponeuroses lead to smaller increases in pennation, consistent with the reduction in changes in rat muscle gearing as it stiffens with age (Holt et al. 2016); (B) the curvatures that develop in contracting fibers are greatest near the aponeuroses, consistent with imaging studies in man (Namburete and Wakeling 2012); and (C) asymmetries in stress distribution across contracting muscles also lead to asymmetries in the transverse deformations of the muscle fibers and muscle belly (Rahemi et al. 2014). These findings are consistent with other experimental studies on muscle belly deformations (Azizi et al. 2008; Siebert et al. 2014b; Wakeling and Randhawa 2014), and provide a framework for understanding how the fundamental contractile properties of the fibers are related to tissue deformations, gearing, and the dynamics of whole muscle contraction. Assumptions, limitations, and future work The purpose of this study was to determine the effect of contractile features that are not accounted for in most Hill-type muscle models on the force and power of muscle during cyclic contractions. As such, we chose to orient the line-of-action of the fibers in these simulations along the longitudinal x-axis of the system (Fig. 1A). This enabled us to directly interpret the effects of the contractile properties on muscle performance because the entirety of the force vector from the fibers acted to accelerate the point mass m of the harmonic oscillator. However, not all muscles have their fibers arranged in parallel, and many have their fibers arranged at an oblique, pennate angle from the line-of-action of the muscle (Wickiewicz et al. 1983). As muscle fibers shorten they must increase in girth to maintain their volume. Such increases in girth allow pennate fibers to shorten and increase their pennation angle so that they can maintain their packing within the muscle belly (Lieber and Ward 2011). When the fibers rotate, their shortening velocity becomes uncoupled from that of the muscle belly in a process termed belly gearing (Wakeling et al. 2011), which lowers the fiber velocity relative to the belly velocity and alters muscle force, work, and power. Thus, architectural features of muscle, including pennation, add a layer of complexity to our understanding of muscle function in vivo. This understanding is further challenged when we consider that the exact nature of how muscles bulge in different directions (Maganaris et al. 1998; Azizi et al. 2008; Randhawa et al. 2013), and how variations in such bulging alters gearing (Azizi et al. 2008) is not fully known. Hence, focusing on parallel-oriented fibers allows us to probe fiber-level questions without being confounded by the additional complexity due to muscle architecture. We also chose to exclude tendinous, or series elastic, elements from the model for similar reasons. Tendons stretch when loaded, and if a tendon is long or compliant, this stretch can uncouple the velocity of the muscle belly from that of the muscle-tendon unit in a process known as tendon gearing (Wakeling et al. 2011). To make more direct interpretations of the effect of muscle (fiber or belly) shortening on the motion of the point mass within the harmonic oscillator, we assumed that tendon stretch was negligible and that the muscle was connected to m via a rigid link. However, it should be acknowledged that where significant tendinous tissue occurs in vivo, the tendon can further modify the behavior of the muscle-tendon-unit by amplifying muscle power output or attenuating muscle power input (Roberts and Azizi 2011). As the primary goal of this study was to examine the relative effect of different contractile phenomena on the work and power production of a muscle during cyclic contractions, the muscle model was designed to generate its maximum power for one of the intermediate contractile regimes with frequency f of 1 Hz, maximum excitation umax of 1, and fast fiber-type properties. Generating maximum power is an example of the muscle acting as a motor, and indeed our test condition operated at 50% of the maximum theoretical average mass-specific power per cycle (calculated for maximal activation and constant velocity) (Weis-Fogh and Alexander 1977). However, muscles can also act as springs, struts, or brakes during locomotion depending on the requirements of the task (Dickinson et al. 2000), and muscles may have specific adaptations to their properties depending on their primary function (Biewener 1998; Mörl et al. 2016). The muscle model in this study was coupled to a spring-mass-damper system that acted as a damped harmonic oscillator to provide an oscillating external load to allow for cyclic contractions of the muscle. Damped harmonic oscillators, by their nature, are only able to store and return or dissipate energy, and as a result, cannot be used to make the muscle behave like a brake. Therefore, the results of this study should be considered in the context of the power producing type of contractions that they were intended to simulate. As a consequence of only examining power-producing contraction cycles where the muscle is primarily active during shortening, force depression had a far greater effect than force enhancement in this study. However, these history-dependent effects are particularly challenging to represent mathematically due to a lack of consensus on the features of, and mechanisms behind, these effects, particularly during cyclic contractions with non-constant activations and length changes. Although force enhancement following active lengthening has been reported at sub-maximal levels of activation during voluntary (Oskouei and Herzog 2005, 2006a, 2006b; Pinniger and Cresswell 2007; Altenburg et al. 2008; Seiberl et al. 2012, 2013) and electrically stimulated (Abbott and Aubert 1952; Meijer 2002) contractions, the relationship between level of activation and magnitude of force enhancement is still poorly characterized. Furthermore, Meijer (2002) found that the magnitude of force enhancement was lower for lower levels of activation, and Oskouei and Herzog (2006b) only observed this phenomenon at sub-maximal activation in a portion of their study sample. Therefore, we chose to use 50% of maximum activation for the given simulation as a threshold value to indicate the start of force enhancement. However, as the excitation is prescribed as a square wave, the activation rises and decays rapidly so this threshold value would have little consequence on the magnitude of force enhancement for a whole cycle. We also chose to model force enhancement using a normalized elastic modulus that varied as a function of the instantaneous length at which this activation threshold was reached (Fig. 3A). To determine this relationship, we assumed that the enhanced muscle force as a function of muscle length (Fig. 3B) was linear and that the modulus depended only on the length where the activation threshold was reached and not on the current muscle length. However, further experimental studies examining force enhancement and depression during sub-maximal cyclic contractions are needed to determine if these assumptions are reasonable. In this study, we chose to model force depression as a function of muscle work, consistent with experimental evidence (De Ruiter et al. 1998; Herzog et al. 2000; Kosterina et al. 2008, 2009) and previous model formulations (McGowan et al. 2010, 2013). However, force depression has been shown to depend on contraction speed (Maréchal and Plaghki 1979; Meijer et al. 1997), and may decay during the period of time when the muscle is shortening (Till et al. 2014). More recently, Fortuna et al. (2018) found the magnitude of isometric, steady-state force depression following a stretch-shortening cycle to be independent of work. Thus, there is still debate as to whether or not force depression is accurately expressed as a function of only mechanical work. To include force depression across muscle model scales it was necessary to non-dimensionalize the force depression effect (Equation 7). Assuming force depression is work-dependent, the force would scale with muscle area, but work would scale with muscle volume in a similar manner to the mass scaling arguments discussed previously. Thus, a large muscle would generate a higher ratio of work to force, and at some limit, the contractile force would be depressed more than the force the muscle could generate. Hence, we assumed that the force depression per sarcomere would be independent of muscle size so that the force depression would result in an equivalent non-dimensional effect across muscle sizes. However, future experimental studies should be conducted to test the validity of this assumption. Hill-type muscle models typically do not use muscle-specific parameters but rather use parameter values taken from different studies on single fibers that vary widely depending on the muscle, species, temperature, preparation method (skinned versus intact), and stimulation protocol. The geometric dimensions and intrinsic properties of the muscle models in this study were chosen to reflect typical or average values reported in the literature so that our results could be more widely applicable than if we were to select values to represent a particular muscle. While we expect that our conclusions would hold if these parameter values were changed, this has not yet been evaluated and so should be considered when interpreting the results of this study. Using computational muscle models in this study allowed us to change features of the muscle that are difficult or impossible to manipulate experimentally. Although doing so provided a means to explore theoretical questions about muscle design and performance across muscles of different sizes under a range of conditions, it also limited our ability to directly validate our models against empirical data. Although we would expect that including more parameters in the muscle model formulation would provide more accurate predictions of in vivo muscle behavior, we have not explicitly tested this assumption. However, changing different model parameters, such as the fiber-type properties, led to changes in the muscle work and power output that were in a direction consistent with experimental evidence. Furthermore, the primary goal of this study was not to improve muscle model predictions, but rather to improve our understanding of the physiological mechanisms that may underlie muscle mechanical behavior during cyclic contractions. Conclusions This study examined a range of contractile properties that may contribute to the work and power output of muscle during cycling contractions. In our simulations, substantial history-dependent effects occurred for the power-producing contractions in these simulations, and were dominated by the force depression effect that occurs during active shortening. Our results show that the mechanical behavior of whole muscle during sub-maximal contractions depends on the fiber-type of the active fibers, even when the inertial effect of inactive muscle tissue is accounted for. There is an inertial cost to moving the tissue mass that results in reduced mass-specific work and power, and this is most pronounced for large muscles and low activations. 3D constraints also affect the mechanics of the muscle across all size scales. As the muscle changes length, the volumetric constraints cause it to undergo transverse deformations, and internal work against the base material stiffness must be done to deform the muscle tissue during contraction. Therefore, internal work must be done both to achieve tissue deformation and allow length changes, and to accelerate the muscle bulk. This study shows that the 3D effects, base material stiffness, and muscle mass all contribute to the dynamics of the muscle tissue, and also that these effects are size-dependent. These properties should thus be considered for a more complete understanding of muscle function during locomotor activities, where muscle activations are often sub-maximal and time-varying, muscles undergo unsteady cycles of length change, and whole muscles may be much larger than the single fibers from which their intrinsic properties have been determined. Acknowledgments The authors thank Dr. Emma Hodson-Tole for insightful discussions. Funding This work was supported by the National Institutes of Health [grant R01AR055648] and National Sciences and Engineering Research Council of Canada (NSERC) Discovery Grants [to both N.N. and J.M.W]. S.D. thanks the support of CONICYT-Chile through Becas Chile. References Abbott BC , Aubert XM. 1952 . The force exerted by active striated muscle during and after change of length . J Physiol 117 : 77 – 86 . Google Scholar CrossRef Search ADS PubMed Altenburg TM , de Ruiter CJ , Verdijk PW , Van Mechelen W , De Haan A. 2008 . Vastus lateralis surface and single motor unit EMG following submaximal shortening and lengthening contractions . Appl Physiol Nutr Metab 33 : 1086 – 95 . Google Scholar CrossRef Search ADS PubMed Altringham JD , Johnston IA. 1990 . Scaling effects on muscle function: power output of isolated fish muscle fibres performing oscillatory work . J Exp Biol 151 : 453 – 67 . Ashton H. 1975 . The effect of increased tissue pressure on blood flow . Clin Orthop Relat Res 113 : 15 – 26 . Google Scholar CrossRef Search ADS Azizi E , Brainerd EL , Roberts TJ. 2008 . Variable gearing in pennate muscles . Proc Natl Acad Sci U S A 105 : 1745 – 50 . Google Scholar CrossRef Search ADS PubMed Azizi E , Deslauriers AR , Holt NC , Eaton CE. 2017 . Resistance to radial expansion limits muscle strain and work . Biomech Model Mechanobiol 16 : 1633 – 43 . Google Scholar CrossRef Search ADS PubMed Baskin RJ , Paolini PJ. 1967 . Volume change and pressure development in muscle during contraction . Am J Physiol 213 : 1025 – 30 . Google Scholar PubMed Biewener AA. 1998 . Muscle function in vivo: a comparison of muscles used for elastic energy savings versus muscles used to generate mechanical power . Am Zool 38 : 703 – 17 . Google Scholar CrossRef Search ADS Biewener AA , Corning WR , Tobalske BW. 1998a . In vivo pectoralis muscle force–length behavior during level flight in pigeons (Columba livia) . J Exp Biol 201 : 3293 – 307 . Biewener AA , Konieczynski DD , Baudinette RV. 1998b . In vivo muscle force–length behavior during steady-speed hopping in tammar wallabies . J Exp Biol 201 : 1681 – 94 . Blake OM , Wakeling JM. 2014 . Early deactivation of slower muscle fibres at high movement frequencies . J Exp Biol 217 : 3528 – 34 . Google Scholar CrossRef Search ADS PubMed Blemker SS , Delp SL. 2005 . Three-dimensional representation of complex muscle architectures and geometries . Ann Biomed Eng 33 : 661 – 73 . Google Scholar CrossRef Search ADS PubMed Boettinger EG. 1957 . The machinery of insect flight. In: Scheer BT , editor. Recent advances in invertebrate physiology . Eugene : University of Oregon Publications . p. 117 – 42 . Cavagna GA , Citterio G. 1974 . Effect of stretching on the elastic characteristics and the contractile component of frog striated muscle . J Exp Biol 239 : 1 – 43 . Daley MA , Biewener AA. 2003 . Muscle force–length dynamics during level versus incline locomotion: a comparison of in vivo performance of two guinea fowl ankle extensors . J Exp Biol 206 : 2941 – 58 . Google Scholar CrossRef Search ADS PubMed Davydov D , Arndt D , Bangerth W , Heister T , Heltai L , Kronbichler M , Maier M , Pelteret JP , Turcksin B , Wells D. 2017 . The deal. II library, version 8.5 . J Numer Math 25 : 137 – 45 . De Ruiter CJ , De Haan AD , Jones DA , Sargeant AJ. 1998 . Shortening‐induced force depression in human adductor pollicis muscle . J Physiol 507 : 583 – 91 . Google Scholar CrossRef Search ADS PubMed Dick TJ , Biewener AA , Wakeling JM. 2017 . Comparison of human gastrocnemius forces predicted by Hill-type muscle models and estimated from ultrasound images . J Exp Biol 220 : 1643 – 53 . Google Scholar CrossRef Search ADS PubMed Dickinson MH , Farley CT , Full RJ , Koehl MAR , Kram R , Lehman S. 2000 . How animals move: an integrative view . Science 288 : 100 – 6 . Google Scholar CrossRef Search ADS PubMed Edman KAP , Elizinga G , Noble MIM. 1982 . Residual force enhancement after stretch of contracting frog single muscle fibers . J Gen Physiol 80 : 769 – 84 . Google Scholar CrossRef Search ADS PubMed Ettema GJ. 2002 . Effects of contraction history on control and stability in explosive actions . J Electromyogr Kinesiol 12 : 455 – 61 . Google Scholar CrossRef Search ADS PubMed Ettema GJ , Meijer K. 2000 . Muscle contraction history: modified Hill versus an exponential decay model . Biol Cybern 83 : 491 – 500 . Google Scholar CrossRef Search ADS PubMed Forcinito M , Epstein M , Herzog W. 1998 . Can a rheological muscle model predict force depression/enhancement? J Biomech 31 : 1093 – 9 . Google Scholar CrossRef Search ADS PubMed Fortuna R , Groeber M , Seiberl W , Power GA , Herzog W. 2017 . Shortening‐induced force depression is modulated in a time‐ and speed‐dependent manner following a stretch-shortening cycle . Physiol Rep published online (doi:10.14814/phy2.13279). Fortuna R , Kirchhuebel H , Seiberl W , Power GA , Herzog W. 2018 . Force depression following a stretch-shortening cycle is independent of stretch peak force and work performed during shortening . Sci Rep published online 8 (doi:10.1038/s41598-018-19657-8). Franklin CE , Johnston IA. 1997 . Muscle power output during escape responses in an Antarctic fish . J Exp Biol 200 : 703 – 12 . Google Scholar PubMed Günther M , Röhrle O , Haeufle DFB , Schmitt S. 2012 . Spreading out muscle mass within a Hill-type model: a computer simulation study. Comput Math Methods Med 2012 : 1 – 13 . Google Scholar CrossRef Search ADS Handsfield GG , Meyer CH , Hart JM , Abel MF , Blemker SS. 2014 . Relationships of 35 lower limb muscles to height and body mass quantified using MRI . J Biomech 47 : 631 – 8 . Google Scholar CrossRef Search ADS PubMed Henneman E , Somjen G , Carpenter DO. 1965 . Excitability and inhibitability of motoneurons of different sizes . J Neurophysiol 28 : 599 – 620 . Google Scholar CrossRef Search ADS PubMed Henneman E , Somjen G , Carpenter DO. 1965 . Functional significance of cell size in spinal motoneurons . J Neurophysiol 28 : 560 – 80 . Google Scholar CrossRef Search ADS PubMed Herzog W. 2014 . The role of titin in eccentric muscle contraction . J Exp Biol 217 : 2825 – 33 . Google Scholar CrossRef Search ADS PubMed Herzog W , Leonard TR. 2000 . The history dependence of force production in mammalian skeletal muscle following stretch-shortening and shortening-stretch cycles . J Biomech 33 : 531 – 42 . Google Scholar CrossRef Search ADS PubMed Herzog W , Leonard TR. 2002 . Force enhancement following stretching of skeletal muscle a new mechanism . J Exp Biol 205 : 1275 – 83 . Google Scholar PubMed Herzog W , Leonard T , Joumaa V , DuVall M , Panchangam A. 2012 . The three filament model of skeletal muscle stability and force production . Mol Cell Biomech 9 : 175 – 91 . Google Scholar PubMed Herzog W , Leonard TR , Wu JZ. 2000 . The relationship between force depression following shortening and mechanical work in skeletal muscle . J Biomech 33 : 659 – 68 . Google Scholar CrossRef Search ADS PubMed Hill AV. First and last experiments in muscle mechanics . UK : Cambridge University Press . p. 52 – 5 . Hisey B , Leonard TR , Herzog W. 2009 . Does residual force enhancement increase with increasing stretch magnitudes? J Biomech 42 : 1488 – 92 . Google Scholar CrossRef Search ADS PubMed Hodson-Tole EF , Wakeling JM. 2009 . Motor unit recruitment for dynamic tasks: current understanding and future directions . J Comp Physiol B 179 : 57 – 66 . Google Scholar CrossRef Search ADS PubMed Holt NC , Danos N , Roberts TJ , Azizi E. 2016 . Stuck in gear: age-related loss of variable gearing in skeletal muscle . J Exp Biol 219 : 998 – 1003 . Google Scholar CrossRef Search ADS PubMed Holt NC , Wakeling JM , Biewener AA. 2014 . The effect of fast and slow motor unit activation on whole-muscle mechanical performance: the size principle may not pose a mechanical paradox . Proc R Soc B 281 : 20140002. Google Scholar CrossRef Search ADS PubMed Huijing PA. 1999 . Muscle as a collagen fiber reinforced composite: a review of force transmission in muscle and whole limb . J Biomech 32 : 329 – 45 . Google Scholar CrossRef Search ADS PubMed Johansson T , Meier P , Blickhan R. 2000 . A finite-element model for the mechanical analysis of skeletal muscles . J Theor Biol 206 : 131 – 49 . Google Scholar CrossRef Search ADS PubMed Josephson RK. 1985 . Mechanical power output from striated muscle during cyclic contraction . J Exp Biol 114 : 493 – 512 . Josephson RK , Edman KAP. 1988 . The consequences of fibre heterogeneity on the force–velocity relation of skeletal muscle . Acta Physiol Scand 132 : 341 – 52 . Google Scholar CrossRef Search ADS PubMed Josephson RK , Stokes DR. 1989 . Strain, muscle length and work output in a crab muscle . J Exp Biol 145 : 45 – 61 . Kosterina N , Westerblad H , Eriksson A. 2009 . Mechanical work as predictor of force enhancement and force depression . J Biomech 42 : 1628 – 34 . Google Scholar CrossRef Search ADS PubMed Kosterina N , Westerblad H , Lännergren J , Eriksson A. 2008 . Muscular force production after concentric contraction . J Biomech 41 : 2422 – 9 . Google Scholar CrossRef Search ADS PubMed Lee SS , Arnold AS , de Boef Miara M , Biewener AA , Wakeling JM. 2013 . Accuracy of gastrocnemius muscles forces in walking and running goats predicted by one-element and two-element Hill-type models . J Biomech 46 : 2288 – 95 . Google Scholar CrossRef Search ADS PubMed Lee HD , Herzog W , Leonard T. 2001 . Effects of cyclic changes in muscle length on force production in in-situ cat soleus . J Biomech 34 : 979 – 87 . Google Scholar CrossRef Search ADS PubMed Lee HD , Suter E , Herzog W. 1999 . Force depression in human quadriceps femoris following voluntary shortening contractions . J Appl Physiol 87 : 1651 – 5 . Google Scholar CrossRef Search ADS PubMed Lieber RL , Ward SR. 2011 . Skeletal muscle design to meet functional demands . Philos Trans R Soc Lond B Biol Sci 366 : 1466 – 76 . Google Scholar CrossRef Search ADS PubMed Loeb GE , Hoffer JA , Pratt CA. 1985 . Activity of spindle afferents from cat anterior thigh muscles. I. Identification and patterns during normal locomotion . J Neurophysiol 54 : 549 – 64 . Google Scholar CrossRef Search ADS PubMed Machin KE , Pringle JW. 1960 . The physiology of insect fibrillar muscle. III. The effect of sinusoidal changes of length on a beetle flight muscle . Proc R Soc B 152 : 311 – 30 . Google Scholar CrossRef Search ADS Maganaris CN , Baltzopoulos V , Sargeant AJ. 1998 . In vivo measurements of the triceps surae complex architecture in man: implications for muscle function . J Physiol 512 : 603 – 14 . Google Scholar CrossRef Search ADS PubMed Maréchal G , Plaghki L. 1979 . The deficit of the isometric tetanic tension redeveloped after a release of frog muscle at a constant velocity . J Gen Physiol 73 : 453 – 67 . Google Scholar CrossRef Search ADS PubMed McGowan CP , Neptune RR , Herzog W. 2010 . A phenomenological model and validation of shortening-induced force depression during muscle contractions . J Biomech 43 : 449 – 54 . Google Scholar CrossRef Search ADS PubMed McGowan CP , Neptune RR , Herzog W. 2013 . A phenomenological muscle model to assess history dependent effects in human movement . J Biomech 46 : 151 – 7 . Google Scholar CrossRef Search ADS PubMed McGuigan MP , Yoo E , Lee DV , Biewener AA. 2009 . Dynamics of goat distal hind limb muscle-tendon function in response to locomotor grade . J Exp Biol 212 : 2092 – 104 . Google Scholar CrossRef Search ADS PubMed Meijer K. 2002 . History dependence of force production in submaximal stimulated rat medial gastrocnemius muscle . J Electromyogr Kinesiol 12 : 463 – 70 . Google Scholar CrossRef Search ADS PubMed Meijer K , Grootenboer HJ , Koopman BF , Huijing PA. 1997 . Fully isometric length–force curves of rat muscle differ from those during and after concentric contractions . J Appl Biomech 13 : 164 – 81 . Google Scholar CrossRef Search ADS Meijer K , Grootenboer HJ , Koopman HF , Van Der Linden BJ , Huijing PA. 1998 . A Hill type model of rat medial gastrocnemius muscle that accounts for shortening history effects . J Biomech 31 : 555 – 63 . Google Scholar CrossRef Search ADS PubMed Mörl F , Siebert T , Häufle D. 2016 . Contraction dynamics and function of the muscle–tendon complex depend on the muscle fibre–tendon length ratio: a simulation study . Biomech Model Mechanobiol 15 : 245 – 58 . Google Scholar CrossRef Search ADS PubMed Namburete AI , Wakeling JM. 2012 . Regional variations in fascicle curvatures within a muscle belly change during contraction . J Biomech 45 : 2835 – 40 . Google Scholar CrossRef Search ADS PubMed Nishikawa KC , Monroy JA , Uyeno TE , Yeo SH , Pai DK , Lindstedt SL. 2012 . Is titin a ‘winding filament’? A new twist on muscle contraction . Proc R Soc B 279 : 981 – 90 . Google Scholar CrossRef Search ADS PubMed Oskouei AE , Herzog W. 2005 . Observations on force enhancement in submaximal voluntary contractions of human adductor pollicis muscle . J Appl Physiol 98 : 2087 – 95 . Google Scholar CrossRef Search ADS PubMed Oskouei AE , Herzog W. 2006 . Force enhancement at different levels of voluntary contraction in human adductor pollicis . Eur J Appl Physiol 97 : 280 – 7 . Google Scholar CrossRef Search ADS PubMed Oskouei AE , Herzog W. 2006 . The dependence of force enhancement on activation in human adductor pollicis . Eur J Appl Physiol 98 : 22 – 9 . Google Scholar CrossRef Search ADS PubMed Perreault EJ , Heckman CJ , Sandercock TG. 2003 . Hill muscle model errors during movement are greatest within the physiologically relevant range of motor unit firing rates . J Biomech 36 : 211 – 8 . Google Scholar CrossRef Search ADS PubMed Pinniger GJ , Cresswell AG. 2007 . Residual force enhancement after lengthening is present during submaximal plantar flexion and dorsiflexion actions in humans . J Appl Physiol 102 : 18 – 25 . Google Scholar CrossRef Search ADS PubMed Rahemi H , Nigam N , Wakeling JM. 2014 . Regionalizing muscle activity causes changes to the magnitude and direction of the force from whole muscles . Front Physiol 5 : 298. Google Scholar CrossRef Search ADS PubMed Rahemi H , Nigam N , Wakeling JM. 2015 . The effect of intramuscular fat on skeletal muscle mechanics: implications for the elderly and obese . J R Soc Interface 12 : 20150365. Google Scholar CrossRef Search ADS PubMed Randhawa A , Jackman M , Wakeling JM. 2013 . Muscle gearing during isotonic and isokinetic movements in the ankle plantarflexors . Eur J Appl Physiol 113 : 437 – 47 . Google Scholar CrossRef Search ADS PubMed Richards CT , Clemente CJ. 2013 . Built for rowing: frog muscle is tuned to limb morphology to power swimming . J R Soc Interface 10 : 20130236. Google Scholar CrossRef Search ADS PubMed Roberts TJ , Azizi E. 2011 . Flexible mechanisms: the diverse roles of biological springs in vertebrate movement . J Exp Biol 214 : 353 – 61 . Google Scholar CrossRef Search ADS PubMed Roberts TJ , Marsh RL , Weyand PG , Taylor CR. 1997 . Muscular force in running turkeys: the economy of minimizing work . Science 275 : 1113 – 5 . Google Scholar CrossRef Search ADS PubMed Rode C , Siebert T , Blickhan R. 2009 . Titin-induced force enhancement and force depression: a ‘sticky-spring’ mechanism in muscle contractions? . J Theor Biol 259 : 350 – 60 . Google Scholar CrossRef Search ADS PubMed Ross SA , Nigam N , Wakeling JM. 2018 . A modelling approach for exploring muscle dynamics during cyclic contractions . PLoS Comput Biol 14 : e1006123. Google Scholar CrossRef Search ADS PubMed Ross SA , Wakeling JM. 2016 . Muscle shortening velocity depends on tissue inertia and level of activation during submaximal contractions . Biol Lett 12 : 20151041. Google Scholar CrossRef Search ADS PubMed Seiberl W , Hahn D , Herzog W , Schwirtz A. 2012 . Feedback controlled force enhancement and activation reduction of voluntarily activated quadriceps femoris during sub-maximal muscle action . J Electromyogr Kinesiol 22 : 117 – 23 . Google Scholar CrossRef Search ADS PubMed Seiberl W , Paternoster F , Achatz F , Schwirtz A , Hahn D. 2013 . On the relevance of residual force enhancement for everyday human movement . J Biomech 46 : 1996 – 2001 . Google Scholar CrossRef Search ADS PubMed Seiberl W , Power GA , Herzog W , Hahn D. 2015 . The stretch‐shortening cycle (SSC) revisited: residual force enhancement contributes to increased performance during fast SSCs of human m. adductor pollicis . Physiol Rep published online (doi:10.14814/phy2.12401). Siebert T , Till O , Blickhan R. 2014a . Work partitioning of transversally loaded muscle: experimentation and simulation . Comp Methods Biomech Biomed Eng 17 : 217 – 29 . Google Scholar CrossRef Search ADS Siebert T , Till O , Stutzig N , Günther M , Blickhan R. 2014b . Muscle force depends on the amount of transversal muscle loading . J Biomech 47 : 1822 – 8 . Google Scholar CrossRef Search ADS Simo JC , Taylor RL , Pister KS. 1985 . Variational and projection methods for the volume constraints in finite deformation elasto-plasticity . Comput Methods Appl Mech Eng 51 : 177 – 208 . Google Scholar CrossRef Search ADS Simo JC , Taylor RL. 1991 . Quasi-incompressible finite elasticity in principal stretches . Comput Methods Appl Mech Eng 85 : 273 – 310 . Google Scholar CrossRef Search ADS Sleboda DA , Roberts TJ. 2017 . Incompressible fluid plays a mechanical role in the development of passive muscle tension . Biol Lett 13 : 20160630. Google Scholar CrossRef Search ADS PubMed Smith LR , Gerace-Fowler L , Lieber RL. 2011 . Muscle extracellular matrix applies a transverse stress on fibers with axial strain . J Biomech 44 : 1618 – 20 . Google Scholar CrossRef Search ADS PubMed Tamura Y , Saito M , Nagato R. 2005 . A new motor model representing the stretch-induced force enhancement and shortening-induced force depression in skeletal muscle . J Biomech 38 : 877 – 84 . Google Scholar CrossRef Search ADS PubMed Till O , Siebert T , Blickhan R. 2014 . Force depression decays during shortening in the medial gastrocnemius of the rat . J Biomech 47 : 1099 – 103 . Google Scholar CrossRef Search ADS PubMed Van Leeuwen JL. 1992 . Muscle function in locomotion. In: Alexander RMcN , editors. Mechanics of animal locomotion . Berlin : Springer . p. 191 – 250 . Van Loocke M , Lyons CG , Simms CK. 2008 . Viscoelastic properties of passive skeletal muscle in compression: stress-relaxation behaviour and constitutive modelling . J Biomech 41 : 1555 – 66 . Google Scholar CrossRef Search ADS PubMed Wakeling JM , Blake OM , Wong I , Rana M , Lee SS. 2011 . Movement mechanics as a determinate of muscle structure, recruitment and coordination . Philos Trans R Soc B 366 : 1554 – 64 . Google Scholar CrossRef Search ADS Wakeling JM , Johnston IA. 1999 . Predicting muscle force generation during fast-starts for the common carp Cyprinus carpio . J Comp Physiol B 169 : 391 – 401 . Google Scholar CrossRef Search ADS Wakeling JM , Lee SS. 2011 . Modelling muscle forces: from scaled fibres to physiological task-groups . Procedia IUTAM 2 : 317 – 26 . Google Scholar CrossRef Search ADS Wakeling JM , Lee SSM , de Boef Miara M , Arnold AS , Biewener AA. 2012 . A muscle’s force depends on the recruitment patterns of its fibers . Ann Biomed Eng 40 : 1708 – 20 . Google Scholar CrossRef Search ADS PubMed Wakeling JM , Randhawa A. 2014 . Transverse strains in muscle fascicles during voluntary contraction: a 2D frequency decomposition of B-mode ultrasound images . Int J Biomed Imaging 2014 : 352910. Google Scholar CrossRef Search ADS PubMed Weis-Fogh T , Alexander RM. 1977 . The sustained power output from striated muscle. In: Pedley TJ , editor. Scale effects in animal locomotion . London : Academic Press . p. 511 – 25 . Weiss JA , Maker BN , Govindjee S. 1996 . Finite element implementation of incompressible, transversely isotropic hyperelasticity . Comput Methods Appl Mech Eng 135 : 107 – 28 . Google Scholar CrossRef Search ADS Wickiewicz TL , Roy RR , Powell PL , Edgerton VR. 1983 . Muscle architecture of the human lower limb . Clin Orthopaed Rel Res 179 : 275 – 83 . Google Scholar CrossRef Search ADS Williams TL. 2010 . A new model for force generation by skeletal muscle, incorporating work-dependent deactivation . J Exp Biol 213 : 643 – 50 . Google Scholar CrossRef Search ADS PubMed Williamson MR , Dial KP , Biewener AA. 2001 . Pectoralis muscle performance during ascending and slow level flight in mallards (Anas platyrhynchos) . J Exp Biol 204 : 495 – 507 . Google Scholar PubMed Wu JZ , Herzog W. 1999 . Modelling concentric contraction of muscle using an improved cross-bridge model . J Biomech 32 : 837 – 48 . Google Scholar CrossRef Search ADS PubMed Yeoh OH. 1993 . Some forms of the strain energy function for rubber . Rubber Chem Technol 66 : 754 – 71 . Google Scholar CrossRef Search ADS Zajac FE. 1989 . Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control . Crit Rev Biomed Eng 17 : 359 – 411 . Google Scholar PubMed Appendix 1 A mixed (P2–P1–P1) Finite Element Method was used to approximate the displacement, dilation, and pressure of the muscle in three dimensions. Piecewise continuous quadratic polynomials (P2) were used to approximate the displacement, whereas piecewise linear polynomials (P1) were used to approximate the pressure and dilation. The deformations of the model were governed by the balance of forces in the system, derived from the system’s total energy. The dilation was added as an unknown using a Lagrange multiplier to represent the internal pressure within the muscle. A three-field formulation was implemented in the library deal.II (Davydov et al. 2017), where the tutorial step-44 was taken as an initial starting point to solve the three resultant dynamic, nonlinear equations. The boundary conditions of the system were set as follows: we imposed zero displacement on some parts on the boundary, and set a traction (superficial stress in the normal direction) free condition on the rest of the boundary. However, a prescribed traction in the x-direction was used on part of this boundary to couple the 3D muscle system with the forces produced by the 1D oscillator. The initial displacement and pressure of the system were set to zero and the initial dilation was set to 1. Prior 3D muscle models using the Finite Element Method have assumed that the inertial effects within muscle are negligible by using quasi-static analysis (e.g., Rahemi et al. 2014). However, to examine the effects of inertia due to tissue mass, we accounted for the muscle’s density and velocity through its kinetic energy. The governing equations were solved using semi-implicit time stepping coupled with a finite element discretization of the space variables. The nonlinear nature of the equations implies that a nonlinear iteration must be applied at every time step to find the subsequent state. As the nonlinear iteration requires an initial guess to start, the converged solution found in the previous time step by the nonlinear solver was used to start the nonlinear iteration in the current state. For the first time step, the initial conditions for the displacement, pressure, and dilation were used as the initial guess. © The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Integrative and Comparative Biology. All rights reserved. For permissions please email: [email protected]. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices)

Journal

Integrative and Comparative Biology – Oxford University Press

Published: Aug 1, 2018

Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Size, History-Dependent, Activation and Three-Dimensional Effects on the Work and Power Produced During Cyclic Muscle Contractions

Size, History-Dependent, Activation and Three-Dimensional Effects on the Work and Power Produced During Cyclic Muscle Contractions

Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Size, History-Dependent, Activation and Three-Dimensional Effects on the Work and Power Produced During Cyclic Muscle Contractions

Size, History-Dependent, Activation and Three-Dimensional Effects on the Work and Power Produced During Cyclic Muscle Contractions

References (114)

Abstract

Journal

Recommended Articles

There are no references for this article.