Koen K. Lemaire

,

Dinant A. Kistemaker

,

Richard T. Jaspers

,

Willem J. van der Laarse,

A. J. ‘Knoek’ van Soest

Willem J. van der Laarse Conceptualization, Formal analysis, Methodology, Resources, Writing – review & editing

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.

When interpreting important features of the human locomotor system, it is often hypothesized that its design (e.g. Lichtwark and Wilson, 2007; Rodman and McHenry, 1980) and concomitant behaviour (e.g. Selinger et al., 2015) are optimized to minimize metabolic energy expenditure. In addition to direct experimentation, hypotheses regarding the role of metabolic energy expenditure in (the control of) human (loco)motion can be tested with the aid of musculoskeletal models (Koelewijn et al., 2019; van den Bogert et al., 2012; Miller et al., 2012; Umberger, 2010; Kistemaker et al., 2010). Typically, musculoskeletal models employed to test such hypotheses contain Hill-type muscle models to represent the actuators. Owing to the original work by A.V. Hill in 1938, metabolic power is usually modelled as the sum of mechanical power and several heat rate terms, which are related to muscle activation, generation of force and shortening/lengthening of the muscle (Lichtwark and Wilson, 2005; Bhargava et al., 2004; Umberger et al., 2003). In this phenomenological approach, these heat terms are parameterized using data obtained from dedicated experiments aimed at a single heat term. The downside of this approach is that it is unclear whether the results of these isolated experiments are representative of contractions in vivo, during which all terms vary simultaneously (Curtin and Barclay, 2023). As such, any number of heat terms may be required to describe a given behaviour. This is a conceptual problem, because the essential feature of muscle tissue is that it can convert metabolic energy into mechanical energy, and that both mechanical behaviour and metabolic energy expenditure result from the same mechanochemical reactions (Smith et al., 2005; Herzog and Ait-Haddou, 2002; Barclay and Curtin, 2023). Incorporating knowledge of the contractile process into the formulation of the muscle model can thus improve the generalizability and structural validity of the model.

The reactions that govern the conversion of metabolic energy into mechanical energy take place during the cross-bridge cycle (Herzog and Ait-Haddou, 2002; Barclay et al., 2010; Huxley, 2000a), and as such cross-bridge cycling determines a large part of the energy required for muscle contraction. In addition to cross-bridge cycling, considerable energy is required in the muscle for the (de)activation process, mainly through the active pumping of intra-cellular Ca²⁺ back into the sarcoplasmic reticulum (Barclay et al., 2007). A mathematical model of the cross-bridge cycling process was first formulated by Huxley (1957). In his scheme, the contractile process is modelled as a time-evolving distribution of the relative fraction of attached cross-bridges over their bond length. In the Huxley model, both muscle force and cross-bridge energy consumption follow from this distribution. In a previous study, we extended the original two-state Huxley model to incorporate parallel and series elasticity, length dependence of contractile element force and activation dynamics, and we showed that this model was capable of capturing the mechanical behaviour of rat soleus muscles during dynamic contractions with constantly varying length and stimulation (Lemaire et al., 2016). Following up on that study, we showed that it is numerically feasible to employ this model in task optimization of a largish-scale musculoskeletal model, and that actuation by muscles represented by Huxley-type cross-bridge models yields mechanical behaviour that is very similar to that obtained using Hill-type muscle models (van Soest et al., 2019). These previous results opened up the possibility of investigating how well metabolic energy expenditure by muscle can be predicted by a Huxley-type muscle–tendon complex (MTC) model. To that end, in an additional study we collected data on the energy consumption of mouse soleus muscle fibre bundles that were subjected to periodic contractions (Lemaire et al., 2019). In that study, the metabolic energy expenditure of the muscle fibre bundles was measured while stimulating the bundles during either the concentric or the eccentric phase. In addition, the metabolic energy expenditure related to cross-bridge cycling was separated from that related to the activation process, by the application of blebbistatin, which blocks the acto-myosin interaction while leaving all other processes intact. In the current study, we investigated the extent to which the previously collected data on muscle mechanical behaviour and metabolic energy expenditure in Lemaire et al. (2019) can be predicted by the Huxley-type MTC model described in Lemaire et al. (2016).

