## Abstract

The purpose of the study was to expand a model of motor unit recruitment and rate coding (30) to simulate the adjustments that occur during a fatiguing contraction. The major new components of the model were the introduction of time-varying parameters for motor unit twitch force, recruitment, discharge rate, and discharge variability, and a control algorithm that estimates the net excitation needed by the motoneuron pool to maintain a prescribed target force. The fatigue-induced changes in motor unit activity in the expanded model are a function of changes in the metabolite concentrations that were computed with a compartment model of the intra- and extracellular spaces. The model was validated by comparing the simulation results with data available from the literature and experimentally recorded in the present study during isometric contractions of the first dorsal interosseus muscle. The output of the model was able to replicate a number of experimental findings, including the time to task failure for a range of target forces, the changes in motor unit discharge rates, the skewness and kurtosis of the interspike interval distributions, discharge variability, and the discharge characteristics of newly recruited motor units. The model output provides an integrative perspective of the adjustments during fatiguing contractions that are difficult to measure experimentally.

- muscle fatigue
- task failure
- first dorsal interosseus
- metabolites

the force exerted by a muscle during a voluntary contraction depends on the extent to which its motor units have been activated. The net force is attributable to concurrent changes in the number of activated motor units and the rates at which they discharge action potentials. The accumulated knowledge on motor unit function enabled Fuglevand et al. (30) to develop a widely used mathematical model of motor unit recruitment and rate coding (3, 41, 46, 56, 64). The model simulates the series of events from the excitation of the motoneurons to the forces exerted at the tendon. The model provides investigators with the opportunity to manipulate selected motor unit properties and to observe the influence on the surface electromyogram (EMG) and net muscle force. Although the model has proven extremely useful in the examination of the EMG-force relation (46, 73), twitch interpolation (41), and force variability (3, 56), it can only be used to simulate brief isometric contractions and not other experimental conditions.

A characteristic feature of physiological models is their ability to synthesize the knowledge acquired from a broad range of experimental conditions, to identify gaps in knowledge, and to determine the relative sensitivity of the output to the input parameters. For example, the original version of the Fuglevand model has triggered research and subsequent refinements in the model to address the influence of motor unit synchronization (72), discharge-rate variability (56), and age-related differences in discharge-rate characteristics (3).

The focus of the present work was to develop a model that can be used to study the adjustments in motor unit activity during fatiguing contractions. Models have been developed to predict the decline in the force-generating capacity of the muscle after exposure to specific electrical stimulation protocols (22, 37, 49) and to characterize populations of motor units in either a state of rest, activity, or fatigue, and to estimate the total muscle force based on the number of active motor units in each state (50, 71). None of these models, however, has simulated the interactions between motoneurons and muscle fibers to predict the change in the force-generating capacity of a muscle during a fatiguing contraction.

The purpose of the study was to expand the Fuglevand model to simulate the adjustments that occur during a fatiguing contraction by introducing time-varying parameters for motor unit twitch force, recruitment, discharge rate, and discharge variability, and a control algorithm that estimates the net excitation needed by the motoneuron pool to maintain a prescribed target force. The model was used to compare the simulation results with data available from the literature and experimentally recorded in this study and to provide examples of simulated physiological parameters related to muscle fatigue that cannot be obtained experimentally.

## METHODS

The motor unit force and motoneuron model of Fuglevand et al. (30) was extended to model the adjustments in motor unit recruitment and rate coding that occur during sustained isometric contractions. The modeled muscle was the first dorsal interosseus (FDI; main index finger abductor), although the parameters can be adjusted to simulate other muscles. The key additions were modulation of motoneuron excitation by a metabolite model of changes in the intra- and extracellular milieu that influenced afferent feedback and drove the adjustments in motor unit activity during the fatiguing contraction. Moreover, a proportional integral derivative (PID) control system was implemented to estimate the excitation level required to maintain a prescribed force profile as the model parameters varied over time (Fig. 1). The model was designed to run at 1,000 Hz.

#### Metabolite model.

Several metabolic substances are known to impact muscle function (1, 2, 12, 28, 62, 69), both individually and via interactions between the metabolites (2, 12). Although metabolite concentrations (MCs) can change in either direction during a sustained contraction, the overall change involves an increase in concentration with respect to baseline conditions. The concentration of K^{+}, for example, decreases in the intracellular environment and increases in the extracellular space (62). Furthermore, the accumulation of some metabolites can influence the synaptic input received by spinal neurons; for example, changes in lactic acid (59) and co-cyclooxygenase (40) can modulate the activity of group III–IV afferents.

The strategy in the present work was to develop a simplified model that represented the average change in MC and could be linked to experimental observations during fatiguing contractions. The model consists of a compartment for the intracellular space of each motor unit (*n* = 120 for the simulations of FDI) and one compartment for the extracellular space (Fig. 2). An assumption with this approach is that the concentration of a substance in each compartment is evenly distributed throughout the entire volume of the compartment instantaneously. Although each motor unit covers only a limited area in the muscle (9), the assumption seems reasonable, given that motor units with similar number of fibers are evenly distributed in the muscle and such motor units produce approximately the same amount of metabolites (9, 38).

