Strain in shock-loaded skeletal muscle and the time scale of muscular wobbling mass dynamics

In terrestrial locomotion, muscles undergo damped oscillations in response to limb impacts with the ground. Muscles are also actuators that generate mechanical power to allow locomotion. The corresponding elementary contractile process is the work stroke of an actin-myosin cross-bridge, which may be forcibly detached by superposed oscillations. By experimentally emulating rat leg impacts, we found that full activity and non-fatigue must meet to possibly prevent forcible cross-bridge detachment. Because submaximal muscle force represents the ordinary locomotor condition, our results show that forcible, eccentric cross-bridge detachment is a common, physiological process even during isometric muscle contractions. We also calculated the stiffnesses of the whole muscle-tendon complex and the fibre material separately, as well as Young’s modulus of the latter: 1.8 MPa and 0.75 MPa for fresh, fully active and passive fibres, respectively. Our inferred Young’s modulus of the tendon-aponeurosis complex suggests that stiffness in series to the fibre material is determined by the elastic properties of the aponeurosis region, rather than the tendon material. Knowing these stiffnesses and the muscle mass, the complex’ eigenfrequency for responses to impacts can be quantified, as well as the size-dependency of this time scale of muscular wobbling mass dynamics.

Isometric force from physiological cross-sectional area (PCSA) The maximum isometric muscle tension σ max is between 2.25·10 5 Nm −2 53 and 3·10 5 Nm −2 54 . If PCSA P CE,0 and the pennation angle α are known, the maximum isometric muscle force F max can be estimated. We found P CE,0 and α values in two cases 55,56 . In case (a) 55 , P CE,0 and α were measured for the gastrocnemius medialis and lateralis muscles separately. The lateralis muscle accounted for 56% of the GAS PCSA (P CE,0 = 111 mm 2 ) and α = 14 • was found in both muscles. In case (b) 56 , measurements were made for only m. gastrocnemius medialis (54 mm 2 , α = 21 • ). Assuming a 44% contribution of this muscle, which is in line with Eng et al. 55 , the overall GAS PCSA in the work of Kadhiresan et al. 56 would then have been P CE,0 = 123 mm 2 . Due to the same boundary condition defining L M T C,90 • in case (i) 55 , we also assumed α = 14 • for our specimens and calculated our mean PCSA (P CE,0 = 112 mm 2 ) as P CE,0 = A CE,0,max cos α Based on the σ max range given above, maximum isometric force can be calculated as F max = σ max · A CE,0 = σ max · P CE,0 · cos α, which would yield values between 24.2 N and 32.3 N in case (i) and between 25.8 N and 34.5 N in case (ii). Compared to these estimations, we found F max = 30 N by extrapolation (Fig. 3).