The Huxley-type MTC model used in the current study was similar to that presented in Lemaire et al. (2016). The experimental data used in this study were previously described in Lemaire et al. (2019). In addition, we used unpublished data that were collected in the same study, consisting of isometric and sinusoidal contractions at various muscle lengths, movement amplitudes and stimulus phases. During these contractions, no oxygen was measured. Parameters governing the mechanical behaviour of the Huxley model were fitted to a subset of the previously collected data, which only concerned measurement of mechanical behaviour. The model predictions of the mechanical behaviour were evaluated by a cross-validation procedure, in which parameters estimated on the data of one animal were used to simulate data from the other animals. Scaling parameters that relate the metabolic energy expenditure associated with cross-bridge cycling and with muscle activation were fitted to the data pertaining to muscle energetics, and this fit was evaluated in terms of its root-mean-square error (RMSE).

In our previous study on the mechanics of rat soleus muscle (Lemaire et al., 2016), the shapes of the rate functions f(x) and g(x) were chosen to result in a good description of the data whilst being computationally inexpensive. In that study, it was possible to choose rate functions arbitrarily, because a given steady-state force–velocity curve can be adequately fitted with various shapes of rate functions. However, even if the force–velocity relationship is similar for two sets of rate functions, the resulting metabolic energy expenditure may not be similar, because force is related to the number of attached cross-bridges and their bond length, whereas metabolic power is related to the rate of cross-bridge unbinding. In the present study, we chose the shapes for the rate functions to be similar to that originally proposed by Huxley (1957), as these resulted in an adequate description of the energetics of frog muscle. However, the original model was not designed to adequately describe the mechanical behaviour of muscle during eccentric contractions, and has been shown to greatly overestimate the force during eccentric contractions (Zahalak, 1981). Therefore, similar to Zahalak (1981), we added an extra detachment rate parameter g₃, in order to capture the mechanics of eccentric contraction. Furthermore, in the Huxley model, the detachment rate parameter g₁ determines the metabolic energy expenditure during a purely isometric contraction. Because we did not have a separate measurement of muscle metabolic energy expenditure during isometric contraction, we fixed this rate parameter to the arbitrary value of 100 Hz. For numerical reasons, the original rate functions were smoothed by a Gaussian sigmoid over a 0.2h wide section at x=0 and x=h. The resulting rate functions are shown in Fig. 1B. Note that in the current study, we normalized the bond length x with the maximum bond length h. The values of the rate parameters reported here pertain to this normalized bond length.

Our formulation of the Huxley MTC model has l_MTC(t) and stim(t) as its inputs, and γ and n as its states, with Eqns 6 and 7 as the corresponding state equations. In order to avoid having to numerically find the root of the force balance Eqn 2 at each time step of a simulation, l_CE was added as a quasi-state variable. As explained in detail in Lemaire et al. (2016), the corresponding state equation can be derived by solving the time derivative of Eqn 2 for ⁠.

The model calibration and evaluation carried out in this study is based on analysis of experimental data on mechanics and energetics of periodic contractions of mouse soleus muscle described in (Lemaire et al., 2019). Below, we will briefly explain the experimental set-up and measurement protocol used in that study.

All experimental procedures were performed in accordance with the regulations of the Animal Experimentation Ethical committee of the Vrije Universiteit Amsterdam (protocol number: DEC FBW1101) and were in accordance with Dutch law. For a complete description of the experimental setup and the experimental conditions, the reader is referred to Lemaire et al. (2019). In brief, 9 C57BL/6J mice were anaesthetized and both soleus muscles were removed from the hindlimbs. A bundle of approximately 70 fibres including their tendon fibres was dissected from one of the muscles and was moved to a jacketed glass chamber containing Tyrode solution equilibrated with a gas mixture containing 95% oxygen and 5% carbon dioxide. The water inside the water jacket was maintained at 32°C. This choice of temperature reflects a trade-off between ecological validity and practical considerations, where we stayed close to the expected physiological temperature, but did not compromise the stability of the preparation. The proximal end of the bundle was attached to the bottom of the chamber. The distal end was attached to a tungsten wire, which left the chamber at the top via a thin capillary and was attached to a servomotor via a force transducer (AE801, SensoNor, Horten, Norway). A custom-made, polarographic oxygen sensor protruded into the chamber through a tightly fitting opening. Through a separate opening, fresh Tyrode solution could be pumped into the chamber. The solution inside the chamber was circulated by a ferritic stainless steel stirrer, which was driven by a magnetic field generated at the bottom of the chamber. The poles of an isolated stimulus generator were attached to stimulus leads placed at the bottom and at the top of the chamber, which were in contact with the Tyrode solution inside. The stimulus isolator generated alternating pulses of 0.4 ms width. Signals from the force transducer, motor position sensor, oxygen electrode and stimulus leads were A/D converted (12 bits) and sampled on a computer at 2 kHz.