The volume of the intracellular compartments (*Eq. 1*) was proportional to the innervation number of the motor units (adjustment of Eq. 13 from Ref. 30):
*i* indicates the motor unit number, V is the volume of the motor unit, V_{O} is the volume offset, which was set to 3, and F_{120}/F_{1} is the ratio of the force for the largest and smallest motor units, which was set to 80. These values produce a range of twitch forces similar to those observed in FDI (15) and in other hand muscles (13, 55); twitch force is assumed to be proportional to motor unit innervation number.

The FDI muscle comprises ∼40,500 fibers (26) and has a muscle cross-sectional area of 1.5 cm^{2} (43) and an average fiber radius of 28 μm^{2} (31). From these data, the extracellular space for this muscle was estimated to be about one-third the volume as the muscle fibers occupied 67% (*Eq. 2*):

The concentrations of metabolites in the intracellular compartment (*Eq. 3*) and extracellular space (*Eq. 4*), which were calculated in 500-ms epochs, were determined with equations commonly used in compartmental models (16):
_{i} is the MC in the *i*th motor unit, MC_{es} is the MC in the extracellular space, and V is volume (arbitrary units) for each compartment. The term MD represents the metabolite diffusion between compartments (positive diffusion implies diffusion from motor unit compartments to the extracellular compartment and vice versa for negative diffusion), MR is the removal rate by blood flow, and MP is the intracellular metabolite production. *Equations 3* and *4* represent the accumulation of metabolites based on the assumptions described in the next few paragraphs.

Metabolite production by each motor unit (*Eq. 5*) was set proportional to the number of action potentials discharged by the motoneuron (1). Thus intracellular metabolite production was related to motor unit volume and the number of discharges (ND) during each epoch:

The transport of metabolites in (*Eq. 6*) and out (*Eq. 7*) of the compartments was modeled as functions of the MCs in the relevant compartments, with the metabolite diffusion also depending on the motor unit volume representing the surface area over which metabolites could diffuse:
*Eq. 6*) should be greater than RC, and that both the intracellular and extracellular MCs should be reduced by ∼90% after 3 min of rest after the contraction. These criteria are met only by assigning DC = 0.01 and RC = 0.04.

The progressive occlusion of blood flow during fatiguing contractions augments the accumulation of metabolites in the extracellular space (18, 20, 60, 63). Blood flow was modeled as a function of the intramuscular pressure (IMP; *Eq. 8*) (20, 63), which depends on both its instantaneous value (IMP_{ins}) and its history (ΔIMP):
_{ins} (*Eq. 9*) denotes a linear relation between force and IMP (18, 60, 63), and force is the average force (%MVC force) produced by all active motor units in the epoch:
*Eq. 9* were determined by linear curve-fitting based on experimental measurements (8, 18, 61).

As contractions sustained at ≤10% MVC do not increase IMP (8), but contractions above ∼15% MVC can result in task failure (8, 52, 65), the increase in IMP was set to occur only during sustained contractions above 15% MVC force, and IMP was set to decrease when force declined below ∼15% MVC force. Complete occlusion occurs at ∼35% MVC force (20, 60, 63), and experimental observations indicate that IMP does not increase at forces above 35–40% MVC force (8). As there are few quantitative data on the rate of change of IMP with contraction level, the parameters that describe the shape of the function ΔIMP were chosen to fit the above observations. The amplitude of the function (defined by the gain of the first hyperbolic tangent functions: 0.12) was determined by matching the simulated rate of change in IMP at 25% MVC to that observed by Crenshaw et al. (18). On this basis, ΔIMP was characterized with the following equation:

The reduction in blood flow during a sustained submaximal contraction was assumed to be sigmoidal, with full occlusion at ∼40 mmHg when force was ∼35% MVC (60). This was achieved by setting the pressure at half occlusion (HO) to 30 mmHg. The following equation met the experimental observations for the sigmoidal decrease in flow with an increase in IMP (*Eq. 11*), where BF attains values between 0 (total blood occlusion) and 1 (no blood occlusion):

#### Estimation of the descending drive.

An algorithm based on the simulated tracking force error (PID) was used to estimated the required descending excitatory drive to the motoneuron pool for a given force profile. The error function, *e*(*t*) (*Eq. 12*), corresponded to the difference between the predefined target force (F_{t}) and the force produced by the model (F_{m}):
*t*); *Eq. 13*]. This was obtained as a linear combination of the error (proportional term), the time course of the error (integral term), and the rate of change of the error (derivative term) (PID), and was updated at a frequency of 3 Hz:
*K*_{p}, *K*_{i}, and *K*_{d} are the gains for the proportional and integral terms, respectively. The values of the gains determine the ability of the model to adjust the simulated force and are thereby responsible for the force fluctuations that arise during fatiguing contractions. The gains were established to produce an initial coefficient of variation for force of ∼3% (21, 45), with most of the power (>80%) in the 0- to 4-Hz band (21, 45) and a coefficient of variation for force of ∼8% at task failure (53). To meet these criteria, *K*_{i} was set to 2.5 × 10^{−4} and K_{p} was set to 2 × 10^{−6}. To achieve the experimentally observed increase in the coefficient of variation for force over time (53), *K*_{d} was not fixed, but increased during the contraction, which impaired the ability of the control algorithm to make the required adjustments. The rate of increase of *K*_{d} (*Eq. 14*) was related to the mean MC (MMC) of all motor unit compartments. The physiological source for this impairment is unknown, but it may be related to a fatigue-induced increase in the sensitivity of the muscle to stretch (17). Due to this possibility, the increase in *K*_{d} was implemented as a function of the mean intracellular MC (*Eq. 14*). The mathematical function describing this relation (exponential) and its parameters were determined by curve-fitting based on the force variability observed experimentally at 20% MVC (53):
*Eq. 15*) was calculated as:

