# Spreading out Muscle Mass within a Hill-Type Model: A Computer Simulation Study

Authors: Michael Günther; Oliver Röhrle; Daniel F. B. Haeufle; Syn Schmitt
Journal: Computational and Mathematical Methods in Medicine
 (2012)

---

## Abstract

It is state of the art that muscle contraction dynamics is adequately
described by a hyperbolic relation between muscle force and
contraction velocity (Hill relation), thereby neglecting muscle
internal mass inertia (first-order dynamics). Accordingly, the
vast majority of modelling approaches also neglect muscle internal
inertia. Assuming that such first-order contraction dynamics yet
interacts with muscle internal mass distribution, this study
investigates two questions: (i) what is the time scale on which the
muscle responds to a force step? (ii) How does this response scale
with muscle design parameters? Thereto, we simulated accelerated
contractions of alternating sequences of Hill-type contractile
elements and point masses. We found that in a typical small muscle the
force levels off after about 0.2 ms, contraction velocity after
about 0.5 ms. In an upscaled version representing bigger mammals'
muscles, the force levels off after about 20 ms, and the
theoretically expected maximum contraction velocity is not
reached. We conclude (i) that it may be indispensable to introduce
second-order contributions into muscle models to understand
high-frequency muscle responses, particularly in bigger
muscles. Additionally, (ii) constructing more elaborate measuring
devices seems to be worthwhile to distinguish viscoelastic and
inertia properties in rapid contractile responses of muscles.


---

## Body

## 1. Introduction

In case the force-velocity relation of a muscle (usually hyperbolic [1]: the Hill relation) is interpreted as a force law and coupled to inertia loads as, for example, in computer models of musculoskeletal multibody systems, it adds first-order dynamics to the systems’ mechanical equations of motion. Usually, the force-velocity relations of a specific isolated muscle preparation (according to [2–4]: earliest known experiments by Jan Swammerdam around 1663; later: e.g., [1, 5–19]) or lumped-muscle assemblies (e.g., [20–22]) are determined through isotonic (e.g., [10, 13, 18]) or isokinetic (e.g., [7, 10, 16, 17]) contractions or by accelerating external inertia loads (e.g., [20, 21]).In the isokinetic condition, a potential external inertia load would be nonaccelerated, whereas the muscle internal velocity distribution could not yet be guaranteed to have cancelled out (maybe internally and locally accelerated). Hypothesising that the muscle consists internally of just two parts arranged in series, the contractile element (CE) and the serial (visco-)elastic element (SEE), the isotonic condition is often chosen to separate an assumed subsequent steady-state response of the CE (and, thus, its isolated properties) from the so-called “initial elastic response” of a muscle. The latter is the finite period (in fibres below a millisecond [13, 23–27]) after a quick release in which the force settles to a new level, as mostly presumed, through mere SEE properties. In analogy to the isokinetic condition, however, the muscle internal force distribution cannot be guaranteed to have cancelled out during an isotonic contraction. That is, both contraction modes are prone to represent in fact nonsteady situations due to potential inertia loads involved, may they come from inertia of the measurement device in series or muscle internal masses themselves.Some studies [28–30] have argued that muscle internal inertia forces can be neglected. All these authors based their argument on assuming that the force within one sarcomere accelerates just one other sarcomere, that is, a very low mass. In real muscle, however, all sarcomeres in parallel (imagine a cut through a muscle’s cross-section) would have to accelerate all masses in series on both sides of the cut. Definitely, force generation of a real muscle cannot work without almost synchronous contraction along its whole mass distribution. Certainly, the force generated by one sarcomere will accelerate the inertia of more than its own and its direct neighbours. And indeed, there are some experimental studies [24–26, 31] in which high-frequency oscillatory responses to step excitations had been found, indicating that muscle inertia may interfere with contraction.As a consequence of this and the fact that finite contraction distances have to be allowed during any contraction experiment, the measured results of both dynamic contraction conditions may delicately depend on muscle internal mass distribution. In fact, all mentioned experimental conditions are meant to finally determine the characteristics of the CE within the muscle. We can conclude that, at least theoretically, the so determined contractile muscle characteristic (force-velocity relation) might not be identified without taking inertia effects into account. Moreover, isolated preparations have been examined solely for invertebrate, arthropod, frog, toad, or small mammal muscles so far, that is, for muscle masses in the range0.0001⋯10g (maximum: 20g in a turkey muscle [32]), thus, much smaller than those of humans or big mammals. Assuming geometric scaling, muscle force would just scale quadratically (cross-sectional area) with body length, whereas muscle mass would scale cubically (volume). Thus, accelerating force per muscle mass should decrease about linearly with body length. Accordingly, it may well be that first-order dynamics does not generally suffice to adequately describe muscle contraction dynamics. Therefore, in our view, at least two questions remain open. (i) Usually, experiments are performed on small muscles, whereas computer models often refer to upscaled muscles. Hence, the first question worth asking seems to be whether and how a mass distribution, scaled up to a design resembling bigger humans’ or mammals’ muscles, would bias the force-velocity relation. (ii) The second question is in regard to all muscle sizes. What is the time scale on which such inertia forces occur, especially compared to the time scale of the “initial elastic response”?To answer these questions, we conducted a computer experiment in which modelling muscle dynamics was purposely reduced to the force-length and force-velocity dependencies of the CE. That is, we neglected any elasticity, whether in parallel or in series, to be able to focus solely on the CE force relationships. Usually, parallel elasticity does not intervene when starting the contraction just below or at about optimal fibre length. Serial elasticity can be neglected in case the contraction does not start in some pre-stretched condition like in, for example, activated isometry. Consequently, we calculated in this study how a muscle, which is modelled on the basis of classical Hill-type CE contraction dynamics, would be biased by and interact with muscle internal second-order dynamics (mass inertia). Thereby, the muscle is assumed to be fully activated initially and subjected to the maximum initial load gradient possible along its length, while being not prestretched: one muscle end is fixed, whereas the other end is released to contract freely.


## 2. Material and Methods

