6139

J Physiol 591.24 (2013) pp 6139–6156

Motor unit recruitment by size does not provide functional advantages for motor performance Jakob L. Dideriksen and Dario Farina Department of Neurorehabilitation Engineering, Bernstein Focus Neurotechnology G¨ottingen, Bernstein Center for Computational Neuroscience, University Medical Center G¨ottingen, Georg-August University, Von Siebold Str. 4, 37075, G¨ottingen, Germany

Key points • Motor unit recruitment according to size (orderly recruitment) is almost universally observed

in human muscles.

• Using genetic optimization methods on neuromuscular models, we explored the functional

The Journal of Physiology

basis for orderly recruitment.

• We show that the key features for optimal muscle performance are the number of motor

neurons innervating the muscle and the distribution of their innervation numbers, and that a random recruitment order is as effective as a recruitment ordered by size. • This implies that size-ordered recruitment of motor neurons universally arises due to its intrinsic linkage to innervation number and not for function optimization. • These findings have implications for our understanding of human neural control of movement.

Abstract It is commonly assumed that the orderly recruitment of motor units by size provides a functional advantage for the performance of movements compared with a random recruitment order. On the other hand, the excitability of a motor neuron depends on its size and this is intrinsically linked to its innervation number. A range of innervation numbers among motor neurons corresponds to a range of sizes and thus to a range of excitabilities ordered by size. Therefore, if the excitation drive is similar among motor neurons, the recruitment by size is inevitably due to the intrinsic properties of motor neurons and may not have arisen to meet functional demands. In this view, we tested the assumption that orderly recruitment is necessarily beneficial by determining if this type of recruitment produces optimal motor output. Using evolutionary algorithms and without any a priori assumptions, the parameters of neuromuscular models were optimized with respect to several criteria for motor performance. Interestingly, the optimized model parameters matched well known neuromuscular properties, but none of the optimization criteria determined a consistent recruitment order by size unless this was imposed by an association between motor neuron size and excitability. Further, when the association between size and excitability was imposed, the resultant model of recruitment did not improve the motor performance with respect to the absence of orderly recruitment. A consistent observation was that optimal solutions for a variety of criteria of motor performance always required a broad range of innervation numbers in the population of motor neurons, skewed towards the small values. These results indicate that orderly recruitment of motor units in itself does not provide substantial functional advantages for motor control. Rather, the reason for its near-universal presence in human movements is that motor functions are optimized by a broad range of innervation numbers. (Resubmitted 19 July 2013; accepted after revision 15 October 2013; first published online 21 October 2013) Corresponding author D. Farina, Medical University G¨ottingen, Von-Siebold-Str. 4,37075 G¨ottingen, Germany G¨ottingen 37075 Germany. Email: [email protected]

 C 2013 The Authors. The Journal of Physiology  C 2013 The Physiological Society

DOI: 10.1113/jphysiol.2013.262477

6140

J. L. Dideriksen and D. Farina

Introduction In 1957 Henneman reported a relation between the motor neuron soma size and its susceptibility to generate action potentials at a given input current (Henneman, 1957). This is predicted by Ohm’s law because the membrane resistance of a neuron depends on soma size, so that a current level induces a greater difference in membrane potential in the smallest neurons (Stein et al. 2005). Therefore, when an increasing current is delivered uniformly to a population of motor neurons, small neurons are recruited before larger ones. The magnitude of the difference in rheobase current between small and large neurons is sufficiently high for the orderly recruitment not to be disrupted by the variability that may occur in input currents across the motor neuron population in most natural conditions (Heckman & Binder, 1993). Therefore, orderly recruitment of motor units has been consistently verified (Cope & Pinter, 1995; Burke, 2002), with few exceptions (Heckman & Enoka, 2012). A relevant exception is that of respiratory muscles whose motor neurons receive a non-uniform synaptic input to accommodate specific mechanical demands (Butler & Gandevia, 2008). The difference in the number of muscle fibres innervated by individual motor neurons within hand muscles can attain an 80-fold range (Enoka & Fuglevand, 2001). Across muscles, such ranges are reflected in a similar or greater variability in the range of motor unit force twitch amplitudes (Monster & Chan, 1977; Calancie & Bawa, 1985). As the size of the soma of a motor neuron is associated to the diameter of its axon (Clamann & Henneman, 1976; Kernell & Zwaagstra, 1981) and the axon diameter determines the degree to which the axon divides into terminal branches (Henneman & Mendell, 1981; Moore & Appenteng, 1991), a large range of innervation numbers (number of muscle fibres innervated by a motor neuron) within a motor neuron population necessarily corresponds to a large range of sizes of motor neuron somas and vice versa. Evolutionary optimization of motor performance may have determined either a recruitment order by size or a broad range in the distribution of innervation numbers; each of these conditions would have forced the presence of the other characteristics. Therefore, the question arises whether the orderly recruitment has been determined by the actual functional need for a progressive recruitment by size or as a consequence of the need for a range in innervation numbers. In the latter case, the orderly recruitment may not per se optimize the motor output but it would be a consequence of an optimization of motor performance by the need of a large range in number of fibres innervated by each motor neuron in a pool. Such a question has never been posed. Rather, it has been classically assumed that recruitment by size corresponds to

J Physiol 591.24

some functional advantages in force production (Kandel et al. 2000; Heckman & Enoka, 2004; Stein et al. 2005). For example, it is an intuitive observation that recruitment by size optimizes the gradation of force by activating first the small motor units that provide small steps in force increase. However, the influence of different recruitment strategies on force generation may be counterintuitive due to the large size of the motor unit population, which is often in the range of hundreds (Heckman & Enoka, 2012). In this study, we challenge the classic view that progressive recruitment ordered by innervation number has functional advantages. For this purpose, we determine if specific optimization criteria in motor output give rise to an orderly recruitment or if such a recruitment order is a physiological consequence rather than a functional outcome. The study was based on parameter optimization of advanced computational models of the neuromuscular system (Fuglevand et al. 1993), based on a genetic algorithm. This algorithm is a biologically inspired optimization scheme implemented to mimic selective survival and mutation in the simulated motor unit populations. The approach consisted of running evolutionary optimizations of a number of parameters in models of motor unit populations and evaluating the generated motor output. In some of the evolutionary optimizations, no link was imposed between motor neuron recruitment threshold and the number of innervated muscle fibres, so that these two parameters could vary and be optimized independently. In these conditions, we determined whether the optimization of motor functions would determine a range of innervation numbers and a specific recruitment order, i.e. an association between recruitment threshold and motor neuron size. If the optimization of functional output parameters does not correspond to an orderly recruitment, but only to a range of innervation numbers, then the recruitment order would not correspond to a functional need but would naturally arise because of the association between innervation number and motor neuron size.

Methods Neuromuscular model

A computational model for motor neuron activity and generation of muscle force (Fuglevand et al. 1993) was adopted. In this model, motor unit recruitment and rate coding is determined by an excitatory drive uniformly delivered to the motor neuron pool, based on three neural and two muscular parameters assigned to each motor unit: recruitment threshold (in arbitrary units of excitation), minimum discharge rate, maximum discharge rate (both in pulses per second; pps), motor unit twitch amplitude (in  C 2013 The Authors. The Journal of Physiology  C 2013 The Physiological Society

Orderly recruitment and motor performance

J Physiol 591.24

arbitrary units of force) and motor unit twitch contraction time (in ms). The assignment of values to the motor unit parameters is described in the next section (Evolutionary algorithm). The neural parameters assigned to each motor unit in the model determine the relation between the input excitation level and its recruitment and discharge pattern, as follows: DRMU (t) = E MU (t) − RTE MU + MDRMU , E MU (t) > RTE MU