For a detailed description of the trials pertaining to measurement of the bundle's metabolic energy expenditure, see Lemaire et al. (2019). In brief, during these trials, the flow of oxygenated Tyrode solution was temporarily stopped, and the bundles were subjected to sinusoidal movements with various amplitudes and frequencies (0.25 mm, 2 Hz; 0.50 mm, 2 Hz; and 0.25 mm, 3 Hz) for a duration of 2 min. Stimulation (5 pulses delivered at 100 Hz) was phased to elicit either concentric or eccentric contraction. For practical reasons, the imposed length changes were in absolute terms. Consequently, the length changes relative to fibre bundle optimum length differed between preparations. Expressed as a percentage of fibre optimal length, the length changes were 4.0±0.5% (range 3.4–4.7%) for the 0.25 mm amplitude condition (see Table 2 for contractile element optimum lengths). After all conditions involving force production were completed, blebbistatin was added to the Tyrode solution. Blebbistatin irreversibly disables the formation of cross-bridges but leaves all other processes unaffected (Kovács et al., 2004; Farman et al., 2008). After addition of blebbistatin, the trials during which oxygen consumption was measured were repeated. The bundle's metabolic energy expenditure in the trials after addition of blebbistatin yielded an estimate of the metabolic energy expenditure that is not associated with cross-bridge cycling. In this study, we assume that this energy expenditure can be ascribed to calcium pumping into the sarcoplasmic reticulum. Between the trials during which oxygen consumption was measured, oxygen content in the chamber was recovered by flushing the chamber with fresh, oxygenated Tyrode solution. During these recovery periods, short-duration contractions were imposed. These included a series of isometric contractions at different bundle lengths, used to characterize the force–length relationship, and a series of isometric contractions at bundle optimum length, with different contraction frequencies (Lemaire et al., 2019). In addition to these isometric contractions, sinusoidal contractions were imposed with various amplitudes (ranging from 0.1 to 0.5 mm), movement frequencies (ranging from 0.5 to 4 Hz) and stimulus trains (ranging from 4 to 8 pulses, always delivered at 100 Hz). Phase of the stimulation with respect to the movement was varied to elicit either concentric or eccentric contraction. Data obtained during the isometric and dynamic trials were used to fit model parameters pertaining to the mechanical behaviour of the Huxley-type MTC model.

In Lemaire et al. (2016), we performed a series of dedicated experiments aimed at estimating the values of the parameters pertaining to the mechanical behaviour of each model element separately. In Lemaire et al. (2016), it was impossible to perform a similar series of dedicated experiments because of the characteristics of the experimental setup. Instead, we varied the amplitude and frequency of the imposed movement to encompass a wide range of contraction conditions. In the current study, all free parameters pertaining to static properties (i.e. SEE, PEE and CE force–length relationships) were fitted to the isometric trials and all other free parameters were fitted simultaneously to the combined data of all dynamic trials. The values of the parameters pertaining to the energetic behaviour were fitted after the parameter values describing the mechanical part were fitted.

For the (static) force–length relationships of CE, SEE and PEE (Eqns 3 and 4), CE-normalized, SEE strain and PEE slack length were taken to be invariant over animals and their values were based on previous literature reports (see Table 1). Maximum CE force, CE optimum length (constrained to be between 0.9 and 1.1 times the fibre length at rest; based on the force–length curve in Lemaire et al., 2019), SEE slack length and PEE shape parameter were taken to be animal specific and were considered free parameters. For each animal separately, these free parameters were fitted to the data of the force–length relationship by finding the values which minimized the sum of squared differences between the measured and the modelled active and passive force–length relationships, using a constrained optimization method based on gradient descent. For ease of interpretation, the root-mean-squared difference between the modelled and the measured forces is expressed as a percentage of the maximal isometric force. The assumptions regarding the PEE force–length properties of the bundles did not have a large influence on our results as passive forces were of the order of 1% of maximum active force.