We simulated linear (onedimensional) dynamic muscle contractions using in-house software “simsys” [33–35] which incorporates the solver “de” [36] for time integration. The solver’s absolute and relative error tolerances were set to 10-12.We approximated the continuous mass distribution along the muscle by modelling (Figure1) a finite number IM+1 of discrete point masses (PMs) in an alternating sequence with a corresponding number of IM contractile elements (CEs) indicated by “i”. To keep the model as simple as possible, we further assumed a homogeneous mass distribution, that is, the same 1/(IM+1) portion of the overall muscle mass ℳ is attributed to each PM. For all models with IM=1,2,4,8,16,32 CEs, we examined just one dynamic load situation: a fully active muscle (all activities qi=1, qi∈[0,1]) contracting concentrically with one muscle end fixed to the world, and the other end is released to contract freely, ignoring gravity. This is the most dynamic, however, non-prestretched situation conceivable, because we chose all IM CEs to be at their (common) optimal length (lCE,i=lCE,opt,i=lCE,opt) initially, not shortening (l˙CE,i=0; the dot “·” denotes the time t derivative) and, therefore, initially pulling with their (common) maximum isometric force (FCE,i=Fmax,i=Fmax). Consequently, the only PM to be accelerated initially was the one at the free end, being exposed to the maximum force gradient possible: FCE,IM=Fmax pulling on the one side, whereas there is no force counteracting on the other side. That is, we computed the decay of this singular force step while going through an intermittent scenario of a force (and acceleration) gradient distributed along the whole muscle length lM=∑i=1IMlCE,i. With this decay, we quantified the time scale on which the muscle attunes to any load difference between its two ends, depending on its CE characteristics and internal inertia.Figure 1
The six models used for simulating accelerated contractions of a muscle with an overall mass of6.5g are arranged from the left to the right. The models differ (i) by having split the muscle mass into an increasing number (IM+1) of point masses (PMs; depicted by circles) connected through a corresponding number IM of CEs (depicted by “brushes”) and (ii) by distributing the optimal muscle length lM,opt homogeneously to the optimal CE lengths lCE,opt,i=lM,opt/IM. The sum of the circle areas, symbolising the muscle mass, is constant across the six models. The uppermost PM is always fixed to the world, whereas all other PMs can be accelerated by their adjacent CEs. Note that there is no gravity. Generally, lM=∑i=1IMlCE,i is the overall muscle length lM.For this loading situation, we consider only two different muscle designs. Basic muscle model parameters are the muscle massℳ, the maximum isometric force Fmax of all the CEs and the overall muscle, and the optimal length lCE,opt of all CEs determining the optimal over-all muscle length lM,opt=∑i=1IMlCE,opt,i=IM·lCE,opt. The first muscle represents an averaged assembly of the four plantar flexors of a piglet [21]: ℳ=6.5g, Fmax=30N, lM,opt=0.015m, the latter two parameters characterising static (isometric) muscle properties. Dynamic concentric contraction properties are characterised by two further parameters Arel,Brel (see next paragraph), fixing the maximum concentric contraction velocity vM,max=0.15m/s and the curvature of an assumed hyperbolic force-velocity relation [1]. The second muscle is a scaled version of this piglet muscle with hundredfold mass (ℳ=650g). This comes from scaling the piglet’s length (to lM,opt=0.15m, vM,max=1.5m/s) as well as the cross-sectional area (meaning Fmax=300N), each by a factor of ten. Thus, both basic contraction parameters consistently scale tenfold, whereas the values of the dynamic parameters Arel,Brel are retained unchanged. The human pectineus muscle [37] or the medial head of the horse triceps muscle [38] would be examples of such a muscle design.Any muscle elasticity (parallel and in series to the CE) is neglected. Also, activation dynamicsqi(t) is not taken into account (qi=1 during complete contractions). Concentric contraction dynamics of the CEs as applied in this study has been completely described elsewhere [21, 39]. Here, we give a short summary. The isometric force function of a CE is modelled as FCE,isom,i(lCE,i,qi)=Fmax·qi·fisom(lCE,i), where fisom(lCE,i)∈[0,1] is constructed from two exponential (bell) curve branches (ascending and descending limb; see parameter values in [21]) steadily and differentiably connected at lCE,opt,i. Besides the isometric force parameter Fmax, the CE’s hyperbolic (Hill-type) force-velocity relation is determined by two further normalised (Hill) parameters Arel,i=Ai/FCE,isom,i, Brel,i=Bi/lCE,opt,i. Ai, Bi are the hyperbola’s asymptotes in nonnormalised units. They are generally, just like FCE,isom,i, rather functions of the CE state variables lCE,i, qi than constant parameters [21, 40]. Together, both asymptotes fix vmax,i and the curvature of the force-velocity relation. Our restriction to identical CEs which contract fully activated around optimal lengths is achieved by globally choosing Arel,i=Arel=0.1 and Brel,i=Brel=1.0(1/s) [21] as constant parameters for each CE. Thus, vmax,i=vmax=Brel/Arel·lCE,opt is also (almost: see Section 3) a fixed parameter in our simulations. Therewith, vM,max=IM·vmax scales directly with lM,opt in our model.Generally, the changedF(x,t)=x¨(x,t)·dℳ in the force F(x,t) acting along the muscle longitudinal axis in x-direction runs in proportion to the mass portion dℳ located between position x and x+dx and to the respective current acceleration x¨(x,t) of this mass portion. Note that F(x,t) means the force acting on one side of each mass portion, that is, represents the force in the corresponding adjacent CE. Our assumption of a homogeneous mass distribution (i.e., ℳ/lM=const) is equivalent to a constant mass increment dℳ=ℳ/lM·dx per change in position dx. Thus, in our model we find
(1)dF(x,t)=x¨(x,t)·ℳ/lM(t)·dx.The situation, in which the muscle is fixed atx=0 and left free at x=lM, is characterised by
(2)x¨(0,t)=0,x¨(xCOM,t)=F(0,t)-F(lM,t)ℳ=F|x=0(t)ℳ,
with the boundary conditions F(lM,t)=0 and F(0,t)=F|x=0(t) applying at any time t. Now, if local accelerations were distributed as homogeneously as mass (implying xCOM(t)=lM(t)/2) along the muscle, that is,
(3)x¨(x,t)=a+b·x,
then, from comparing (3) with the two equations (2), we could deduce a=0 and b=2F|x=0(t)/(ℳ·lM(t)). For homogeneously distributed accelerations we, thus, find
(4)x¨|x=lM(t)=x¨(lM,t)=2F|x=0(t)ℳ
as the end point acceleration. Once the acceleration distribution is known, the corresponding force distribution can be calculated according to (1).WithF|x=0(t) being the force difference between the resting (here: at x=0) and the most accelerated (here: at x=lM(t)) muscle part, and with x¨|x=lM(t) being the difference in acceleration between both of these muscle parts, we may define
(5)μeff(t)=F|x=0(t)x¨|x=lM(t)
as the “effective mass” of a muscle for which inertia, forces, and kinematics are distributed along its longitudinal axis. According to (4),
(6)μeff=ℳ2
would be an indicator of accelerations increasing linearly from the fixed to the free muscle end.


## 3. Results

The maximum contraction velocityvmax is generally a function of of activity q [21, 40] and length lCE [21, 40–43]: vmax(lCE,q)=Brel(lCE,q)/Arel(lCE,q)·lCE,opt. Essentially, in our case of a fully active muscle (q=1) always at lCE≤lCE,opt, this means vmax(lCE)=Brel,0/Arel,0·lCE,opt·fisom(lCE) [21, 40] with Arel,0=Arel(lCE=lCE,opt,q=1), Brel,0=Brel(lCE=lCE,opt,q=1) [21]. However, the time periods analysed here are short enough for the CEs not to shorten by more than 0.25·lCE,opt. Thus, vmax is well approximated by its absolute maximum value vmax,0=Brel,0/Arel,0·lCE,opt (see Section 3.3 below: vmax>0.95·vmax,0).


### 3.1. Smaller Muscle