(1)

where DR denotes instantaneous discharge rate, E the instantaneous excitation level, RTE the recruitment threshold, MDR the minimum discharge rate and MU the motor unit identification number. E was modelled as the sum of the excitation offset and zero-mean Gaussian noise with a standard deviation scaled to get realistic motor unit discharge rate variability across all contraction levels (Dideriksen et al. 2010). The gain of the excitation rate coding relation (equation 2 in Fuglevand et al. 1993) was set to 1 for all motor units. The motor unit twitch force was modelled as a critically dampened second-order system with a gain dependent on the discharge rate. Based on the complete discharge pattern of all active motor units, the force generated by the whole muscle could be computed. In this way, the model characterized the behaviour of a population of motor units and the resulting motor output. Figure 1 shows graphical examples of model parameter distributions. Evolutionary algorithm

In the genetic (or evolutionary) algorithm, a population is a set of models, each with specific values of the parameters. A population of 50 models, corresponding to 50 sets of parameters that characterized each model, was defined. At the first iteration of the algorithm, for each model, all neuromuscular parameters were set to random values within broad ranges. Moreover, the parameters could vary in a full independent way with respect to each other, without any assigned association between them. The innervation number of each motor unit was expressed as a percentage of the total number of fibres in the muscle and its limits for each motor unit were from the percentage corresponding to one single innervated fibre to 100% (full possible range of conditions). In this way, a muscle could comprise, as limit conditions, one single motor unit innervating all its muscle fibres, or a number of units equal to the number of muscle fibres, each innervating a single fibre. The motor unit recruitment thresholds were assigned in the range 0.1–50 excitation units (eu) and the motor unit twitch contraction times in the range 20–140 ms. The ranges of minimum and maximum motor unit discharge rates were between 1–20 pps and 20–100 pps, without imposed associations  C 2013 The Authors. The Journal of Physiology  C 2013 The Physiological Society

6141

between these values. For example, valid parameters for a motor unit (or even for all motor units in a model) could have been, as limit conditions, a minimum and maximum rate both equal to 20 pps. Furthermore, the number of muscle fibres and the number of motor units in the muscle were randomly assigned in the ranges 1000–100,000 and 1–500, respectively. The peak twitch force of each muscle fibre was assumed to be the same for all fibres and was set to 1 mN. The twitch force for a motor unit was the sum of twitches of the individual muscle fibres. By using the broad ranges for parameter values reported above, we did not make any assumptions on the distributions or relations between any of the model parameters, which could vary broadly within large limits and independently of each other. Any model with parameters fixed within the ranges described was a valid member of the model population. Assigning specific values to all model parameters determined a specific member (model) of the population. The first population in the evolutionary optimization process was made of 50 models, each with a different random selection of the parameters. These models could generate muscle force with characteristics according to the chosen parameters and thus different for each member (model) of the population. The genetic optimization approach consisted in iteratively generating populations of models with ‘better characteristics’. The characteristics for optimizing models in a population had to be chosen according to known desired motor functions. Once these characteristics have been decided, the genetic approach generates populations, by selective survival of the fittest members, which progressively contain models that fit better the optimization criteria. The criteria for optimization were chosen based on well known important properties of force generation and are described in the following. The performance (fitness) of each model in the population was assessed using a fitness function (f total ), defined as the linear combination of five criteria (f ), inspired by characteristics of the ideal controller-motor system proposed by Heckman & Enoka (2004). The five criteria were precision of force, rate of force development, fatigue resistance, force generating capacity and muscle compactness. For each criterion a range of realistic values representing the motor output quality (from poor to optimal performance) was linearly related to fitness function, so a value between 0 (lowest fitness) and 1 (highest fitness) was obtained. The total fitness was expressed as: f total =

5 

wc · f c

(2)

c=1

where w denotes the weight assigned for each criterion (c). Changing w for each criterion allows to assign more

6142

J. L. Dideriksen and D. Farina

relative importance to certain characteristics with respect to others. The first fitness criterion (f 1 ) was related to force precision and was defined as the linear combination of a rating

A

J Physiol 591.24

of the force grading accuracy and the steadiness of force. The force grading accuracy was quantified by imposing an 8 s excitation trajectory increasing linearly from 0 to 90% of the maximum excitation level. The linearity of the

B MU1

MU1

MU2

MU2 MU181

MU368

RT

MDR

MU1

28

19

31

1.4%

30

74

MU2

48

11

25

7.4%

21

...

...

...

...

...

...

...

...

0.4%

119

MU368

19

4

76

1.2%

48

RT

MDR

PDR

IN

CT

MU1

22

19

31

1.4%

32

MU2

14

7

80

4.0%

...

...

...

...

MU181

29

9

59

PDR

IN

CT

Σ=38,929

Σ=89,310

C

D MU1

MU1

MU2 MU2

MU390

MU181

RT

MDR

PDR

IN

CT

MU1

2

9

66

2.5%

74

MU2

48

15

25

5.3%

...

...

...

...

MU390

19

4

71

RT

MDR

PDR

MU1

2

6

66

2.5%

74

21

MU2

14

11

80

4.0%

61

...

...

...

...

...

...

...

...

1.2%

48

MU181

11

9

59

0.2%

124

Σ=42,180

IN

CT

Σ=89,310

Figure 1. The generation of new motor unit populations in the evolutionary optimization process A and C, represent the two most functionally effective motor unit populations of the total number of 50 populations in each optimization. The three representative neurons from each population receive synaptic input (from left-hand side) and their axons innervate the muscle. The soma size indicates its recruitment threshold. The minimum motor unit discharge rate is represented as a 500 ms spike train above the axon, whereas the motor unit twitch is shown next to the muscle (duration: 500 ms). The twitch amplitude reflects the innervation number. As shown in the tables below the graphical representation of the neuromuscular system in each panel, each population consisted of 181 and 390 motor units respectively with random values assigned to the parameters of the neuron: recruitment threshold (RT, in arbitrary units of excitation), minimum discharge rate (MDR, in pps), peak discharge rate (PDR, in pps), innervation number [IN, expressed in percentage of the total number of fibres (bottom of this column)], and twitch contraction time (CT). B and D, two new motor unit populations defined as random combinations of the properties of the single motor units of the populations in (A) and (C), as indicated by the different shades of grey. For example, the new population in (B) adopts the properties of the first motor unit from the population in (A) and from the second and last motor unit from the population in (C). In the new populations, a subset of the parameter values (bold and underlined) are exposed to random variations (mutations). For example, the recruitment threshold of the first motor unit of the population in (B) (whose values were adopted from the population in A), was changed from 22 to 28.  C 2013 The Authors. The Journal of Physiology  C 2013 The Physiological Society

Orderly recruitment and motor performance

J Physiol 591.24