The initial and dynamic strain contributions
Fibre material strain can be divided into two contributions: 'initial strain' due to the initial, isometric force and additional 'dynamic strain' due to the impact response. We measured the latter, and three examples are plotted in Fig. 4. For the initial strain, two limit cases of locating elasticity as a potential cause are conceivable: strain is either (i) solely located in the myosin heads or (ii) distributed across all sarcomere structures that act in series, that is, heads and actin plus myosin filament backbones. By further assuming in the case (i) that decreasing initial force was solely due to the reduced number of attached crossbridges, the initial strain of the fibre material would always the same, regardless of the force -that is, any single cross-bridge generates its maximum force (about 0.4% 30,35,48,49,50 , see bottom, horizontal, blue, dashed line in Fig. 5). In any other case, the most simple model relating initial fibre force to initial strain would be a linear stiffness (see bottom, sloped, red, solid line in Fig. 5 with zero strain at zero force and -0.4% at maximum isometric force), which is our case (ii). The overall maximum fibre strain during impact can then be calculated by the sum of the estimated initial strain and the measured dynamic strain. For one trial (no. 1 of 'd4s2'), initial and dynamic strain contributions are visualised in Fig. 5 by vertical arrows -a blue, dashed arrow (i) and a red, solid arrow (ii) for the two initial strain cases as well as a black, solid arrow on top of them for the measured dynamic strain. The sloped, thick, solid line represents the best linear fit of measured dynamic strains in all trials at initial forces higher than 12 N; the thinner linear line represents extrapolations to zero force and zero strain. This force boundary value at 12 N comes from a stiffness analysis of the MTC (Fig. 6).
Calculating dynamic muscle force change and the elements' stiffnesses Beyond describing CE-internal kinematics by strain , we used Newton's second law ∆F = M · (a COM + g) to calculate the dynamic force difference ∆F between MTC ends in response to the impact. The symbols are: muscle mass M, COM acceleration a COM , and gravitational acceleration g. Because CE mass is much higher than TAC mass, ∆F also well approximates the force difference between CE ends. On average, maximum values of ∆F are 0.36 ± 0.07 N, and they are practically independent of fatigue (see left ordinate in Fig. 7). Since the dynamic force difference (∆F (t)) as well as COM deflection relative to the frame ∆L M T C (t) = y COM (t) − y f rame (t) − (y COM (t = 0) − y f rame (t = 0)) and CE-internal elongation ∆L CE (t) are all known as functions of time t in each single trial, the MTC (K M T C ) and CE (K CE ) stiffnesses can be calculated based on all data points in the respective force-elongation relations ∆F (∆L i (t)) between TD and the instant of a COM ,max , with i = M T C (examples in Fig. S2) From this, the stiffness of the 'tendon-aponeurosis complex' (TAC) was inferred based on the model idea that the MTC consists of a mass that is suspended to the frame (or bone) by two compliant elements arranged in series (CE and TAC) given the overall MTC and the local CE stiffness values. We examined whether our calculated CE and MTC stiffness values may depend monotonically on isometric force F by using the simplest approach, that is, we fit the data with a linear function. The corresponding straight-line fits through all data points above 12 N are plotted in Fig. S3. We extracted this force boundary value at 12 N by separating the calculated K M T C (F ) and K CE (F ) data points into two ranges of approximately constant values (low forces) and assumed linear change with force F (high forces: straight lines), and searching for the boundary value between them with the best continuity of the inferred K T AC (F ) curve (see black, solid line). Except for the cases where either (i) the slopes of K M T C (F ) and K CE (F ) are zero or (ii) both stiffness values have their x-intercept at zero force, the K T AC (F ) curve will be a non-linear function. The inferred curve represents an almost constant value around K T AC (F ) ≈ 5,600 N m −1 down to ca. F ≈ 18 N, where the K CE (F ) curve starts to undergo the K T AC (F ) curve. Non-linear deviation of inferred K T AC (F ) from near-constancy increases with a further decrease in force down to the boundary F = 12 N. Yet, three points should be noted: (i) K T AC (F ) curvature may be a technical artefact of the fact that K CE (F ) and K M T C (F ) themselves do not exactly trend linearly with physiological reality; (ii) the independent inferences from below and above F = 12 N are almost the same at the boundary itself, namely, K T AC (F ) ≈ 6,700 N m −1 ; (iii) the uncertainty of the values K T AC (F ) = 5,600-6,700 N m −1 is the same as in other experiments on solely the m. gastrocnemius medialis head 42 , and the predicted absolute values are slightly lower than found in two other analyses based on Hill-type models 6, 42 . Here, K T AC is equivalent to K SEC = 9000±1,400 N m −1 (pers. comm. to T. Siebert) 42 and K SEC = 8,000 N m −1 6 there.