In Figure2, the simulation results of a contraction with one end fixed are shown for the smaller muscle (ℳ=6.5g, Fmax=30N, lM,opt=0.015m). After about 0.5ms, the PM at the free end of the muscle (and, thus, the whole muscle) has reached vmax (Figure 2(a)). Then, the muscle has contracted by about 0.01·lM,opt. Modelling just one CE (and, thus, one PM free to move: IM=1) already approximates the rise in contraction velocity of a continuous mass distribution. An analytic solution for the force step response of one CE accelerating the overall muscle mass (see Appendix A) confirms our numerical results. The typical time constant for this muscle, when approaching vmax, is 3.6·10-4s. The initial response at v=0 is faster by two orders of magnitude.Computer simulation of accelerated contractions of six muscle models with the same overall mass (6.5g), maximum isometric force (30N), optimal length (0.015m), and maximum contraction velocity (0.15m/s). The muscles were fixed at one end and always fully active (q=1). Models differ just with respect to the number of accelerated discrete point masses approximating a continuous distribution of muscle mass. The point masses were connected by an equal number IM (see insets in Figure 2) of contractile elements. Their respective optimal lengths lCE,opt,i were chosen equal to the optimal muscle length lM,opt divided by IM. The graph depicts the velocity of the point mass at the free end (a), the force at the fixed end (b), and the effective mass (c), that is, the ratio of force at the fixed end and acceleration of the point mass at the free end, versus time. The effective mass to be expected for an exactly linear acceleration distribution along the muscle would be half of the muscle mass (MassM=ℳ; (6): μeff=ℳ/2=3.25g). The analytic solution for one CE accelerating one point mass predicts (see Appendix A) a typical time of 3.6·10-4s for this muscle to approach vmax.
(a)
(b)
(c)The mass distribution is well approximated by modelling few tens of PMs as the difference between the time courses of theIM=16 and IM=32 models is not discernible any more in Figure 2. Due to the hyperbolic force-velocity relation the force exerted at the fixed end drops even sharper (Figure 2(b)) than velocity rises. The IM=1 model is enough to predict the dynamics of this force decay. The ratio between the fixed end force and the acceleration of the PM at the free end (i.e., their differences between both ends, resp.), which defines the effective mass μeff (5), reaches its stationary value after about the very time in which the force has dropped to almost zero (0.2ms: Figure 2(c)).Not surprisingly, we always get exactlyμeff=ℳ/2 (6) in case of IM=1 because (3) is consistently fulfilled throughout. In case of IM>1, the initial value is μeff=ℳ/IM (i.e., limIM→∞(μeff)=0) because maximum isometric force acts initially on both sides of each PM, except for the very PM at the free end. With increasing IM the stationary effective mass value converges to a little more than μeff=ℳ/2. Note that μeff=ℳ/2 is predicted in (4) exactly (6), if the PM accelerations increase linearly with the distance x from the fixed to the free end (3). A stationary effective mass nearby μeff=ℳ/2 indicates such an approximately linear acceleration (and, thus, force) distribution already after about 0.2ms. Subsequently, the muscle still accelerates (see time course of the velocity), albeit slightly and at low forces.


### 3.2. Scaling

Here, we condense the scaling characteristics of the muscle model, going without explicit plots. In our one-dimensional model, which assumes a homogeneous mass distribution along the muscle, the parametersArel,Brel of the hyperbolic force-velocity relation fully determine the dimensionless time courses of forces and kinematic variables (length, velocity, acceleration). As a consequence thereof, the effective mass scales in proportion to force and time, respectively, and reciprocally to length (5). When varying solely one single muscle parameter, we find the following strictly linear scaling:(i)
forces just scale withFmax;(ii)
kinematic variables just scale withlM,opt;(iii)
time (i.e., length of time responses) scales withℳ/Fmax and with lM,opt (see also Appendix A).As we see in the next paragraph, time does still scale approximately with fixedℳ·lM,opt/Fmax; however, the dimensionless time courses of all variables are biased. Note that we express scaling with the maximum isometric force Fmax instead of the current isometric force FCE,isom because (i) Fmax is a measure of cross-sectional area in muscles (maximum stress is nearly a constant across muscle dimensions), and (ii) our simulations are performed with a fully activated muscle around the optimal length.


### 3.3. Bigger Muscle

The course of the force (Figure3(b)) in a muscle representing an upscaled version (ℳ=650g, Fmax=300N, lM,opt=0.15m) of the one in Section 3.1 proves that the time scale is hundredfold, as expected due to tenfold ℳ/Fmax and ten-fold lM,opt. Yet, the courses of the velocity (Figure 3(a)) and effective mass (Figure 3(c)) also demonstrate that contraction kinematics is biased. The effective mass does not saturate, as in the smaller muscle (Figure 2(c)), before end point acceleration has reached zero. This occurs at about 35ms for IM=1 and at about 30ms for IM=32 (approximating continuous mass distribution IM→∞), leading to a pole in the effective mass (Figure 3(c)). This also means that, due to mass inertia, the force in the (outermost: i=IM) CE acting on the PM at the free end vanishes earlier than that of the opposite (i=1) CE at the fixed end. As a consequence, in the instant when the PM at the free end is not accelerated any more (see the pole in effective mass), since the CE force has dropped to zero (thus, the outermost CE has reached vmax), the inner elements (closer to the fixed end) shorten still slightly below vmax. Note that the muscle has contracted by about 0.24·lM,opt after 30ms (compared to the smaller muscle: there, vM,max is reached already after about 0.5ms with about 0.01·lM,opt change in length; see Section 3.1). That is, at this instant vM,max has additionally dropped by (but not more than: due to broad ascending branch of force-length relation [21]) about five percent below its absolute maximum value vM,max,0=Brel/Arel·lM,opt due to its length dependency. This latter effect superposes with the fact that the inner elements still contract at submaximal velocity at this instant. As a result, the theoretical absolute maximum vM,max,0 is never reached.Computer simulation of accelerated contractions of six muscle models with the same overall mass (650g), maximum isometric force (300N), optimal length (0.15m), and maximum contraction velocity (1.5m/s). The muscles were fixed at one end and always fully active (q=1). Models differ just with respect to the number of accelerated discrete point masses approximating a continuous distribution of muscle mass. The point masses were connected by an equal number IM (see insets in Figure 3) of contractile elements. Their respective optimal lengths lCE,opt,i were chosen equal to the optimal muscle length lM,opt divided by IM. The graph depicts velocity (a), force (b), and effective mass (c) versus time. The effective mass to be expected for an exactly linear acceleration distribution along the muscle would be half of the muscle mass (MassM=ℳ; (6): μeff=ℳ/2=325g). The analytic solution for one CE accelerating one point mass predicts (see Appendix A) a typical time of 3.6·10-2s for this muscle to approach vmax. Note the hundredfold muscle mass, tenfold maximum isometric force and optimal length, respectively, and hundredfold time scale as compared to the results presented in Figure 2.
(a)
(b)
(c)In our model, muscle velocity decreases after that instant of zero force at the free end, occurring at30⋯35ms. The corresponding energy dissipation originates from those CEs that contract beyond their vmax,i(lCE,i) (due to kinetic energy of the adjacent PM). In this loading situation, we allowed such low compressive forces (FCE,i<0) that would arise from steadily extrapolating the hyperbolic force-velocity relation to vCE,i>vmax,i(lCE,i). Physiologically, some low compressive forces should in fact occur in real muscles as these are always somehow constrained in transverse direction. Whether this steady extrapolation is backed physiologically remains open owing to a lack of literature data. Due to such compressive forces, we find a maximum in the free end’s velocity, which is lower than the absolute maximum value vM,max,0 that would be theoretically expected from pure CE properties. Even without compressive forces, the theoretical vM,max,0 would never be reached. This is due to the fact that mere interaction of mass inertia and force-velocity relation delay the increase in contraction velocity long enough for part of the CEs to shorten into the ascending branch of their force-length dependency.