force response was estimated as the root mean square error (RMSE) from a linear function from 0 to 90% of maximum force. A perfectly linear force response (RMSE = 0) was considered as having optimal force grading accuracy whereas RMSE higher than 20% maximum voluntary contraction (MVC) was given the lowest fitness (Heckman & Enoka, 2004). In additional simulations, the optimal excitation–force relation was defined as an exponential function, so the gain of this relation depended on the contraction level, with smaller gain at lower forces. This relation corresponds to a more precise grading of low forces and a rapid grading of high forces. Two exponential functions were applied, with the ratio between the gains for forces 0–10% MVC and for forces 80–90% MVC equal to 0.4 and 0.2. For example, a ratio of 0.4 indicated that the gain in the relation excitation–force for low forces was 40% of the gain for high forces. Values of RMSE within the interval 0–20% MVC were linearly related to the fitness value. The level of force steadiness depended on the coefficient of variation (CoV) for force generated by a constant level of synaptic input to the motor neurons. This input level was maintained for 5 s and corresponded to 5% MVC force. This particular contraction level was selected as force steadiness is the lowest at low forces (45 ms) at a contraction level of 20% MVC, as these units exhibit little decline in twitch force in sustained contractions (Burke et al. 1973) and have a smaller energy consumption (Bolstad & Ersland, 1978) than faster contracting motor units (fitness range: 0.6–1). It is important to note that this separation in slow and fast twitch motor units based on contraction time was arbitrary as the population of units is continuous and not binomial (Enoka, 2012). Nevertheless, pilot tests showed clearly that the results change only minimally (without any implications for the general conclusions) by other arbitrary divisions of motor units, including the possibility of adding further classes of units, such as that of fast and fatigue-resistant motor units (results not shown). With normalization, the fitness criterion f 3 was defined as: f 3 = 2.5 ·

F CT>45ms − 1.5 F total

(5)

The fourth fitness criterion (f 4 ) was related to the muscle force generating capacity and was linearly associated to the absolute MVC force (fitness range: 0–1000 N): f4=

F MVC 1000 N

(6)

The last fitness criterion (f 5 ) favoured muscle compactness as it was assumed that excessive muscle size and weight would be a disadvantage for the overall body energy expenditure. Therefore, this criterion was simply inversely related to the muscle cross-sectional area (A; fitness range: 50–5000 mm2 ), assuming a specific muscle tension of 0.2 N mm−2 (Powell et al. 1984). With normalization, the criterion was defined as follows: f 5 = A · 2 · 10−4 + 1

(7)

The global fitness function was defined as in eqn (2) as the weighted sum of the five criteria listed above. The weights were varied in the different simulations to reach solutions corresponding to models optimized with emphasis on different motor performance characteristics. For example, a weighted average f total for which w1 = 0 corresponded to a model with characteristics tuned to optimize all criteria with the exception of the force generation accuracy. Such optimization would lead to models able to generate high and fast but not accurate forces. Starting from the initial population of 50 models, with parameters randomly selected, the 10 population members (models) that determined the highest value of the overall fitness function [eqn (2)] were selected to ‘reproduce’ and to generate 10 new members of a new population (models), replacing the 10 members with the lowest fitness. These 10 most fit members

J. L. Dideriksen and D. Farina

Simulation strategies

Two main simulation strategies were applied. In the first strategy, the five weights of the global fitness function (eqn (2)) were chosen to reflect the functional requirements of widely investigated muscles of the hand and leg. This allowed to compare the optimal solutions generated by the genetic algorithm with known ranges of physiologic parameters. In the second strategy, the five weights in the global fitness function were systematically varied independently to cover a greater class of optimal solutions with a variety of combinations of weights. Both simulation strategies were performed with and without a fixed association between the innervation number and recruitment threshold of the motor unit and. When this association was not imposed, innervation numbers and motor neuron sizes were independently adjusted to

J Physiol 591.24

A

29

…

60

82

30

…

64

64

24 …

24 …

…

80 …

80 …

18

32

…

79

79

26

…

67

67

28

28

…

61

81

1

2

20 30

Population fitness

26

B

61

62

Iteration Force (% MVC)

were divided randomly into five pairs of two members (parent-members), whose neuromuscular parameter values were used to generate new sets of parameters characterizing two new members. The two new members were generated by randomly mixing the model parameters of the two parent-members. First, the number of muscle fibres and motor neurons of the new members were set to the values of one of the parent-members (randomly choosing between the parents). Next, the new members adopted the values of recruitment threshold, innervation number, twitch contraction time, and minimum and peak discharge rate for each motor unit of one or the other parent-member. The innervation numbers were normalized so their sum was equal to the total number of muscle fibres. To mimic mutation, a random value (normal distribution, zero mean, standard deviation equal to 0.2 × parameter range) was added to 40% (randomly selected) of all values of the neuromuscular parameters of the new members. The addition of mutation finalized the generation of a new population of 50 models. Figure 2 schematically explains this evolutionary optimization method. The above process of model generation was repeated by creating iteratively a new population of 50 models until the global fitness function of all the models in the population had converged to stable values, so that further iterations would not change the members of the population. This was considered as the optimal population. The results are described based on the optimal population, so that all results are reported for all the 50 models in the final, optimal population. If the converged fitness value was below 75% of the maximum value, however, the process was repeated, so that all optimized populations presented in the results had a minimum fitness (averaged over all 50 models) always at least equal to 75% of the best possible fitness.

5.5

5

4.5 Force (% MVC)

6144

3

2

4

CoV=2.1% CoV=0.7% 5

Time (s)

100 80 60 40 20 0

0

0.05

0.1 Time (s)

t=119 ms t=54 ms 0.15 0.2

Figure 2. Evolutionary optimization scheme of the neuromuscular models A, each square represents a motor unit population (members). In each column (representing each iteration of the optimization process), six of the 50 members are shown. The number within each square represents its estimated fitness value in per cent. Crossed squares indicate excluded members after each iteration (lowest numbers) and the arrows show how members are transferred to next iteration, including the generation of new members (two arrows merging from small circles). For example, in the first iteration (first column), the first and fourth member were excluded (as 18% and 20% were the lowest fitness values), to be replaced by ‘children’ of the second and sixth member (with highest fitness values; 28% and 30%) in the next iteration (second column). The average fitness of all members is shown on the graph below, illustrating how it increases progressively throughout the optimization process until convergence after 62 iterations. In this example, fitness criteria f 1 (force precision) and f 2 (rate of force development) were assigned the highest weights. B, examples of the performance (force steadiness and maximum rate of force generation) of two selected motor unit populations from (A). The grey lines indicates a motor unit population that was excluded in the first iteration (grey square, first line of first column with fitness 20%) and the black lines the best motor unit population of the last iteration (black square, first line of fourth column with fitness 82%). As indicated by the values of CoV of force and time to reach 50% MVC force, the optimized member provided the steadiest force and the fastest force generation. CoV, coefficient of variation; MVC, maximum voluntary contraction.

 C 2013 The Authors. The Journal of Physiology  C 2013 The Physiological Society

J Physiol 591.24

Orderly recruitment and motor performance