The parameter values governing the metabolic energy expenditure (c₁–c₃ in Eqn 13) were estimated for each animal separately, according to the following procedure. For each trial in which metabolic energy expenditure was measured, a simulation was done only of the first contraction in that trial, using the parameter values obtained previously for that animal, and using the bundle length and stimulation as imposed during the experiment. From these simulations, the time-averaged values for the power terms associated with cross-bridge cycling and calcium pumping (see Eqn 13) were obtained. The experimentally observed force in each trial declined throughout the course of the trial. However, the shape of the workloops was similar between subsequent contractions within one trial (see fig. 8 in Lemaire et al., 2019). Because the steady-state value of n(x,t) scales with the relative isometric force, we modelled this within-trial decline in force as a change in the maximum CE force over time, and we scaled the integrals in Eqn 13 with the ratio of the average force obtained during the whole trial to the average force during the first contraction. In the simulations of the trials that were conducted after blebbistatin was added, the integrals in Eqn 11 were set to zero. Equating the modelled P_met,tot to the measured metabolic power results in one equation with three unknowns c₁–c₃, for each trial. Combining all trials per animal results in an overdetermined system of equations, which is linear in the unknown terms c₁–c₃. First, c₃ was solved separately, by only considering the trials that were conducted after blebbistatin was added. For the remaining set of equations, the now known term c₃·γ was subtracted from the measured metabolic energy expenditure and the remaining parameters c₁ and c₂ were solved in a least squares sense using the Moore–Penrose pseudo inverse.

The model was evaluated in terms of its prediction of force by a cross-validation procedure. For each animal, simulations of the data of all other 8 animals were made using the optimal parameter set obtained for this animal, and the difference between these simulation results and the experimental data was taken as a measure of the model's generalizability. For the metabolic energy expenditure, we compared the prediction error with the variation occurring within the experimental data.

Results are presented as means±s.d. of the sample, with N=9, unless mentioned otherwise. All cases in which a typical example is shown or referred to pertain to the same animal for which typical examples were also shown in Lemaire et al. (2019), unless mentioned otherwise.

For the fibre bundle of each animal, the force–length relationship could be adequately fitted. The normalized RMSE between the modelled and the measured forces was small (2.3±0.8%), and the fitted parameters were similar for all animals (Table 2). Because of the small between-animal variation in fibre optimum length, we do not expect that the between-animal differences in relative length changes imposed in the dynamic conditions had a substantial effect on other aspects of the model predictions or fits. A typical example of a modelled and measured force–length relationship, for both active and passive conditions, is shown in Fig. 2.

The dynamic trials could also be adequately described by the model, as evidenced by the small normalized RMSE between modelled and measured forces (5.1±1.6%, Table 3). Typical examples of simulation results for various contraction types and the corresponding data are shown in Fig. 3. We will discuss the modelled metabolic power traces in Fig. 3 in the next section. Note that the RMSE values of the force traces (range 2.7–4.2% of maximal isometric force) are similar to the mean values for all trials listed in Table 3. Given the similarity between modelled and measured force traces in Fig. 3, this indicates good correspondence between modelled and measured forces throughout. Nevertheless, two peculiarities can be noted. First, immediately after the initial stimulus pulse, the predicted forces were lower than the measured values (Fig. 3). This finding was consistent across animals. The experimentally measured force response to stimulation is highly non-linear in this regard; the initial pulse generates a far greater change in force than subsequent pulses. A potential way to incorporate this behaviour in future modelling efforts would be to make the activation dynamics dependent on the recent history regarding stimulation.

Table 3.

Fitted model parameters and RMSE pertaining to dynamic behaviour

f₁ (Hz)g₂ (Hz)g₃ (Hz)τ_act (s)τ_deact (s)kRMSE (%⁠)