#### Motor unit discharge rate.

The discharge rates of motor units recruited from the onset of a sustained contraction decline as the contraction progresses, whereas those recruited during the contraction tend to increase discharge rate (6, 15, 19, 25, 35, 47, 57, 58). The instantaneous discharge rate (DR) of each motor unit (*Eq. 16*) was computed as follows (adapted from *Eq. 2* in Ref. 30):
*Eq. 17*), which was determined by the sum of the descending drive (DD), the excitatory afferent input (EAI), and the inhibitory afferent input (IAI):
*Eq. 18*; Eq. 1 in Ref. 30) varied across the motor unit pool:
*a* = ln(30)/*n*(30). The minimum discharge rate (MDR; *Eq. 19*) of each motor unit was defined as:
_{1} is the minimum discharge of the smallest motor unit, which was set to 6.2 pulses/s (pps), and MDRD is the difference in MDR between the smallest and the largest motor unit and was set at 10 pps (3). The peak discharge rate (PDR) was defined similarly, with PDR_{1} = 15.6 pps and PDRD = 15.6 pps (3). The excitation level required to produce the PDR (E_{max}; *Eq. 20*; Eq. 8 in Ref. 30) was defined as:
*Eqs. 21* and *22*) received by each motoneuron had to be maximal (E_{max}) to reach its PDR (30):

As the DD (*Eq. 23*) accounts for ∼70% of the instantaneous discharge rate (IDR) at all levels of excitability (51), the remaining 30% was assumed to be produced by the EAI (*Eq. 24*):
*k* (*Eq. 25*) reduces the EI so the resulting IDR is reduced to 70%:
*k* and substituting the EI at the two levels of IDR with *Eq. 2* in Fuglevand et al. (30), *k* (*Eq. 26*) was expressed as:

The decline in discharge rate appears to be correlated with the initial discharge rate of the motor unit and thus to motor unit size (19, 35). In contractions sustained to task failure, the discharge rate of high-threshold motor units can decrease by ∼50% (6, 19), whereas that for low-threshold motor units declines by ∼35% (25). According to these experimental observations, it was assumed that the decline in EAI and the increase in IAI could maximally cause a reduction in discharge rate of ∼35 and ∼50% for the low- and high-threshold motor units, respectively. Although changes in the intrinsic excitability of the motoneurons were not explicitly included in the model, the net effect is incorporated in the inhibitory inputs and EIs.

The change in EAI was based on the change in the amplitude of the Hoffmann reflex (HR) during fatiguing contractions. Duchateau et al. (23) found that the amplitude of the H-reflex exhibited an exponential decline of ∼40% relative to the initial amplitude at all contraction levels, whereas the rate of decline depended on the contraction level. Based on these results, the EAI (ΔHR) during a sustained contraction {when force [F(*t*)] > 0} was mathematically described as an exponential function with a time constant of 3 × 10^{6}/F(*t*); the value was chosen to match the rates of decline observed by Duchateau et al. (22) at 25 and 50% MVC force. The rate of change during the decline was expressed as the first derivative of the exponential function (*Eq. 27*). In recovery [when F(*t*) = 0], the EAI returned linearly to the initial value (equal to an increase of 40%) within 5 min (3 × 10^{5} ms) (23), which is reflected in *Eq. 28*:

The initial value for ΔHR was set to 1, as the initial amplitude of the H reflex was not depressed, and the EAI (*Eq. 29*) estimated from *Eq. 24* was adjusted according to the sum of all of the changes in H-reflex amplitude:

The sum of all ΔHR was limited to an interval between 0.6 and 1, which meant that the reduction in EAI could maximally account for 12% of the decline in net synaptic input, that is, a decline in 40% in the EAI (23) that constitutes 30% of the net synaptic input (51). The remaining reduction in discharge rate was ascribed to an increase in IAI (*Eq. 30*), which, therefore, could reach a maximum (Inh max) of ∼23 and ∼38% for low- and high-threshold motor units, respectively:
*Eq. 2* in Fuglevand et al. (30), the desired magnitude of the decline in discharge rate due to IAI described by *Eq. 30* could be converted into values representing excitation for each motor unit (Exc inh max; *Eq. 31*):

Given that a major contributor to the inhibitory input is the chemically sensitive group III-IV afferents (7, 33), the magnitude of IAI was defined by the MC_{es} at each point in time, which assumes that the interstitium is the site at which metabolic accumulation influences these afferents. The sigmoidal relation for the magnitude of the inhibitory input (IIM; *Eq. 32*), which ranged from 0 (no inhibitory input) to 1 (maximal inhibitory input), was obtained as:
_{es} at one-half of maximal IIM. The value selected for HM and the sigmoidal shape for the relation were determined by matching the rate of the decline in discharge rate during an MVC to that observed experimentally (6). MC_{es} was calculated with *Eq. 4*.

As Exc inh max represents the residual excitation level for each motor unit at maximal inhibition, the difference between E_{max} and this value represented the maximal inhibitory input. Based on these relations, the IAI (*Eq. 33*) was calculated as:

#### Discharge rate variability.

Variability in the interval between successive action potentials (interspike interval) is due to fluctuations in the membrane potential during the afterhyperpolarization phase of the motoneuron (14). In the absence of fatigue, the interspike intervals are most variable when the net EI to the motor unit is close to its recruitment threshold, and the variability decreases when the net synaptic input increases above this level (54). The variability increases during sustained contractions, presumably due to the increase in synaptic input received by the motoneuron pool (5, 25, 35, 57, 58).

The variability in interspike interval was implemented by adding synaptic noise (SN; *Eq. 35*) to the NSI (*Eq. 34*) received by the motoneuron pool, as described in *Eq. 18*:
*N*(0, 1/6) is a normally distributed stochastic variable with zero mean and standard deviation 1/6. These parameters were chosen so that the coefficient of variation for interspike interval was ∼10% when the excitation level for a motoneuron was above its recruitment threshold, and ∼30% when it was close to its recruitment threshold in the nonfatigued state (56). As the amount of synaptic noise (*Eq. 36*) is related to the total synaptic input (4), the contributions of DD, EAI, and IAI to the synaptic noise was related to their magnitude:

#### Force.

The amplitude and time course of the force produced by each motor unit can change during sustained contractions (29, 67). These changes were implemented by introducing separate gain functions for twitch contraction time (*T*) and twitch amplitude (*P*) in the equation describing the motor unit twitch force (*Eq. 37*; from Ref. 30):
*Eq. 38*): gain_{max} and Δgain_{mc}. The term gain_{max} (*Eq. 39*) described the maximal value for the gain for each motor unit, thereby ensuring that relaxation time varied across motor units, as experimentally observed (29, 67). The term Δgain (*Eq. 40*) corresponded to the time-varying term of *Eq. 39* and ranged from 0 to 1, depending on the MC, thereby defining the magnitude of the gain at each instant in time. The changes in contraction time during a sustained contraction were assumed to be sigmoidal, as described in (10).
_{ref} (set to 1,150) indicates the MC at which the relaxation time and the tetanic amplitude reached the value indicated in Fig. 3. The value of MC_{ref} was determined by matching the simulated times to task failure across contraction levels to those observed experimentally (see Fig. 5).

As T gain (*Eq. 38*) also influences twitch amplitude, a correction factor (CF; *Eqs. 41* and *42*) was necessary to eliminate this effect:

The changes achieved with the two gain functions were modified to reflect the fatigue-induced changes in the tetanic force observed experimentally. The amplitude of the decline in tetanic force was represented with a sigmoidal function and to reach zero for all motor units (27), albeit at a faster rate for high-threshold motor units (29, 67). The adjustments in motor unit force (*Eq. 43*) were implemented with the function *h* (*Eq. 44*), which defined the MC at which the sigmoidal curve reached one-half of its maximal decline. *Equation 44* enabled the decline in motor unit force at a given intracellular MC to be greater for high-threshold motor units and was implemented by adjusting the twitch forces:
*s*(*i*), *t*(*i*), and *u*(*i*) were fit to make the decrease in motor unit tetanic force match the distribution described in Fig. 3 and are expressed as:
*Eqs. 38* and *43*) on the force produced by three motor units when activated to discharge 13 action potentials at 40 pps in a 330-ms interval is depicted in Fig. 4, which corresponds to a classic protocol used to quantify the fatigability of motor units (29, 67).

#### Model output.

The output of the model was compared with published data in the literature and with the results of an experiment. The purpose of the experiment was to determine the change in the coefficient of variation for the force exerted by the index finger and the frequency content of force during submaximal contractions sustained to task failure. Eight healthy men (age, mean ± SD, 23.4 ± 4.3 yr) volunteered to participate in three experimental sessions, which were separated by at least 24 h. The subjects provided informed consent before participating in the study, and the local ethics committee approved the procedures.

The subjects' left forearm was secured in a brace with the wrist joint in a neutral position and the palm facing downwards. The middle and ring fingers were secured with padded straps, and the index finger was free to rotate about its metacarpophalangeal joint. A force transducer (Interface MB-25) was aligned with the most distal phalanx of the index finger and used to measure abduction force exerted by the index finger. The force signal was amplified and displayed as visual feedback for the subject (Miso II, LISiN-OT Bioelettronica, Torino, Italy). Each session comprised the measurement of the MVC force (three trials separated by 2 s of rest) and one contraction of the FDI sustained at a target force of 30, 50, or 70% MVC force until task failure. The order of the contraction levels for the three sessions was randomly assigned for each subject.

The coefficient of variation for force as well as the frequency content of the force fluctuations in three bands (0–4, 4–8, 8–12 Hz) were calculated in 10, 7, or 4 epochs (for the 30, 50, and 70% MVC contractions, respectively). The values of the first and final epoch of these parameters were compared using two-way ANOVAs with the factors time (beginning and end) and contraction force.

## RESULTS

The results comprise the adjustments exhibited by selected model components during a simulated fatiguing contraction. The similarities between the function of the model and some of the classic experimental findings on fatiguing contractions, especially experimental observations on FDI and other hand muscles, were compared when possible. The model parameters were adjusted in some simulations to match experimental observations, whereas other simulations represent model predictions that were not obtained by tuning model parameters, as specified in each case. Table 1 summarizes the model parameter values used for all simulations.

#### Time to task failure.

Figure 5 depicts the time to task failure obtained with the model for target forces in the range 20–70% MVC. The simulated times to task failure were fit to the experimental data by tuning the parameter MC_{ref} in *Eqs. 40* and *41*. In both the simulations and the experiment, task failure was defined as the instant when the force dropped below 90% of the target force for 2–3 s. This criterion was adopted from one of the previously published studies used for comparison in Fig. 5 (32) and is similar to the criterion of a deviation of the exerted force from the target force by 1.5% of the target force for 3 s (52). Although a less stringent criterion was used in the other study to which the present data were compared (52, 74), the endurance time for the contraction sustained at 50% MVC force was similar to that observed in the present experiment.

The time to task failure for the model declined with an increase in target force. The simulated data were compared with that reported in the literature for sustained isometric contractions of the FDI (32, 52, 74), as well as those obtained in the present experiment (Fig. 5). In addition, the influence of central fatigue (65) on the time to task failure was estimated by reducing the maximal DD by 35%, which is a value that has been reported for FDI (24). The reduction in maximal DD by 35% had only a minor effect on the time to task failure; it declined by 12–18 s across the simulated target forces (dashed line in Fig. 5).

The simulated relation between time to task failure (ttf) and target force in the absence of the reduction in maximal DD (*Eq. 45*) was approximated with an exponential function (root mean square error: 16.1):

Time to task failure was moderately influenced by changes in twitch characteristics and afferent feedback in the model. When the rates of change in twitch amplitude and contraction time due to fatigue were varied in the model by changing the parameter MC_{ref} within ±25%, the simulated time to task failure at 20% MVC force was reduced from 532 s at MC_{ref} to 407 s at −25% and prolonged to 682 s at +25%. These values are within the confidence intervals for experimental studies (32, 52). Similarly, the rates of change of the EAI and IAI and the rate of change of IMP had only a minor influence on the time to failure. Variation in the rate of change of the IAI (by means of the parameter HM; *Eq. 32*) by ±25%, for example, reduced time to failure to 517 s at −25% and prolonged it by 547 s at +25%. Furthermore, variations in the volume of the extracellular compartment also had only a modest influence on the time to task failure; an increase in the volume by 25% prolonged the time to 573 s, whereas a decrease in the volume by 25% reduced the time to 498 s at 20% MVC force.

#### Force variability.

The amplitude and frequency content of the simulated force fluctuations changed with the duration of the sustained contraction. The coefficient of variation for force during a simulated contraction sustained at 30% MVC force began at 3.6% and increased to 7.3% at task failure (Fig. 6). The magnitude of the force variability was tuned to match the experimental results of Maton and Gamet (53) (*Eq. 14*), and not those presented in Fig. 6. The coefficient of variation for force measured at the three contraction levels in the experiment ranged from 3.4 to 4.0% at the beginning of the contraction and increased to final values of 11.9 to 15.0% at task failure. There was a significant difference from the start to task failure for all contraction levels in the experimental data (*P* < 0.001), and no significant difference across the three contraction levels (*P* = 0.40). The coefficient of variation for the simulated force increased only for contractions below 40% MVC force, whereas it remained at ∼4% throughout the contraction for target forces >40% MVC force.

The relative power (%total power) in the lowest frequency band (0–4 Hz) during the experiment ranged from 62.5 to 65.9% at the beginning of the contraction and from 78.3 to 86.4% at the end (Fig. 6). There was less power in the other two bands, which ranged from 3.6–13.6 to 4.4–8.0% for the 4- to 8-Hz band and from 6.1–7.0 to 2.4–3.5% for the 8- to 12-Hz band at the beginning and end of the contraction, respectively. There was a significant difference in relative power from the start to task failure for all bands and contraction levels in the experimental data (*P* < 0.001) and no significant difference among the three contraction levels. The general trends in these experimental values were similar to those obtained with the model, especially for the coefficient of variation and the power in the low-frequency band. The simulations produced the following ranges for relative powers: 57.2–67.3% (beginning) and 59.5–85.6% (task failure) for the 0- to 4-Hz band, 10.7–16.4 and 5.4–12.5% for the 4- to 8-Hz band, and 3.7–13.6 and 3.8–8.6% for the 8- to 12-Hz band. These distributions of energy represent independent predictions of the model.

Varying the rate of change in the derivate gain as a function of the MC (*Eq. 14*) by ±25% did influence the coefficient of variation for force. Compared with the default condition in which the coefficient of variation for force increased to 7.3% at task failure, a decrease in the rate of change by 25% during a contraction sustained at 30% MVC force produced a final coefficient of variation for force of 4.0%, whereas an increase in the rate of change by 25% increased the value to 8.5%.

#### Contribution from individual motor units to the muscle force.

The degree to which each motor unit contributed to the force produced by the muscle varied with fatigue, as smaller motor units are more resistant to fatigue. Fig. 7 depicts the cumulative contribution of the motor units to the muscle force in three conditions: during an MVC without fatigue, at the task failure for a contraction sustained at 60% MVC force, and at task failure for a contraction sustained at 20% MVC force. The 100% value on the *y*-axis of Fig. 7, therefore, corresponds to muscle forces of 100, 60, and 20% MVC force. There was a small difference between the contribution of motor units to the total force for the MVC and at task failure for the 60% MVC contraction. In contrast, the distribution at task failure differed markedly for the 20% MVC contraction, when the 80 smallest motor units contributed for ∼50% of the total force compared with ∼30% in the two other conditions.

#### MC and force loss.

As an example of the capacity of the model to replicate experimental findings, Fig. 8 depicts the relation between the total MC in all compartments and the decline in force during a simulated MVC. The simulated data were compared with the increase in H_{2}PO_{4}^{−} concentration reported by Cady et al. (11) during three 15-s MVCs. As the simulated MC was reported in arbitrary units, the H_{2}PO_{4}^{−} concentration was normalized to the simulated MC at maximal force loss (normalization ratio: 1:25); thus the comparison is relevant only for relative changes. The association between the decline in MVC force and the simulated increase in MC exhibited a similar trend as the experimental measurement of the increase in H_{2}PO_{4}^{−} concentration. As the accumulation of metabolites was not tuned to these data, this result represents an independent prediction of the model.

#### Discharge rate.

Figure 9 depicts the change in mean discharge rate for three motor units during simulated ramp-and-hold contractions (3-s slopes, 10-s hold, 4-s rest) for comparison with the experimental recordings of Carpentier et al. (15). The ramp-and-hold contractions were repeated at a rate of 3/min until task failure (∼700 s). The low- and medium-threshold motor units exhibited an approximately exponential decline in discharge rate. The discharge rate of the high-threshold motor unit (no. 112) is shown for the time interval when it began to discharge action potentials repetitively. The mean discharge rate of *motor units 1–96* decreased during the simulated contraction, whereas it increased for *motor units 97–120*. *Motor unit 107* was the first among those not active at the beginning of the contraction to be recruited, which occurred at 220 s after the start of the contraction. As the modulation in discharge rate was not imposed in the simulations, the results represent an independent prediction of the model.

#### Discharge rate variability.

Figure 10 shows the mean discharge variability for selected groups of motor units during a simulated contraction at 35% MVC force sustained to task failure. Each value for the coefficient of variation was calculated from at least 50 interspike intervals. Most of the motor units recruited from the beginning of the contraction (nos. 1–100) showed a slow increase in the coefficient of variation for interspike interval throughout the contraction. Three motor units (nos. 101–103) recruited at the onset of the contraction, however, had more variable values during the first half of the simulated contraction because they were activated at levels just above recruitment threshold. At task failure, the coefficient of variation for all motor units converged to the 10–15% range. The initial variability in discharge rate was tuned to match experimental data (*Eq. 35*).

Figure 11 shows the skewness (*A*) and kurtosis (*B*) for the distributions of interspike intervals for selected groups of motor units during simulated contractions sustained at 60% MVC force to task failure (95 s). The interspike intervals for motor units with excitation levels close to recruitment threshold were more positively skewed and had a smaller kurtosis than those with excitation levels above recruitment threshold. The dominant feature of the distributions was the marked kurtosis, which indicated that the interspike intervals were more similar to the central value than might be expected based on a normal distribution. Furthermore, the skewness of the motor units with an initial excitation level close to recruitment threshold declined to a value similar to that for the smaller motor units as the DD increase near task failure. The magnitude of the simulated force did not influence the values of skewness and kurtosis. The model parameters were not tuned for skewness and kurtosis of the interspike interval distributions, which thus represent an independent prediction of the model.

#### Recruitment threshold.

Figure 12 indicates the relation between the motor unit recruitment thresholds before and after a simulated fatiguing contraction in which the force was sustained at either 25 or 50% MVC force to task failure. Time to failure was 397 s for the lower target force and 136 s for the greater target force. Except for motor units with very low recruitment thresholds, recruitment thresholds were greater before the fatiguing contractions, which is consistent with the values measured after the fatiguing contraction, providing an index of the change in the force contributed to the contraction by motor units with lower recruitment thresholds. Experimental results from Garland et al. (35) are shown for comparison. Figure 13 shows the percent change in recruitment threshold as a function of the initial value after sustaining a simulated contraction at 50% MVC force to failure; the experimental results of Carpentier et al. (15) are shown for comparison. The simulated recruitment thresholds were obtained without specific tuning of parameters in the model.

## DISCUSSION

The study describes the expansion of a model of motor unit recruitment and rate coding to simulate the adjustments that occur during sustained isometric contractions. The major new components of the model include a compartment model to account for the distribution and transport of metabolites, formal descriptions of the rates of change in motor unit activity during sustained contractions, and a variable activation of the motoneuron pool to sustain the target force. The model has been evaluated by comparison with experimental data in representative simulations.

#### PID algorithm.

In a model that includes time-varying changes in its parameters, it was necessary to include an algorithm that allows the tracking of a target force profile by continuously estimating the necessary level of excitation for the task, given the instantaneous values of the model parameters. This was accomplished with a PID controller. The algorithm, however, was not intended to represent the many peripheral and central processes by which the central nervous system (CNS) controls motor output (44). Rather, the aim was to identify the input required by the motoneurons to sustain a target force, given the other model parameters. The PID controller was not intended to represent the function of the CNS, but instead it estimated the descending input required by the motoneurons, given the set of muscle and afferent input parameters to achieve the target force. Due to the complexity of the model, it was not possible to determine the input needed by the motoneurons with a closed analytic formulation.

Despite the absence of explicit pathways to the motoneuron pool in the model, there is some correspondence between the control algorithm and the inputs received by the motoneurons. For example, the integral term of the PID algorithm reduces the low-frequency error in force tracking as occurs when visual feedback of target force is provided to the subject (4, 57). Similarly, the discharge rates of Golgi tendon organs and muscle spindles, which signal the instantaneous force and rate of change in length of a muscle (34), are quantified by the proportional and derivative terms of the PID algorithm. Furthermore, the likely increase in the gain of muscle spindles during fatiguing contractions (17) is represented by the increase in the derivate gain in the model. Despite these associations, the PID controller in the model, nonetheless, only estimates the input needed to perform the task, which is the only variable that can be controlled in this system, and it cannot be derived analytically.

#### Metabolite model.

Several fatigue-induced changes in motor unit activity are driven by the simulated MC. Although a number of metabolites can influence muscle function (1, 2, 12, 28, 62, 69), the exact influence of each compound is unknown, and some are interrelated (2, 12). As there are insufficient experimental data in the literature on the influence of the various metabolites and their interactions on specific muscle properties, modeling the role of individual metabolites is not currently possible. As an alternative approach, a single metabolite parameter was chosen and tuned to have similar recovery dynamics to those indirectly observed experimentally (70). It was not possible with the current data to describe more complex dynamics at an integrative level; however, the model can be revised as more data become available. As the simulated concentration refers to a compound parameter, its change during the simulated contractions cannot be validated directly, although it exhibited similar dynamics as H_{2}PO_{4}^{−} during an MVC (11). Instead, the validity of the metabolite model can be inferred by its ability to elicit adjustments in neuromuscular activity to match experimental findings. Accordingly, a number of simulations were performed to confirm that the output of the expanded model could match experimental data on time to task failure and changes in the recruitment thresholds and discharge rates of motor units.

The compartment structure of the metabolite model made it possible to describe changes in properties of quiescent motor units by diffusion of metabolites in the extracellular space, consistent with the findings of Kostyukov et al. (48) and Gazzoni et al. (36). Removal of the metabolites was modeled as a function of the extracellular concentration, a parameter describing the RC and blood flow. The RC was set at 4% of the metabolites in the extracellular space within a period of 500 ms. Blood flow was assigned a normalized value between 0 and 1 that depended on the level of the IMP. It is recognized, however, that the modeling of blood flow was a simplification of the actual dynamics, as it is the pressure gradient between the increases in blood pressure and IMP that influences blood flow. Little is known about the pressure gradient during sustained contractions; thus the modeling of blood flow includes some assumptions that are not directly based on experimental data. As mentioned in the Introduction, however, the model design provides a basis for identifying the experimental data needed to increase the knowledge of the complex interactions that occur during fatiguing contractions.

#### Force variability.

The simulated increase in force fluctuations was achieved by increasing the derivative gain in the PID control algorithm throughout the contraction and thereby enhancing the tendency for the control algorithm to overcompensate for error in the force. The derivative gain allows slow changes in the error and thereby favors changes in the low-frequency content of the force signal. When the values of the three gains in the controller were held constant, the force fluctuations decreased or remained stable during simulations at all contraction levels (results not shown), contrary to experimental observations (35, 53, 57, 58). Despite not modeling the actual physiological controller, the model predicts that the task of maintaining a constant force becomes easier (less variability) with an increase in the number of motor unit discharges and a lengthening of the twitches. Nonetheless, force fluctuations do actually increase in experimental conditions, despite the more favorable conditions for maintaining constant force, which indicates that the ability of the CNS to maintain a constant force declines over time. Accordingly, it has been shown that reaction time increases and accuracy decreases in performing dynamic tasks after a fatiguing protocol (42, 66), which suggests impaired CNS control in muscle fatigue. The physiological adjustments responsible for the decrease in force control are presumably associated with the change in the derivative gain of the simplified controller.

The degree of variability in the simulated force and its frequency content corresponded to the data reported in the literature, both in the absence of fatigue (21, 45, 68) and during fatiguing contractions (Fig. 6) at low target forces (≤40% MVC). In contrast with experimental data, however, the oscillations in simulated force did not increase with target force. The simulation results are explained by the observation that the metabolite accumulation at task failure is higher at low-contraction forces than at higher forces and by the relation between the gain of the PID control algorithm and the metabolite accumulation. As the derivative gain of the PID control algorithm is exponentially related to the MMC (*Eq. 14*), the ability of the algorithm to control and to sustain a constant force is compromised most at low forces. This discrepancy may indicate that the reduction in the ability to control force during fatiguing contractions is not exclusively related to metabolite accumulation, which is another example of the performance of the model and provides predictions that need to be examined experimentally. Another possible source for the decline in force control is signal-dependent noise in the DD to the motoneuron pool (39).

#### Model fitting and predictions.

Many parameters of the model, but not all of them, were tuned to match the experimental data described in the results. For example, the relation between the rate of change in the inhibitory afferent feedback and the MC (*Eq. 32*) was chosen to describe the experimentally observed rates of change in motor unit discharge rates (6, 15). However, no model parameters were imposed to induce the recruitment of additional motor units and their characteristic initial increase in discharge rate, as shown, for example, by Carpentier et al. (15) (Fig. 9). Similarly, the rates of change in the force capacity of each motor unit (*Eqs. 38* and *43*) were adjusted so that the simulated time to task failure was within the range of experimentally observed values (Fig. 5). Changes to a number of other parameters within their physiological ranges (±25%), however, did not cause the simulated times to task failure to exceed the confidence intervals reported in the literature (see results). Additionally, the parameters in the function that defined the magnitude of the synaptic noise (*Eq. 35*) were chosen so that the variability in discharge times was similar to experimental values (25, 35, 53, 56–58); nonetheless, the interspike intervals exhibited similar characteristics to those observed by Matthews (54) (Fig. 11), despite no parameters in the model being specified to promote these trends. Finally, changes in the recruitment thresholds of the motor units during the fatiguing contractions, which reflect the combined influence of several model functions, displayed similar trends to those reported by Garland et al. (35) and Carpentier et al. (15) (Figs. 12 and 13).

As the model was able to predict several known characteristics of muscle fatigue, the results suggest that the model captured at least some of the major physiological adjustments that occur during fatiguing. Nonetheless, the association does not constitute a validation of the model, as the simulations were limited to representative conditions. Rather, the study proposes a model that integrates current knowledge on motor unit physiology during fatiguing and identifies critical gaps in this knowledge.

#### Application of the model.

The model provides a tool for extracting information about parameters that cannot be measured experimentally. It is not possible, for example, to measure the force contributed by each motor unit to the net muscle force, but this was possible with the model (Fig. 7). Furthermore, a sensitivity analysis of the model parameters can identify experimental studies that are needed to improve our knowledge on the physiology of muscle fatigue. For example, the current simulation results suggest that variability across subjects in the time to task failure is likely due to differences in the rate of change in the force capacity of muscle fibers, as relatively large changes in the rates of change of the inhibitory afferent feedback and IMP did not have significant influences on the time to task failure. The simulations also showed that a significant reduction in the DD had a relatively minor effect on the time to task failure (Fig. 5), which was due to the exponential increase in excitatory drive (as estimated by *Eq. 14*) at all contraction forces up to maximal value at task failure.

An important feature of physiological models is their promotion of research to fill the gaps in the mathematical descriptions of physiological functions. For example, the current model underscores the need for more work on identifying the factors that are responsible for the decline in force control during fatiguing contractions: the model predicts that the task is easier when the peripheral properties change, as in fatigue conditions, which suggests that a central impairment is responsible for the decline in performance. Similarly, more data are needed on metabolite accumulation and its influence on muscle performance, so that the compound metabolite in the current version of the model can be replaced by the dynamics of actual physiological substances that are responsible for the adjustments.

The model is designed explicitly to reflect the adjustments during fatiguing contractions performed with the FDI muscle. Model parameters that would vary across muscles include the number of motor units, the range of PDRs (Eq. 5 in Ref. 30), the innervations numbers of the motor units (*Eq. 1*), the number and proportion of muscle fiber types (reflected in *Eqs. 38* and *43*) (46), and the muscle force at which blood flow is occluded (*Eq. 11*) (20). Changes in these parameters will surely alter the simulated results.

In conclusion, the present study describes an expanded model of motor unit recruitment and rate coding (30) that includes the relevant central and peripheral adjustments during fatiguing contractions. Many of the adjustments were implemented as a function of a compound metabolite model and a PID control algorithm that provided the net excitation needed to achieve the prescribed target force. The model provided values for time to task failure and motor unit characteristics comparable to experimental findings over a broad range of target forces.

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the author(s).

- Copyright © 2010 the American Physiological Society

## REFERENCES

- 1.↵
- 2.↵
- 3.↵
- 4.↵
- 5.↵
- 6.↵
- 7.↵
- 8.↵
- 9.↵
- 10.↵
- 11.↵
- 12.↵
- 13.↵
- 14.↵
- 15.↵
- 16.↵
- 17.↵
- 18.↵
- 19.↵
- 20.↵
- 21.↵
- 22.↵
- 23.↵
- 24.↵
- 25.↵
- 26.↵
- 27.↵
- 28.↵
- 29.↵
- 30.↵
- 31.↵
- 32.↵
- 33.↵
- 34.↵
- 35.↵
- 36.↵
- 37.↵
- 38.↵
- 39.↵
- 40.↵
- 41.↵
- 42.↵
- 43.↵
- 44.↵
- 45.↵
- 46.↵
- 47.↵
- 48.↵
- 49.↵
- 50.↵
- 51.↵
- 52.↵
- 53.↵
- 54.↵
- 55.↵
- 56.↵
- 57.↵
- 58.↵
- 59.↵
- 60.↵
- 61.↵
- 62.↵
- 63.↵
- 64.↵
- 65.↵
- 66.↵
- 67.↵
- 68.↵
- 69.↵
- 70.↵
- 71.↵
- 72.↵
- 73.↵
- 74.↵