understand which of these parameters separately had an influence on fitness function. To determine the five weights in eqn (2) that would reflect typical muscles of the hand and the leg, we referred to long-term surface EMG recordings from the first dorsal interosseous and the vastus lateralis muscles of healthy subjects (Kern et al. 2001). In that study, EMG was recorded continuously during daily life activities (outside the laboratory) using portable devices during two 10 h periods in 14 healthy adults. In these periods, the subjects were instructed to follow a normal daily routine. It was found that the average EMG burst amplitude was lowest for the first dorsal interosseous, indicating a higher demand for force precision for this muscle [this led us to choose w 1 = 3.5 for the hand and w 1 = 2.5 for the leg in eqn (2)], whereas the significantly higher occurrence of near-maximal EMG bursts combined with a short average burst duration for the vastus lateralis indicated a greater priority for contraction speed and maximum force for the leg muscle [which was translated in w 2 = 2 for the hand and w 2 = 3.5 for the leg in eqn (2)]. The demand for fatigue resistance was assumed slightly higher for the hand muscle, as the integral of the EMG was highest for the first dorsal interosseous (thus, hand w 3 = 3 and leg w 3 = 2.5). Finally, it was assumed that maximal force has an higher priority for the leg, and muscle compactness is functionally more important for the hand muscle, due to the relatively smaller bones of the hands and more proximal location (hand: w 4 = 0.5, w 5 = 1; leg: w 4 = 1, w 5 = 0.5). The set of weights determined in this way should not be considered as reflecting the demands of specific muscles, but only as separating the functional demands of representative muscles in different parts of the body. Moreover, the weights are relevant in relative terms, thus dividing all the weights for both muscles by a constant value would not change the results. The simulations for each of the two representative muscles were repeated four times with and without the fixed association between innervation number and recruitment threshold. In the second simulation set, the weights w 1 –w 3 corresponding to the criteria f 1 –f 3 were varied independently across three values [0, 0.5, 1] (26 combinations, when excluding the combination with all three weights set to 0). Each combination was repeated twice, yielding a total of 52 simulations. Data from pilot simulations, in which the weights of the criteria f 4 and f 5 were varied in a similar manner (results not shown), indicated that these two criteria only significantly influenced the number of muscle fibres but not the number of motor neurons or other neuromuscular parameters. Thus, the values for f 4 and f 5 were fixed to 0.25 for all simulations in the second paradigm. The second simulation set thus explored the effect of varying the relative importance of the first three criteria over their entire range. In the additional set of simulations in which  C 2013 The Authors. The Journal of Physiology  C 2013 The Physiological Society

6145

the optimal force grading sensitivity [eqn (3)] was defined according to exponential functions, w 1 was set to 1 and w 2 –w 5 to 0.25. Simulations with the linear and each of the two exponential functions were repeated twice, without a fixed association between innervation number and recruitment threshold.

Analysis and statistics

For the second set of simulations, the best linear fit between the innervation number and recruitment threshold, twitch contraction time as well as minimum and peak discharge rate were estimated for the motor neurons of each model in the last population (optimized). To ensure that the linear regression between innervation number and recruitment threshold was not biased in cases with a low number of large motor units (high degree of skewness in the innervation numbers), the difference between the average recruitment thresholds of the motor units with innervation numbers below and above 50% of the maximum innervation number, was computed. Furthermore, the proportion of motor units with the low innervation number (less than 50% of the maximum innervation number within the motor unit population) and the innervation ratio (the ratio between the lowest and highest innervation numbers within the motor unit population) were calculated for each model of the optimized population. The average values of these parameters for the 50 models within each of the 52 applied combinations of the fitness criterion weights were obtained (yielding a total of 52 values for each parameter). The results are presented as a function of each of the three fitness criteria whose relative importance was varied, averaging over the variability due to the other two criteria. In this way, three mean values (and standard deviations) were obtained for each fitness criteria, representing all simulations in which that criteria was assigned the value 0, 0.5 or 1, respectively. This allowed for direct comparison of the influence of the weights of each of the three fitness criteria on the optimized model parameters and the correlation between them. Student’s t test was used to determine if the correlation values were different from zero or different across simulation conditions. The main outcome variable was the correlation between innervation number and recruitment threshold when the association between innervation number and recruitment threshold was not imposed. The strength of this correlation indeed indicated the functional relevance of an orderly recruitment. Absence of correlation would indicate that orderly recruitment does not result in functional advantages with respect to other random recruitment orders. In this case, orderly recruitment would not have emerged for functional needs but as an indirect consequence of a range in recruitment numbers, due to

6146

J. L. Dideriksen and D. Farina

the physiological association between innervation number and motor neuron size and between motor neuron size and excitability. A second main outcome of the simulations was the comparison between the overall fitness level with and without a fixed association between motor neuron recruitment threshold and innervation number, i.e. with and without orderly recruitment. A paired t test was employed to test the difference. A negligible difference between such fitness levels would indicate that orderly recruitment is not necessary for function optimization. Results In the first set of simulations, the weights of the optimization fitness criteria (i.e. force precision, rate of force development, fatigue resistance, maximum force and muscle compactness) were set to reflect the functional demands of typical muscles of the hand and leg. These simulations were performed to verify the ability of the evolutionary optimization approach to capture known differences in the neuromuscular properties across muscles. Simulations were performed with and without a fixed association between the innervation number and recruitment threshold of the motor units. Figure 3 depicts innervation numbers from one simulation for the leg muscle without this fixed association. The distributions of innervation numbers covered a large range (average

Proportion of motor unit population (%)

70 60 50 40 30 20 10 0

0

0.5

1.5 1 Innervation number (%)

2

Figure 3. Optimized motor unit innervation numbers for one of the four simulations representing the leg muscle performed without a fixed association between innervation number and recruitment threshold The optimized distributions of innervation numbers (expressed as a percentage of the maximum innervation number) were skewed toward low innervation numbers. The average fitness for all 50 members at the final iteration was 75.9 ± 1.8%.

J Physiol 591.24

values: hand, 83-fold; leg, 262-fold) and were skewed towards lower values, as observed experimentally (Enoka & Fuglevand, 2001). Table 1 shows the differences in neuromuscular properties for the two simulated conditions representing the leg and hand muscle with and without the imposed association between innervation number and recruitment threshold. The simulated number of motor units was in the range of those observed in many muscles (Heckman & Enoka, 2004). As experimentally observed, the average innervation number was highest for the leg muscle (Heckman & Enoka, 2004). Furthermore, the results reflected the experimental observation that larger proximal muscles of the extremities tend to have a larger proportion of motor units with fast contraction times compared to smaller hand muscles (Enoka & Fuglevand, 2001).The ranges of innervation numbers were large (average values: hand, 56-fold; leg, 267-fold) and were consistently skewed towards low values (Fig. 3). There was no correlation between innervation number and recruitment threshold when this relation was not imposed, whereas the twitch contraction times (hand, 77.7 ± 34.4 ms; leg, 72.1 ± 35.7 ms) were inversely correlated with the motor unit innervation number (Stephens & Usherwood, 1977). This correlation tended to be weaker for the hand muscle (not significant; P = 0.2), which may explain why this relation is rarely found experimentally for this type of muscles (Fuglevand, 2011). In the second set of simulations, the weights for the three fitness criteria related to force precision, maximum rate of force development, and fatigue resistance were systematically varied to investigate their relative influence on the neuromuscular parameters. As in the first set of simulations, these simulations were performed with and without a fixed association between motor unit innervation number and recruitment threshold. Figure 4 shows examples of representative optimized model performance in the second set of simulations without this association. Here, examples of motor unit populations are shown when optimized exclusively with respect to one of the criteria f 1 –f 3 (force precision, rate of force development and fatigue resistance). A motor unit population optimized mainly for force precision can generate highly steady forces (black line in Fig. 4A) but it can only provide relatively slow rapid contractions (black line in Fig. 4B). Similarly, a motor unit population optimized for rapid force generation (grey line, Fig. 4B) may have poor fatigue resistance (grey bar, Fig. 4C). Furthermore, differences in the functional demands involved large differences in the neuromuscular properties, such as number of motor units, which was high (479) when optimized for force precision, but very low (∼20) when optimized exclusively with respect to the other criteria.  C 2013 The Authors. The Journal of Physiology  C 2013 The Physiological Society