## 3.1. Smaller Muscle

In Figure2, the simulation results of a contraction with one end fixed are shown for the smaller muscle (ℳ=6.5g, Fmax=30N, lM,opt=0.015m). After about 0.5ms, the PM at the free end of the muscle (and, thus, the whole muscle) has reached vmax (Figure 2(a)). Then, the muscle has contracted by about 0.01·lM,opt. Modelling just one CE (and, thus, one PM free to move: IM=1) already approximates the rise in contraction velocity of a continuous mass distribution. An analytic solution for the force step response of one CE accelerating the overall muscle mass (see Appendix A) confirms our numerical results. The typical time constant for this muscle, when approaching vmax, is 3.6·10-4s. The initial response at v=0 is faster by two orders of magnitude.Computer simulation of accelerated contractions of six muscle models with the same overall mass (6.5g), maximum isometric force (30N), optimal length (0.015m), and maximum contraction velocity (0.15m/s). The muscles were fixed at one end and always fully active (q=1). Models differ just with respect to the number of accelerated discrete point masses approximating a continuous distribution of muscle mass. The point masses were connected by an equal number IM (see insets in Figure 2) of contractile elements. Their respective optimal lengths lCE,opt,i were chosen equal to the optimal muscle length lM,opt divided by IM. The graph depicts the velocity of the point mass at the free end (a), the force at the fixed end (b), and the effective mass (c), that is, the ratio of force at the fixed end and acceleration of the point mass at the free end, versus time. The effective mass to be expected for an exactly linear acceleration distribution along the muscle would be half of the muscle mass (MassM=ℳ; (6): μeff=ℳ/2=3.25g). The analytic solution for one CE accelerating one point mass predicts (see Appendix A) a typical time of 3.6·10-4s for this muscle to approach vmax.
(a)
(b)
(c)The mass distribution is well approximated by modelling few tens of PMs as the difference between the time courses of theIM=16 and IM=32 models is not discernible any more in Figure 2. Due to the hyperbolic force-velocity relation the force exerted at the fixed end drops even sharper (Figure 2(b)) than velocity rises. The IM=1 model is enough to predict the dynamics of this force decay. The ratio between the fixed end force and the acceleration of the PM at the free end (i.e., their differences between both ends, resp.), which defines the effective mass μeff (5), reaches its stationary value after about the very time in which the force has dropped to almost zero (0.2ms: Figure 2(c)).Not surprisingly, we always get exactlyμeff=ℳ/2 (6) in case of IM=1 because (3) is consistently fulfilled throughout. In case of IM>1, the initial value is μeff=ℳ/IM (i.e., limIM→∞(μeff)=0) because maximum isometric force acts initially on both sides of each PM, except for the very PM at the free end. With increasing IM the stationary effective mass value converges to a little more than μeff=ℳ/2. Note that μeff=ℳ/2 is predicted in (4) exactly (6), if the PM accelerations increase linearly with the distance x from the fixed to the free end (3). A stationary effective mass nearby μeff=ℳ/2 indicates such an approximately linear acceleration (and, thus, force) distribution already after about 0.2ms. Subsequently, the muscle still accelerates (see time course of the velocity), albeit slightly and at low forces.


## 3.2. Scaling

Here, we condense the scaling characteristics of the muscle model, going without explicit plots. In our one-dimensional model, which assumes a homogeneous mass distribution along the muscle, the parametersArel,Brel of the hyperbolic force-velocity relation fully determine the dimensionless time courses of forces and kinematic variables (length, velocity, acceleration). As a consequence thereof, the effective mass scales in proportion to force and time, respectively, and reciprocally to length (5). When varying solely one single muscle parameter, we find the following strictly linear scaling:(i)
forces just scale withFmax;(ii)
kinematic variables just scale withlM,opt;(iii)
time (i.e., length of time responses) scales withℳ/Fmax and with lM,opt (see also Appendix A).As we see in the next paragraph, time does still scale approximately with fixedℳ·lM,opt/Fmax; however, the dimensionless time courses of all variables are biased. Note that we express scaling with the maximum isometric force Fmax instead of the current isometric force FCE,isom because (i) Fmax is a measure of cross-sectional area in muscles (maximum stress is nearly a constant across muscle dimensions), and (ii) our simulations are performed with a fully activated muscle around the optimal length.


## 3.3. Bigger Muscle

The course of the force (Figure3(b)) in a muscle representing an upscaled version (ℳ=650g, Fmax=300N, lM,opt=0.15m) of the one in Section 3.1 proves that the time scale is hundredfold, as expected due to tenfold ℳ/Fmax and ten-fold lM,opt. Yet, the courses of the velocity (Figure 3(a)) and effective mass (Figure 3(c)) also demonstrate that contraction kinematics is biased. The effective mass does not saturate, as in the smaller muscle (Figure 2(c)), before end point acceleration has reached zero. This occurs at about 35ms for IM=1 and at about 30ms for IM=32 (approximating continuous mass distribution IM→∞), leading to a pole in the effective mass (Figure 3(c)). This also means that, due to mass inertia, the force in the (outermost: i=IM) CE acting on the PM at the free end vanishes earlier than that of the opposite (i=1) CE at the fixed end. As a consequence, in the instant when the PM at the free end is not accelerated any more (see the pole in effective mass), since the CE force has dropped to zero (thus, the outermost CE has reached vmax), the inner elements (closer to the fixed end) shorten still slightly below vmax. Note that the muscle has contracted by about 0.24·lM,opt after 30ms (compared to the smaller muscle: there, vM,max is reached already after about 0.5ms with about 0.01·lM,opt change in length; see Section 3.1). That is, at this instant vM,max has additionally dropped by (but not more than: due to broad ascending branch of force-length relation [21]) about five percent below its absolute maximum value vM,max,0=Brel/Arel·lM,opt due to its length dependency. This latter effect superposes with the fact that the inner elements still contract at submaximal velocity at this instant. As a result, the theoretical absolute maximum vM,max,0 is never reached.Computer simulation of accelerated contractions of six muscle models with the same overall mass (650g), maximum isometric force (300N), optimal length (0.15m), and maximum contraction velocity (1.5m/s). The muscles were fixed at one end and always fully active (q=1). Models differ just with respect to the number of accelerated discrete point masses approximating a continuous distribution of muscle mass. The point masses were connected by an equal number IM (see insets in Figure 3) of contractile elements. Their respective optimal lengths lCE,opt,i were chosen equal to the optimal muscle length lM,opt divided by IM. The graph depicts velocity (a), force (b), and effective mass (c) versus time. The effective mass to be expected for an exactly linear acceleration distribution along the muscle would be half of the muscle mass (MassM=ℳ; (6): μeff=ℳ/2=325g). The analytic solution for one CE accelerating one point mass predicts (see Appendix A) a typical time of 3.6·10-2s for this muscle to approach vmax. Note the hundredfold muscle mass, tenfold maximum isometric force and optimal length, respectively, and hundredfold time scale as compared to the results presented in Figure 2.
(a)
(b)
(c)In our model, muscle velocity decreases after that instant of zero force at the free end, occurring at30⋯35ms. The corresponding energy dissipation originates from those CEs that contract beyond their vmax,i(lCE,i) (due to kinetic energy of the adjacent PM). In this loading situation, we allowed such low compressive forces (FCE,i<0) that would arise from steadily extrapolating the hyperbolic force-velocity relation to vCE,i>vmax,i(lCE,i). Physiologically, some low compressive forces should in fact occur in real muscles as these are always somehow constrained in transverse direction. Whether this steady extrapolation is backed physiologically remains open owing to a lack of literature data. Due to such compressive forces, we find a maximum in the free end’s velocity, which is lower than the absolute maximum value vM,max,0 that would be theoretically expected from pure CE properties. Even without compressive forces, the theoretical vM,max,0 would never be reached. This is due to the fact that mere interaction of mass inertia and force-velocity relation delay the increase in contraction velocity long enough for part of the CEs to shorten into the ascending branch of their force-length dependency.