3.4e+02	1.6e+04	5.3e+02	0.036	0.079	0.26	4.3
1.4e+03	1.9e+04	3.9e+03	0.090	0.053	0.15	5.5
5.2e+02	7.8e+04	4.8e+02	0.034	0.083	0.40	3.2
7.2e+02	1.3e+04	4.6e+02	0.028	0.070	0.34	6.2
1.2e+03	1.6e+03	1.3e+02	0.022	0.066	0.79	4.8
5.3e+02	1.4e+05	6.2e+02	0.024	0.083	0.41	8.7
7.8e+02	1.1e+04	4.3e+02	0.024	0.075	0.35	5.1
6.4e+02	1.0e+04	4.9e+02	0.033	0.064	0.35	4.1
6.7e+02	1.2e+04	4.0e+02	0.032	0.063	0.40	4.0

f₁ (Hz)g₂ (Hz)g₃ (Hz)τ_act (s)τ_deact (s)kRMSE (%⁠)

3.4e+02	1.6e+04	5.3e+02	0.036	0.079	0.26	4.3
1.4e+03	1.9e+04	3.9e+03	0.090	0.053	0.15	5.5
5.2e+02	7.8e+04	4.8e+02	0.034	0.083	0.40	3.2
7.2e+02	1.3e+04	4.6e+02	0.028	0.070	0.34	6.2
1.2e+03	1.6e+03	1.3e+02	0.022	0.066	0.79	4.8
5.3e+02	1.4e+05	6.2e+02	0.024	0.083	0.41	8.7
7.8e+02	1.1e+04	4.3e+02	0.024	0.075	0.35	5.1
6.4e+02	1.0e+04	4.9e+02	0.033	0.064	0.35	4.1
6.7e+02	1.2e+04	4.0e+02	0.032	0.063	0.40	4.0

Each row in the table corresponds to a single animal. The RMSE in the last column is the value of RMSE_rel in Eqn 14 for the dataset of each animal.

The second peculiarity concerns the variation in estimated parameter values between animals. Because the animals have the same genetic background, it might be expected that morphology-independent parameters such as time constants and rate constants would be similar among animals. However, the fitted parameter values differ substantially between animals (Table 3, albeit not by more than one order of magnitude). The time constants in the experimental data are expected to vary with fibre-type distribution, which varied in the experimental data. The bundles had an average of 37±12% type I fibres, ranging from 22% to 53% (Lemaire et al., 2019). Part of the variation in fitted parameter values may thus be related to variations in fibre-type distribution, which were not accounted for in the current model. In future work, introducing a fibre-type distribution dependency of the model's time and rate constants may be expected to reduce the variance in their values. This is desirable, as ideally one (dimensionless) set of parameter values suffices to adequately model the behaviour of each animal.

Contrary to our previous publication (Lemaire et al., 2016), we did not have an independent set of data to validate the mechanical behaviour of the model in the current study. To nonetheless explore the generalizability and to cross-validate the predictions of the mechanical behaviour, we performed simulations of the experimental data pertaining to all animals, using the optimal parameter values fitted to each animal's experimental data. This resulted in 9×9 comparisons of simulated and measured data, for each of which the cost RMSE_rel (Eqn 14) was computed. The results of these comparisons in terms of c are listed in Table 4. In Table 4, the values in each row pertain to the set of parameter values fitted to the data for a particular animal (corresponding to the parameters listed in each row of Table 3), and the values in each column pertain to simulations of the experimental trials of a particular animal. For example, the value in row 2, column 3 in Table 4 indicates the value of RMSE_rel computed from simulations of the data pertaining to animal 3, using parameter values fitted to the data of animal 2. The values on the diagonal are thus identical to the values of RMSE_rel listed in Table 3. The cost values in Table 4 can be interpreted as a measure of the quality of the description of the data. The mean normalized RMSE of each set of parameter values of the data of the other 8 animals was 11.4% (s.d. 2.8%). The model predictions are thus reasonably robust to changes in parameter values.

In Table 5, we show the coefficients c₁–c₃ of Eqn 13; the resulting relative RMSE between the predicted and the measured average metabolic power for the pre-blebbistatin trials; and the relative RSME value for the comparison of predicted and measured forces of these trials. This is shown for each animal in a separate row. As can be appreciated from the last column of Table 5, the relative RMSE of the force predictions is similarly low for these trials as for the data on which the mechanical behaviour was fitted. Moreover, the relative RMSE of the metabolic energy prediction appears to be unrelated to the relative RMSE of the force prediction. Regarding the relative RMSE of the energy predictions of the post-blebbistatin trials (on which only the value of c₃ was fitted, and not listed in Table 5), we note that that is completely determined by the within-animal variation of the empirical data. This is because in the model, we do not distinguish between concentric and eccentric post-blebbistatin trials, as we did not find evidence for such a difference previously (Lemaire et al., 2019). We thus used two post-blebbistatin trials per animal to fit c₃. The mean RMSE of the post-blebbistatin metabolic energy expenditure prediction was 27% (s.d. 14%) of the measured value, or 1.7 W kg⁻¹ (s.d. 0.7 W kg⁻¹). Thus, the relatively large RMSE reflects the between-trial, within-animal variation, rather than indicating a poor model fit.