Orderly recruitment and motor performance

J Physiol 591.24

6147

Table 1. Optimal neuromuscular parameters for hand and leg muscles With fixed association between recruitment threshold and innervation number

Number of motor neurons Average innervation number Proportion of fast-twitch units (%) Correlation (innervation number – recruitment threshold) Correlation (innervation number – contraction time)

Hand

Leg

Hand

Leg

341.4 ± 162.4 40.8 ± 48.2 16.7 ± 1.0

352.6 ± 112.5 314.8 ± 126.2 21.7 ± 6.7

285.4 ± 232.7 68.5 ± 105.0 17.9 ± 3.1

257.3 ± 108.6 504.0 ± 319.4 26.0 ± 6.5

1 ± 0

1 ± 0

0.05 ± 0.21

0.01 ± 0.11

−0.14 ± 0.15

−0.25 ± 0.09

−0.11 ± 0.16

−0.27 ± 0.15

The examples in Fig. 4 are purely theoretical as in each case two of the optimization criteria are given weight zero, i.e. they have no relevance in the final optimal model. None the less, they indicate the working principle of the optimization scheme that identifies A

solutions that are optimal according to the importance (weight) given to each criteria. In the more general case, the weights of the criteria all differ from zero and their relative values determine the final optimal solution.

B

C 100

MU: 21±34 IN ratio: 0.3±0.3 % Low IN MUs: 87.1±7.5% r (IN-CT): -0.44±0.14 r (IN-RT): -0.44±0.05

MU: 479±46 IN ratio: 5.1±2.1 % Low IN MUs: 89.2±7.9% r (IN-CT): 0.04±0.04 r (IN-RT): -0.01±0.04

5.5

MU: 19±52 IN ratio: 0.6±0.5 % Low IN MUs: 86.1±7.1% r (IN-CT): 0.15±0.19 r (IN-RT): 0.05±0.45

5

4.5 3 Time (s)

4

Fatigue-resistant force (%)

Force (% MVC)

Force (% MVC)

100

2

Without fixed association between recruitment threshold and innervation number

50

0 0

100 200 Time (ms)

300

50

0

Figure 4. Performance of three MU populations, each optimized with respect to one fitness criterion, and without an imposed association between MU IN and RT Black lines/bar represents one MU population optimized with respect to force precision (w 1 = 1, w 2 = 0, w 3 = 0), whereas grey lines/bar represents optimized populations with respect to rate of force development (w 1 = 0, w 2 = 1, w3 = 0) and black dashed lines represent optimizations with respect to fatigue resistance (w 1 = 0, w 2 = 0, w 3 = 1). A, performance of the population optimized with respect to force precision (black) is highlighted (bold line) and the average parameters of all 50 MU populations are reported (black box). Similarly, the performance of the two other populations are highlighted (with bold line in B and C) and their parameters (average ± S.D.) reported in (B) and (C). The level of force steadiness (A) was highest for the MU population optimized for this feature (black line). This was achieved with a high number of MUs (black box in A). Similarly the other MU populations (grey and black, dashed) were visibly better with respect to rate of force development (B) and fatigue resistance (C), and both had very low number of MUs (grey box in B and black, dashed-line box in C). IN ratio, lowest innervation numbers expressed as a percentage of the highest INs in the MU population; low IN MUs, the proportion of MUs with INs less than half of the maximum IN in the MU population; MU, number of motor units; MVC, maximum voluntary contraction; r, correlation coefficient between IN and MU twitch contraction time (CT) and recruitment threshold (RT), respectively.  C 2013 The Authors. The Journal of Physiology  C 2013 The Physiological Society

6148

J. L. Dideriksen and D. Farina

J Physiol 591.24

B Normalized recruitment threshold

A 0.7

1

Force precision: w1 = 0: r=0.16±0.21 * w1 = 0.5: r=0.04±0.05 w1 = 1: r=0.01±0.06

0.8 0.6

0.6

0.4 0.2

0.5 0

0.5

0 0

1

0.5

1

D Normalized recruitment threshold

C 0.7

1

Rate of force development w2 = 0: r=0.05±0.11 w2 = 0.5: r=0.11±0.11 w2 = 1: r=0.05±0.10

0.8 0.6

0.6

w1 = 0 w2 = 0.5 w3 = 0.5 r = 0.93 RMSE = 0.09

0.4 0.2 0.5 0

0.5

0

1

E

0

0.5

1

Normalized recruitment threshold

F 0.7

Fatigue resistance w3 = 0: r=0.01±0.07 w3 = 0.5: r=0.12±0.11 * w3 = 1: r=0.06±0.13

1 0.8 0.6

0.6

w1 = 0.5 w2 = 0.5 w3 = 0.5 r = 0.12 RMSE = 0.29

0.4 0.2 0.5 0

0.5 1 Normalized innervation number

0

0

0.5 1 Normalized innervation number

Figure 5. The relation between innervation number and recruitment threshold in the optimized models Recruitment thresholds were expressed in units of excitation, normalized to the maximum recruitment threshold of the motor unit population. A, each of the three lines represents the average of simulations in which the weight of the force precision criteria was assigned either 0 (black), 0.5 (dark grey) or 1 (light grey). The lines represent the average regression lines for all optimized models with this fitness weight. The mean ± S.D. correlation coefficient of each line is reported (∗ correlations significantly different from zero). The average RMSE was between 0.28 and 0.31. A, strongest relation between innervation number and recruitment threshold (r = 0.16 ± 0.21) occurred at w 1 = 0. Similarly, (C) and (E) show the average regression lines for rate of force development and fatigue resistance, respectively. B, all average regression lines for the 52 simulations performed with different combinations of fitness weights. In this way, the three lines in (A) represent the average of three subgroups of the lines in (B), divided according to the value assigned to w 1 (and similarly for w 2 and w 3 in C and E). D, data from the motor unit population (of the 50 populations optimized within each setting of the weights of the fitness criteria) with the highest correlation coefficient of the optimization setting that yielded the highest average correlation (w 1 = 0, w 2 = 0.5, w 3 = 0.5). Each dot represents one motor unit and the line represents the best linear fit. Maximum innervation number was 3555. A clear relation is seen in this case (r = 0.93, P < 0.001), but this model contained only 10 motor units. F, in contrast, represents the motor unit population with the largest correlation between innervation number and recruitment threshold for the populations containing more than 150 motor units (w 1 = 0.5, w 2 = 0.5, w 3 = 0.5). The maximum innervation number was 2080. In this case, the correlation was lower, but still significant (r = 0.12, P = 0.04). RMSE, root mean square error.

 C 2013 The Authors. The Journal of Physiology  C 2013 The Physiological Society