Young's moduli -a comparison to literature
Ai,0 is defined as a local material property that can be found from known macroscopic properties of a finite mass portion i, including stiffness K i , reference length L i,0 , and cross-sectional area A i,0 . With the chosen reference length L CE,0 = 0.0182 m (Table 1) in the centre of the muscle belly, A CE,0 = 1.09·10 −4 m 2 the maximum estimate of the anatomical crosssectional area (ACSA) within L CE,0 (Table 1), and K CE ≈ 9050 N m −1 the stiffness of fully active, fresh muscle fibres, we calculated the lower limit of their Young's modulus to be E CE ≈ 1.5 MPa. Due to low pennation angles α ≈ 20 • at optimal fibre length 55, 56, 47 , our ACSA and physiological cross-sectional area (PCSA) values found in literature 55,56 are practically identical. For fixed A i,0 and L i,0 , Young's modulus scales linearly with K i . Thus, in strongly fatigued and passive muscle fibres with K CE ≈ 3700 N m −1 , we found E CE ≈ 0.6 MPa as a minimum, which is ca. 40% of the stiffness and modulus values for fresh, fully active fibres. Due to using an upper limit for ACSA, these values of Young's modulus represent lower limits. Using A CE,0 = 0.73·10 −4 m 2 estimated at its upper (proximal) end, which is the minimum ACSA value within L CE,0 , we estimated an upper limit of Young's modulus of E CE ≈ 2.3 MPa for fresh and E CE ≈ 0.9 MPa for fatigued/passive muscle fibres. Assuming volume constancy, homogeneity, and uniformity for fibres and the belly's 5, 57, 58, 59, 52 , Young's modulus equals the shear modulus. With this, our results match others' data 60 well; this is demonstrated by the local shear values that were found at the Z-disc and M-line locations: 2.6 MPa and 1.3 MPa, respectively, in a fully active fibre as well as 1.7 MPa and 1.2 MPa, respectively, in a passive fibre.
TAC stiffness K T AC according to Eq. (7) is calculated (inferred) from K T AC contains all compliances in series to the fibre material, and is equivalent to what is often referred to as 'serial elastic element' (SEE) or 'serial elastic component' (SEC) in muscle modelling. The basic notion is described in three papers 61,62,63 . TAC stiffness is dominated by aponeurosis rather than tendon properties in rat GAS. The line of reasoning is concise. First, to check whether TAC stiffness was dominated by tendon properties, we calculated the tendon stiffness from anatomical data L tendon,0 = 0.012 m, A tendon,0 = 2.8·10 −6 m 2 , and a wellestablished literature value for Young's modulus of mammalian tendon material E tendon = 1.5 GPa 33, 34, 64, 65, 66 as K tendon = 350 kN m −1 . This K tendon value would be about 50 times higher than our upper limit of inferred TAC stiffness values K T AC = 6,700 N m −1 . Conversely, an upper limit estimate of TAC's Young's modulus can be calculated from assuming (a) the highest stiffness value K T AC = 6,700 N m −1 inferred from our analysis, (b) the smallest cross-sectional area value A T AC,0 = A tendon,0 = 2.8·10 −6 m 2 measured for rat GAS tendon, and (c) the length L T AC,0 = L M T C,0 − L CE,0 = 0.027 m that bridges the distance between the fibre material (analysed within L CE,0 ) and the frame. The corresponding theoretical, upper limit E T AC = 64 MPa is 23 times lower than mammalian tendons' Young's modulus. Alternatively, the lower limit E T AC = 2.1 MPa can be calculated from using (a) the lowest TAC stiffness value K T AC = 5,600 N m −1 occurring in our analysis and (b) the minimum fibre ACSA value A T AC,0 = A CE,0 = 0.73·10 −4 m 2 at the boundary between fibre material (CE) and TAC. As K CE results from the serial arrangement of all local elasticities within the cross-sectional areas along the finite length L CE,0 , the most likely value of inferred E CE should be approximately the arithmetic mean value of the upper and lower limits at each force, that is, E CE = 1.9 MPa for fresh and E CE = 0.75 MPa for fatigued/passive muscle. The local net material property E T AC in a belly's aponeurosis region may well increase continuously -when approaching the aponeurosis-tendon junction -from E T AC = E CE = 0.75 . . . 1.9 MPa at the fibreaponeurosis boundary to E T AC = E tendon = 1.5 GPa at the aponeurosis-tendon junction. Besides the fibre activity, this continuous change will depend on the geometrical arrangement of the aponeurosis and the fibres as well as on the shifting fibre-to-aponeurosis ratios of material contributions within the belly and on the belly surface. Average values of TAC's Young's modulus in the range E T AC = 2.1-64 MPa supports this idea.