## 4. Discussion

### 4.1. Inertia in Finite Element Models

Besides modelling skeletal muscles using Hill-type models, there exist many other modelling approaches. For this study, Hodgkin-Huxley models aiming to investigate electrophysiological cellular properties [44, 45], structural three-dimensional models using the governing equations of finite elasticity to describe muscle mechanics [46–50] and combinations thereof [51, 52] seem to be of primary interest. Hodgkin-Huxley-like models typically focus only on the electrophysiological properties and hence ignore muscle mass entirely. In the case of structural models, mass is included within the governing equations of finite elasticity through mass density, representing weight forces and inertia terms. Except for [47, 53], inertia terms are typically ignored. In most cases, quasi-static approaches are adopted, which simplify, from a computational point of view, the governing equations of finite elasticity. Ignoring muscle inertia and appealing to quasi-static approaches might be a valid assumption when considering isometric contractions or slow motions. However, as shown here, there exist contraction scenarios in which it is no longer justifiable to ignore muscle inertia in structural models.


### 4.2. Contraction Dynamics of Second Order?

Muscle inertia must be taken into account in specific modes of contraction, particularly in responses to impact loads [54–60]. In such a case, whole muscle contraction dynamics is of second order. In concentric contractions, muscle masses do not seem to play a relevant role because the characteristic of contraction dynamics was originally found to be a force-velocity relation, namely, the hyperbolic Hill relation [1]. The relation definitely depends explicitly on the external load (force), whereas it is explicitly independent of internal muscle inertia. However, force-velocity relations have ever been determined in preparations of isolated small muscles of a few grams rather than in isolated muscles of big mammals. So far, in these preparations of isolated small muscles inertia forces have not been found to bias contraction dynamics, at least not in muscle states at a few milliseconds after release.In contrast, our simulation results predict that maximum concentric contraction velocity, being one parameter in the Hill relation, cannot be directly measured in whole muscle preparations of big mammals. The other way round, if the hyperbolic force-velocity relation really mapped the concentric muscular contraction dynamics microscopically [61, 62] and macroscopically [1, 63], then a big muscle would never reach the corresponding maximum contraction velocity at zero load (Figure 3(a)). In big muscles, therefore, muscle internal inertia should have to be factored in, even during concentric contractions, when striving to identify the contractile properties only. Consequently, a first-order force-velocity relation would not be expected to fully describe the macroscopic dynamics of larger muscles.


### 4.3. Rapid Responses in Isotonic after Load Experiments: Just Elastic?

Contractions of frog fibre bundles withℳ=0.06g, Fmax=0.3N, lCE,opt=3.0cm have been examined in a former study [13]. The muscle force per mass ratio between our small (piglet-like) muscle and a frog fibre bundle, which has a hundredth of our muscle’s cross-sectional area, is about the same. Consequently, time characteristics should approximately be equal. The authors [13] stated: “We conclude that the change of velocity follows the change in force very quickly probably in less than 1 msec and certainly in less than 6 msec.” According to our simulations which include mass inertia, vmax is reached after about 0.5ms. Thus, the question arises whether the time for levelling off at the new force in isotonic quick-release contractions is really just due to serial (visco-)elasticity rather than a mixture of (visco-)elastic and inertia properties.In other experiments on single frog fibres [64], again with a comparable force per mass ratio, the measuring device had an internal elasticity (in series with the muscle) with an unloaded eigenfrequency of about 200Hz, corresponding to a time period of 5ms. As an “initial elastic response” within about 1ms can be estimated from their plots on contracting muscle, we would conclude that the muscle stiffness must have been much higher than the device stiffness. Without any device stiffness the usually proposed “initial elastic response” should have occurred even faster than within 1ms. If this was true, the rapid response to force steps would even more likely be a superposition of (visco-)elastic and inertia effects.


### 4.4. Experimental Implications for Studying Inertia in Muscle Preparations

In case one aims to examine inertia effects in the rapid responses of small muscle preparations (as frog, toad, piglet) experimentally, measuring devices of very low inertia, high stiffness, and at least10-6s time resolution were indispensable. A force step should be favoured as the experimental protocol, rather than a length or velocity step. This is, because the interpretation of the experimental results of kinematic disturbances depends on assuming unknown material characteristics, whereas a force step is a disturbance provoking the contractile kinematics in a natural cause-effect chain, that is, not interfering with the material characteristics during the (ideally infinitely short) step. As to our knowledge, there has been one modern study [27] fulfilling already three of the experimental requirements: it used a force-step protocol for frog muscles with time resolution of 10-6s. Unfortunately, the resonance frequency of the apparatus was not more than 5·104Hz. Eventually, they applied force steps at a length of about 10-4s. To resolve inertia effects, a setup tuned for step lengths of at least 10-6s would be required (ideally 10-7s), while choosing for a time resolution of at least the step length, accordingly. As an alternative, longer fibres could be used. Additionally, any tendinous material in series to the fibres should be removed as good as possible, and the stiffness of the apparatus should always be clearly higher than that of any remaining tendon material.


### 4.5. Muscle-Inertia Interaction in Nonstationary Conditions