Orderly recruitment and motor performance

J Physiol 591.24

Does the size principle have functional determinants?

Figure 5 depicts the correlation between innervation number and recruitment threshold (in the simulations without the fixed association between these) when varying the value assigned to each of the three fitness criteria (force precision: Fig. 5A; rate of force development:

60 40 20

C

RMSE = 2.4% 1

2

3

100

Force (% MVC)

80 60 40 20

E

0 0

RMSE = 3.8% 1

2

3

100

Force (% MVC)

80 60 40 20 0 0

RMSE = 3.6% 1

2

3

Time (s)

D Normalized recruitment threshold

Force (% MVC)

80

Normalized recruitment threshold

B

100

0 0

Fig. 5C; fatigue resistance: Fig. 5E). The three lines in each of these figures represent the average best linear fit for all simulations performed with the weight set to 0, 0.5 and 1 for that particular criterion. There was a general tendency towards a positive correlation, but on average this was only significant in two cases (w 1 = 0 and w 3 = 0.5). This finding was confirmed by the 1 0.8 0.6 0.4 0.2 0

r = -0.16±0.05 0

0.5

0

0.5

1

1 0.8 0.6 0.4 0.2 0

F Normalized recruitment threshold

A

6149

r = 0.10±0.05 1

1 0.8 0.6 0.4 0.2 0

r = 0.17±0.06 0 0.5 Normalized innervation number

1

Figure 6. The relation between innervation number and recruitment threshold with different optimization criteria for force grading sensitivity Here, the models were optimized with high emphasis on force precision (w 1 = 1, w 2−5 = 0.25). A, C and E, representative force trace in response to a linear increase in excitation (thick line) and the optimal relation (dashed line), which was either linear (A), or exponential (with a ratio of 0.4, B, and 0.2, C, between the gains for forces 0–10% MVC and 80–90% MVC) (see text for details). B, D and F, corresponding average linear relations (lines) between normalized values of innervation number and recruitment threshold across the two repetitions of each setting. The grey dots represent the single motor unit innervation number and recruitment threshold for one model realization in each setting with the highest correlation coefficient: (B) r = 0.06; (D) r = 0.22; (F) r = 0.29. MVC, maximum voluntary contraction; RMSE, root mean square error.  C 2013 The Authors. The Journal of Physiology  C 2013 The Physiological Society

6150

J. L. Dideriksen and D. Farina

difference between the average recruitment thresholds of the motor units with low and high innervation numbers (below and above 50% of maximum innervation number), which was on average 5.5 ± 15.9% (range: −26.4–43.4%). Even in cases with significant linear correlation, the correlation was much lower than observed experimentally (Milner-Brown et al. 1973). The few cases where a strong correlation arose (Fig. 5D reports one of these cases) were associated with low numbers of motor units. In contrast, all simulations that involved a more realistic number of motor units (>150; Heckman & Enoka, 2004), strong associations never emerged (Fig. 5F reports a typical association between innervation number and threshold in realistic conditions). Hence, the general lack of the orderly recruitment observed in the simulations representing the hand and leg muscles (Table 1) was here confirmed in the more general case. Figure 6 shows the correlation between innervation number and recruitment threshold when the models were optimized with high emphasis on force precision (w 1 = 1, w 2−5 = 0.25) and with exponential excitation–force relations. With exponential relations, the association between innervation number and recruitment threshold remained weak, with correlation coefficient values always below 0.3. For example, in the simulations where the gain was five times higher for the range 80–90% MVC than from 0 to 10% MVC (Fig. 6E), the average correlation coefficient was 0.17 ± 0.06, whereas a negative correlation was present in the linear case (Fig. 6A). The average fitness value of all simulations (n = 52 with and without the imposed association) was 79.2 ± 3.9% and 78.2 ± 5.2% of the maximum fitness for the conditions with (corresponding to orderly recruitment) and without the fixed association between innervation number and recruitment threshold (not significantly different; P = 0.07). Figure 7 shows the similarity of the functional quality of models optimized with (Fig. 7B, D, F and H) and without (Fig. 7A, C, E and G) the imposed orderly recruitment. In this example, the models were optimized with highest weights assigned to force precision and rate of force development (w 1 = 1, w 2 = 1, w 3 = 0.5). Only modest, non-systematic differences in the force steadiness, force grading accuracy and time to reach 50% MVC force were observed between the case of orderly recruitment (obtained by fixing a priori the association between innervation number and excitability) (Fig. 7B, D, F and H) and random recruitment (Fig. 7A, C, E and G). These results indicate that the functional quality of the simulated muscles was not substantially better with orderly recruitment according to motor unit size. There was no significant difference in the optimal number of motor units (in most cases 200–300; Fig. 8A–C) for the different weights of the fitness criteria when simulations were performed with and without a fixed association between the motor unit recruitment threshold and innervation number. Whereas

J Physiol 591.24

the differences in the proportion of motor units with low innervation number (Fig. 8F) and range of innervation numbers (Fig. 8I) were significant across the two conditions, they clearly exhibited similar trends (>80% motor units with low innervation number and a 50–200-fold difference in highest and lowest innervation number). These observations indicate that the similar average fitness values across the simulations with and without orderly motor unit recruitment according to size were not achieved by a radically different configuration of other neuromuscular parameters. Overall, orderly recruitment only emerged in a few non-physiological cases (number of motor units 80%) for all combinations of fitness criteria (Fig. 8D–F). The same was the case for the range of innervation numbers (Fig. 8G–I). Rate of force development was the parameter with the highest influence on the distribution of innervation numbers, increasing the skewness and range when assigned a high weight (Fig. 8E and H). Figure 9 depicts the correlation between innervation number and motor unit contraction time (in the simulations without a fixed association between innervation number and recruitment threshold) as a function of the value assigned to each of the three fitness criteria (force precision: Fig. 9A; rate of force development: Fig. 9C; fatigue resistance: Fig. 9E). The three lines in each of these figures represent the average best linear fit for all simulations performed with the weight set to 0, 0.5 and 1 for that particular criterion. In most cases, an inverse association is present, as experimentally observed (Stephens & Usherwood, 1977), with the rate of force development being the most influential determinant of its strength. When the weight of this criterion was set to 0, a positive correlation was present (Fig. 9C). The strongest  C 2013 The Authors. The Journal of Physiology  C 2013 The Physiological Society

Orderly recruitment and motor performance

J Physiol 591.24

6151

B

Force (% MVC)

A 60

60

40

40

20

20

0

Force (% MVC)

C

500 ms 0

1 Time (s)

2

5.1

5

5

4.9

E Force (% MVC)

2

500 ms

CoV = 0.52±0.10% 2

1

1.5 Time (s)

0

RMSE = 7.6±2.3% 3 1.5 Time (s)

0

T = 64.2±12.0 ms 50 100 Time (ms)

F 100

80

80

60

60

40

40

20 0

G

RMSE = 7.1±2.6% 1.5 3 Time (s)