For the pre-blebbistatin trials, the RMSE of the metabolic energy prediction differed substantially between animals (Table 5, column 4). To examine these differences more closely, we plotted the predicted metabolic energy expenditure against the measured values, separated by condition (symbols) and animal (colours; Fig. 4, large symbols indicate the across-animal average per condition, with the error bars indicating across-animal variation of the within-animal means). Note that the average fit for the single eccentric conditions was near perfect, which was expected because there was one free parameter associated with the metabolic energy expenditure during eccentric contractions. Further note that the animals for which the RSME values were large (higher than 15%; animals 3, 4, 5, 8 and 9) had a large within-animal variation in the empirical data. For instance, in animal 5 (pink), both the concentric 0.50 mm, 2 Hz (circle) and the 0.25 mm, 3 Hz (square) condition had a lower metabolic energy expenditure than the concentric 0.25 mm, 2 Hz condition (downward triangle). This result goes against expectations and against the mean trend in the experimental data. In animal 3 (bright green), the empirical results do agree with the mean trend in the data, but there is a large variation (horizontal distance) in repeated measures, which is not reflected by a change in model predictions, as the inputs to the model are the same for repeated trials. In other words, for a given condition and set of parameter values, the model will always yield the same result, whereas the experimental data might differ between successive trials. Any within-animal variation in the empirical trials will therefore directly appear in the RMS difference between model predictions and empirical data. In the animals for which the RMSE between the predicted and measured metabolic energy expenditure was small (lower than 15%; animals 1, 2, 6 and 7), the empirical data more closely followed expected trends (e.g. animal 2, red) and the variation in repeated measures was much smaller (e.g. animal 1, blue). Across animals, the RMSE was 20.3% (12.6%) or 4.3 W kg⁻¹ (2.6 W kg⁻¹) muscle. To qualitatively assess this prediction, it merits a comparison to the variation, expressed as the RMS difference from the mean value of the trials for which we have repeated measures. These were the 0.25 mm, 2 Hz concentric condition (N=7) and the post-blebbistatin (eccentric and concentric) (N=9) conditions. The RMS difference values for these conditions were 1.8 W kg⁻¹ (s.d. 1.7 W kg⁻¹) and 1.2 W kg⁻¹ (0.7 W kg⁻¹), respectively. This indicates that almost half of the differences between the model predictions and the experimental data might be explained by within-animal variation. In the Supplementary Materials and Methods and Fig. S1, we further investigate the cause of the experimental within-animal variation and conclude that this is at least not due to cell necrosis during the course of the experiment. All in all, given the substantial variation in the experimental data, we consider the model predictions of metabolic energy expenditure to be in reasonable agreement with the experimental results.

To explore the generalizability of these results to other contraction conditions, we simulated some of the trials that were not used in the fitting procedure, and qualitatively assess the simulation results below. In Fig. 3, examples of the predicted metabolic power are shown (bottom graph in each panel), for trials that were not included in the dataset used during parameter fitting of the model of metabolic energy expenditure. Note that no experimental data on the metabolic energy expenditure are shown here, as we did not measure oxygen consumption during these particular trials. Further note that for all contractions, the (de)activation energy expenditure is only dependent on the number and frequency of stimulus pulses. In the isometric contractions (Fig. 3A,B), total metabolic power scales with the force–time integral, as expected. Furthermore, the ratio between (de)activation and total metabolic work was in the expected range, for both low (37%) and high (30%) forces, even though no data from isometric contractions were available for fitting of the parameter values. For the concentric and eccentric contractions, the mechanical efficiency was 0.15 and −0.13, which is similar to previously reported values for mouse soleus muscle (Barclay and Curtin, 2023). These results are by no means conclusive, but they seem at least qualitatively to be in agreement with previous findings, and we thus expect that the predictions of metabolic energy expenditure will generalize well to other contraction conditions. This means that the current modelling effort is suitable for describing the energetics of muscle contraction.