If the hyperbolic force-velocity relationF(v) is linearised at any operating point, an analytic solution can be derived for the dynamic response of the muscle model (including an accelerated mass) to a step in the external force (see Appendix A). This solution predicts an aperiodic exponential decay of the force step towards a limit velocity. Furthermore, limit velocity and time-constant τ of the decay vary along the operating points F(v). In fact, τ is simply the reciprocal of the current slope of the force-velocity hyperbola times the accelerated mass. Therefore, depending on the curvature, the time-constant τ of the decay can vary one to two orders of magnitude across all possible contraction states of a muscle. Thus, adaptations to force steps happen quickly for contractions with high force and slow velocity (i.e., fastest around the isometric state), while it happens much slower at low forces and high velocities. In the example calculated in Appendix A (parameters of our smaller muscle from Section 3.1), we find τ=3.0·10-6 at v=0, and τ=3.6·10-4 at v=-vmax.This step response is induced by the interaction of inertia and force-velocity relation. In a technical sense, this characteristic represents a low-pass filter with a state-dependent cut-off frequency. Mechanically spoken, the force-velocity relation is equivalent to a damper with a state-dependent damping coefficient. This is in accordance with a recently proposed macroscopic model for the CE [65]. In that model, the CE consists of three simple mechanical elements: an active energy source (AE), a parallel damper (PDE), and a serial element (SE). There, the PDE is assumed to be a viscous damper, with the damping coefficient depending linearly on the CE force (high damping for high CE forces). With this assumption, the model predicts the operating points on the hyperbolic force-velocity relation from the force equilibrium of the three elements. Obviously, the differential equation defining a hyperbola (see (A.8) in Appendix A) is mathematically equivalent to a manifold of force equilibria of the CE elements such as proposed in [65]. All in all, it seems the muscle’s hyperbolic force-velocity relation can be viewed from at least three perspectives, which are equivalent. (i) It can be considered a linearly force-dependent damping coefficient for external inertia loads. (ii) It performs as a low-pass filter for whole muscle contraction exposed to external steps in force. And (iii) it represents a potential to compensate (homogenise) muscle internal force and velocity distributions.In first attempts to implement the above-mentioned CE as a mechanical assembly [66, 67], the SE was realised with a linear mechanical spring introducing additional elasticity (in series to possible SEE elasticity representing tendon and aponeurosis). Such a CE-internal elasticity from the SE superposes oscillations to the above-described aperiodic exponential decay. Alternatively, this CE model allows also to detail the assumed SE properties as viscoelastic or purely damping. Each would be expected to reduce any oscillatory behaviour. Both in technical implementations as proposed above and in experimental preparations of physiological muscle, these effects should be as much visible and, thus, verifiable. Such experiments could, therefore, help to further reveal the internal structure of assemblies of active muscle fibres and help to push forward CE models proposed to represent them. Heat measurements on muscles are another means to check for the validity of different CE models with the same force-velocity relation but different internal designs (element properties) [65]. We would expect that heat measurements should also be prone to inertia effects, may they come from muscle internal masses or potential inertia of the measuring device [68].As a first approximation we can usually consider muscle massℳ to be scaled geometrically by either changing its cross-sectional area S or its length as represented by its optimal fibre length lM,opt. This implies that both mass density ρ and maximum stress σmax are about invariant parameters of muscles. Thus, in a real muscle its mass is not an independent parameter as in our model, and the rule to scale contraction times in proportion to ℳ·lM,opt/Fmax as expressed in Section 3.1 and in (A.11) and (A.12) can be condensed to (A.13). We should be aware of the fact that (A.11) and (A.12) do apply to any contracting muscle, may it accelerate just its internal mass (as analysed numerically in this paper) or an external mass that is usually much higher than its own mass. In contrast, (A.13) applies to the specific case of a muscle contracting just against its internal mass under the above-mentioned approximative assumptions of constant density and maximum stress. Thus, (A.13) tells us that the contraction time for a muscle exposed to a step in force, rather than analysed when accelerating an external load, scales in fact with the square of its fibre length. We would, therefore, expect that the effect of not reaching the theoretical maximum contraction velocity vmax (Figure 3(a)) would occur in any muscle preparation with fibre lengths that exceed a value somewhere between 1.5cm and 15cm, regardless of muscle mass. In the end, it seems that it is not the value of the muscle mass itself but the fibre length that scales the response times to force steps (even quadratically, see (A.13)). The effect occurs because contraction times for the whole muscle increase due to the force gradient decreasing with muscle extension, whereas the outermost muscle portions (released immediately in a force step) are always accelerated to their local vmax (and already decelerated consecutively due to local compression) much faster on an invariant time scale. For a 15cm fibre, a muscle is predicted to reach just 93% of vmax.We dare this preliminary prediction subject to future theoretical analyses explicitly including muscle internal elasticity. Serial elasticity might not really help to further speed whole muscle acceleration up. Parallel elasticity, however, might be effective because it adds accelerative forces to the force-velocity relation. Yet, interaction with damping like the force-velocity relation itself and with the mass distribution in series might reduce its accelerating effect down to a secondary level. Parallel elasticity is also a challenge to experimenters because it usually biases contraction measurements as soon as they are performed above or even at optimal fibre lengths. To make our theoretical analysis of the force-velocity relation contributing to accelerated contractions as comparable as possible with experiments, we chose initial conditions at which parallel elasticity is (almost) negligible, that is, at optimal lengths. As a result, our bigger muscle has already contracted by about0.24·lM,opt in the instant when it starts to decelerate its all over contraction (see first paragraph of Section 3.3). Because our model force-length relation has an unusually broad ascending branch [21] the isometric force for a CE shortened by considerable 24% has only decreased a little to about 95% of its absolute maximum. Therefore, it seems that reaching a contraction velocity of just 0.93·vmax is due to a mixture of local compressive forces and the force-length relation in our model simulations. Additional parallel elasticity would only have the potential to compensate for the latter relation. We would, therefore, expect that the phenomenon of not reaching theoretical vmax in big muscles should also occur when starting at initial lengths above the optimal fibre length.Finally, we would like to point to another interesting aspect. The propagation velocity of muscle excitation on its surface is about4m/s in vertebrate cross-striated muscles [69]. In the smaller muscle, this means that an excitation spike would need about an order of magnitude more time to spread along the complete length of a fibre than to reach vmax after a maximum force step down to zero: for 1.5cm fibres we calculated τv=-vmax=3.6·10-4s, and we find 0.015m/(4(m/s))=3.8·10-3s. As τv=-vmax scales quadratically with the fibre length (see Appendix A) we calculate, however, τv=-vmax=3.6·10-2s in the bigger muscle with a tenfold fibre length. In 15cm fibres the estimated propagation time of 0.15m/(4(m/s))=3.8·10-2s from one end to the other would, thus, already slightly exceed τv=-vmax. Accordingly, when thinking of the muscle response to a twitch as a response to an internal step in force, the propagation velocity, rather than inertia effects, may in fact limit the twitch time for a rise in force up to about 7cm fibre length (assumed the only end-plate is located in the middle of the fibre). Yet, even in shorter fibres it may be that inertia, rather than excitation propagation, could limit twitch response times because τv=-vmax increases linearly with muscle stress (see (A.13)), and twitch stress should be clearly lower than maximum stress during a fused tetanus. Additionally, any twitch time limitation by excitation propagation would be brought nearer to inertia limitation if a couple of end-plates were distributed along a fibre. All in all, excitation propagation contributing to the electromechanical delay restricts response times when a rapid rise to a higher force level, that is, to an increased activity is required [39]. In contrast, inertia effects, rather than excitation propagation, should be a time limiting factor in physiological loading situations when rapid changes in external forces occur at some given activity level at which Ca2+ is sufficiently available anyway along the whole fibre.


## 4.1. Inertia in Finite Element Models