20 0

H

80

80

60

60

40

40

20 0

1

4.9 CoV = 0.55±0.09% 2 1.5 Time (s)

100

0

0

Time (s)

5.1

1

Force (% MVC)

0

D

20 0

T = 68.6±11.6 ms 50 100 Time (ms)

0

Figure 7. The performance of models optimized with (B, D, F, H) and without (A, C, E, G) the imposed association between innervation number and recruitment threshold The weights of the optimization criteria were: w 1 = 1, w 2 = 1, w 3 = 0.5 (highest priority of force precision and rate of force development). A and B, representative force ramp with the recruitment threshold (grey circles) and twitch properties (grey box) of five motor units. A, motor units with high twitch amplitudes are recruited throughout the contraction, and not only at high contraction levels as in (B). Note the scale of force amplitude is different in the grey box compared to that in the force ramp. In the other panels the performance of the 10 members from each population with the highest fitness are shown concerning the criteria reflecting force precision (C–F) and rate of force development (G, H). The differences between the performance of the two conditions were small. For example, the force steadiness at 5% MVC (C, D) was slightly higher in the case with the imposed association (C), whereas the opposite trend was present for the force grading accuracy (E, F). CoV, coefficient of variation; MVC, maximum voluntary contraction; RMSE, root mean square error.

 C 2013 The Authors. The Journal of Physiology  C 2013 The Physiological Society

6152

J. L. Dideriksen and D. Farina

correlations occurred in cases with very few motor units present (n = 9; Fig. 9D), but the inverse correlation was still present in more realistic conditions (Fig. 9F). The average minimum and peak discharge rates were 10.5 ± 5.1 pps and 59.7 ± 23.0 pps, respectively across all simulations. The average peak discharge rate was substantially higher (71.1 ± 6.7 pps) in the simulations with the values of the weights set to focus on rate of force development (w 1 = 0, w 2 = 1, w 3 = 0). The same trend has been observed experimentally in subjects after ballistic contraction training, although here higher peak discharge rates were observed (probably as phenomenons such as doublet and triplet discharges are not included in the applied model) (Van Cutsem et al. 1998). The average minimum discharge rate and its tendency towards a slightly positive correlation

A

for all weights of all fitness criteria with innervation number (statistically significant in four of nine conditions) corresponded well to experimental observations (Moritz et al. 2005). Furthermore, the correlations between peak discharge rates and innervation number were not strongly influenced by variations in fitness weights (average correlation coefficients below 0.2; statistically significant in five of nine cases), indicating that the various associations between peak discharge rates and motor unit size, experimentally observed across muscles (Gydikov & Kosarov, 1974; De Luca et al. 1982) may not have a large functional impact. To summarize the results, optimization of fitness criteria related to functional tasks produced a large range of innervation numbers (Fig. 3) but not the orderly recruitment of motor units (when the recruitment

B Number of motor units

J Physiol 591.24

C

500

500

500

400

400

400

300

300

300

200

200

200

100

100

100

0

0

0

Low-IN motor units (%)

D

E

100

100

95

95

95

90

90

90

85

85

85

80

80

80

75

75

75

70

70

70

G Innervation ratio (%)

F

100

H

I

10

10

10

8

8

8

6

6

6

4

4

4

2

2

2

0

0

0

-2

-2

-2 0

1 0.5 Force precision weight

0

0.5 1 Rate of force development weight

0

0.5 1 Fatigue resistance weight

Figure 8. Optimized motor unit model parameters expressed as a function of the weight assigned to several fitness criteria with (black lines) and without (grey lines) an imposed association between innervation number and recruitment threshold The parameters depicted are: the optimal number of motor units (A–C), the proportion of motor units with innervation numbers below 50% of the maximum value within that population (D–F) and the ratio between the innervation number of the 5% of the motor units with the lowest and highest innervation numbers respectively (G–I). Each point represents the mean (±S.D.) value over all the simulations in which the weight of one fitness criterion (indicated at the bottom of each column) was assigned one specific value. A, for example, the left-most point represents the 16 simulations in which w 1 (weight of force precision criterion) was assigned the value 0.  C 2013 The Authors. The Journal of Physiology  C 2013 The Physiological Society

Orderly recruitment and motor performance

J Physiol 591.24

A

B

Contraction time (ms)

110

55

0 0

6153

140

70 Force precision: w1 = 0: r=-0.31±0.17 * w1 = 0.5: r=-0.09±0.06 * w1 = 1: r=-0.07±0.06 * 0.5

0

1

0

0.5

1

D Contraction time (ms)

C

110

140

55

70

0

Rate of force development w2 = 0: r=0.10±0.09 * w2 = 0.5: r=-0.28±0.11 * w2 = 1: r=-0.32±0.10 * 0

0.5

E

F

Contraction time (ms)

110

0

1

0.5

140

w1 = 1 w2 = 1 w3 = 1 r = -0.36 RMSE = 32.1

70

55

0

0

1

w1 = 0 w2 = 1 w3 = 1 r = -0.95 RMSE = 12.8

Fatigue resistance w2 = 0: r=-0.19±0.07 * w2 = 0.5: r=-0.15±0.10 * w2 = 1: r=-0.15±0.11 * 0

0.5

1

Normalized innervation number

0

0

0.5

1

Normalized innervation number

Figure 9. The relation between innervation number and motor unit twitch contraction times in the optimized models with an imposed association between motor unit innervation number and recruitment threshold A, each of the three lines represents the average regression lines calculated from all simulations in which the weight of the force precision criteria was assigned either 0 (black), 0.5 (dark grey) or 1 (light grey). The mean ± S.D. correlation coefficient of each line is reported (∗ correlations significantly different from zero). The average RMSE was between 29.5 and 33.7. In A, the strongest relation between innervation number and recruitment threshold (r = −0.31 ± 0.17) occurred at w 1 = 0. Similarly, C and E show the average regression lines for rate of force development and fatigue resistance, respectively. B, all average regression lines for the 52 simulations performed with different combinations of the fitness weights. In this way, the three lines in (A) represents the average of three subgroups of the lines in (B), divided according to the value assigned to w 1 (and similarly for w 2 and w 3 in C and E). D, data from the motor unit population (of the 50 populations optimized within each setting) with the highest correlation coefficient of optimization setting that yielded the highest negative average correlation (w 1 = 0, w 2 = 1, w 3 = 1). Each dot represents one motor unit and the line represents the best linear fit. Maximum innervation number was 14,662. A strong negative relation is seen in this case (r = −0.95, P < 0.001), but this model contained only nine motor units. In contrast, (F) represents the motor unit population with the highest negative correlation between innervation number and recruitment threshold for the populations containing more than 150 motor units (w 1 = 1, w 2 = 1, w 3 = 1). The maximum innervation number was 2080. In this case, a strong correlation was still present (r = −0.36, P = 0.008). RMSE, root mean square error.

 C 2013 The Authors. The Journal of Physiology  C 2013 The Physiological Society

6154

J. L. Dideriksen and D. Farina

threshold was not directly associated to the innervation number) (Fig. 5). Imposing the association between innervation number and motor neuron recruitment threshold, which also implied an association between size and recruitment threshold and thus an orderly recruitment, did not have a substantial influence on fitness levels (e.g. Fig. 7). Overall, these simulations indicate that orderly recruitment is not related to the functional advantages of such recruitment order, but rather it indirectly emerges from the advantage of a large range of innervation numbers across motor neurons. Discussion This study examined the important principle in motor physiology that orderly recruitment of motor units confers a functional advantage during force generation. This principle is valid at the level of intrinsic excitability of individual motor neurons and has been verified in a variety of experimental conditions during natural movements, with few exceptions. Here we argue against the assumption that orderly recruitment in itself has important functional implications as often assumed (Cope & Pinter, 1995; Kandel et al. 2000; Heckman & Enoka, 2004; Stein et al. 2005). The approach we followed for verifying the functional implications of orderly recruitment is based on the concurrent optimization of the number of muscle fibres and motor units, recruitment order, innervation numbers, values and range of discharge rates, and twitch contraction times, without any a priori assumption about the values, distributions or associations between these parameters, which were freely modified in the optimization scheme within broad ranges of values. The optimization function was a combination of force steadiness, speed and fatigue resistance, as well as absolute force producing capacity and muscle size. The evolutionary approach of optimization allowed convergence to optimal solutions in limited computational time, despite the very large size of the input parameter space. This optimization method was employed as a full evaluation of the entire model parameter space for optimization was not possible. The main result of the simulations is the universal importance of a large range of innervation numbers skewed towards low values. Conversely, orderly recruitment never strongly emerged without links to other parameters, and when it was imposed, it did not influence the task performance, according to the optimization fitness criteria. Thus, an orderly recruitment did not provide substantial benefits in task optimization with respect to a random recruitment order whereas the fitness criteria were optimized only when the innervation numbers across motor neurons had a broad (and skewed) distribution. These results indicate that the reason for the almost universal presence of orderly recruitment of motor