As with the parameters pertaining to the mechanical behaviour, one would expect the parameter values c₁–c₃ to be similar between animals because of their similar genetic background. However, again, there is a fairly large variation in these parameter values (Table 5). This variation is due to both between-animal variation in the empirical data and between-animal variation in the parameters underlying the mechanical behaviour (see above). Similar to the mechanical behaviour, part of the variation in the empirical data may be due to variation in fibre-type distribution. Previously, we found that the percentage of type I fibres was correlated to the mechanical efficiency (Lemaire et al., 2019). Including a fibre-type dependency in the model is therefore expected to reduce the between-animal variation in the values of c₁–c₃. In the current study, we did not do so, because of the risk of over-fitting parameters within one animal. If a fibre-type dependency is added to the current model in future studies, we recommend using a single set of rate constants and muscle activation time constants across animals, where these parameters are scaled by the ratio of type I fibres in such a way that a higher percentage of type I fibres leads to slower dynamics. In this way, fibre-type dependency would affect both the mechanical behaviour and metabolic energy expenditure in an internally consistent manner.

In the general discussion below, we will suggest potential strategies to further improve the prediction of metabolic energy expenditure in this dataset.

The aim of the present study was to investigate how well the data on muscle metabolic energy expenditure presented in Lemaire et al. (2019) can be predicted by the Huxley-type muscle–tendon complex model described in Lemaire et al. (2016). Similar to our previous work, the parameter fitting resulted in good correspondence between predicted and measured forces of the fitted trials. Furthermore, the force predictions generalized well across animals. For the metabolic trials (the data of which were not used to fit the parameters pertaining to mechanical behaviour), the prediction of the force traces was also good. Overall, the prediction of metabolic energy expenditure was reasonably good, even though there were large between- and within-animal variations. Below, we discuss possibilities for improving upon the current study, and future directions.

Similar to Lemaire et al. (2016), we chose to split the parameter-fitting procedure into consecutive stages, as opposed to lumping all free parameters into one large-scale optimization. The reasoning behind this choice was that separate parts of the model can be parameterized in isolation using dedicated experiments. In general, partitioning the parameter estimation procedure is conceptually attractive because it allows us to test and evaluate contained parts of the model separately, and it results in a closer correspondence to observed (muscle) physiology. This may be expected to increase a model's generalizability to conditions other than the ones the parameters were fitted to.

As mentioned before, our method of parameter estimation proved successful in terms of predicting mechanical muscle behaviour. Yet, in contrast to the approach in Lemaire et al. (2016), both the parameters pertaining to the force–velocity relationship (rate parameters) and the parameters governing the activation dynamics were fitted simultaneously in the second stage, as we did not have data from dedicated force–velocity experiments. This introduced the possibility of redundancy in the parameter space, as the activation dynamics and the muscle contractile dynamics constitute two dynamical systems in series, for which different parameter sets can result in similar dynamic behaviour [i.e. similar moments of the distribution n(x,t)]. An example where this may have happened is in animal 2, which had a large value for f₁, combined with a large value for τ_act and a low value for κ (row 2 in Table 3). In terms of the overall dynamics of the model, the net effect of these parameter values may well be ‘average’ behaviour. Considering that the large variation in parameter values observed in this study did not lead to large deviations in the quality of the fit (Table 3), it indeed seems that there is redundancy in the parameter space. This redundancy may have had consequences for the prediction of metabolic energy expenditure, which depends on both the activation dynamics and the cross-bridge distribution (Eqn 13). It is thus conceivable that a different combination of rate parameter values and activation dynamics parameter values would have resulted in a least squares solution for the metabolic parameters that has a slightly lower prediction error, while not affecting the description of the mechanical behaviour to any meaningful extent. We chose not to pursue this further in the current study, as we took as a starting point the description of the mechanical behaviour, which can only get worse if an additional term is added to the cost function.