Besides modelling skeletal muscles using Hill-type models, there exist many other modelling approaches. For this study, Hodgkin-Huxley models aiming to investigate electrophysiological cellular properties [44, 45], structural three-dimensional models using the governing equations of finite elasticity to describe muscle mechanics [46–50] and combinations thereof [51, 52] seem to be of primary interest. Hodgkin-Huxley-like models typically focus only on the electrophysiological properties and hence ignore muscle mass entirely. In the case of structural models, mass is included within the governing equations of finite elasticity through mass density, representing weight forces and inertia terms. Except for [47, 53], inertia terms are typically ignored. In most cases, quasi-static approaches are adopted, which simplify, from a computational point of view, the governing equations of finite elasticity. Ignoring muscle inertia and appealing to quasi-static approaches might be a valid assumption when considering isometric contractions or slow motions. However, as shown here, there exist contraction scenarios in which it is no longer justifiable to ignore muscle inertia in structural models.


## 4.2. Contraction Dynamics of Second Order?

Muscle inertia must be taken into account in specific modes of contraction, particularly in responses to impact loads [54–60]. In such a case, whole muscle contraction dynamics is of second order. In concentric contractions, muscle masses do not seem to play a relevant role because the characteristic of contraction dynamics was originally found to be a force-velocity relation, namely, the hyperbolic Hill relation [1]. The relation definitely depends explicitly on the external load (force), whereas it is explicitly independent of internal muscle inertia. However, force-velocity relations have ever been determined in preparations of isolated small muscles of a few grams rather than in isolated muscles of big mammals. So far, in these preparations of isolated small muscles inertia forces have not been found to bias contraction dynamics, at least not in muscle states at a few milliseconds after release.In contrast, our simulation results predict that maximum concentric contraction velocity, being one parameter in the Hill relation, cannot be directly measured in whole muscle preparations of big mammals. The other way round, if the hyperbolic force-velocity relation really mapped the concentric muscular contraction dynamics microscopically [61, 62] and macroscopically [1, 63], then a big muscle would never reach the corresponding maximum contraction velocity at zero load (Figure 3(a)). In big muscles, therefore, muscle internal inertia should have to be factored in, even during concentric contractions, when striving to identify the contractile properties only. Consequently, a first-order force-velocity relation would not be expected to fully describe the macroscopic dynamics of larger muscles.


## 4.3. Rapid Responses in Isotonic after Load Experiments: Just Elastic?

Contractions of frog fibre bundles withℳ=0.06g, Fmax=0.3N, lCE,opt=3.0cm have been examined in a former study [13]. The muscle force per mass ratio between our small (piglet-like) muscle and a frog fibre bundle, which has a hundredth of our muscle’s cross-sectional area, is about the same. Consequently, time characteristics should approximately be equal. The authors [13] stated: “We conclude that the change of velocity follows the change in force very quickly probably in less than 1 msec and certainly in less than 6 msec.” According to our simulations which include mass inertia, vmax is reached after about 0.5ms. Thus, the question arises whether the time for levelling off at the new force in isotonic quick-release contractions is really just due to serial (visco-)elasticity rather than a mixture of (visco-)elastic and inertia properties.In other experiments on single frog fibres [64], again with a comparable force per mass ratio, the measuring device had an internal elasticity (in series with the muscle) with an unloaded eigenfrequency of about 200Hz, corresponding to a time period of 5ms. As an “initial elastic response” within about 1ms can be estimated from their plots on contracting muscle, we would conclude that the muscle stiffness must have been much higher than the device stiffness. Without any device stiffness the usually proposed “initial elastic response” should have occurred even faster than within 1ms. If this was true, the rapid response to force steps would even more likely be a superposition of (visco-)elastic and inertia effects.


## 4.4. Experimental Implications for Studying Inertia in Muscle Preparations

In case one aims to examine inertia effects in the rapid responses of small muscle preparations (as frog, toad, piglet) experimentally, measuring devices of very low inertia, high stiffness, and at least10-6s time resolution were indispensable. A force step should be favoured as the experimental protocol, rather than a length or velocity step. This is, because the interpretation of the experimental results of kinematic disturbances depends on assuming unknown material characteristics, whereas a force step is a disturbance provoking the contractile kinematics in a natural cause-effect chain, that is, not interfering with the material characteristics during the (ideally infinitely short) step. As to our knowledge, there has been one modern study [27] fulfilling already three of the experimental requirements: it used a force-step protocol for frog muscles with time resolution of 10-6s. Unfortunately, the resonance frequency of the apparatus was not more than 5·104Hz. Eventually, they applied force steps at a length of about 10-4s. To resolve inertia effects, a setup tuned for step lengths of at least 10-6s would be required (ideally 10-7s), while choosing for a time resolution of at least the step length, accordingly. As an alternative, longer fibres could be used. Additionally, any tendinous material in series to the fibres should be removed as good as possible, and the stiffness of the apparatus should always be clearly higher than that of any remaining tendon material.


## 4.5. Muscle-Inertia Interaction in Nonstationary Conditions