J Physiol 591.24

units in force control is probably indirect and does not itself offer a direct functional advantage. Although no a priori assumptions were imposed on the distribution of the neuromuscular properties, the evolutionary optimization approach had a remarkable ability to generate common experimental observations across muscles, for example the near-exponential distribution of motor unit innervation numbers (Figs 3 and 8) (Enoka & Fuglevand, 2001) and the inverse association between innervation number and twitch contraction time (Fig. 8) (Stephens & Usherwood, 1977). Furthermore, the simulations predicted known differences across muscles, such as average innervation number (Heckman & Enoka, 2004) and fibre-type composition (Table 1) (Enoka & Fuglevand, 2001), when the corresponding weights for the fitness function contributions were tuned according to the demands of different muscles. None the less, some other model parameters were less sensitive to optimization, such as the number of motor units. Although the number of motor units was indeed within the range experimentally observed for many muscles, none of the optimization criteria could explain the functional advantage for a muscle to consist of more than approximately 350 motor units (Fig. 8A–C), contrary to the case of large muscles (Heckman & Enoka, 2004). One possible explanation for this result may be muscle compartmentalization, which implies that some muscles consist of several subunits, each receiving separate synaptic input (English et al. 1993). Presumably this phenomenon arises to improve the functional efficiency of the muscle and has been reported for arm (ter Haar Romeny et al. 1984), neck (Schomacher et al. 2012) and respiratory (Butler & Gandevia, 2008) muscles. In such cases, the optimization criteria would apply to each subunit rather than to the entire muscle, so that the total number of motor units in the entire muscle would be much greater than in current simulations. As in any optimization process, the selection of fitness criteria is critical for the outcome. Although based on intuitive and generally accepted views on optimal motor systems (Heckman & Enoka, 2004), the validity of the exact mathematical formulations is uncertain. For this reason, the complexity of the algorithm was minimized by defining each criterion as simply as possible, using low numbers of simulations and linear relations between model output and estimated fitness. In this way, some particular aspects of specific functional advantages may have been overlooked, as only the core of the motor system functionality were the primary focus. For example, whereas the division of motor unit fatigability into two distinct groups probably oversimplifies the complex relation between muscle contractile properties and fatigability in humans (Thomas et al. 1991) and ignores the fact that muscle properties change continuously across motor units rather than in distinct classes (Enoka, 2012), the contraction  C 2013 The Authors. The Journal of Physiology  C 2013 The Physiological Society

J Physiol 591.24

Orderly recruitment and motor performance

time is generally regarded to be a good predictor of motor unit fatigue resistance (Carpentier et al. 2001; Heckman & Enoka, 2004). Such simplifications of the fitness criteria also implied that the focus of the results was on general trends (e.g. whole population correlation coefficients) rather than on details of specific single motor unit properties. The main conclusion of the present study may appear partly in conflict with previous findings showing mathematically that orderly recruitment optimizes the force accuracy of the motor output (Hatze, 1979). The large number of motor units in a muscle explains this apparent contradiction (see below). Moreover, these studies exclusively focused on the influence of recruitment order on one specific feature of motor function, whereas, as emphasized by Henneman and Mendell, all physiological and functional features must be taken into account as some may be mutually exclusive (Henneman & Mendell, 1981). For example, whereas a high proportion of slow twitch muscle fibres favours resistance to fatigue, it may compromise the ability to generate rapid forceful contractions. The approach we used in the present study considers several realistic demands concurrently for finding a solution as an optimal compromise among such different demands, according to the relative importance of the demands that could be changed across simulations. Furthermore, the current results do not imply that an orderly recruitment according to size in itself is necessarily worse than other recruitment modalities, but rather that there are configurations of neuromuscular properties that would provide the same advantages of orderly recruitment at least as effectively, so that there would not be an evolutionary drive for ordered recruitment strategies. Indeed, in the simulated evolutionary approach, the orderly recruitment only emerged in few extreme simulation conditions (Fig. 5D), so that optimization of motor function was generally found with a random recruitment configuration (as depicted in Fig. 5F). Therefore, if the size and excitability of motor neurons were not related to the innervation number, an orderly recruitment would not emerge. The reported results may be counterintuitive. For example, it seems evident that small contributions to small forces are better for a smooth gradation of force. However, this reasoning depends on the number of contributions. Even at very low forces, the number of active motor units is very large if the force contributions of each unit are small enough (Enoka & Fuglevand, 2001), so that even ‘large’ units provide a relatively small contribution to the total force. Smoothness or precision in force generation therefore depends on the large number of motor units in the muscle, which thus fractions the force contributions in many relatively small parts, rather than to a precise recruitment of the smallest contributions. Accordingly, we have recently shown that force steadiness is not influenced  C 2013 The Authors. The Journal of Physiology  C 2013 The Physiological Society

6155

by synaptic noise (random oscillations in membrane potential) when a critical number of motor units are active, and that this critical number depends on contractile motor unit properties (Dideriksen et al. 2012). This result is because the contribution of each motor unit to the total force is sufficiently small (independently on the motor unit size) with respect to the total muscle force, so that it has a minimal impact on force steadiness. This applies also to the motor units with high innervation numbers, especially if the majority of motor units innervate relatively few muscle fibres. In the optimization functions, the model properties almost always converged to reflect these characteristics by means of a total number of motor units in the range 200–300 (Fig. 8A–C), with a skewed distribution of innervation numbers towards a greater proportion of small motor units (Fig. 8D–I). Interestingly, in the few simulations where the number of motor units was small (