It should be noted that the two-state cross-bridge model as currently presented is a minimal version of the cross-bridge model, which closely follows the original paper by Huxley (1957). For one thing, the contribution of titin to the contractile process is not modelled. Titin is a molecular spring that plays a role predominantly during active eccentric contractions. Through calcium-mediated modulation of its stiffness, titin is hypothesized to explain, amongst others, residual (isometric) force enhancement after stretch (Herzog, 2018). Implementation of the mechanisms in which titin contributes to muscle contraction in a model requires the inclusion of multiple additional states, equations and corresponding parameters (e.g. Millard et al., 2024; Sampaio de Oliveira and Uchida, 2024). Implementing such changes in the current model would most likely improve the predictions of the experimental data, albeit at the cost of introducing additional free parameters and complexity to the model. In the current work, we did not explicitly model the role of myosin-binding protein-C (MyBP-C). MyBP-C is a family of accessory proteins of striated muscles that contributes to the assembly and stabilization of thick filaments, and regulates the formation of actomyosin cross-bridges, via interactions with both the myosin and actin filaments (Ackermann and Kontrogianni-Konstantopoulos, 2011; Ackermann et al., 2015). Three distinct MyBP-C isoforms have been characterized, in cardiac muscle and in slow and fast skeletal muscle. Detailed discussion on how to model the action of MyBP-C in muscle contraction are beyond the scope of this work. However, we do note that any model that incorporates fibre-type dependency in its dynamics would capture some of the effects of different MyBP-C isoforms.

In addition to data from studies on muscle structure, there are behavioural data that were not taken into account in the current modelling effort. There is strong evidence from temperature jump and rapid length change studies that the cross-bridge cycle includes at least one more weakly attached state (for reviews, see Huxley, 2000b; Ranatunga, 2010; Barclay and Curtin, 2023). Including such additional states may improve the overall fit of the data. However, we note that (i) this would again involve adding additional parameters with the associated risk of overfitting and (ii) the additional states were needed to describe the transients due to temperature jumps and rapid length changes; it is unclear whether these are necessary to describe the experiments considered here. Recently, Seow and Seow (2022) suggested that variation in the curvature of the steady-state force–velocity relationship can be interpreted as changes in the velocity dependence of the rate of cross-bridge detachment. In a similar vein, we could consider making the rate of cross-bridge detachment contraction velocity dependent in the current model. Also, we could consider improving upon the model of activation dynamics. Lastly, our assumptions regarding the metabolic energy requirement for the calcium pumps are rather simplistic, and not reflective of the current state of calcium pump models (e.g. Baylor and Hollingworth, 2007; Tran et al., 2009). Incorporating these more realistic models of calcium transport is likely to improve the general validity of the model. In the current study, this may not have had a large effect as the variation in contraction conditions was in terms of contraction velocity and activation phase, and not in terms of activation duration. It would be interesting to investigate the above modifications and additions to the model in future work.

In the above, we argued that different combinations of rate parameter values and activation dynamics parameter values may lead to similar mechanical behaviour but different metabolic energy expenditure. This argument can be extended to different shapes of rate functions, because there is both a many-to-one mapping of distributions n(x) to force and (a different) many-to-one mapping of distributions n(x) to metabolic power. Different rate function shapes may thus lead to similar mechanical behaviour but different metabolic power. If so, the latter could be exploited to further reduce the prediction error of the metabolic energy expenditure. As the amount of data for this study is limited, we did not explore these options further, to prevent over fitting.

Finally, given the relatively large variation in the experimental data, it would be helpful if more high-quality data on muscle energetics becomes available. In particular, data from experiments on mammalian muscle at the appropriate temperature and in which the contributions of calcium pumping and cross-bridge cycling to metabolic energy expenditure are separated would be valuable. Such data should be combined with available knowledge on the details of the muscle contraction process, to further develop and validate structural models of the mechanics and energetics of muscle contraction.

As far as the authors are aware, this study is the first attempt to evaluate a Huxley-type MTC model in terms of its simultaneous predictions of mechanical behaviour and metabolic energy expenditure. In our view, the current results are promising; directions for improvement of the analysis and model in future work are provided in the above. In combination with our previous work on the feasibility of implementing a Huxley-type MTC model in musculoskeletal modelling (van Soest et al., 2019), this study represents an important step towards better prediction of whole-body energy expenditure, using the Huxley cross-bridge model. If this proves successful, this will aid investigations into the role of energy expenditure in human motor control, and the design of assistive devices in sports and health applications.

© 2025. Published by The Company of Biologists

2025