If the hyperbolic force-velocity relationF(v) is linearised at any operating point, an analytic solution can be derived for the dynamic response of the muscle model (including an accelerated mass) to a step in the external force (see Appendix A). This solution predicts an aperiodic exponential decay of the force step towards a limit velocity. Furthermore, limit velocity and time-constant τ of the decay vary along the operating points F(v). In fact, τ is simply the reciprocal of the current slope of the force-velocity hyperbola times the accelerated mass. Therefore, depending on the curvature, the time-constant τ of the decay can vary one to two orders of magnitude across all possible contraction states of a muscle. Thus, adaptations to force steps happen quickly for contractions with high force and slow velocity (i.e., fastest around the isometric state), while it happens much slower at low forces and high velocities. In the example calculated in Appendix A (parameters of our smaller muscle from Section 3.1), we find τ=3.0·10-6 at v=0, and τ=3.6·10-4 at v=-vmax.This step response is induced by the interaction of inertia and force-velocity relation. In a technical sense, this characteristic represents a low-pass filter with a state-dependent cut-off frequency. Mechanically spoken, the force-velocity relation is equivalent to a damper with a state-dependent damping coefficient. This is in accordance with a recently proposed macroscopic model for the CE [65]. In that model, the CE consists of three simple mechanical elements: an active energy source (AE), a parallel damper (PDE), and a serial element (SE). There, the PDE is assumed to be a viscous damper, with the damping coefficient depending linearly on the CE force (high damping for high CE forces). With this assumption, the model predicts the operating points on the hyperbolic force-velocity relation from the force equilibrium of the three elements. Obviously, the differential equation defining a hyperbola (see (A.8) in Appendix A) is mathematically equivalent to a manifold of force equilibria of the CE elements such as proposed in [65]. All in all, it seems the muscle’s hyperbolic force-velocity relation can be viewed from at least three perspectives, which are equivalent. (i) It can be considered a linearly force-dependent damping coefficient for external inertia loads. (ii) It performs as a low-pass filter for whole muscle contraction exposed to external steps in force. And (iii) it represents a potential to compensate (homogenise) muscle internal force and velocity distributions.In first attempts to implement the above-mentioned CE as a mechanical assembly [66, 67], the SE was realised with a linear mechanical spring introducing additional elasticity (in series to possible SEE elasticity representing tendon and aponeurosis). Such a CE-internal elasticity from the SE superposes oscillations to the above-described aperiodic exponential decay. Alternatively, this CE model allows also to detail the assumed SE properties as viscoelastic or purely damping. Each would be expected to reduce any oscillatory behaviour. Both in technical implementations as proposed above and in experimental preparations of physiological muscle, these effects should be as much visible and, thus, verifiable. Such experiments could, therefore, help to further reveal the internal structure of assemblies of active muscle fibres and help to push forward CE models proposed to represent them. Heat measurements on muscles are another means to check for the validity of different CE models with the same force-velocity relation but different internal designs (element properties) [65]. We would expect that heat measurements should also be prone to inertia effects, may they come from muscle internal masses or potential inertia of the measuring device [68].As a first approximation we can usually consider muscle massℳ to be scaled geometrically by either changing its cross-sectional area S or its length as represented by its optimal fibre length lM,opt. This implies that both mass density ρ and maximum stress σmax are about invariant parameters of muscles. Thus, in a real muscle its mass is not an independent parameter as in our model, and the rule to scale contraction times in proportion to ℳ·lM,opt/Fmax as expressed in Section 3.1 and in (A.11) and (A.12) can be condensed to (A.13). We should be aware of the fact that (A.11) and (A.12) do apply to any contracting muscle, may it accelerate just its internal mass (as analysed numerically in this paper) or an external mass that is usually much higher than its own mass. In contrast, (A.13) applies to the specific case of a muscle contracting just against its internal mass under the above-mentioned approximative assumptions of constant density and maximum stress. Thus, (A.13) tells us that the contraction time for a muscle exposed to a step in force, rather than analysed when accelerating an external load, scales in fact with the square of its fibre length. We would, therefore, expect that the effect of not reaching the theoretical maximum contraction velocity vmax (Figure 3(a)) would occur in any muscle preparation with fibre lengths that exceed a value somewhere between 1.5cm and 15cm, regardless of muscle mass. In the end, it seems that it is not the value of the muscle mass itself but the fibre length that scales the response times to force steps (even quadratically, see (A.13)). The effect occurs because contraction times for the whole muscle increase due to the force gradient decreasing with muscle extension, whereas the outermost muscle portions (released immediately in a force step) are always accelerated to their local vmax (and already decelerated consecutively due to local compression) much faster on an invariant time scale. For a 15cm fibre, a muscle is predicted to reach just 93% of vmax.We dare this preliminary prediction subject to future theoretical analyses explicitly including muscle internal elasticity. Serial elasticity might not really help to further speed whole muscle acceleration up. Parallel elasticity, however, might be effective because it adds accelerative forces to the force-velocity relation. Yet, interaction with damping like the force-velocity relation itself and with the mass distribution in series might reduce its accelerating effect down to a secondary level. Parallel elasticity is also a challenge to experimenters because it usually biases contraction measurements as soon as they are performed above or even at optimal fibre lengths. To make our theoretical analysis of the force-velocity relation contributing to accelerated contractions as comparable as possible with experiments, we chose initial conditions at which parallel elasticity is (almost) negligible, that is, at optimal lengths. As a result, our bigger muscle has already contracted by about0.24·lM,opt in the instant when it starts to decelerate its all over contraction (see first paragraph of Section 3.3). Because our model force-length relation has an unusually broad ascending branch [21] the isometric force for a CE shortened by considerable 24% has only decreased a little to about 95% of its absolute maximum. Therefore, it seems that reaching a contraction velocity of just 0.93·vmax is due to a mixture of local compressive forces and the force-length relation in our model simulations. Additional parallel elasticity would only have the potential to compensate for the latter relation. We would, therefore, expect that the phenomenon of not reaching theoretical vmax in big muscles should also occur when starting at initial lengths above the optimal fibre length.Finally, we would like to point to another interesting aspect. The propagation velocity of muscle excitation on its surface is about4m/s in vertebrate cross-striated muscles [69]. In the smaller muscle, this means that an excitation spike would need about an order of magnitude more time to spread along the complete length of a fibre than to reach vmax after a maximum force step down to zero: for 1.5cm fibres we calculated τv=-vmax=3.6·10-4s, and we find 0.015m/(4(m/s))=3.8·10-3s. As τv=-vmax scales quadratically with the fibre length (see Appendix A) we calculate, however, τv=-vmax=3.6·10-2s in the bigger muscle with a tenfold fibre length. In 15cm fibres the estimated propagation time of 0.15m/(4(m/s))=3.8·10-2s from one end to the other would, thus, already slightly exceed τv=-vmax. Accordingly, when thinking of the muscle response to a twitch as a response to an internal step in force, the propagation velocity, rather than inertia effects, may in fact limit the twitch time for a rise in force up to about 7cm fibre length (assumed the only end-plate is located in the middle of the fibre). Yet, even in shorter fibres it may be that inertia, rather than excitation propagation, could limit twitch response times because τv=-vmax increases linearly with muscle stress (see (A.13)), and twitch stress should be clearly lower than maximum stress during a fused tetanus. Additionally, any twitch time limitation by excitation propagation would be brought nearer to inertia limitation if a couple of end-plates were distributed along a fibre. All in all, excitation propagation contributing to the electromechanical delay restricts response times when a rapid rise to a higher force level, that is, to an increased activity is required [39]. In contrast, inertia effects, rather than excitation propagation, should be a time limiting factor in physiological loading situations when rapid changes in external forces occur at some given activity level at which Ca2+ is sufficiently available anyway along the whole fibre.


## 5. Conclusions

We expect that a significant progress in the time resolution of muscle experiments and their sensitivity with respect to measured mechanical variables and heat should allow to verify the predictions made by such models of the contractile muscular machinery that include the muscle mass distribution. The vast majority of current models has neglected muscle inertia and is, thus, not appropriate to understand high-frequency oscillations that have already been seen during some experiments [21, 56, 59, 60, 68, 70–74], let alone even potentially higher frequencies of muscle internal mass repositioning [24–26, 31]. The time scales, on which highly dynamic responses occur when muscle fibre contractions are disturbed by rapid (ideally infinitely steep) force steps, are very likely not to be explainable without muscle inertia. Such responses are eventually expected to depend on three basic mechanical properties: elasticity, viscosity or damping (leading to a release in heat), and inertia. In the model as presented here, we neglected explicit modelling of elasticity (and, therefore, oscillations), but we rather concentrated on the prediction of the time scale on which (exponentially damped) decays of disturbances would occur in case just the well-known hyperbolic Hill relation between force and velocity of active muscle fibres was present. Of course, these time scales would be expected to be also depending on fibre internal elasticity. This is even more the case, because the Hill relation itself has been predicted to possibly emerge from serial (visco-)elasticity that interacts with damping properties [65–67]. Therefore, to explain the superposition of (i) time constants for (exponentially damped) decays of disturbances, (ii) oscillations from the interaction of elasticity and muscle plus external inertia, and (iii) current heat production, thus, to understand muscle design by examining its mechanical and thermodynamical coupling with the environment, more focus on muscle inertia as part of muscle models seems to be indicated. In big muscles, even the asymptotic performance when approaching high velocities at low loads might depend on muscle inertia.



---

*Source: 848630-2012-11-22.xml*