A reflection on muscle-internal mechanics during impact
In the subsection 'Calculating dynamic muscle force change and the elements' stiffnesses', we explained how we used Newton's second law to calculate -from the kinematics of the centre of the distributed muscle masses (COM) estimated by marker tracking -the dynamic force difference ∆F between muscle ends (origin and insertion) that must have acted to accelerate the COM. As additional information, we measured the initial, isometric force F with a force sensor in series to the muscle (plotted in Fig. 3 as a function of time after muscle extraction). If the muscle masses are suspended to a rigid construction (i.e., the frame), one can calculate a finite distance between the COM and the suspension positions that approximates an anatomical length, namely, the muscle-tendon complex (MTC) length. Our mechanical analysis of a real muscle implies at least two further prerequisites: (i) we wanted to formulate a mechanical model of the MTC to potentially subdivide it into further sub-elements, and (ii) the MTC is suspended to a single rigid construction at both ends. Both prerequisites are implied in our calculation of the MTC length, and (ii) means that it does not matter whether MTC length changes are calculated with respect to the upper or lower suspension point.
Our initial idea was that the minimum mechanical model description starts with just one net suspending structure for the COM: a spring acting between COM and one suspension point; note that any length change calculated this way represents the length change of all MTC material. Consistent with this, the calculated force difference ∆F , which is proportional to the COM acceleration, is interpreted as the net force change that acts on the MTC spring while it changes its length by ∆L M T C . Likewise, a net MTC stiffness can be calculated. We further know that (iii) few percent of the MTC masses are located in the tendons, which justifies estimating COM acceleration from kinematic information of fibre material solely. We then can use this to make a first MTC length subdivision within the model, namely, into the length of mass distribution itself (termed CE) and the remaining material in series. Thus, the COM of the (masses in the) CE are accelerated in space as calculated above from ∆F , and we determine an additional length change ∆L CE of this CE part superposed to the COM movement. This CE length change between the upper and lower ends of the CE region subtracted from all MTC length change is the length change ∆L T AC of all MTC parts in series to the CE region, which we indicate by TAC.
Using these arguments, we can now explain in more detail our model idea of two deforming (length changing) structures in series that suspend the muscle mass to the skeleton or frame, respectively: while the CE can locally change length in any instant due to finite fibre stiffness, which can easily induce, e.g., double fibre stress (and thus mean force F in a cross-sectional area) by doubling strain as compared to initial, isometric strain, the force difference ∆F between both ends of the MTC (and thus CE) provides its acceleration in space. Note that we measured reliable, static (isometric) force values only before touchdown. ∆F should also be reflected by the sum of both length changes in the upper and lower anatomical structures in series to CE, that is, the upper and lower TAC parts, respectively.
In a nutshell, F and ∆F are, for the time being, two independent mechanical variables that are probably locally and continuously distributed during wave propagation. The variable F represents the concept of 'joint force' that causes linear momentum transmission between two adjacent mass distributions (compare, e.g., the study 51 [appendix 2]) that can be calculated according to Newton's third law (action=reaction) in any conceived area of cutting through the muscle. The variable ∆F represents the difference in external forces that accelerates one selected mass distribution localised between its spatial boundaries ('joints').
In a real muscle, the stress or force, respectively, and their gradients are expected to vary temporally and spatially during wave propagation. The local force F is the integral of stress over a finite cross-sectional area, and the force difference ∆F is an integral of a force gradient distribution over a finite distance along the gradient. Adding the idea of elastic cross-bridges, changes in strain and force F are proportionally connected, which applies to both temporal and spatial variations -that is, rates and gradients, respectively -provided there is no forcible detachment. For example, a local gradient in strain generates a local force gradient and thus finite force differences ∆F for finite distances along the gradients. As we could only spatially resolve the net (mean) strain in the whole fibre material region of ca. 18 mm using our present pioneering methodological approach, our strain signals (plotted in Figs. 4 and 5) represent the mean strain values along this CE region. With higher spatial resolution, local strain and thus stress or force F values, respectively, are probably distributed around their mean values. . Vertical, grey lines indicate time instants of TD (t = 0 s) and maximum dynamic strain amplitude CE ,max during shock response t > 0 s. The bottom trial 'd5s201' was excluded from our analysis because its earliest extreme in strain (at ca. 9 ms) was a minimum rather than a maximum as in the top trial 'd4s100' (at ca. 4.5 ms).
are also plotted using three parameters k, b, d to fit a linear function F (L i ,L i ) to measured data (here: L i = L M T C and F = ∆F ). Accordingly, parameter k is a direct estimation of the MTC stiffness K M T C (slope of dashed line) that results from the serial arrangement of muscle fibre (CE: K CE ) and tendon/aponeurosis (TAC: K T AC ) material. For comparison, the effect on estimating K M T C when neglecting the damping parameter d is shown (slope of '2 parameter fit').  (Fig. S2) to their measured force-length data and identifying slope parameter k as the respective element's stiffness. Through each of the two data clouds, a linear fit K = a · F + b to all data points was calculated for F > 12 N. Additionally, the following are also plotted: (i) circles represent the median value for trials with active muscles in the region F < 12 N, and (ii) a cross represent the median value for all trials with passive muscles. To illustrate the uncertainty of the linear approximations of the stiffness dependencies on force (fatigue), the 'undisturbed linear fits' (light blue, thick, dashed line in A and red, thick, dashed line in B), we recalculated the linear fits of Eq. (12) using two additional worst case (wc) data points that were hypothetically placed at the maximum observed force (plus and minus observed, maximum deviation from each undisturbed linear fit in A and B: blue and green squares, respectively). The corresponding linear fits, each including the respective wc data point, are plotted as dashed blue and green lines, respectively, in A and B. Additionally, the 95% confidence intervals for the undisturbed linear fits in A and B are plotted with red, thin, solid lines. Note that, although MTC stiffness does only weakly depend on force, the slope of the linear fit remains positive even in the wc of the hypothetically low stiffness value at maximum force.