Accuracy of experimentally estimated muscle properties: Evaluation and improvement using a newly developed toolbox
Abstract
The mechanical behaviour of a muscle-tendon complex depends on properties such as the force-length relationships, the force-velocity relationship, and the excitation dynamics. Quick-release and step-ramp experiments are commonly used to estimate these properties. The accuracy of these methods is unclear, as the actual values of these properties are unknown in experiments on real muscle. We conducted a modelling study using a Hill-type muscle-tendon complex model with literature-derived parameter values and simulated quick-release, step-ramp, and isometric experiments. From the simulated experiments, we assessed how accurately the model’s parameter values could be retrieved. Using a method traditionally used in literature, the series elastic element stiffness was underestimated by ~35%, due to the incorrect assumption that muscle fibres do not shorten during quick releases. Consequently, this yielded an overestimation of the excitation dynamics activation time constants of ~20%. We developed an improved method that accounted for muscle fibre length shortening during quick releases. Using our improved method, all parameter values closely matched their actual values. A sensitivity analysis showed that the most critical parameters were robust to perturbations in experimental data. Lastly, we compared Hill-type MTC model predictions against in situ data from three rat m. gastrocnemius medialis muscles. Predictions based on parameters from the improved method showed closer agreement than those based on the traditional method — both for quick-release, step-ramp, and isometric experiments, as well as for independent stretch-shortening cycles. In conclusion, the improved method enables more accurate estimates of muscle-tendon complex properties, addressing limitations of the traditionally used method.
1 Introduction
In the early 20th century, A.V. Hill laid the foundation of muscle mechanics by recognising that muscles contain not only muscle fibres (the contractile element; CE) but also connective tissue, both parallel to the muscle fibre (the parallel elastic element; PEE) and in series with the muscle fibre (the tendon or serial elastic element; SEE). Together these three elements form the muscle-tendon-complex (MTC). The contraction dynamics describe how CE force emerges from the interplay between the force-length properties of CE, PEE and SEE, and the force-velocity properties of CE. In addition, CE force depends on the active state (the relative amount of \(Ca^{2+}\) bound to troponin C, Ebashi and Endo 1968). Active state, in turn, depends on activation via the excitation dynamics. Consequently, the mechanical behaviour of MTC arises from the intricate interplay between the contraction and the excitation dynamics.
In his famous (1938) paper A.V. Hill performed so-called ‘isotonic quick-release experiments’ using the now-renowned Levin-Wyman lever (Levin and Wyman 1927) to estimate the CE force-velocity relationship. In these experiments, MTC was attached to a lever and maximally stimulated. When SEE force plateaued (i.e., MTC was isometrically delivering force), a magnet was released from the lever causing a sudden drop in the force acting on the MTC. As SEE force remained constant right after this drop in SEE force, SEE length remained constant and therefore all MTC length changes were attributed to CE length changes. The CE velocity corresponding to the constant SEE force after the drop was then calculated as the maximum rate of change of MTC length over time. Consequently, each isotonic quick-release experiment contributed a single data point of the CE force-velocity relationship. Nowadays, servomotors are typically used in experiments on isolated MTCs to control MTC length (changes). In contrast to the force-controlled experiments with the Levin-Wyman lever, servomotor control the MTC length over time (see Figure S1). In order to estimate the CE force-velocity relationship using a servomotor, Cecchi, Colomo, and Lombardi (1978) introduced ‘step-ramp experiments’. In these experiments, MTC is kept isometrically under maximal stimulation until SEE force plateaus. Then, a rapid shortening of MTC length is imposed (the ‘step’), immediately followed by a constant-velocity shortening of the MTC (the ‘ramp’). This experimental protocol results in a brief plateau in SEE force just after the ramp. During this brief period, SEE length remains constant and therefore CE velocity equals MTC velocity. Thus, step-ramp experiments with servomotors serve as a valid alternative for isotonic quick-release experiments to estimate the CE force-velocity relationship.
In addition to estimating the CE force-velocity relationship using isotonic quick-release experiments, the research group of A.V. Hill also used these experiments to approximate SEE stiffness. The underlying idea was that the drop in MTC force was so fast that, due to contraction dynamics, CE had no time to shorten and all MTC shortening could thus be directly attributed to SEE shortening. Consequently, this experiment allowed to directly relate SEE length changes to changes in SEE force and therefore to estimate SEE stiffness. As mentioned above, in isotonic quick-release experiments, a sudden drop in force acting on the MTC causes a rapid change in MTC length. While in these lever-based experiments MTC force is controlled, in experiments using a servomotor MTC length is controlled. To estimate SEE stiffness with a servomotor setup, the opposite approach is taken: a rapid change in MTC length is imposed by the motor, resulting in a rapid change in SEE force. Somewhat confusingly, these servomotor-based experiments are also referred to as ‘quick-release experiments’, even though their approach differs fundamentally from that of an isotonic quick-release experiment. In the remainder of this paper, the term quick-release experiments will exclusively refer to servomotor-based length-controlled quick-release experiments. In quick-release experiments using a servomotor, the drop in MTC length takes several milliseconds (e.g., approximately 10 ms with the commonly used Aurora Scientific’s 300C series motor). As CE shortens within this short time period, the assumption that MTC length changes are exclusively due to SEE length changes may not be entirely valid. This raises the question to what extent SEE stiffness can be accurately estimated based on quick-release experiments using servomotors.
A.V. Hill already acknowledged in (1950) that CE shortens during the brief period in which SEE force decreases – even for his isotonic quick-release experiments using the Levin-Wyman lever – and introduced a method to correct for this CE shortening. In this method, the CE shortening that occurred due to the quick-release was estimated based on the CE force-velocity relationship of the MTC investigated. Thus, additional experimental work – such as performing step-ramp experiments and estimating the CE force-velocity relationship – are required to correct for CE shortening during the quick-release. Zandwijk et al. (1997) and Lemaire et al. (2016) developed similar methods to correct for CE shortening. As expected, Zandwijk et al. (1997) showed that this method results in higher estimates of SEE stiffness than with the method without correcting for CE shortening. Nevertheless, the accuracy of SEE stiffness estimates remains uncertain, as the actual SEE stiffness is unknown in experiments, making it infeasible to compare the estimated value against the actual one. Consequently, it remains unclear whether the additional complexity and experimental work required for the improved method are truly justified in comparison with the simpler traditional method, especially when the only interest is SEE stiffness.
In experimental studies in which MTC properties are estimated, SEE stiffness is typically estimated first in order to discriminate between SEE stiffness and the other properties, such as those of the CE force-length relationship. In experiments, CE length is unknown and therefore the properties of the CE force-length relationship cannot be directly estimated based on experimental data. Typically, the CE force-length properties are estimated based on the measured MTC force-length relationship (e.g., Blümel et al. 2012; Lemaire et al. 2016). However, different combinations of CE and SEE properties can yield almost identical MTC force-length relationships. Consequently, inaccuracies in SEE stiffness estimates can propagate to errors in the estimation of the CE force-length properties. While predictions under isometric conditions may remain reasonably accurate, these errors can lead to substantially different mechanical behaviour during dynamic contractions such as stretch-shortening cycles (see Figure 1). As such, incorrect estimation of SEE stiffness might (partially) explain difference between predictions derived by a Hill-type MTC model and experimental results. In sum, the question arises to which extent the estimates of other MTC properties are affected by inaccurate SEE stiffness estimates and whether more accurate SEE stiffness estimates improves predictions derived by a Hill-type MTC model.
The potential effect of error propagation also extends to the estimation of the excitation dynamics properties. In the absence of histochemical data, excitation dynamics properties are commonly estimated by fitting predicted SEE force of the Hill-type MTC model to experimental data (e.g., Blümel et al. 2012; Lemaire et al. 2016). As SEE force over time obviously depend on the contraction dynamics properties (e.g., the force-length relationships and the force-velocity relationship), the estimation of the excitation dynamics properties is not independent of the estimation of the contraction dynamics properties. Consequently, inaccuracies in SEE stiffness estimates can propagate to inaccuracies in estimated excitation dynamics properties. This interdependence raises an additional question of what experimental data should be used in order to minimise this potential error propagation.
The primary objective of this study was to systematically evaluate the accuracy of estimating contraction and excitation dynamics properties on the basis of data collected in commonly used experiments. Specifically, we aimed to compare the results of a method that does not correct for CE shortening in quick-release experiments (the ‘traditional method’) with a method that does correct for CE shortening (the ‘improved method’). As explained earlier, the limitation of relying exclusively on an experimental approach is that the actual properties are unknown, which renders a reliable assessment of the accuracy of the estimated properties not feasible. To circumvent this problem, we conducted a modelling study using a Hill-type MTC model. First, we obtained three different sets of parameter values of a Hill-type MTC model from existing literature. Using these parameter values, we simulated data of quick-release, step-ramp and isometric protocols mimicking experiments on isolated MTCs using servomotors. We then investigated how accurately we could retrieve the model’s parameter values (with respect to their ‘actual’ value). In addition, we employed a comprehensive sensitivity analysis to assess the sensitivity of parameter values against perturbations in experimental data and to examine the interdependency of the estimated properties. The secondary objective of this study was to evaluate whether predictions from a Hill-type MTC model using parameter values estimated with the improved method better matched experimental data than predictions using parameter values estimated with the traditional method. To this end, we used in situ data from experiments on three rat m. gastrocnemius medialis. We estimated the contraction and excitation dynamics properties using both the traditional and the improved method, and compared the resulting model predictions — based on each parameter set — to experimental force data. Overall, this study provides insights into the accuracy of existing commonly used methods for contraction and excitation dynamics property estimation and their influence on predictions of mechanical behaviour derived by a Hill-type MTC model. The approach that we designed in this study for estimating contraction and excitation dynamics properties based on data of quick-release, step-ramp and isometric experiments using servomotors, is made available as an open-source toolbox (https://github.com/edwinreuvers/mp-estimator).
2 Methods
The primary objective of this study was to assess the accuracy of estimating contraction and excitation dynamics properties. For this purpose, we employed a Hill-type MTC model to simulate MTC mechanical behaviour, as detailed in Section 2.1. Using parameter values obtained from existing literature, we simulated data of quick-release, step-ramp, and isometric experiments, mimicking experiments on isolated MTCs with servomotors, as described in Section 2.2.1. We then estimated the contraction and excitation dynamics parameter values using two methods: one with (the ‘traditional method’) and one without (the ‘improved method’) a correction for CE shortening due to the quick-release (see Section 2.3). Additionally, we conducted comprehensive sensitivity analyses to assess the accuracy of the improved method, as outlined in Section 2.4. The secondary objective was to evaluate predictions derived by a Hill-type MTC model using two sets of parameter values: one estimated by the traditional method, the other by the improved method. For this purpose, we used in situ data (see Section 2.2.2) from three rat m. gastrocnemius medialis. Predictions were assessed against quick-release, step-ramp, and isometric experiments, as well as independently measured stretch-shortening cycles. This allowed us to investigate whether the improved method enhanced the predictions derived by a Hill-type MTC model.
2.1 Muscle model
2.1.1 Contraction dynamics
We employed a Hill-type MTC model consisting of a contractile element (CE) and a parallel elastic element (PEE), which were both in series with a serial elastic element (SEE), as depicted in Figure 2. CE represented the contractile element of the muscle fibres, while PEE and SEE represented the tissues arranged in parallel and in series with the muscle fibres, respectively. Given the negligible small mass of MTC compared to the forces typically delivered by CE, PEE, and SEE, we simplified the model by neglecting the second-order dynamics of the system:
\[ \begin{aligned} L_{MTC} &= L_{CE} + L_{SEE} + \\ F_{SEE} &= F_{CE} + F_{PEE} \end{aligned} \tag{1}\]
\(L_{MTC}\), \(L_{CE}\), \(L_{SEE}\) and \(L_{PEE}\) denote MTC, CE, PEE and SEE length respectively, while \(F_{SEE}\), \(F_{CE}\), and \(F_{PEE}\) denote the CE, PEE and SEE force. PEE and SEE were assumed be purely elastic and were modelled as quadratic springs (Zajac 1989,):
\[ F_{EE}= \begin{cases} k_{SEE} \cdot (l_{EE}-L_{EE}^{0})^2, & \text{if}\ l_{EE} \ge L_{EE}^{0} \\ 0, & \text{otherwise} \end{cases} \tag{2}\]
\(L_{EE}\) denotes the length of either PEE or SEE and \(L_{EE}^{0}\) the slack length of either PEE or SEE. \(k_{EE}\) denotes a parameter that scales the stiffness of either PEE or SEE.
Isometric CE force depends on CE length (Blix 1892, 1893; A. V. Hill 1925). We simplified the classic sarcomere force-length relationship (Gordon, Huxley, and Julian 1966; Walker and Schrodt 1974,) to a second order polynomial, which describes the force-length relationship reasonably well (Bobbert, Ettema, and Huijing 1990; Woittiez et al. 1984):
\[ F_{CE}^{isom,rel} = \begin{cases} \frac{-1}{w^{2}} \cdot (L_{CE}^{rel}-1)^2+1, & \text{if} \ F_{CE}^{isom,rel} > 0.1 \\ s_{exp} \cdot e^{k_{exp} \cdot L_{CE}^{rel}}, & \text{otherwise} \end{cases} \tag{3}\] \(F_{CE}^{isom,rel}\) denotes the CE isometric force normalised by \(F_{CE}^{max}\), \(w\) determines the width of the CE force-length relationship (for \(F_{CE}^{isom,rel} > 0.1\)) and \(L_{CE}^{rel}\) denotes CE length normalised by CE optimum length (\(L_{CE}^{opt}\), the CE length at \(F_{CE}^{isom,rel}=1\)). Exponential tails were added to the isometric CE force-length relationship such that it had a continuous first derivative with respect to \(L_{CE}^{rel}\) (see Figure 2). The parameter values of \(s_{exp}\) and \(k_{exp}\) were determined separately for the ascending and descending limb of the CE force-length relationship.
The concentric CE force-velocity was modelled according to Archibald Vivian Hill (1938). The original description was adjusted as suggested by Soest and Bobbert (1993) to account for situations where the active state (\(q\), the relative amount of \(Ca^{2+}\) bound to troponin C, Ebashi and Endo 1968) is not maximal:
\[ F_{CE} = \frac{q \cdot (F_{CE}^{isom,rel} \cdot F_{CE}^{max} \cdot b + a \cdot v_{CE})}{b-v_{CE}} \tag{4}\]
\(F_{CE}\) denotes the CE force and \(v_{CE}\) denotes the CE velocity. \(a\) is a parameter that defines the curvature of the concentric CE force-velocity relationship and defines together with \(b\) and \(F_{CE}^{max}\) the maximum shortening velocity (\(v_{CE}^{max}\), \(v_{CE}^{max} = -b/a \cdot F_{CE}^{max}\)). The maximum shortening velocity is scaled by \(F_{CE}^{isom,rel}\) due to the formulation of the CE force-velocity relationship (see Equation 4). However, the maximum shortening velocity has been reported to be more or less constant above optimum CE length (Gordon, Huxley, and Julian 1966; Stern 1974). Based on this observation, the parameter \(a\) was scaled by \(F_{CE}^{isom,rel}\) above optimum CE length (\(L_{CE}^{rel} > 1\)) to make the maximum shortening velocity independent of CE length above optimum CE length. The parameter \(b\) was scaled with \(b_{scale}\) for low levels of active state in order to make the maximal contraction velocity dependent on the active state (Petrofsky and Phillips 1981). The original description of this scaling by Soest and Bobbert (1993) was somewhat reformulated, such that \(b_{scale}\) was continuously differentiable with respect to \(q\):
\[ \begin{aligned} b_{scale} &= \frac{1}{1 + e^{-b_{shape} \cdot (q-q_{sc})}} \\[0.5em] \text{with} \quad q_{sc} &= \frac{log(\frac{1}{b_{scale}^{min}-1}+q_0 \cdot b_{shape})}{b_{shape}} \end{aligned} \tag{5}\]
\(q_0\) denotes the minimum value of \(q\), \(b_{scale}^{min}\) denotes the minimal scale factor of \(b\) (i.e., the minimal value of \(b_{scale}\)) and \(b_{shape}\) defines the steepness of the relation between \(b_{scale}\) and \(q\). The value of \(b_{shape}\) was obtained by minimising the least square error between the formulation of Soest and Bobbert (1993) and the reformulated version of the relation between \(b_{scale}\) and \(q\).
The eccentric part of the CE force-velocity relationship was modelled as a slanted hyperbola function (Van Soest, Gföhler, and Casius 2005) of the following form: \[ F_{CE}^{rel} = \frac{c_4 - (c_3 + v_{CE}^{rel} \cdot (c_1 + c_2 \cdot v_{CE}^{rel})}{c_3 \cdot v_{CE}^{rel}} \tag{6}\]
\(c_1\), \(c_2\), \(c_3\) and \(c_4\) denote parameters that define the shape of the slanted hyperbola and depend on \(F_{CE}^{isom,rel}\) and \(q\). These parameters were chosen such that (1) the concentric and eccentric curves were continuous in the point \(v_{CE}^{rel}=0\); (2) the ratio between the eccentric and concentric derivatives of \(\frac{dF_{CE}^{rel}}{dv_{CE}^{rel}}\) in the point \(v_{CE}^{rel}=0\) was defined by \(r_{slope}\); (3) the oblique asymptote of the eccentric curve was defined by \(F_{asymp}\) and (4) the slope of the slanted asymptote was defined by \(r_{as}\).
2.1.2 Excitation dynamics
The excitation dynamics were modelled in two steps. The first step, to which we refer as the calcium dynamics, related the rate of change in normalised free \(Ca^{2+}\) concentration between the myofilaments (\(\gamma\)) to normalised \(CE\) stimulation (\(STIM\)) and the normalised free \(Ca^{2+}\) concentration between the myofilaments itself (Hatze 1981, pp 31-42):
\[ \dot{\gamma}= \begin{cases} \frac{STIM \cdot (1 - \gamma_0) - \gamma + \gamma_0}{\tau_{act}}, & \text{if} \ STIM \ge \gamma \\ \frac{STIM \cdot (1 - \gamma_0) - \gamma + \gamma_0}{\tau_{deact}}, & \text{otherwise} \end{cases} \tag{7}\]
\(\gamma_0\) denotes the minimum value of \(\gamma\) and had a arbitrary small value such that Equation 8 was always solvable. \(\tau_{act}\) and \(\tau_{deact}\) are both time constants of the first-order differential equations describing the calcium dynamics. The second step involved the relation between active state (\(q\)) and \(\gamma\), to which we refer as the \(q-[Ca^{2+}]\) relation. It is well-known that muscle fibres become increasingly sensitive to \([Ca^{2+}]\) as their length increases (Kistemaker, Van Soest, and Bobbert 2005; Rack and Westbury 1969; Stephenson and Williams 1982). Consequently, \(q\) was also made dependent on \(L_{CE}^{rel}\). The \(q-[Ca^{2+}]\) relation was modelled according to Hatze (1981, pp 31-42), but was mathematically reformulated to ensure that the parameter values were physiologically meaningful and to facilitate its application for Optimal Control (e.g., Kistemaker et al. 2023):
\[ \begin{aligned} q &= q_0 + \frac{1-q_0}{1+(kCa \cdot \gamma)^{A_{act}} \cdot e^{a_{act} \cdot B_{act}}} \\[0.5em] \text{with} \quad A_{act} &= \log_{10}(e^{-a_{act}}) \\ \text{and} \quad B_{act} &= b_{act,1} + b_{act,2} \cdot L_{CE}^{rel} + b_{act,3} \cdot {L_{CE}^{rel}}^2 \end{aligned} \tag{8}\]
\(kCa\) relates \(\gamma\) to the actual \(Ca^{2+}\) concentration between the myofilaments and \(B_{act}\) denotes the \(pCa^{2+}\) level at which \(q=0.5\). \(B_{act}\) depended on \(L_{CE}^{rel}\) (Equation 8), while parameter \(a_{act}\) determined the steepness of the relation between \(\gamma\) and \(q\).
2.2 Simulated and in situ data
2.2.1 Simulated data
To simulate quick-release, step-ramp, and isometric experiments, we used three parameter sets of a Hill-type MTC model from existing literature for three rat m. gastrocnemius medialis (GM1, GM2, and GM3; see Table S1).
Quick-release experiments Each quick-release experiment consisted of an isometric phase until SEE force plateaued, followed by a rapid (step) change (10 ms) in MTC length and then followed by another isometric phase (see Figure 3A). CE stimulation was maximal (i.e., \(STIM = 1\)) during the first isometric phase and continued at this maximal level until it was switched off shortly after the step change in MTC length occurred. The step change in MTC length was 0.2 mm and was chosen such that the resulting change in SEE force was about 5-10% of \(F_{CE}^{max}\). To prevent irreversible damage to muscle fibres, the maximal MTC length used in experiments is typically not far above the MTC length that yields maximal isometric SEE force (\(L_{MTC}^{opt}\)). Accordingly, we simulated quick-release experiments at various initial MTC lengths, ranging from a very short MTC length (i.e., a length yielding very low isometric SEE force) and at every 1 mm increment up to a MTC length that was maximal 3 mm above \(L_{MTC}^{opt}\).
Step-ramp experiments Each step-ramp experiment consisted of an isometric phase slightly above (\(\pm\) 0.5 mm) \(L_{MTC}^{opt}\), followed by a rapid (step) change (10 ms) in MTC length and then a constant MTC velocity ramp (see Figure 3B). CE stimulation was maximal (i.e., \(STIM = 1\)) during the first isometric phase and continued at this maximal level until it was switched off 0.1 s after the step change in MTC length occurred. The change in MTC length of this step and the constant MTC velocity of the ramp were chosen such that this resulted in a more or less a constant SEE force (and therefore constant CE force) at the beginning of the ramp. We simulated 9 step-ramp experiments with different combinations of step sizes and constant MTC velocity ramps, to cover a substantial part of the concentric CE force-velocity relationship.
Isometric experiments Isometric experiments were simulated for a combination of different MTC lengths (about 0, -2, -4 and -6 mm below \(L_{MTC}^{opt}\)) and stimulation durations (35, 65 and 95 ms). CE stimulation was maximal (i.e., \(STIM = 1\)) during the indicated stimulation duration and fully off elsewhere (see Figure 3C). The combination of four different MTC lengths and three different stimulation durations resulted in twelve different isometric experiments.
2.2.2 In situ data
We also used data from an in situ experiment on isolated rat m. gastrocnemius medialis, conducted on three male Wistar rats. The experiment included quick-release, step-ramp and isometric experiments, which were similar to those described in Section 2.2.1. For each rat, we performed between 11 and 13 quick-release experiments, 10 and 12 step-ramp experiments and 12 isometric experiments. Below, we provide a brief description of the experimental procedure, which is fully detailed in (Reuvers et al. 2025).
In the experiment, approved by the Committee on the Ethics of Animal Experimentation at the Vrije Universiteit (Permit Number: FBW- AVD11200202114471), rats were first anesthetised with urethane. The hindlimb was then shaved, and the overlying skin and m. biceps femoris were removed. The medial and lateral part of m. gastrocnemius were carefully separated from their surrounding tissue and exposed as much as possible. The rats were placed in the experimental setup with the hindlimb, femur and foot fully fixed. The distal end of the calcaneal tendon was attached to a servomotor (Aurora 309C, Aurora Scientific, Aurora, Canada) using Kevlar thread and aligned to ensure that m. gastrocnemius medialis pulled in its natural direction. All nerves not innervating m. gastrocnemius medialis were severed. M. gastrocnemius medialis was stimulated via a cuff-electrode placed on the sciatic nerve, with the proximal nerves crushed to prevent spinal reflexes. This experimental setup allowed for precise control of m. gastrocnemius medialis length changes and stimulation as well as accurate measurement of m. gastrocnemius medialis force.
2.3 Parameter value estimation procedure
The general procedure to estimate the parameter values involved minimising the sum of squared differences between the data and the values based on the estimated parameter values. To distinguish between these, we used subscripts: for example, the value of the SEE force before the quick-release of the data is indicated by \(F_{SEE}^{QR-,data}\) while the value based on the estimated parameter values is indicated by \(F_{SEE}^{QR-,est}\).
The data of the quick-release experiments were used to estimate the parameter values of the CE, PEE and SEE force-length relationships. The data of the step-ramp experiments were used to estimate the parameter values of the CE force-velocity relationship. The data of the isometric experiments were used to estimate the parameter values of the excitation dynamics. Below we provide an overview on the methods used to estimate the parameter values. We also offer an open-source toolbox that automates the estimation of contraction and excitation dynamics parameter values.
2.3.1 CE, SEE and PEE force-length parameter value estimation
As explained in the introduction, different combinations of contraction dynamics parameters can yield almost identical mechanical behaviour under isometric conditions (see also Figure 1). For this reason, it is important to estimate SEE stiffness first, in order to discriminate between SEE stiffness on the one hand, and \(L_{CE}^{opt}\) and \(L_{SEE}^0\) on the other. In our approach, the parameter that scales SEE stiffness (\(k_{SEE}\)) is therefore estimated first. This is followed by the estimation of the PEE parameters, as \(k_{SEE}\) is also required for their estimation. Finally, the parameters \(F_{CE}^{max}\), \(L_{CE}^{opt}\), and \(L_{SEE}^0\) are estimated.
Estimation of SEE stiffness. The first step was to estimate SEE stiffness. The model SEE force depends on the parameter values of \(k_{SEE}\) and \(L_{SEE}^0\), and obviously on SEE length. The problem, however, is that SEE length is unknown in experiments. This makes it challenging to estimate the parameter values concerning the SEE force-length relationship. Fortunately, quick-release experiments provide a way out. During quick-release experiments with a servomotor, the motor quickly shortens MTC length such that there is a rapid decline in SEE force (the ‘quick-release’). Due to the shortness of this timeframe, CE shortening is minimal such that almost all MTC shortening is taken up by SEE. Often, in experimental studies, it is assumed that all MTC shortening can be attributed to SEE shortening. In reality, this assumption leads to an overestimation of SEE shortening, as CE also shortens during this short timeframe. We introduced a method to correct for this overestimation of SEE shortening (see Section 2.3.4) and estimated the parameter values considering both methods: without and with correcting for CE shortening during the quick-release.
Obtaining the SEE length change due to the quick-release (\(\Delta L_{SEE}^{QR}\)) is a crucial step. This is because SEE length after the quick-release (\(L_{SEE}^{QR+}\)) can then be expressed as SEE length before the quick-release (\(L_{SEE}^{QR-}\)) plus the change in SEE length due to the quick-release: \(L_{SEE}^{QR+} = L_{SEE}^{QR-} + \Delta L_{SEE}^{QR}\). This allowed us to rewrite Equation 2 to estimate SEE force immediately before (\(F_{SEE}^{QR-,est}\)) and after (\(F_{SEE}^{QR+,est}\)) the quick-release. This yielded:
\[ \begin{split} F_{SEE}^{QR-,est} &= k_{SEE} \cdot (L_{SEE}^{QR-} - L_{SEE}^0)^2 \\ F_{SEE}^{QR+,est} &= k_{SEE} \cdot (L_{SEE}^{QR-} + \Delta L_{SEE}^{QR} - L_{SEE}^0)^2 \end{split} \tag{9}\]
Now, there are two unknown parameters (i.e., \(k_{SEE}\) and \(L_{SEE}^0\)), while also SEE length before the quick-release is also unknown for each quick-release experiment. Consequently, there are more unknowns than equations. To address this issue, we replaced \(L_{SEE}^{QR-} - L_{SEE}^0\) with a temporary parameter \(c_{SEE}\). This yielded:
\[ \begin{split} F_{SEE}^{QR-,est} &= k_{SEE} \cdot (c_{SEE})^2 \\ F_{SEE}^{QR+,est} &= k_{SEE} \cdot (c_{SEE} + \Delta L_{SEE}^{QR})^2 \end{split} \tag{10}\]
As \(k_{SEE}\) scales SEE stiffness, its value is constant across all quick-release experiments. In contrast, \(c_{SEE}\) depends on the parameter \(L_{SEE}^0\) and SEE length before the quick-release, which differs among each quick-release experiment. Consequently, \(c_{SEE}\) should be determined for each quick-release experiment individually. Visually, \(c_{SEE}\) determines the shift along the x-axis on the SEE force-length relationship, while \(k_{SEE}\) scales SEE stiffness and consequently SEE force (see Figure 4A). The value of \(k_{SEE}\) (and \(c_{SEE}\) for each quick-release experiment) were found by minimising the sum of squared differences between the data and estimated SEE force before (\(F_{SEE}^{QR-}\)) and after (\(F_{SEE}^{QR+}\)) the quick-release.
Estimation of PEE parameter values. The second step was to estimate PEE stiffness. The model PEE force depends on the parameter values of \(k_{PEE}\) and \(L_{PEE}^0\), and obviously on PEE length. In experiments, SEE force is measured, while PEE force is required to estimate PEE stiffness. As such, experimental data should be used in which CE force is negligible such that PEE force is approximately equal to SEE force. At the beginning of the quick-release experiments, the active state is so low that CE force is negligible. We selected an interval of 10 ms in which the SEE force was minimal (\(F_{SEE}^{min}\), orange dot in Figure 3A).
In experiments, PEE length is unknown. This problem can be addressed as follows. First, PEE length equals MTC length minus SEE length: \(L_{PEE} = L_{MTC} - L_{SEE}\) (Equation 1). Second, SEE length is the sum of SEE slack length and SEE elongation (\(E_{SEE}\): \(L_{SEE} = L_{SEE}^0 + E_{SEE}\)). This allowed us to rewrite Equation 2 to estimate PEE force (\(F_{PEE}^{est}\)), yielding:
\[ F_{PEE}^{est} = k_{PEE} \cdot (L_{MTC} - - L_{SEE}^0 - E_{SEE} - L_{PEE}^0)^2 \tag{11}\]
Now, there are two unknown parameter values (i.e., \(k_{PEE}\) and \(L_{PEE}^0\)), while also SEE elongation is unknown for every quick-release experiment. Consequently, there are more unknowns than equations. To address this, we replaced \(L_{PEE}^0 + L_{SEE}^0\) with a temporary parameter \(c_{PEE}\). Subsequently, SEE elongation was computed as \(\sqrt{\frac{F_{SEE}^{MIN,data}}{k_{SEE}}}\), given the estimated SEE stiffness scaling factor in the previous step. This yielded:
\[ F_{PEE}^{est} = k_{PEE} \cdot (L_{MTC} - \sqrt{\frac{F_{SEE}^{MIN,data}}{k_{SEE}}} - c_{PEE})^2 \tag{12}\]
As \(k_{PEE}\) only scales PEE stiffness, its value is constant across all quick-release experiments. Similarly, \(c_{PEE}\) is the sum of \(L_{SEE}^0\) and \(L_{PEE}^0\) and should therefore be constant for each quick-release experiment. The values of \(k_{PEE}\) and \(c_{PEE}\) were then computed by minimising the sum of squared differences between \(F_{PEE}^{est}\) and \(F_{SEE}^{min,data}\).
Estimation of \(F_{CE}^{max}\), \(L_{CE}^{opt}\) and \(L_{SEE}^0\). The third step was to estimate \(F_{CE}^{max}\), \(L_{CE}^{opt}\) and \(L_{SEE}^0\). Before the quick-release, MTC is isometrically delivering force. We used the SEE force at this instant (\(F_{SEE}^{QR-,data}\), black dot in Figure 3A) and the corresponding MTC length to obtain the MTC force-length relationship from the data. We then estimated maximal isometric CE force (\(F_{CE}^{max}\)), CE optimum length (\(L_{CE}^{opt}\)) and SEE slack length (\(L_{SEE}^0\)) by minimising the sum of squared differences between the data and estimated MTC force-length relationship. Lastly, \(L_{PEE}^0\) was computed by subtracting \(L_{SEE}^0\) from \(c_{PEE}\).
2.3.2 CE force-velocity parameter value estimation
The values of the CE force–velocity relationship parameters \(a\) and \(b\) were estimated from data obtained during the plateau phase of SEE force in step-ramp experiments. This approach was chosen because SEE length is constant when SEE force is constant. Consequently, CE velocity equals MTC velocity under these conditions. Following this argument, we first identified a 10 ms interval in which SEE force changed the least. Second, we derived CE length and CE force as functions of time using Equation 1 and Equation 2. Third, we calculated the CE velocity as the time-derivative of CE length. Finally, we averaged CE force and CE velocity over the 10 ms interval. This procedure yielded the data of the CE force-velocity relationship.
Now, the model CE force-velocity relationship can be fit to that of the data to obtain value of \(a\) and \(b\). However, employing a method that simply minimises the sum of ordinary least squares differences would not be appropriate because in experiments there is uncertainty in both measured CE force and CE velocity. Instead, we used a total least squares method that minimised the distance between the modelled CE force (\(F_{CE}^{est}\)) and CE velocity (\(v_{CE}^{est}\)) (now called: model values) and those of the data (see also Figure 5). To do this, we had to find the nearest model values to the CE force of the data (\(F_{CE}^{data}\)) and CE velocity of the data (\(v_{CE}^{data}\)) (now called: datapoints). We used the following observations: 1) the model value has to satisfy Equation 4 and 2) the derivative of \(\frac{dF_{CE}}{dv_{CE}}\) of the model value should be perpendicular to the line from the datapoint to the model value (see Figure 5). Hence, the model value could be found by solving Equation 4 and Equation 13:
\[ \frac{dv_{CE}^{est}}{dF_{CE}^{est}} \cdot (v_{CE}^{est} - v_{CE}^{data}) = F_{CE}^{data} - F_{CE}^{est} \tag{13}\]
The following cost-function was then minimised to find the parameter values of \(a\) and \(b\):
\[ J = \sum_{i=1}^{n} {\biggl(\frac{F_{CE}^{data} - F_{CE}^{est}}{c_1}\biggr)}^2 + {\biggl(\frac{v_{CE}^{data} - v_{CE}^{est}}{c_2}\biggr)}^2 \]
\(c_1\) and \(c_2\) denote scaling factors such that both terms of the cost-function are more or less equally weighted, which was done by setting \(c_1\) equal to the maximal range in CE force of the data and setting \(c_2\) equal to the maximal range of CE velocity of the data.
2.3.3 Excitation dynamics parameter value estimation
We estimated the parameters values of \(\tau_{act}\) and \(\tau_{deact}\) based on isometric experiments. The choice to estimate only these parameters of the excitation dynamics was made because it is generally challenging to discriminate between parameter values of both the calcium dynamics (Equation 7) and the \(q-[Ca^{2+}]\) relation (Equation 8) based on mechanical measurements outlined above. The parameter values of the \(q-[Ca^{2+}]\) relation were set to those used in simulating the isometric experiments. Consequently, these parameters matched their actual values for the simulated data, but obviously were suboptimal for the in situ data. To estimate the parameter values of \(\tau_{act}\) and \(\tau_{deact}\), we minimised the sum of squared differences between the SEE force of the data and the estimated SEE force based on the parameter values. For this, an interval of the data was used starting from maximal CE stimulation to 0.1 s after CE stimulation ceased off.
2.3.4 Parameter estimation: traditional method & improved method
To estimate the SEE stiffness it is generally assumed that the quick-release is so fast that all MTC shortening is taken up by SEE. As A.V. Hill already acknowledged in (1950), this is not the case in reality as CE also shortens. Although this CE shortening might be small, SEE is in general stiff and therefore small errors in SEE length may result in large errors in predicted SEE forces and thus in estimated SEE stiffness. We used two different methods to estimate the parameter values: 1) the traditional method and the improved method. In the traditional method, we followed the procedure outlined in section Section 2.3.1 and Section 2.3.2. In the improved method, we also followed the procedure outlined in Section 2.3.1 and Section 2.3.2 after which we incorporated an additional procedure to correct for CE shortening due to the quick-release. This additional procedure was as follows: 1) We calculated CE length just before the (step) change in MTC length occurred using Equation 1, 2, 3 and 8, assuming that CE velocity was 0 and that \(\gamma\) equalled 1. 2) We performed a (short) simulation from the time just before the (step) change in MTC length to just after the (step) change in MTC length. 3) We computed the change in CE length over this time interval as the average CE length slightly before and slightly after the (step) change in MTC length. 4) We subtracted the change in CE length from the change in MTC length to obtain the change in SEE length. In addition to correcting for CE length change during the quick-release, we incorporated a step to obtain better estimates of the actual PEE force at the time instance of minimum SEE force (\(F_{SEE}^{min,data}\); see Figure 3). We computed CE force at this time instance using Equation 3 and 8, assuming that \(\gamma\) equalled its minimum value. This CE force was then subtracted from \(F_{SEE}^{min,data}\).
The improved method leads to a different (i.e., higher) SEE stiffness, which affect the estimation of the parameter values of the CE, PEE and SEE force-length relationships (see Section 2.3.1). These estimated parameter values were then used to again estimate the parameter values of the CE force-velocity relationship (see section Section 2.3.2). As the parameter values of the CE force-velocity relationship slightly changed, the estimate of CE length change due to the quick-release also changed and therefore estimated SEE stiffness changed. We found that this process always converged to a stable solution and therefore we used an iterative process until the change in all parameter values was less than 0.1% (see Figure S2).
2.4 Sensitivity analysis
2.4.1 Monte Carlo simulations
In experiments on isolated MTC, the MTC force-length relationship can shift due to irreversible damage of SEE (Aubert, Roquet, and Van der Elst 1951). We observed this phenomenon in an experiment involving isolated rat m. gastrocnemius medialis, where the MTC force-length relationship shifted by about 1 mm (~8% \(L_{CE}^{opt}\)) over a period of approximately 8 hours (Reuvers et al., unpublished observation). To assess the influence of these shifts on the parameter value estimation, we performed Monte Carlo simulations. We assumed that the observed shifts of the MTC force-length relationship were caused solely by a decrease in SEE stiffness. Accordingly, we decreased SEE stiffness to cause random shifts of the MTC force-length relationship between 0 and 1 mm. This way, we simulated data as if they were collected at random intervals throughout an 8-hour period. Importantly, changes in SEE stiffness affect the step size and constant MTC velocity ramp required for a more or less a constant SEE force at the beginning of the ramp in step-ramp experiments. Therefore, we adjusted the step size and constant MTC velocity ramp for each simulated step-ramp trial individually. Parameter values of each MTC were then estimated using the ‘perturbed’ data using the improved method (see Section 2.3.4). This process was repeated 50 times for each MTC, allowing us to evaluate both the mean and the spread of the estimated parameter values.
2.4.2 Interdependency of parameter values
The parameter estimation procedure detailed above follows a fixed order in which certain parameters are estimated before others. For instance, \(k_{SEE}\) is estimated first and therefore affects the estimation of \(L_{SEE}^0\) and \(L_{CE}^{opt}\). This sequence creates an interdependency among the parameter values, which has been identified as a potential cause for substantial variation in parameter values among individuals, even within the same muscle type (Lemaire et al. 2016). To investigate this interdependency, we systematically changed parameter values of our model. Specifically, we increased and decreased each parameter value by 5% relative to the values obtained with the improved method. After this, we re-estimated all other parameter values. This approach allowed us to assess the influence of changes in one parameter on the estimation of others.
3 Results
3.1 Evaluating parameter value estimation accuracy - simulated data
Code
# %% Load packages
import os, pickle, sys
import numpy as np
import pandas as pd
from pathlib import Path
# Set directories
cwd = Path.cwd()
baseDir = cwd.parent
dataDir = baseDir / 'data'
funcDir = baseDir / 'analysis' / 'functions'
sys.path.append(str(funcDir))
# Custom functions
import hillmodel, stats
# %% Extract parameters
muscles = ['GMs1','GMs2','GMs3']
param_keys = ['a', 'b', 'kpee', 'ksee', 'fmax', 'lce_opt', 'lpee0', 'lsee0', 'tact', 'tdeact']
orParms, tmParms, imParms = [], [], []
for mus in muscles:
orPar = pickle.load(open(os.path.join(dataDir,mus,'parameters',mus+'_OR.pkl'), 'rb'))
tmPar = pickle.load(open(os.path.join(dataDir,mus,'parameters',mus+'_TM.pkl'), 'rb'))[0]
imPar = pickle.load(open(os.path.join(dataDir,mus,'parameters',mus+'_IM.pkl'), 'rb'))[0]
parms = [orPar[k] for k in param_keys]
orParms.append(parms)
parms = [tmPar[k] for k in param_keys]
tmParms.append(parms)
parms = [imPar[k] for k in param_keys]
imParms.append(parms)
# Store in pandas dataframe
df_or = pd.DataFrame(list(zip(*orParms)), columns=['GM1', 'GM2', 'GM3'], index=param_keys) # df with actual/real parameter values
df_tm = pd.DataFrame(list(zip(*tmParms)), columns=['GM1', 'GM2', 'GM3'], index=param_keys) # df with estimated parameter values using the traditional method
df_im = pd.DataFrame(list(zip(*imParms)), columns=['GM1', 'GM2', 'GM3'], index=param_keys) # df with estimated parameter values using the improved method
# %% Compute percentage differences from real/actual and estimated parameter values
# Positive: Estimated is that % higher than Real
# Negative: Estimated is that % lower than Real
df_tm_p = stats.pdiff(df_tm, df_or)
df_im_p = stats.pdiff(df_im, df_or)
# %% Variables related to overestimation of length changes due to QR
deltaLsee_tm = np.full((len(muscles),10), np.nan)
deltaLsee_real = np.full((len(muscles),10), np.nan)
for iMus,mus in enumerate(muscles):
orPar = pickle.load(open(os.path.join(dataDir,mus,'parameters',mus+'_OR.pkl'), 'rb'))
tmPar,dataQRtm = pickle.load(open(os.path.join(dataDir,mus,'parameters',mus+'_TM.pkl'), 'rb'))[0:2]
fseeQRpre = dataQRtm['fseeQRpre']
fseeQRpst = dataQRtm['fseeQRpst']
lseeQRpre = (fseeQRpre/orPar['ksee'])**0.5 + orPar['lsee0']
lseeQRpst = (fseeQRpst/orPar['ksee'])**0.5 + orPar['lsee0']
# Compute change in SEE length due to QR
deltaLsee_tm[iMus] = dataQRtm['lseeQRpre']-dataQRtm['lseeQRpst'] # [m] estimated SEE length change due to QR
deltaLsee_real[iMus] = lseeQRpre-lseeQRpst # [m] real/actual SEE length change due to QR
# Compute percentage difference from real to estimated SEE length change due to QR
# Positive: TM is that % higher than Real
# Negative: TM is that % lower than Real
p_deltaLsee = stats.pdiff(deltaLsee_tm, deltaLsee_real) # [%] percentage difference
p_deltaLseeAvg = np.mean(p_deltaLsee,1) # [%] average per muscle
p_deltaLseeStd = np.std(p_deltaLsee,1) # [%] std per muscle
# Compute difference between real and estimated CE length change due to QR
deltaLce = deltaLsee_real-deltaLsee_tm # [m]
deltaLceAvg = np.mean(deltaLce,1) # [m] average per muscle
deltaLceStd = np.std(deltaLce,1) # [m] std per muscle
# %% 3.1.1: Traditional vs. Improved method
# Extract all contraction dynamics parameters except ksee
p_TM_con = df_tm_p.loc[['a','b', 'kpee', 'fmax', 'lce_opt', 'lpee0', 'lsee0']]
# Avg. difference: from real to estimated SEE length change
p_TM_deltaLsee_avgall = f'{p_deltaLseeAvg.mean():0.0f}' # [%] avg. percentage difference
# Avg. difference over all muscles: from real to estimated SEE stiffness
p_TM_ksee_avg = f'{df_tm_p.loc['ksee'].mean():0.0f}' # [%] avg. percentage difference
# Avg. difference over all muscles: from real to estimated tact
p_TM_tact_avg = f'{df_tm_p.loc['tact'].mean():0.0f}%' # [%] avg. percentage difference
# Max. difference over all muscles: from real to estimated tact
p_TM_con_max = f'{p_TM_con.abs().max().max():0.1f}' # [%] max. percentage difference
# Max. difference over all muscles: from real to estimated tact
p_IM_max = f'{df_im_p.abs().max().max():0.0f}' # [%] max. percentage difference
# %% 3.1.2: Estimation of SEE stiffness.
# Difference between SEE elongation of estimated parameter value and real/actual value
# Positive: estimated value is higher than real/actual value
d_esee = (df_tm.loc['fmax']/df_tm.loc['ksee'])**0.5 - (df_or.loc['fmax']/df_or.loc['ksee'])**0.5 # [m]
# Correlation coefficients
r_see = np.full(len(muscles), np.nan)
r_pee = np.full(len(muscles), np.nan)
r_mtc = np.full(len(muscles), np.nan)
for iMus,mus in enumerate(muscles):
parFile = os.path.join(dataDir,mus,'parameters',mus+'_TM.pkl')
muspar, dataQR, dataSR, dataACTout = pickle.load(open(parFile, 'rb'))
# SEE
lseeData = np.hstack((dataQR['lseeQRpre'],dataQR['lseeQRpst']))
fseeData = np.hstack((dataQR['fseeQRpre'],dataQR['fseeQRpst']))
fseeMdl = hillmodel.LEE2Force(lseeData,0,muspar)[0]
r_see[iMus] = np.corrcoef(fseeData, fseeMdl)[0,1]
# PEE
fpeeMdl = hillmodel.LEE2Force(0,dataQR['lpeeQR'],muspar)[1]
r_pee[iMus] = np.corrcoef(dataQR['fpeeQR'], fpeeMdl)[0,1]
# MTC
fseeMdl = hillmodel.ForceEQ(dataQR['lmtcQRpre'],1,muspar)[0]
r_mtc[iMus] = np.corrcoef(dataQR['fseeQRpre'], fseeMdl)[0,1]
r_all = np.hstack((r_see,r_pee,r_mtc))
# Percentage difference from real to estimated parameter
p_TM_ksee = [f'{v:0.0f}' for v in df_tm_p.loc['ksee']] # [%] per muscle
# Avg&Std per muscle: differences from real to estimated parameters
p_TM_deltaLsee_avgmus = [f'{v:0.0f}' for v in p_deltaLseeAvg] # [%]
p_TM_deltaLsee_stdmus = [f'{v:0.0f}' for v in p_deltaLseeStd] # [%]
# Avg&Std per muscle: CE shortening during QR
p_TM_deltaLce_avgmus = [f'{v*1e6:0.0f}' for v in deltaLceAvg] # [um] avg. per muscle
p_TM_deltaLce_stdmus = [f'{v*1e6:0.0f}' for v in deltaLceStd] # [um] std per muscle
# Difference from real to estimate in SEE elongation @ fmax
d_TM_esee = [f'{v*1e3:0.2f}' for v in d_esee] # [mm]
# Minimum correlation coefficient
r_TM_min = f'{np.min(np.floor(r_all*100)/100):0.2f}' # []
# %% 3.1.2: Estimation of PEE parameter values
nPEEdata = np.full(len(muscles), np.nan)
for iMus,mus in enumerate(muscles):
muspar, dataQR, dataSR, dataACT = pickle.load(open(os.path.join(dataDir,mus,'parameters',mus+'_TM.pkl'), 'rb'))
lpeeQR = dataQR['lpeeQR'] # [m] PEE length of data
fpeeQR = dataQR['fpeeQR'] # [N] PEE force of data
fpeeModel = hillmodel.LEE2Force(0,lpeeQR,muspar)[1] # [N] PEE force estimated
# Compute how many data points available
nPEEdata[iMus] = np.sum(fpeeModel>0)
# Dislay percentage differences from real to estimated parameters
p_TM_kpee = [f'{v:0.1f}' for v in df_tm_p.loc['kpee']] # [%] percentage difference
p_TM_lpee0 = [f'{v:0.1f}' for v in df_tm_p.loc['lpee0']] # [%] percentage difference
# Display amount of datapoints available
n_PEE_data = [f'{v:0.0f}' for v in nPEEdata] # [ ]
# %% 3.1.2: Estimation of Fcemax, Lceopt, Lsee0
# Dislay percentage differences from real to estimated parameters
p_TM_fl_max = f'{df_tm_p.loc[['fmax','lce_opt','lsee0'],:].abs().max().max():0.1f}' # [%] maximal percentage difference
# %% 3.1.2: Estimation of CE force-velocity parameter values
# Dislay percentage differences from real to estimated parameters
p_TM_a_max = f'{df_tm_p.loc['a',:].abs().max():0.1f}' # [%] maximal percentage difference
p_TM_b_max = f'{df_tm_p.loc['b',:].abs().max():0.1f}' # [%] maximal percentage difference
# %% 3.1.2: Estimation of excitation dynamics parameter values
p_tact_QR = np.full(len(muscles), np.nan)
p_tdeact_QR = np.full(len(muscles), np.nan)
p_tact_SR = np.full(len(muscles), np.nan)
p_tdeact_SR = np.full(len(muscles), np.nan)
for iMus, mus in enumerate(muscles):
orPar = pickle.load(open(os.path.join(dataDir,mus,'parameters',mus+'_OR.pkl'), 'rb'))
qrPar = pickle.load(open(os.path.join(dataDir,mus,'parameters',mus+'_TM_QR.pkl'), 'rb'))[0]
srPar = pickle.load(open(os.path.join(dataDir,mus,'parameters',mus+'_TM_SR.pkl'), 'rb'))[0]
# Compute percentage difference from real to estimated parameter
# Positive: estimated is that % higher than Real
# Negative: estimated is that % lower than Real
p_tact_QR[iMus] = stats.pdiff(qrPar['tact'],orPar['tact'])
p_tact_SR[iMus] = stats.pdiff(srPar['tact'],orPar['tact'])
p_tdeact_QR[iMus] = stats.pdiff(qrPar['tdeact'],orPar['tdeact'])
p_tdeact_SR[iMus] = stats.pdiff(srPar['tdeact'],orPar['tdeact'])
# Dislay percentage differences from real to estimated parameters by using ISOM data
p_TM_tact = [f'{v:0.0f}' for v in df_tm_p.loc['tact']] # [%] percentage difference
p_TM_tdeact = [f'{v:0.1f}' for v in df_tm_p.loc['tdeact']] # [%] percentage difference
# Display percentage difference from real to estimated parameters by using
# either QR or SR data
p_tact_QR_avg = f'{p_tact_QR.mean():0.0f}' # [%] percentage difference when using QR data
p_tact_SR_avg = f'{p_tact_SR.mean():0.0f}' # [%] percentage difference when using SR data
p_tdeact_QR_avg = f'{p_tdeact_QR.mean():0.0f}' # [%] percentage difference when using QR data
p_tdeact_SR_avg = f'{p_tdeact_SR.mean():0.0f}' # [%] percentage difference when using SR data
# For tdeact: differences are similar so show avg of the two
p_tdeact_QRSR = np.concatenate([p_tdeact_QR, p_tdeact_SR])
p_tdeact_QRSR_avg = f'{p_tdeact_QRSR.mean():0.0f}' # [%] percentage difference when using QR/SR data3.1.1 Traditional method versus Improved method
The traditional method did not correct for CE shortening during quick-release experiments. Consequently, SEE shortening due to the quick-release was overestimated by 24% on average. This overestimation of SEE shortening caused a decrease in \(k_{SEE}\) (and thus SEE stiffness) of 34% on average. The estimation of the activation time constant was most affected by the underestimated SEE stiffness, resulting in an underestimation of 19%% on average. All other estimated parameter values were within 8.7% of their actual values (Table 1). These findings show the influence of CE shortening — even within a very brief 10 ms interval — on the parameter value estimation of both the contraction dynamics and excitation dynamics.
In contrast, the improved method corrected for CE shortening during the quick-release experiments. The correction substantially reduced the overestimation of SEE shortening, yielding an accurate estimate of \(k_{SEE}\). As a result, all estimated parameter values deviated by no more than 3% from their actual values (Table 1). These results demonstrate that accounting for CE shortening - even over a brief 10 ms interval - substantially improves the accuracy of the parameter value estimation.
For interested readers, we discuss below the specific factors contributing to the differences between the actual contraction and excitation dynamics parameter values and those estimated based on the traditional method.
3.1.2 Understanding errors in parameter value estimation
Estimation of SEE stiffness. The SEE stiffness parameter (\(k_{SEE}\)) was underestimated by 39%, 27% and 37% for GM1, GM2 and GM3, respectively (see Table 1). This overestimation directly resulted from the overestimation of SEE shortening due to the quick-release, which was overestimated by 28±2%, 17±2% and 26±2% for GM1, GM2 and GM3, respectively, averaged across all quick-release experiments. In absolute terms, CE shortened only 43±3, 29±2 and 41±2 µm over the 11 ms interval between the time points at which SEE force and SEE length were sampled (i.e., immediately before and after the quick-release). The amount of CE shortening was smallest in GM2 because GM2 was a slower muscle than GM1 and GM3. Consequently, the overestimation of SEE shortening was also smallest in GM2, which in turn led to the smallest underestimation of \(k_{SEE}\). All in all, this shows that even very small amounts of CE shortening has substantial influence on the estimation of SEE stiffness.
Since SEE was very stiff, the overestimation in SEE elongation at maximal isometric CE force was only 0.50, 0.32 and 0.49 mm for GM1, GM2 and GM3, respectively. This is an important finding because the SEE stiffness is the first parameter in the estimation process and therefore affects all subsequent estimated parameter values. Therefore, since the overestimation of SEE elongation (in mm) at maximal isometric CE force was only small, the impact on the estimation of the other contraction dynamics parameter values was also small.
Lastly, it should be noted that the correlation coefficient between the data SEE force‑length relationship and the one based from the estimated parameters is not an adequate measure of how accurately SEE stiffness is captured. Obviously, a correlation coefficient is not sensitive to the overestimation of SEE shortening due to the quick-release. Consequently, high correlation coefficients (\(R^2 >\) 0.99) can still be observed even when SEE stiffness is substantially underestimated. To further investigate this issue, we simulated the quick-release experiments after obtaining all contraction and excitation dynamics parameter values. The simulations clearly showed a slower rise in SEE force after the quick-release in comparison with the experimental data (Figure 6A), indicating that SEE stiffness was underestimated. These results suggest that the rise in SEE force is a better measure of SEE stiffness estimation accuracy.
3.2 Sensitivity analysis
Code
# %% Load packages
import os, pickle, sys
import pandas as pd
from pathlib import Path
# Set directories
cwd = Path.cwd()
baseDir = cwd.parent
dataDir = baseDir / 'data'
funcDir = baseDir / 'analysis' / 'functions'
sys.path.append(str(funcDir))
# Custom functions
import stats
# %% 3.2.1: Monte Carlo simulations
muscles = ['GMs1', 'GMs2', 'GMs3']
# Parameters to track
params = ['a', 'b', 'fmax', 'kpee', 'ksee', 'lce_opt', 'lpee0', 'lsee0', 'tact', 'tdeact']
# Store Monte-Carlo results
results = []
for mus in muscles:
par_file = os.path.join(dataDir, mus, 'parameters', f'{mus}_OR.pkl')
or_par = pickle.load(open(par_file, 'rb'))
eseeMaxOR = (or_par['fmax']/or_par['ksee'])**0.5
for iMC in range(1, 51):
mc_file = os.path.join(dataDir, mus, 'parameters','mc', f'{mus}_MC{iMC:02d}.pkl')
mc_par = pickle.load(open(mc_file, 'rb'))[0]
eseeMaxMC = (mc_par['fmax']/mc_par['ksee'])**0.5
entry = {'muscle': mus, 'MC': iMC, 'd_esee': (eseeMaxMC-eseeMaxOR)}
for param in params:
entry[param] = stats.pdiff(mc_par[param], or_par[param])
results.append(entry)
# Convert to pandas DataFrame
df = pd.DataFrame(results)
# Compute summary statistics
df_mean = df.groupby('muscle').mean().T
df_std = df.groupby('muscle').std().T
# Extract all contraction dynamics parameters except ksee
p_MC_con_avg = df_mean.loc[['a','b', 'kpee', 'fmax', 'lce_opt', 'lpee0', 'lsee0']]
p_MC_con_std = df_std.loc[['a', 'b', 'fmax', 'kpee', 'lce_opt', 'lpee0', 'lsee0']]
# Percentage difference over all muscles: from real to estimated ksee
p_MC_ksee_avg = f'{df_mean.loc['ksee'].mean():0.0f}' # [%] percentage difference
# Absolute difference over all muscles of SEE elongation @ Fcemax
d_MC_esee_avg = f'{df_mean.loc['d_esee'].mean()*1e3:0.1f}' # [mm]
# Max. avg. percentage difference over all muscles: from real to estimated value of contraction dyn. parms
p_MC_con_maxavg = f'{p_MC_con_avg.abs().max().max():0.1f}'
# Max. std. percentage difference over all muscles: from real to estimated value of contraction dyn. parms
p_MC_con_maxstd = f'{p_MC_con_std.abs().max().max():0.0f}'
# Max. std percentage difference over all muscles: from real to estimated value of kpee
p_MC_kpee_maxstd = f'{df_std.loc['kpee'].max():0.0f}'
# Max. avg percentage difference over all muscles: from real to estimated value of tact
p_MC_tact_maxavg = f'{df_mean.loc['tact'].max():0.1f}'
# Avg. std percentage difference over all muscles: from real to estimated value of tact
p_MC_tact_avgstd = f'{df_std.loc['tact'].mean():0.0f}'3.2.1 Monte Carlo simulations
Monte Carlo simulations were used to examine how shifts in the MTC force–length relationship caused by a decrease in SEE stiffness (e.g., due to irreversible SEE damage in experiments) affect the accuracy of the estimated contraction and excitation dynamics parameter values. The induced shifts of the MTC force-length relationship were between 0–1 mm, and therefore it was no surprise that SEE stiffness decreased by 36% on average — which corresponds to a .5 mm shift.
All other estimated contraction dynamics parameter values were within 6.3% on average of their actual value (Table 1). The variance in the estimated contraction dynamics parameter values was below 31% for all parameters except for the PEE stiffness scaling parameter. The PEE stiffness scaling parameter showed substantially higher variance (with a standard deviation up to 31%), because it was based on four or fewer quick-release trials, making it more sensitive to perturbations in the data. Regarding the excitation dynamics, the influence on the average time constants was minimal (within 0.7%), but affected the variance in the estimated values of the activation time constant (with a standard deviation of 4% on average). These findings indicate that an average decrease in SEE stiffness of 36% has a much smaller effect on the estimated parameter values of the contraction and excitation dynamics, even at the level of an individual muscle.
3.2.2 Interdependency of parameter values
We investigated the interdependency of the estimated parameter values by adjusting each parameter value by 5% and re-estimating all other parameter values. \(F_{CE}^{max}\) and \(L_{SEE}^0\) were the parameters that had most effect on the estimates of the others.
First, \(F_{CE}^{max}\) affected the estimation of parameter \(a\) and \(b\) of the CE force-velocity relationship. Underestimating \(F_{CE}^{max}\) leads to an underestimation of CE force at 0 velocity. Due to the formulation of the CE force-velocity relationship, the curve, by definition, crosses the point at \(v_{CE}=0\) and \(F_{CE}=F_{CE}^{max}\). Consequently, the best fit to the data with an underestimated \(F_{CE}^{max}\) is a much flatter CE force-velocity relationship, which is realised by an increase in \(a\) and a decrease in \(b\). This flatter CE force-velocity relationship affects the estimation of all other contraction dynamics parameter values because the CE force-velocity relationship is used to estimate the CE shortening due to the quick-release (see Section 2.3.4). As a result, \(k_{SEE}\) and \(L_{SEE}^0\) decreased due to underestimating \(F_{CE}^{max}\), whereas \(L_{CE}^{opt}\) increased. In summary, \(F_{CE}^{max}\) substantially influenced the estimation of all contraction dynamics parameter values mainly by its effect on the CE force-velocity relationship.
Second, \(L_{SEE}^0\) affected the estimation of \(L_{CE}^{opt}\) and \(F_{CE}^{max}\). As explained earlier, when \(L_{SEE}^0\) is overestimated, \(L_{CE}^{opt}\) is underestimated, such that \(L_{MTC}^{opt}\) remains more or less unaffected. Underestimating \(L_{CE}^{opt}\) causes a narrower MTC force-length relationship, leading to an overestimation of \(F_{CE}^{max}\) to preserve a good fit between the data and the model. As previously discussed, overestimating \(F_{CE}^{max}\), in turn, resulted in an underestimation of parameters \(a\) and \(b\). Thus, \(L_{SEE}^0\) substantially influenced the estimation of all contraction dynamics parameter values mainly by its effect on the MTC force-length relationship.
Taken together, the interdependence of \(F_{CE}^{max}\) and \(L_{SEE}^0\) with other parameter values underscores the need for precise estimation of these key parameters. In this regard, it is reassuring that these two parameters were found to be robust for perturbations in the experimental data.
3.3 Evaluating model predictions - in situ data
Code
# %% Load packages
import os, pickle, sys, glob
import numpy as np
import pandas as pd
from pathlib import Path
# Directories
cwd = Path.cwd()
baseDir = cwd.parent
dataDir = baseDir / 'data'
funcDir = baseDir / 'analysis' / 'functions'
sys.path.append(str(funcDir))
# Custom functions
import helpers, stats
# %% Extract parameters
muscles = ['GMe1', 'GMe2', 'GMe3']
param_keys = ['a', 'b', 'kpee', 'ksee', 'fmax', 'lce_opt', 'lpee0', 'lsee0', 'tact', 'tdeact']
orParms, tmParms, imParms = [], [], []
for mus in muscles:
tmPar, dataQRtm = pickle.load(open(os.path.join(dataDir,mus,'parameters',mus+'_TM.pkl'), 'rb'))[0:2]
imPar, dataQRim = pickle.load(open(os.path.join(dataDir,mus,'parameters',mus+'_IM.pkl'), 'rb'))[0:2]
parms = [tmPar[k] for k in param_keys]
tmParms.append(parms)
parms = [imPar[k] for k in param_keys]
imParms.append(parms)
# Store in pandas dataframe
df_tm = pd.DataFrame(list(zip(*tmParms)), columns=['GM1', 'GM2', 'GM3'], index=param_keys) # df with estimated parameter values using the traditional method
df_im = pd.DataFrame(list(zip(*imParms)), columns=['GM1', 'GM2', 'GM3'], index=param_keys) # df with estimated parameter values using the improved method
# %% Compute percentage differences from real/actual and estimated parameter values
# Positive: Imporoved is that % higher than Traditional
# Negative: Imporoved is that % lower than Traditional
df_p = stats.pdiff(df_im, df_tm)
# %% Compute RMSDs
par_models = ['TM', 'IM']
experiments = ['QR', 'SR', 'ISOM', 'SSC']
muscles = ['GMe1', 'GMe2', 'GMe3']
avg_mean_all = []
columns = [f'{model}_{exp}' for model in par_models for exp in experiments]
rows = muscles + ['Avg ± Std']
df = pd.DataFrame(index=rows, columns=columns, dtype=str)
for exp in experiments:
for model in par_models:
all_rmsd = []
for mus in muscles:
# Load parameter
parFile = os.path.join(dataDir, mus, 'parameters', f'{mus}_{model}.pkl')
muspar = pickle.load(open(parFile, 'rb'))[0]
dataDirExp = os.path.join(dataDir, mus, 'dataExp', exp)
rrunDirExp = os.path.join(dataDir, mus, 'simsExp', model, exp)
dataFiles = sorted(glob.glob(os.path.join(dataDirExp, '*')))
rrunFiles = sorted(glob.glob(os.path.join(rrunDirExp, '*')))
rms_list = []
for dataFile, rrunFile in zip(dataFiles, rrunFiles):
dataFilename = os.path.basename(dataFile)[:-4]
rrunFilename = os.path.basename(rrunFile)[:-4]
if dataFilename[:-3] != rrunFilename[:-3]:
raise ValueError(f"File mismatch: {dataFilename} vs {rrunFilename}")
dataData = pd.read_csv(dataFile).T.to_numpy()[0:4]
rrunData = pd.read_csv(rrunFile).T.to_numpy()
time1, _, stim1, fsee1 = dataData
time2, _, _, fsee2 = rrunData
tStimOn, tStimOff = helpers.get_stim(time1, stim1)[1:]
iStart = int(np.argmin(abs(time1 - tStimOn[0])))
if exp == 'ISOM':
iStop = int(np.argmin(abs(time1 - 0.1 - tStimOff[0])))
elif exp in ['QR', 'SR']:
iStop = int(np.argmin(abs(time1 - tStimOff[0])))
else: # SSC
iStop = int(np.argmin(abs(time1 - tStimOff[-1])))
rms = stats.rmse(fsee1[iStart:iStop], fsee2[iStart:iStop]) / muspar['fmax'] * 100
rms_list.append(rms)
M = np.mean(rms_list)
S = np.std(rms_list)
df.loc[mus, f'{model}_{exp}'] = f'{M:.1f} ± {S:.1f}'
all_rmsd += rms_list # just mean, for averaging across muscles
# Average over muscles for this (exp, model)
avg_M = np.mean(all_rmsd)
avg_S = np.std(all_rmsd)
df.loc['Avg ± Std', f'{model}_{exp}'] = f'{avg_M:.1f} ± {avg_S:.1f}'
avg_mean_all.append(avg_M)
avg_mean_all = np.array(avg_mean_all)
#%% 3.3: Evaluating model predictions
# Avg. difference: from traditional to improved method: SEE stiffness
p_IS_ksee_avg = f'{df_p.loc['ksee'].mean():0.0f}' # [%] avg. percentage difference
# Avg. difference: from traditional to improved method: tact
p_IS_tact_avg = f'{df_p.loc['tact'].mean():0.0f}' # [%] avg. percentage difference
# Avg. difference: from traditional to improved method: tdeact
p_IS_tdeact_avg = f'{df_p.loc['tdeact'].mean():0.0f}' # [%] avg. percentage difference
# Avg. difference: from traditional to improved method: lce_opt
p_IS_lceopt_avg = f'{df_p.loc['lce_opt'].mean():0.0f}' # [%] avg. percentage difference
# Take the mean over all each sort of experiment
p_IS_rmsd = f'{stats.pdiff(avg_mean_all[1::2],avg_mean_all[0::2]).mean():2.0f}'We estimated contraction and excitation dynamics parameter values from in situ data of rat m. gastrocnemius medialis (Table S2). The improved method yielded a 67% higher SEE stiffness compared to the traditional method. As explained earlier, higher SEE stiffness also affects the estimation of all other parameter values. The most noticeable changes were longer time constants for both the activation (\(\tau_{act}\)) and deactivation (\(\tau_{deact}\)) dynamics, with increases of 61% and 16% on average, respectively. Another noticeable change was in CE optimum length, which was about 11% higher using the improved method compared to the traditional method.
To assess model predictions, we re-ran quick-release, step-ramp and isometric experiments, as well as stretch-shortening cycles with experimentally measured MTC length and CE stimulation over time as inputs to the model. We compared the model predictions using parameter sets obtained from both methods. For all experiments and rats, the average difference between experimental data and model predictions was -30% smaller on average using parameters from the improved method (Figure 8; Table S3). Overall, the improved method substantially enhanced the predictions derived by the Hill-type MTC model compared to the traditional method.
4 Discussion
The aim of this study was to evaluate the accuracy of estimating contraction and excitation dynamics parameter values on the basis of data collected in commonly used experimental protocols. In real experiments, the actual parameter values are unknown, making it impossible to directly assess estimation accuracy. In this study, we took a different approach: using a Hill-type MTC model with parameter values obtained from literature we generated synthetic ‘data’ by simulating quick-release, step-ramp, and isometric experiments. We then estimated all important contraction and excitation dynamics parameter values by using the synthetic data. Since the actual parameter values were known in this case, we could assess how accurately they were retrieved. Two different estimation methods were compared. The first was the traditional method, commonly used in the literature, which does not account for muscle fibre shortening during quick-release experiments. The second was an improved method that includes a correction for muscle fibre shortening during the quick release. Both methods developed in this study—designed to estimate contraction and excitation dynamics parameters from quick-release, step-ramp, and isometric experiments using servomotors—are made available as an open-source toolbox. In the remainder of this paper, we will discuss 1) the difference between the two methods; 2) the robustness of the improved method and 3) the implications for muscle modelling.
In quick-release experiments using servomotors, muscle fibres shorten, even within the short duration of the release. Obviously, the extent of muscle fibre shortening depends on the duration of the quick-release. Consequently, the longer the quick-release, the more important it becomes to account for muscle fibre shortening when estimating SEE stiffness (see Figure S3). For a typical quick-release duration of 10 ms, we found that muscle fibre shortening was nearly 25% of MTC shortening (and thus 33% of SEE shortening), resulting in an underestimation of SEE stiffness by approximately 35%. This finding is important, especially in studies aiming to estimate SEE energy storage from SEE force or length, since energy storage estimates scale linearly with SEE stiffness. Although SEE was substantially underestimated with the traditional method, this had minimal effect on the estimated parameter values of the PEE and CE force-length relationships and the CE force-velocity relationship, but impacted the estimated excitation dynamics parameter values. The improved method, in turn, yielded SEE stiffness values close to the actual ones and substantially enhanced the estimates of the other parameter values.
Our sensitivity analysis revealed substantial interdependencies between certain parameters. Particularly, we observed high sensitivity of estimated parameter values to variations in the maximal isometric CE force (\(F_{CE}^{max}\)) and the SEE slack length (\(L_{SEE}^0\)) parameter values. Despite this sensitivity, it is important to note that these parameters were accurately estimated even from perturbed experimental data. This was true in general: contraction and excitation dynamics parameter values were quite robust against shifts in the MTC force-length relationship of the simulated data according to the Monte Carlo simulations. These findings indicate the robustness of the improved method.
We applied both the traditional and the improved method to in situ data from quick-release, step-ramp, and isometric experiments on rat m. gastrocnemius medialis. The improved method yielded substantially higher estimates of SEE stiffness compared to the traditional method (67% on average). Using parameter sets from both methods, we simulated the quick-release, step-ramp, and isometric experiments, as well as independent stretch-shortening cycles, all with experimentally obtained MTC length and CE stimulation over time as input to the model. With parameter estimates obtained with the traditional method, the root mean squared difference between predicted and experimentally measured force was about 5.6% of maximal CE force, averaged across all rats and experiments; with parameter estimates obtained with the improved method this was only 3.9%. The small difference between predicted and experimentally observed force over time with parameters from the improved method suggests that muscle force can be accurately predicted with a Hill-type MTC model across a wide variety of contractions. This is a remarkable finding for two reasons. First, the Hill-type MTC model is obviously a simplification of real muscle, abstracting the muscle belly, aponeurosis, and tendon into the CE, SEE, and PEE components, and thus does not account for certain complexities (e.g., non-constant pennation angle, muscle inhomogeneities, etc.). Second, the Hill-type MTC model used in this study consists of 26 parameters, of which only 10 values were estimated. This means that all other parameters were left suboptimal. Despite these simplifications, the model was able to accurately predict muscle behaviour, demonstrating that a relatively simple model with a limited number of estimated parameter values can still capture important aspects of muscle dynamics.
In conclusion, our results demonstrate that the improved parameter estimation method provides accurate and robust estimates of contraction and excitation dynamics properties, offering a reliable approach for muscle property estimation in both biomechanical and muscle physiology research. The approach that we designed in this study is made available as an open-source toolbox (https://github.com/edwinreuvers/mp-estimator).
Acknowledgments
The authors thank Maarten F. Bobbert for helpful comments on a draft of the paper. The authors also thank Koen K. Lemaire for insightful discussions.
Funding
This work was funded by The Dutch Research Council (NWO) [21728 to D.A.K.].
Data and resource availability
All data, code, and materials used in this study are openly available:
GitHub repository: All raw data, processed data, and analysis code are hosted on GitHub at https://github.com/edwinreuvers/acc-mp-estimation.
Reproducible analysis website: Full analysis pipeline — including data, analysis code, and figure/table generation — is available at https://edwinreuvers.github.io/publications/acc-mp-estimation.
Glossary
| Symbol | Description |
|---|---|
| \(a\) | Hill constant |
| \(b\) | Hill constant |
| \(b_{min}^{scale}\) | Minimal scale factor of \(b\) (i.e., the minimal value of \(b_{scale}\)) |
| \(b_{shape}\) | Determines the steepness of the relation between \(b_{scale}\) and \(q\) |
| \(F_{asymp}\) | Oblique asymptote of the eccentric part of the CE force-velocity relationship |
| \(F_{CE}^{max}\) | Maximum isometric CE force |
| \(k_{PEE}\) | Scales the PEE stiffness |
| \(k_{SEE}\) | Scales the SEE stiffness |
| \(L_{CE}^{opt}\) | CE length at which CE can produce maximal isometric force |
| \(L_{PEE}^0\) | PEE slack length |
| \(L_{SEE}^0\) | SEE slack length |
| \(r_{as}\) | Slope of the slanted asymptote of the eccentric part of the CE force-velocity relationship |
| \(r_{slope}\) | Ratio between the eccentric and concentric derivatives of \(\frac{dF_{CE}^{rel}}{dv_{CE}^{rel}}\) @ \(v_{CE}=0\) |
| \(w\) | Determines the width of the CE force-length relationship |
| \(a_{act}\) | Determines the steepness of the relation between \(\gamma\) and \(q\) |
| \(b_{act,1}\) | Determines together with \(b_{act,1}\) and \(b_{act,3}\) the \(pCa^{2+}\) level at which \(q=0.5\) |
| \(b_{act,2}\) | Determines together with \(b_{act,2}\) and \(b_{act,3}\) the \(pCa^{2+}\) level at which \(q=0.5\) |
| \(b_{act,3}\) | Determines together with \(b_{act,1}\) and \(b_{act,2}\) the \(pCa^{2+}\) level at which \(q=0.5\) |
| \(k_{Ca}\) | Relates \(\gamma\) to the actual \(Ca^{2+}\) concentration |
| \(\gamma_0\) | Minimal value of \(\gamma\) |
| \(q_0\) | Minimal value \(q\) |
| \(\tau_{act}\) | Activation time constant of the calcium dynamics |
| \(\tau_{deact}\) | Deactivation time constant of the calcium dynamics |
References
Aubert, X., M. L. Roquet, and J. Van der Elst. 1951. “The Tension-Length Diagram of the Frog’s Sartorius Muscle.” Archives Internationales de Physiologie 59 (2): 239–41. https://doi.org/10.3109/13813455109145002.
Blix, Magnus. 1892. “Die Länge Und Die Spannung Des Muskels.” Skandinavisches Archiv Für Physiologie 3 (1): 295–318. https://doi.org/10.1111/j.1748-1716.1892.tb00660.x.
———. 1893. “Die Länge Und Die Spannung Des Muskels.” Skandinavisches Archiv Für Physiologie 4 (1): 399–409. https://doi.org/10.1111/j.1748-1716.1893.tb00749.x.
Blümel, Marcus, Scott L. Hooper, Christoph Guschlbauerc, William E. White, and Ansgar Büschges. 2012. “Determining All Parameters Necessary to Build Hill-Type Muscle Models from Experiments on Single Muscles.” Biological Cybernetics 106 (10): 543–58. https://doi.org/10.1007/s00422-012-0531-5.
Bobbert, Maarten F., Gertjan C. Ettema, and Peter A. Huijing. 1990. “The Force-Length Relationship of a Muscle-Tendon Complex: Experimental Results and Model Calculations.” European Journal of Applied Physiology and Occupational Physiology 61 (3): 323–29. https://doi.org/10.1007/BF00357621.
Bortolotto, Susan K., Maria Cellini, D. George Stephenson, and Gabriela M. M. Stephenson. 2000. “MHC Isoform Composition and Ca2+- or Sr2+-Activation Properties of Rat Skeletal Muscle Fibers.” American Journal of Physiology-Cell Physiology 279 (5): C1564–77. https://doi.org/10.1152/ajpcell.2000.279.5.C1564.
Burkholder, Thomas J., and Richard L. Lieber. 2001. “Sarcomere Length Operating Range of Vertebrate Muscles During Movement.” Journal of Experimental Biology 204 (9): 1529–36. https://doi.org/10.1242/jeb.204.9.1529.
Cecchi, G, F Colomo, and V Lombardi. 1978. “Force-Velocity Relation in Normal and Nitrate-Treated Frog Single Muscle Fibres During Rise of Tension in an Isometric Tetanus.” The Journal of Physiology 285 (1): 257–73. https://doi.org/10.1113/jphysiol.1978.sp012570.
Ebashi, S., and M. Endo. 1968. “Calcium and Muscle Contraction.” Progress in Biophysics and Molecular Biology 18 (January): 123–83. https://doi.org/10.1016/0079-6107(68)90023-0.
Gordon, A. M., A. F. Huxley, and F. J. Julian. 1966. “The Variation in Isometric Tension with Sarcomere Length in Vertebrate Muscle Fibres.” The Journal of Physiology 184 (1): 170–92. https://doi.org/10.1113/jphysiol.1966.sp007909.
Hatze, Herbert. 1981. Myocybernetic Control Models of Skeletal Muscle. Pretoria: University of South Africa.
Hill, A. V. 1925. “Length of Muscle, and the Heat and Tension Developed in an Isometric Contraction.” The Journal of Physiology 60 (4): 237–63. https://doi.org/10.1113/jphysiol.1925.sp002242.
Hill, Archibald Vivian. 1938. “The Heat of Shortening and the Dynamic Constants of Muscle.” Proceedings of the Royal Society of London. B. Biological Sciences 126 (843): 136–95. https://doi.org/10.1098/rspb.1938.0050.
———. 1950. “The Series Elastic Componet of Muscle.” Proceedings of the Royal Society of London. B. Biological Sciences 137 (887): 273–80. https://doi.org/10.1098/rspb.1950.0035.
Katz, Bernhard. 1939. “The Relation Between Force and Speed in Muscular Contraction.” The Journal of Physiology 96 (1): 45–64. https://doi.org/10.1113/jphysiol.1939.sp003756.
Kistemaker, Dinant A., Robert M. Terwiel, Edwin D. H. M. Reuvers, and Maarten F. Bobbert. 2023. “Limiting Radial Pedal Forces Greatly Reduces Maximal Power Output and Efficiency in Sprint Cycling; an Optimal Control Study.” Journal of Applied Physiology 134 (4): 980–91. https://doi.org/10.1152/japplphysiol.00733.2021.
Kistemaker, Dinant A., Arthur (Knoek) J. Van Soest, and Maarten F. Bobbert. 2005. “Length-Dependent [Ca2+] Sensitivity Adds Stiffness to Muscle.” Journal of Biomechanics 38 (9): 1816–21. https://doi.org/10.1016/j.jbiomech.2004.08.025.
Lemaire, Koen K., Guus C. Baan, Richard T. Jaspers, and A. J. ‘Knoek’van Soest. 2016. “Comparison of the Validity of Hill and Huxley Muscle–Tendon Complex Models Using Experimental Data Obtained from Rat m. Soleus in Situ.” Journal of Experimental Biology 219 (14): 2228–28. https://doi.org/10.1242/jeb.144394.
Levin, A., and J. Wyman. 1927. “The Viscous Elastic Properties of Muscle.” Proceedings of the Royal Society of London. Series B, Containing Papers of a Biological Character 101 (709): 218–43. https://doi.org/10.1098/rspb.1927.0014.
Petrofsky, Jerrold Scott, and Chandler Allen Phillips. 1981. “The Influence of Temperature, Initial Length and Electrical Activity on the Force-Velocity Relationship of the Medial Gastrocnemius Muscle of the Cat.” Journal of Biomechanics 14 (5): 297–306. https://doi.org/10.1016/0021-9290(81)90039-7.
Rack, P. M. H., and D. R. Westbury. 1969. “The Effects of Length and Stimulus Rate on Tension in the Isometric Cat Soleus Muscle.” The Journal of Physiology 204 (2): 443–60. https://doi.org/10.1113/jphysiol.1969.sp008923.
Reuvers, Edwin D. H. M, Huub Maas, Wendy Noort, Maarten F. Bobbert, and Dinant A. Kistemaker. 2025. “Maximising Average Mechanical Power Output During Stretch-Shortening Cycles of Rat Medial Gastrocnemius Muscle.” biorXiv. https://doi.org/10.1101/2025.09.25.678481.
Rijkelijkhuizen, J. M., C. J. de Ruiter, P. A. Huijing, and A. de Haan. 2003. “Force/Velocity Curves of Fast Oxidative and Fast Glycolytic Parts of Rat Medial Gastrocnemius Muscle Vary for Concentric but Not Eccentric Activity.” Pflügers Archiv 446 (4): 497–503. https://doi.org/10.1007/s00424-003-1052-9.
Soest, Arthur J. van, and Maarten F. Bobbert. 1993. “The Contribution of Muscle Properties in the Control of Explosive Movements.” Biological Cybernetics 69 (3): 195–204. https://doi.org/10.1007/BF00198959.
Stephenson, D. G., and D. A. Williams. 1982. “Effects of Sarcomere Length on the Force—pCa Relation in Fast- and Slow-Twitch Skinned Muscle Fibres from the Rat.” The Journal of Physiology 333 (1): 637–53. https://doi.org/10.1113/jphysiol.1982.sp014473.
Stern, Jack T. 1974. “Computer Modelling of Gross Muscle Dynamics.” Journal of Biomechanics 7 (5): 411–28. https://doi.org/10.1016/0021-9290(74)90004-9.
Van Soest, A. J., Margit Gföhler, and L. J. Richard Casius. 2005. “Consequences of Ankle Joint Fixation on FES Cycling Power Output: A Simulation Study.” Medicine & Science in Sports & Exercise 37 (5). https://journals.lww.com/acsm-msse/fulltext/2005/05000/consequences_of_ankle_joint_fixation_on_fes.14.aspx.
Walker, Sheppard M., and G. Randolph Schrodt. 1974. “I Segment Lengths and Thin Filament Periods in Skeletal Muscle Fibers of the Rhesus Monkey and the Human.” The Anatomical Record 178 (1): 63–81. https://doi.org/10.1002/ar.1091780107.
Woittiez, R. D., P. A. Huijing, H. B. K. Boom, and R. H. Rozendal. 1984. “A Three-Dimensional Muscle Model: A Quantified Relation Between Form and Function of Skeletal Muscles.” Journal of Morphology 182 (1): 95–113. https://doi.org/10.1002/jmor.1051820107.
Zajac, F. E. 1989. “Muscle and Tendon: Properties, Models, Scaling, and Application to Biomechanics and Motor Control.” Critical Reviews in Biomedical Engineering 17 (4): 359–411.
Zandwijk, Jan Peter van, Guus C. Baan, Maarten F. Bobbert, and Peter A. Huijing. 1997. “Evaluation of a Self-Consistent Method for Calculating Muscle Parameters from a Set of Isokinetic Releases.” Biological Cybernetics 77 (4): 277–81. https://doi.org/10.1007/s004220050388.
Supplementary material
4.1 Supplemental figures
