Mathematical relationships between spinal motoneuron properties

  1. Arnault H Caillet
  2. Andrew TM Phillips
  3. Dario Farina
  4. Luca Modenese  Is a corresponding author
  1. Department of Civil and Environmental Engineering, Imperial College London, United Kingdom
  2. Department of Bioengineering, Imperial College London, United Kingdom
  3. Graduate School of Biomedical Engineering, University of New South Wales, Australia

Abstract

Our understanding of the behaviour of spinal alpha-motoneurons (MNs) in mammals partly relies on our knowledge of the relationships between MN membrane properties, such as MN size, resistance, rheobase, capacitance, time constant, axonal conduction velocity, and afterhyperpolarization duration. We reprocessed the data from 40 experimental studies in adult cat, rat, and mouse MN preparations to empirically derive a set of quantitative mathematical relationships between these MN electrophysiological and anatomical properties. This validated mathematical framework, which supports past findings that the MN membrane properties are all related to each other and clarifies the nature of their associations, is besides consistent with the Henneman’s size principle and Rall’s cable theory. The derived mathematical relationships provide a convenient tool for neuroscientists and experimenters to complete experimental datasets, explore the relationships between pairs of MN properties never concurrently observed in previous experiments, or investigate inter-mammalian-species variations in MN membrane properties. Using this mathematical framework, modellers can build profiles of inter-consistent MN-specific properties to scale pools of MN models, with consequences on the accuracy and the interpretability of the simulations.

Editor's evaluation

This article will be of interest to scientists studying the contribution of motoneuron behaviour to motor control as well as anyone interested in the relation between neuron morphology, intrinsic properties, and neuron behaviour. The authors have distilled decades of research on motoneuron properties into a set of mathematical relationships that can guide both experimentalists and modellers interested in developing realistic models of populations of motoneurons. In fact, Caillet et al. present a data-driven regression analysis to infer the relationships between morphological and electrophysiological measures from spinal motor neurons in different animal species. Finally, the authors emphasize the value of this approach, but also carefully consider its limitations, including inter-study variability and limited sample sizes in the experimental datasets used to derive the relationships between multiple intrinsic properties.

https://doi.org/10.7554/eLife.76489.sa0

eLife digest

Muscles receive their instructions through electrical signals carried by tens or hundreds of cells connected to the command centers of the body. These ‘alpha-motoneurons’ have various sizes and electrical characteristics which affect how they transmit signals. Previous experiments have shown that these properties are linked; for instance, larger motoneurons transfer electrical signals more quickly. The exact nature of the mathematical relationships between these characteristics, however, remains unclear. This limits our understanding of the behaviour of motoneurons from experimental data.

To identify the equations linking eight motoneuron properties, Caillet et al. analysed published datasets from experimental studies on cat motoneurons. This approach uncovered simple mathematical associations: in fact, only one characteristic needs to be measured experimentally to calculate all the other properties. The relationships identified were also consistent with previously accepted approaches for modelling motoneuron activity. Caillet et al. then validated this mathematical framework with data from studies on rodents, showing that some of the equations hold true for different mammals.

This work offers a quick and easy way for researchers to calculate the characteristics of a motoneuron based on a single observation. This will allow non-measured properties to be added to experimental datasets, and it could help to uncover the diversity of motoneurons at work within a population.

Introduction

Assessing the morphological and electrophysiological properties of individual spinal alpha-motoneurons (MNs) is crucial for understanding their recruitment and discharge behaviour, and exploring the neuromechanical interplay that controls voluntary motion in mammals. As measured in individual spinal MNs in experimental studies and reported in review papers (Henneman, 1981; Burke, 1981; Binder et al., 1996; Powers and Binder, 2001; Olinder et al., 2006; Heckman and Enoka, 2012), significant correlations exist in mammalian MN pools between the MN morphological and/or electrophysiological properties reported in Table 1. For example, the soma size (Dsoma) and the current threshold for spike initiation (Ith) of an MN are both positively correlated to the axonal conduction velocity (ACV) and afterhyperpolarization duration (AHP). Also, Ith decreases with increasing input resistance (R) and AHP, R is negatively correlated to ACV and Dsoma, and ACV varies inversely with AHP. Besides supporting and sequentially leading to extensions of the Henneman’s size principle (Henneman, 1957; Wuerker et al., 1965; Henneman et al., 1965a; Henneman et al., 1965b; Henneman et al., 1974; Henneman, 1981; Henneman, 1985), these empirical associations find strong consistency with Rall’s theoretical approach of representing the soma and the dendritic trees of MNs as an equivalent cylinder (Rall, 1977; Rall et al., 1992; Powers and Binder, 2001). In such case, MNs behave like resistive-capacitive (RC) circuits, as successfully simulated with computational RC models (Izhikevich, 2004; Dong et al., 2011; Negro et al., 2016; Teeter et al., 2018), that rely on further property associations, such as the MN membrane time constant τ precisely equalling the product of the MN membrane-specific capacitance Cm and resistivity Rm (τ=RmCm).

Table 1
The motoneuron (MN) and muscle unit (mU) properties investigated in this study with their notations and SI base units.

SMN is the size of the MN. As reproduced in Table 2, the MN size SMN is adequately described by measures of the MN surface area Sneuron and the soma diameter Dsoma. R and Rm define the MN-specific electrical resistance properties of the MN and set the value of the MN-specific current threshold Ith (Binder et al., 1996; Powers and Binder, 2001; Heckman and Enoka, 2012). C and Cm (constant among the MN pool) define the capacitance properties of the MN and contribute to the definition of the MN membrane time constant τ (Gustafsson and Pinter, 1984a; Zengel et al., 1985). Vth is the amplitude of the membrane voltage depolarization threshold relative to resting state required to elicit an action potential. Ith is the corresponding electrical current causing a membrane depolarization of Vth. AHP duration is defined in most studies as the duration between the action potential onset and the time at which the MN membrane potential meets the resting state after being hyperpolarized. ACV is the axonal conduction velocity of the elicited action potentials on the MN membrane. SmU is the size of the mU. As indicated in Table 2, the mU size SmU is adequately described by measures of (1) the sum of the cross-sectional areas (CSAs) of the fibres composing the mU CSAtot, (2) the mean fibre CSA CSAmean, (3) the innervation ratio IR, that is, the number of innervated fibres constituting the mU, and (4) the mU tetanic force Ftet. Ftw is the MU twitch force.

PropertiesNotationUnit
MN propertiesSize:
Neuron surface area
Soma diameter
SMN
Sneuron
Dsoma
[m2]
[m]
ResistanceRΩ
Specific resistance per unit areaRm[Ωm2]
CapacitanceC[F]
Specific capacitance per unit areaCm[Fm2]
Time constantτ[s]
Rheobase (current recruitment threshold)Ith[A]
Voltage thresholdVth[V]
Afterhyperpolarization durationAHP[s]
Axonal conduction velocityACV[ms1]
mU propertiesSize:SmU
Total fibre cross-sectional areaCSAtot[m2]
Mean fibre cross-sectional areaCSAmean[m2]
Innervation ratioIR[]
Tetanic forceFtet[N]
Twitch forceFtw[N]
Table 2
Measurable indices of motoneuron (MN) and muscle unit (mU) sizes in mammals.

SMN and SmU are conceptual parameters which are adequately described by the measurable and linearly inter-related quantities reported in this table.

MN size (SMN)mU size (SmU)
SneuronFtet
DsomaIR
CSAmean
CSAtot

Yet, because the empirical correlations between MN properties were obtained from scattered data from individual experimental studies, the quantitative mathematical associations between the MN properties reported in Table 1, beyond their aforementioned global relative variations, remain unclear. For example, the negative correlation between R and ACV was described in the literature with linear (Fleshman et al., 1981), exponential (Burke, 1968) or power relationships (Kernell, 1966; Kernell and Zwaagstra, 1981; Ulfhake and Kellerth, 1984; Gustafsson and Pinter, 1984a), while Zwaagstra and Kernell, 1980 reported for this association a slope twice greater than Gustafsson and Pinter, 1984a in the double-logarithmic space. This makes it difficult to reconcile the conclusions of multiple empirical studies that investigated different property associations. For example, Burke, 1968 and Zwaagstra and Kernell, 1980, respectively, proposed exponential R-ACV and linear ACV-AHP associations, suggesting a hybrid exponential relationship between R and AHP, that is different from the power R-AHP relationship directly reported by Ulfhake and Kellerth, 1984. These divergences between studies are a major limitation for our understanding of the associations and the distribution of MN properties in an MN pool, which cannot be directly investigated experimentally either, as measuring multiple MN morphometric properties in vitro and electrophysiological properties in vivo for a large sample of MNs is challenging.

Here, we reprocessed and merged the published data from 19 available experimental studies in adult cat preparations to derive and validate a unique set of mathematical power relationships between all pairs of the MN morphometric and electrophysiological cat properties listed in Table 1. The significance of these quantitative relationships, which are consistent with the conclusions from Rall’s cable theory, supports the notion that the properties reported in Table 1 are all associated to each other. The uniqueness of these mathematical relationships tackles the aforementioned inter-study variability in describing the data and clarifies our understanding of the quantitative association between spinal MN properties, including for pairs of properties never simultaneously measured in experiments. These relationships also provide a convenient mathematical framework for modellers for the derivation of appropriate and consistent MN profiles of MN-specific morphometric and electrophysiological properties for the realistic scaling of pools of computational models of MNs, improving the interpretability of model predictions. Using the relationships, experimenters can readily complete their datasets by deriving MN-specific values, which are representative of the literature, for the MN properties that were not experimentally measured.

After deriving the mathematical framework from cat data, we then demonstrate that the normalized cat relationships apply to other mammals by validating them with data from nine adult rat and mouse electrophysiological studies in vivo. This approach helps better understanding the inter-mammalian-species variations in MN properties. Finally, using additional correlations obtained between some MN and muscle unit (mU) properties from 14 mammal studies, we discuss the empirical relationships obtained between MN properties in the context of the Henneman’s size principle.

Methods

Selected studies reporting processable adult mammal data

Identification of the selected studies

To optimize inter-study consistency, we only selected studies that concurrently measured at least two of the morphometric and/or electrophysiological properties reported in Table 1 in individual spinal alpha-MNs of healthy adult cats, rats, and mice. To build an extensive set of relevant studies so that the mathematical relationships derived in this study describe a maximum of the data published in the literature, the output of the systematic analysis provided in Highlander et al., 2020 was screened and used for a further search by reference lists. Among the retrieved studies, larval and postnatal specimens were disregarded because of the pronounced age-related variance (Highlander et al., 2020) in the mean values of the electrophysiological properties listed in Table 1, which corroborates the non-extrapolability of neonatal neuronal circuitry to older ages (Song et al., 2006; Carp et al., 2008; Nakanishi and Whelan, 2010; Mitra and Brownstone, 2012). Due to the nonlinear inter-species scalability of the spinal alpha-MN electrophysiological and morphometric properties (Manuel et al., 2019) and a pronounced variance in the inter-species mean values (Highlander et al., 2020), non-mammalian species were also disregarded, while the retained data from cats, rats, and mice were processed separately. Highlander et al., 2020 also report an ‘unexplained’ variance in the mean electrophysiological property values reported between studies investigating specimens of the same species, sex, and age, which may be explained by the differences between in vivo and in vitro protocols (Carp et al., 2008). For this reason, only studies measuring the electrophysiological properties in vivo were considered, while most in vitro studies had already been disregarded at this stage as dominantly performed on neonatal specimens due to experimental constraints (Manuel et al., 2009; Mitra and Brownstone, 2012). As all but two of the remaining selected studies focussed on lumbosacral MNs, the final set of considered publications was constrained to MNs innervating hindlimb muscles to improve inter-study consistency. From a preliminary screening of the final set of selected studies, it was finally found that relatively more cat studies (19) were obtained than rat (9) and mouse (4) studies, while the cat studies also investigated more pairs of morphometric and electrophysiological properties. The mathematical relationships sought in this article between MN properties were thus derived from the cat data reported in the selected studies, and then validated for extrapolation to rat and mouse data.

Selected studies providing cat data: Experimental approaches

The 19 selected studies focussing on cats that were included in this work were published between 1966 and 2001. They applied similar experimental protocols to measure the morphometric and electrophysiological properties reported in Table 1. All the selected studies performing morphometric measurements injected in vitro the recorded MNs intracellularly with horseradish peroxidase (HRP) tracer by applying a continuous train of anodal current pulses. The spinal cord was eventually sliced frozen or at room temperature with a microtome or a vibratome and the MN morphometry was investigated. All morphometric measurements were manually performed, and the MN compartments (soma, dendrites, axon) approximated by simple geometrical shapes, yielding some experimental limitations discussed later in ‘Methods’ and Appendix 1.

To measure the electrophysiological properties, animals were anaesthetized, immobilized, and kept under artificial breathing. Hindlimb muscles of interest were dissected free, and their nerves were mounted onto stimulating electrodes, while maintaining body temperature between 35 and 38°C. After a laminectomy was performed over the lumbosacral region of the spinal cord, the MNs were identified by antidromic invasion following electrical stimulation of the corresponding muscle nerves. All selected cat studies reported the use of electrodes having stable resistances and being able to pass currents up to 10 nA without evident polarization or rectification. Some aspects of the experimental protocol however diverged between the selected studies, such as the age and the sex of the group of adult cats, the size population of cats and recorded MNs, the muscles innervated by the recorded MNs, the means of anaesthesia, the level of oxygen, and whether the spinal dorsal roots were severed (Kernell, 1966; Burke, 1968; Barrett and Crill, 1974; Fleshman et al., 1981; Gustafsson and Pinter, 1984a; Gustafsson and Pinter, 1984b; Zengel et al., 1985; Foehring et al., 1987), the complete surrounding hindlimbs were denervated (Burke, 1968; Kernell and Zwaagstra, 1981), or the spinal cord was transected (Burke, 1968; Krawitz et al., 2001).

In all the selected cat studies, ACV was calculated as the ratio of the conduction distance and the antidromic spike latency between the stimulation site and the ventral root entry. The studies however did not report how the nerve length was estimated. AHP was measured using brief suprathreshold antidromic stimulations as the time interval between the spike onset to when the voltage deflection of the hyperpolarizing phase returned to the prespike membrane potential. In all studies, Ith was obtained as the minimal intensity of a 50–100-ms-long depolarizing rectangular current pulse required to produce regular discharges. Ith was obtained with a trial-and-error approach, slowly increasing the intensity of the pulse. The selected studies did not report any relevant source of inaccuracy in measuring ACV, AHP, and Ith. The studies measured R using the spike height method (Frank and Fuortes, 1956). In brief, small (1–5 nA) depolarizing and/or hyperpolarizing current pulses (Iin) of 15–100 ms duration were injected through a Wheatstone bridge circuit in the intracellular microelectrode amplifier, and the subsequent change in the amplitude of the membrane voltage potential (Vm) was reported in steady-state Vm-Iin plots. In all studies, R was obtained by calculating the slope of the linear part of the Vm-Iin plots near resting membrane potential by analogy with Ohm’s law. It is worth noting for inter-study variability that depolarizing pulses return slightly higher R values than hyperpolarizing currents (Sasaki, 1991). To measure the membrane time constant τ, all the selected studies analysed the transient voltage responses (Vm) to weak (1–12 nA), brief (0.5 ms), or long (15–100 ms) hyperpolarizing current pulses. Both brief- and long-pulse approaches were reported to provide similar results for the estimation of τ (Ulfhake and Kellerth, 1984). Considering MNs as equivalent isopotential cables, the membrane time constant τ was identified in all studies as the longest time constant when modelling the measured voltage transient response to a current pulse as the sum of weighted time-dependent exponential terms (Rall, 1969; Rall, 1977; Powers and Binder, 2001). Semi-logarithmic plots of the time history of the voltage transient Vm or of its time derivative dVmdt were drawn, and a straight line was fit by eye to the linear tail of the resulting plot (Fleshman et al., 1988), the slope of which was τ. A graphical ‘peeling’ process (Rall, 1969) was undertaken to recover the first equalizing time constant τ1, required to estimate the electronic length (L) of Rall’s equivalent cylinder representation of the MN and the MN membrane capacitance with Rall’s equations (Rall, 1969; Rall, 1977; Gustafsson and Pinter, 1984b; Powers and Binder, 2001):

Lπττ1-112
C=τRLtanhL

These electrophysiological measurements were reported to be subject to three main experimental sources of inaccuracy. First, a variable membrane ‘leak’ conductance, which can be estimated from the ratio of two parameters obtained from the ‘peeling’ process (Gustafsson and Pinter, 1984b), arises from the imperfect seal around the recording micropipette (Ulfhake and Kellerth, 1984; Gustafsson and Pinter, 1984b; Pinter and Vanden Noven, 1989). As reviewed in Powers and Binder, 2001 and Olinder et al., 2006, this ‘leak’ can affect the measurements of all the properties, notably underestimating the values of R and τ and overestimating those of C. However, this ‘leak’ is probably not of major significance in the cells of large spike amplitudes (Gustafsson and Pinter, 1984b), on which most of the selected studies focussed. Importantly, the membrane ‘leak’ is also not expected to affect the relations between the parameters (Gustafsson and Pinter, 1984b).

Secondly, some nonlinearities (Ito and Oshima, 1965; Burke and ten Bruggencate, 1971; Ulfhake and Kellerth, 1984) in the membrane voltage response to input current steps arise near threshold because of the contribution of voltage-activated membrane conductances to the measured voltage decay (Fleshman et al., 1988). This contradicts the initial assumption of the MN membrane remaining passive to input current steps (Rall, 1969; Burke and ten Bruggencate, 1971). Because of this issue, the ends of the Vm-Iin plot are curvilinear and can affect the estimation of R, while the transient voltage response to an input current pulse decays faster than exponentially (it undershoots) at the termination of the current injection, which makes it impossible to plot the entire course of lnVt and challenges the graphical procedure taken to estimate τ (Fleshman et al., 1988). Therefore, some of the selected studies discarded all the MNs that displayed obvious large nonlinearity in the semi-logarithmic V–t or dVdt -t plots (Burke and ten Bruggencate, 1971). The membrane nonlinearities can be corrected by adding a constant voltage to the entire trace (Fleshman et al., 1988) or by considering the three time-constant model of Ito and Oshima, 1965 to approximate and subtract to the voltage signal the membrane potential change produced by the current steps. However, none of the selected studies performed such corrections, potentially yielding a systematic underestimation of the values of R and τ (Gustafsson and Pinter, 1984b; Zengel et al., 1985; Powers and Binder, 2001). Yet, this systematic error for R, which should not contribute to inter-study variability, is expected to remain low because the selected studies used input current pulses of low strength (Ulfhake and Kellerth, 1984) and measured R from the linear part of the I–V plots near the resting membrane potential where no membrane nonlinearity is expected to occur. This systematic error may also not affect the distribution of recorded R values, as displayed in Zengel et al., 1985. Zengel et al., 1985 besides corrected for the membrane nonlinearities with the three time-constant model approach of Ito and Oshima, 1965 when measuring τ and returned values of τ around 40% higher than the other selected studies. However, this systematic discrepancy is expected to disappear with the normalization of the datasets described in the following, as the reported normalized distributions of τ are very similar between the selected studies (see density histogram for the τ;R dataset in Appendix 1—figure 1). Thirdly, the highest source of inaccuracy and inter-study variability arises from the subjective fit ‘by eye’ of a straight line to the transient voltage when estimating τ (Pinter and Vanden Noven, 1989).

Overall, all the selected studies used similar experimental approaches to measure ACV, AHP, Ith, R, C, and τ, and little sources of inaccuracy and inter-study variability were identified.

Relationships between MN properties

For convenience, in the following, the notation A;B refers to the pair of morphometric or electrophysiological MN properties A and B, with A,BSMN;ACV;AHP;R;Ith;C; τ, defined in Table 1. The selected studies generally provided clouds of data points for pairs A;B of concurrently measured MN properties through scatter graphs. These plots were manually digitized using the online tool WebPlotDigitizer (Ankit, 2020). When a study did not provide such processable data – most reported the mean ± sd property values of the cohort of measured MNs – the corresponding author was contacted to obtain the raw data of the measured MN properties, following approval of data sharing. Upon reception of data, datasets of all possible A;B pairs were created and included into the study.

Normalized space and choice of regression type

The sets of data retrieved from different cat studies for each property pair A;B were merged into a ‘global’ dataset dedicated to that property pair. The A;B data was also normalized for each study and transformed as a percentage of the maximum property value measured in the same study and normalized ‘global’ datasets were similarly created. A least-squares linear regression analysis was performed for the lnA-lnB transformation of each global dataset yielding relationships of the type lnA=alnB+k, which were then converted into power relationships of the type A=kBa (eq. 1). The adequacy of these global power trendlines and the statistical significance of the correlations were assessed with the coefficient of determination r2 (squared value of the Pearson’s correlation coefficient) and a threshold of 0.01 on the p-value of the regression analysis, respectively. For each A;B pair, the normalized global datasets systematically returned both a higher r2 and a lower p-value than the datasets of absolute A;B values, in agreement with the ‘unexplained’ inter-experimental variance reported in Highlander et al., 2020. It was therefore decided to investigate the A;B pairs of MN properties in the normalized space using the normalized global datasets to improve the cross-study analysis. For this preliminary analysis and the rest of the study, power regressions (eq. 1) were preferred to linear (A=k+aB) or exponential (A=keaB) regressions to maintain consistency with the mathematical structure of the equations from Rall’s cable theory (Rall, 1957; Rall, 1959; Rall, 1960; Powers and Binder, 2001) and of other well-defined relationships, such as R=kIth1. Also, a least-squares linear regression analysis was preliminary performed for each normalized experimental dataset and for its lnA-lnB and lnA-B transformations, yielding linear, power, and exponential fits to the data. The power regressions overall returned r2>0.5 for more experimental datasets than the linear and exponential fits (see Appendix 1—table 1). To avoid bias, power regressions were the only type of fit used in this study, despite a few datasets being more accurately fitted by linear or exponential regressions (the difference in the quality of the fit was very small in all cases). Other regression types, such as polynomial fits, were not justified by previous findings in the literature.

Global datasets and data variability

The data variability between the studies constituting the same global dataset A;B was assessed by analysing the normalized distributions of the properties A and B with four metrics, which were the range of the measured data, the mean of the distribution, the coefficient of variation (CoV) calculated as the standard deviation (sd) divided by the mean of the distribution, and the ratio MeMd of the median and the mean of the distribution. The inter-study data variability in the A;B global dataset was assessed by computing the meang±sdg across studies of each of the four metrics. Low variability between the data distributions from different studies was concluded if the normalized distributions (in percentage maximum value) returned similar values for the range, mean, CoV, and MeMd metrics, that is, when sdg<10% and sdgmeang<0.15 were obtained for all four metrics. In such case, the data distributions would respectively span over similar bandwidth length of the MU pool, be centred around similar mean values, and display similar data dispersion and similar skewness. The relative size of the experimental datasets reported by the selected studies constituting the same global datasets was also compared to assess their relative impact on the regression curves fitted to the global datasets. Then, the inter-study variability in associating the distributions of properties A and B was assessed by computing the 95% confidence interval of the linear model fitted to the global dataset A;B in the log–log space, which yields the power regression A=kBa. Low inter-study variability was considered if the value of a varied less than 0.4 in the confidence interval.

The variability of the data distribution of a property A between multiple global datasets A;B was then assessed by computing the same four metrics as previously (range, mean, CoV, MeMd) for each global dataset, and then computing the meanG±sdG of each of the four metrics across the global datasets A;B. Low variability between the global datasets in the distribution of property A was considered if sdG<10% and sdGmeanG<0.15 were obtained for all four metrics. If the inter-study and inter-global dataset variability was low, the global datasets were created and processed to derive the mathematical relationships between MN properties according to the procedure described below.

Size-dependent normalized relationships

From inspection of the considered cat studies, most of the investigated MN property pairs comprised either a direct measurement of MN size, noted as SMN in this study, or another variable well-admitted to be strongly associated to size, such as ACV, AHP, or R. Accordingly, and consistently with Henneman’s size principle, we identified SMN as the reference MN property with respect to which relationships with the electrophysiological MN properties in Table 1 were investigated. To integrate the data from all global normalized datasets into the final relationships, the MN properties in Table 1 were processed in a step-by-step manner in the order ACV, AHP, R, Ith,C,τ, as reproduced in Figure 1 for two arbitrary properties, to seek a ‘final’ power relationship of the type eq.1 between each of them and SMN. For each electrophysiological property A and each A;B dataset, we considered two cases. If B=SMN, the global dataset was not further processed as the electrophysiological property A was already related to MN size SMN from experimental measurements. If BSMN, the global A;B dataset was transformed into a new A;SMN dataset by converting the discrete values of B with the trendline regression SMN=kdBd, which was obtained at a previous step of the data processing. With this dual approach, as many MN size-dependent A;SMN intermediary datasets as available A;B global datasets were obtained for each property A. These size-dependent datasets were merged into a ‘final’ A;SMN dataset to which a least-squares linear regression analysis was performed for the lnA-lnSMN transformation, yielding relationships of the type lnA=clnSMN+kc, which were converted into the power relationships:

A=kcSMNc
Detailed example for the process adopted to successively create the two final {R;SMN} and {Ith;SMN} datasets (right-side thick-solid contour rectangular boxes).

These final datasets were obtained from respectively three and three normalized global datasets of experimental data obtained from the literature (dashed-contour grey-filled boxes) {{R;ACV}, {R;AHP}, {R;SMN}} and {{Ith;R}, {Ith;AHP}, {Ith;ACV}}. The {R;ACV} and {R;AHP} datasets were first transformed (⊗ symbol) into two intermediary {R;SMN} datasets (dotted-contour boxes) by converting the ACV and AHP values to equivalent SMN values with two ‘inverse’ SMN=fACV(ACV) and SMN=fAHP(AHP) power relationships (oval boxes with triple dots), which had been previously obtained from two unshown steps that had yielded the final {ACV;SMN} and {AHP;SMN} datasets. The two intermediary {R;SMN} datasets were merged with the remaining global {R;SMN} dataset to yield the final {R;SMN} dataset, to which a power relationship of the form R=kSMNc was fitted. If r2>0.3 and p<0.01, an ‘inverse’ power relationship SMN=fR(R) (oval box) was further fitted to this final dataset. In a similar approach, the three normalized global datasets {Ith;R}, {Ith;AHP}, and {Ith;ACV} were transformed with the three ‘inverse’ relationships into intermediary {Ith;SMN} datasets, which were merged to yield the final {Ith;SMN} dataset. An ‘inverse’ SMN=fIth(Ith) power relationship was further derived to be used in the creation of the final {C;SMN} and {τ;SMN} datasets in the next taken steps.

The adequacy of these final power trendlines and the statistical significance of the correlations were assessed identically to the A;B relationships derived directly from normalized experimental data, using the coefficient of determination r2 (squared value of the Pearson’s correlation coefficient) and a threshold of 0.01 on the p-value of the regression analysis, respectively. If r2>0.3 (i.e., r > 0.55) and p<0.01, another power trendline SMN=kdAd was fitted to the final dataset to describe the inverse relationship for SMN;A and used in the processing of the next-in-line property.

Normalized mathematical relationships between electrophysiological properties

When a final relationship with SMN was obtained for two MN properties A and B by the procedure described above, that is, A=kcASMNcA and B=kcBSMNcB, one of the expressions was mathematically inverted and a third empirical relationship was derived for the property pair A;B :

A=kcAkcBcAcBBcAcB=keBe

This procedure was applied to all possible A;B pairs in Table 1.

Validation of the normalized relationships

The normalized relationships were validated using a standard fivefold cross-validation procedure (Chollet, 2021). The data in each A;B global normalized dataset was initially randomized and therefore made independent from the studies constituting it. Each shuffled global dataset was then split into five non-overlapping complementary partitions, each containing 20% of the data. Four partitions including 80% of the data were taken to constitute a training set from which the normalized relationships were built as described previously. The latter equations were validated against the last data partition, which includes the remaining 20% of the data and which is called test set in the following. To perform this validation, the final size-related relationships A=kcSMNc and relationships between electrophysiological properties A=keBe obtained with the training set were applied to the SMN or B data in each test set, yielding predicted values A. It was then assessed to what extent the mathematical relationships predicted the test data A by calculating the normalized maximum error (nME), the normalized root mean squared error (nRMSE) and the coefficient of determination rpred2 between predicted and control values of A in each test set. This process was repeated for a total of five times by permutating the five data partitions and creating five different pairs of one training set and one test set. The nME, nRMSE, and rpred2 values were finally averaged across the five permutations. For each global dataset, the average rpred2 was compared to the rexp2 obtained from the power trendlines least-squared fitted to the control data in the log–log space. Once validated with this fivefold cross-validation method, the final normalized size-dependent relationships A=kcSMNc and relationships between electrophysiological properties A=keBe were computed for the complete global datasets.

Scaling of the normalized relationships

The A=kcSMNc and A=keBe normalized relationships were finally scaled to typical cat property values in three steps. First, it was assessed from the literature the fold range over which each MN property was to vary. SMN varied over a qSE -fold range, taken as the average across cat studies of the ratios of minimum and maximum values measured for SMN :

qSE=1Nstudiesi=1NstudiesSMNimaxSMNimin

An experimental ratio qAE was similarly obtained for each electrophysiological property A. As A=kcSMNc, any electrophysiological property A could theoretically vary over a qAT=qSEc -fold range, which was compared to qAE for consistency in the results. Then, a theoretical range Amin; Amax of values were derived for each electrophysiological property A. Defining AminE and AmaxE as the respective average of the minimum and maximum A-values across studies, it was enforced AmaxAmin=qAT and Amin+Amax2=AminE+AmaxE2 so that the Amin; Amax theoretical ranges were consistent with the normalized size-dependent relationships while best reproducing the experimental data from cat literature. A theoretical range for SMNSMNmin; SMNmax was similarly built over the qSE -fold range previously derived. Finally, the intercept kc in the size-dependent relationships A=kcSMNc was scaled as

{ kc=Amin(SMN)minc  ;  if c>0 kc=Amax(SMN)minc  ;  if c<0

A similar approach was used to mathematically derive the intercept ke and scale the relationships A=keBe between electrophysiological properties.

Extrapolation to rats and mice

It was finally assessed whether the A=kcSMNc and A=keBe scaled relationships derived from cat data accurately predicted rat and mouse quantities reported in other experimental studies. These mathematical relationships were applied to the SMN or B data in each A;B global dataset of absolute values obtained for rats and mice independently. The accuracy in predicting the quantity A was assessed for each available dataset with the same three validation metrics used to validate the normalized power trend lines: the normalized maximum error (nME), the normalized root-mean-square error (nRMSE) and the coefficient of determination r2 between predicted and experimental values of A.

Definitions of MN size SMN and mU size SmU

MN size SMN

The selected studies that performed both electrophysiological and morphometric measurements on the same MNs dominantly measured the MN soma diameter Dsoma as an index of MN size. Therefore, we chose SMN=Dsoma in the described methodology. In the literature on spinal MNs in adult mammals, the size of an MN SMN is also related to the measures of the somal cross-sectional area CSAsoma (Gao et al., 2009; Friese et al., 2009; Deardorff et al., 2013; Mierzejewska-Krzyżowska et al., 2014; Dukkipati et al., 2018), soma surface area Ssoma (Ulfhake and Kellerth, 1984; Brandenburg et al., 2020), axonal diameter Daxon (Cullheim, 1978), individual dendrite diameter Ddendrite (Amendola and Durand, 2008; Carrascal et al., 2009; Mantilla et al., 2018), individual dendritic surface area Sdendrite (Li et al., 2005; Obregón et al., 2009; Carrascal et al., 2009; Filipchuk and Durand, 2012; Kanjhan et al., 2015), total dendritic surface area (Ulfhake and Cullheim, 1988; Amendola and Durand, 2008; Brandenburg et al., 2020), or total MN surface area Sneuron (Burke et al., 1982) defined as the summed soma and dendritic surface areas.

In the references (Zwaagstra and Kernell, 1981; Ulfhake and Kellerth, 1981), a linear correlation between Dsoma and the average diameter of the stem dendrites Ddendritemean has been reported (r > 0.62, population size in {40; 82} cells). Moreover, a linear or quasi-linear correlation has been found (Cullheim et al., 1987; Ulfhake and Cullheim, 1988; Moschovakis et al., 1991; Prakash et al., 2000; Li et al., 2005; Obregón et al., 2009; Mantilla et al., 2018) between the stem dendrite diameter Ddendrite and the membrane surface area of the corresponding dendritic tree Sdendrite (r > 0.78, population size in 33;342 dendrites). Therefore, from these studies, we can assume an approximate linear correlation between Dsoma and the average dendritic surface area Sdendritemean. Moreover, the number of dendritic trees Ndendrites per cell increases with increasing soma surface Ssoma (Brandenburg et al., 2018), and a linear correlation between Dsoma or Sdendrites with Ndendrites has been observed (Zwaagstra and Kernell, 1981; Ulfhake and Cullheim, 1988; Amendola and Durand, 2008) (r > 0.40, population size in {14; 32; 87} cells). It was therefore assumed that Dsoma and total dendritic surface area Sdendritetot are linearly correlated. This assumption/approximation is consistent other conclusions of a linear correlations between Dsoma and Sdendritetot (Amendola and Durand, 2008; Brandenburg et al., 2020). Then, according to typical values of Dsoma, Ssoma, and Sdendritestot obtained from the studies previously cited, yielding SdendritestotSsoma, it was also assumed Sdendritestot Sneuron. It is thus concluded that Dsoma is linearly related to total neuron membrane area Sneuron, consistently with results by Burke et al., 1982 (r=0.61, 57 cells).

For the above reasons and assumptions, the mathematical relationships A=kcSMNc derived previously, with SMN=Dsoma, were extrapolated with a gain to relationships between A and SMN=Sneuron, notably permitting the definition of surface specific resistance Rm and capacitance Cm. Following the same method as described above, theoretical ranges Sneuronmin; Sneuronmax were derived from additional morphometric studies on adult cat spinal alpha-MNs. The new relationships were extrapolated as A=kcDsomac=SneuronminDneuronminSneuronc.

It must however be highlighted that the conclusion of a linear correlation between Dsoma and Sneuron, while plausible, is crude. Morphometric measurements of individual MNs are indeed difficult and suffer many limitations. MN staining, slice preparations, and MN reconstruction and identification are complex experimental procedures requiring a large amount of work; the results from most of the cited studies thus rely on relatively small pools of investigated MNs. Most cited studies did not account for tissue shrinkage after dehydration, while some may have failed to assess the full dendritic trees from MN staining techniques (Brandenburg et al., 2020), thus underestimating dendritic membrane measurements. Moreover, before automated tools for image segmentation, landmark mapping, and surface tracking (Amendola and Durand, 2008; Obregón et al., 2009; Mierzejewska-Krzyżowska et al., 2014) were available, morphometric measurements were performed from the manual reproduction of the cell outline under microscope, yielding important operator errors reported to be of the order of ~0.5 µm, that is, ~20% stem diameter. In these studies, morphological quantities were besides obtained from geometrical approximations of the 2D MN shape, such as modelling the individual dendritic branches as one (Cullheim et al., 1987) or more (Prakash et al., 2000) equivalent cylinders for the derivation of Sdendrite. Morphometric measurement approaches also varied between studies; for example, oval or circle shapes were best fitted by sight onto the soma outline and equivalent Dsoma quantities were derived either as the mean of measured maximum and minimum oval diameters or as the equivalent diameter of the fitted circle. In these studies, CSAsoma and Ssoma were therefore derived by classical equations of circle and spheric surface areas, which contradicts the results by Mierzejewska-Krzyżowska et al., 2014 obtained from surface segmentation of a linear relationship between Dsoma and CSAsoma (r=0.94, 527 MNs). Finally, no correlation between soma size and Ndendrites was found (Ulfhake and Kellerth, 1981; Cullheim et al., 1987). Therefore, the linear correlation between Dsoma and Sneuron is plausible but crude and requires awareness of several important experimental limitations and inaccuracies. Conclusions and predictions involving Sneuron should be treated with care in this study as these morphometric measurements lack the precision of the measures performed for the electrophysiological properties listed in Table 1.

Muscle unit size SmU

To enable future comparisons between MN and mU properties, we here assess potential indices of the mU size SmU suggested in the literature. The size of an mU (SmU) can be defined as the sum CSAtot of the CSAs of the innervated fibres composing the mU. CSAtot depends on the mU innervation ratio (IR) and on the mean CSA (CSAmean) of the innervated fibres: SmU=CSAtot=IRCSAmean. CSAtot was measured in a few studies on cat and rat muscles, either indirectly by histochemical fibre profiling (Burke and Tsairis, 1973; Dum and Kennedy, 1980; Burke, 1981), or directly by glycogen depletion, periodic acid Schiff (PAS) staining and fibre counting (Burke et al., 1982; Rafuse et al., 1997). The mU tetanic force Ftet is however more commonly measured in animals. As the fibre mean specific force σ is considered constant among the mUs of one muscle in animals (Lucas et al., 1987; Enoka, 1995), the popular equation Ftet= σIRCSAmean (Burke, 1981; Enoka, 1995) returns a linear correlation between Ftet and IRCSAmean=SmU in mammals. Experimental results (Burke and Tsairis, 1973; Bodine et al., 1987; Chamberlain and Lewis, 1989; Kanda and Hashizume, 1992; Rafuse et al., 1997; Hegedus et al., 2007) further show a linear correlation between Ftet and IR and between Ftet and CSAmean. Consequently, Ftet, IR, and CSAmean are measurable, consistent, and linearly related measures of SmU in animals, as summarized in Table 2.

Relationships between MN and mU properties

To assess whether the size-dependent relationships A=kcSMNc derived in this study were in accordance with Henneman’s size principle of motor unit recruitment, we identified a set of experimental studies that concurrently measured an MN property BMN (Table 1) and an mU property AmU for the same MU. The normalized global datasets obtained for the pairs AmU;BMN were fitted with power trendlines, as previously described for MN properties, yielding AmU=kBMNb relationships. Using both the definition of SmU (Table 2) and the A=kcSMNc relationships derived previously, the AmU=kBMNb relationships were then mathematically transformed into SmU=kSMNc relationships. If all c-values obtained from different AmU;BMN pairs were of the same sign, it was concluded that mU and MN sizes were monotonically related. Considering the limited data available obtained from the literature, data obtained for MNs innervating different hindlimb muscles were merged. Also, cat, rat, and mouse studies were processed independently but the resulting c-values were compared without regards to the related species.

Results

We respectively identified 19, 6, and 4 studies respecting our desired criteria on cats, rats, and mice that reported processable experimental data for the morphometric and electrophysiological properties listed in Table 1. Additional publications including some from the past 10 years were identified but could not be included in this work as no processable data could be recovered. An additional 14 studies were found to perform concurrent MN and mU measurements on individual motor units. From the selected cat studies, the 17 pairs of MN properties and the 8 pairs of one MN and one mU property represented in the bubble diagram of Figure 2A were investigated.

Experimental (A) and unknown (B) relations between motoneuron (MN) and muscle unit (mU) properties.

(A) Bubble diagram representing the pairs of MN and/or mU properties that could be investigated in this study from the results provided by the 40 studies identified in our web search. MN and mU properties are represented by circle and square bubbles, respectively. Relationships between MN properties are represented by coloured connecting lines; the colours red, blue, green, yellow, and purple are consistent with the order ACV, AHP, R, Ith,C,τ in which the pairs were investigated (see Table 3 for mathematical relationships). Relationships between one MN and one mU property are represented by black dashed lines. (B) Bubble diagram representing the mathematical relationships proposed in this study between pairs of MN properties for which no concurrent experimental data has been measured to date.

Table 3
Fitted experimental data of pairs of motoneuron (MN) properties and subsequent normalized final size-related relationships.

For information, the r2, p-value, and the equation A=kaBa are reported for each fitted global dataset. The normalized MN-size dependent relationships A=kcSMNc are mathematically derived from the transformation of the global datasets and from the power trendline fitting of the final datasets (N data points) as described in ‘Methods’. The minimum and maximum values of ka, kc, a, and c defining the 95% confidence interval of the regression are also reported in parenthesis for each global and final dataset. The r2 values reported in this table are consistent with the r2 values obtained when directly fitting the normalized experimental datasets with power regressions (see Appendix 1—table 1).

MN propertyA=kaBa(normalized global datasets)A=kcSMNc(final MN-size dependent datasets)
ARelationshipkaar2p-valueReference studieskccr2p-valueN points
ACVACV=kaSMNa4.0
(2.5; 6.4)
0.7
(0.6; 0.8)
0.58< 10-5Cullheim, 1978; Kernell and Zwaagstra, 1981; Burke et al., 19824.0
(2.5; 6.4)
0.7
(0.6; 0.8)
0.58< 10-5109
AHPAHP=kaSMNa6.1 · 103
(1.2 · 103; 3.2 · 104)
−1.2
(−1.6; −0.8)
0.34< 10-5Zwaagstra and Kernell, 19802.5 · 104
(1.2 · 104; 5.0 · 104)
−1.5
(−1.7; −1.3)
0.41< 10-5492
AHP=kaACVa1.5 · 104
(7.4 · 103; 2.9 · 104)
−1.4
(−1.5; −1.2)
0.41< 10-5Eccles et al., 1958a; Zwaagstra and Kernell, 1980; Gustafsson and Pinter, 1984b; Foehring et al., 1987
RR=kaSMNa1.5 · 105
(2.7 · 104; 7.9 · 105)
−2.1
(−2.5; −1.7)
0.61< 10-5Kernell and Zwaagstra, 1981; Burke et al., 19829.6 · 105
(4.1 · 105; 2.3 · 106)
−2.4
(−2.6; −2.2)
0.37< 10-5745
R=kaACVa6.3 · 105
(1.9 · 105; 2.1 · 106)
−2.3
(−2.6; −2.0)
0.38< 10-5Kernell, 1966; Burke, 1968; Barrett and Crill, 1974; Kernell and Zwaagstra, 1981; Fleshman et al., 1981; Gustafsson and Pinter, 1984b; Sasaki, 1991
R=kaAHPa6.2 · 10−1
(4.1 · 10−1; 9.2 · 10−1)
1.1
(0.9; 1.2)
0.65< 10-5Gustafsson, 1979; Gustafsson and Pinter, 1984b; Foehring et al., 1987; Pinter and Vanden Noven, 1989; Sasaki, 1991
IthIth=kaRa1.1 · 103
(0.8 · 103; 1.3 · 103)
−1.0
(−1.1; −0.9)
0.37< 10-5Kernell, 1966; Fleshman et al., 1981; Gustafsson and Pinter, 1984a; Zengel et al., 1985; Munson et al., 1986; Foehring et al., 1987; Krawitz et al., 20019.0 · 10−4
(4.7 · 10−4; 1.7 · 10−3)
2.5
(2.4; 2.7)
0.37< 10-5722
Ith=kaACVa3.2 · 10−6
(1.3 · 10−7; 8.2 · 10−5)
3.7
(3.0; 4.4)
0.37< 10-5Kernell and Monster, 1981; Gustafsson and Pinter, 1984a
Ith=kaAHPa2.5 · 104
(1.3 · 104; 4.8· 104)
−1.7
(−1.9; −1.6)
0.60< 10-5Gustafsson and Pinter, 1984a
CC=kaRa2.4 · 102
(2.0 · 102; 3.9 · 102)
−0.4
(−0.4; −0.3)
0.57< 10-5Gustafsson and Pinter, 1984b1.2
(0.7; 2.0)
1.0
(0.9; 1.2)
0.28< 10-5444
C=kaItha2.9 · 101
(2.4 · 101; 3.5 · 101)
0.3
(0.2; 0.3)
0.51< 10-5Gustafsson and Pinter, 1984a
C=kaAHPa2.8 · 102
(1.8 · 102; 4.4 · 102)
−0.4
(−0.5; −0.3)
0.24< 10-5Gustafsson and Pinter, 1984b
C=kaACVa2.5
(0.7; 8.4)
0.8
(0.5; 1.0)
0.17< 10-5Gustafsson and Pinter, 1984b
ττ=kaRa8.7
(7.2; 10.6)
0.5
(0.4; 0.6)
0.52< 10-5Burke and ten Bruggencate, 1971; Barrett and Crill, 1974; Gustafsson, 1979; Gustafsson and Pinter, 1984b; Zengel et al., 1985; Pinter and Vanden Noven, 1989; Sasaki, 19912.6 · 104
(1.5 · 104; 4.5 · 104)
−1.5
(−1.6; −1.4)
0.46< 10-5649
τ=kaAHPa2.2
(1.3; 3.5)
0.8
(0.7; 1.0)
0.63< 10-5Gustafsson and Pinter, 1984b
τ=kaItha2.3 · 102
(1.9 · 102; 2.7 · 102)
−0.4
(−0.5; −0.3)
0.72< 10-5Gustafsson and Pinter, 1984a
τ=kaACVa1.2 · 104
(2.2 · 103; 6.6 · 104)
−1.3
(−1.7; −0.9)
0.30< 10-5Gustafsson and Pinter, 1984b

Relationships between MN properties

Global datasets and data variability

The experimental data retrieved from the selected studies for the 17 pairs of MN properties drawn in Figure 2A were merged into 17 normalized global datasets plotted in Figure 3. Eight global datasets ACV;SMN, AHP;ACV, R;SMN, R;ACV, R;AHP, Ith;R, Ith;ACV, τ;R included the data from two or more experimental studies.

Normalized global datasets.

These were obtained from the 19 studies reporting cat data that measured and investigated the 17 pairs of motoneuron (MN) properties reported in Figure 2A. For each {A;B} pair, the property A is read on the y-axis and B on the x-axis. For information, power trendlines A=kBa (red dotted curves) are fitted to the data of each dataset and reported in Table 3. The 95% confidence interval of the regression is also displayed for each dataset (green dotted lines). The studies are identified with the following symbols: • (Gustafsson, 1979; Gustafsson and Pinter, 1984a; Gustafsson and Pinter, 1984b), ○ (Munson et al., 1986), ▲ (Zengel et al., 1985), ∆ (Foehring et al., 1987), ■ (Cullheim, 1978), □ (Burke, 1968; Burke and ten Bruggencate, 1971; Burke et al., 1982), ◆(Krawitz et al., 2001), ◇ (Fleshman et al., 1981), + (Eccles et al., 1958b), ☓ (Kernell, 1966; Kernell and Zwaagstra, 1981; Kernell and Monster, 1981), - (Zwaagstra and Kernell, 1980), — (Sasaki, 1991), ✶ (Pinter and Vanden Noven, 1989). The axes are given in % of the maximum retrieved values in the studies consistently with ‘Methods’ section.

These studies reported in each global dataset similar normalized property distributions, as visually displayed in the frequency histograms in Appendix 1—figure 1. This was confirmed by the meang±sdg calculations described in ‘Methods’ of the four metrics range, mean, CoV, and MeMd displayed as error bars in Appendix 1—figure 2. Out of the 64 meang +/- sdg calculations (4 metrics for 16 property distributions of the 8 global datasets), 55 (86%) returned sdg<10% and sdgmeang<0.15. Otherwise, the distributions of only three, two, and four properties showed sensibly higher inter-study variability in the global datasets for the range, mean, and CoV metrics, respectively, however, still verifying sdg<15% and sdgmeang<0.33. Additional details are provided in Appendix 1. Then, as displayed in Figure 3 and reported in Table 3, power trendlines and their 95% confidence interval were fitted to the global datasets. All 17 trendlines were statistically significant and adequately represented by a power model A=kaBa (p-value<10-5 and r20.34;0.72 for 15 datasets; p-value<10-5 and r20.24;0.17 for the C;AHP and {C;ACV} datasets, respectively). As displayed in green dotted lines in Figure 3, the confidence interval remained narrow for all datasets with the value of power a varying less than 0.4, except for the Ith;ACV dataset, suggesting a low inter-study variability in associating the distributions of two properties. As displayed in Appendix 1—figure 3, the selected studies however reported datasets of different sizes. Yet, out of the 35 experimental datasets constituting the 8 global datasets identified previously, only 1 and 7 experimental datasets were identified to respectively constitute more than 50% and less than 10% of the global dataset they were included in. Therefore, despite a general imbalance towards the over-appearance of the work from Gustafsson and Pinter, 1984a; Gustafsson and Pinter, 1984b and a tendency towards overlooking some small datasets (Kernell, 1966; Burke and ten Bruggencate, 1971; Sasaki, 1991), the remaining studies all played a significant role in the procedure to constitute the final datasets. Therefore, the inter-study data variability was globally low, and the data distributions reported in the experimental studies were confidently merged into the global datasets plotted in Figure 3. It is worth noting that three research groups provided 40, 27, and 15% of the 2717 data points describing the SMN, ACV, AHP, R, Ith, and τ property distributions in this study, while the C-related data was provided by one only research group, leading to potential group-specific methodological bias.

The merged property distributions obtained in the 17 global datasets showed low variability between global datasets, as displayed in Appendix 1—figure 4. Indeed, sdG<10% and sdGmeanG<0.15 were obtained for all four metrics for all the investigated properties across the global datasets they appeared in. Therefore, the property distributions were similar enough between global datasets to apply the process displayed in Figure 1 and derive the final size-dependent datasets.

Size-dependent normalized relationships

The 17 normalized global datasets were processed according to the procedure described in Figure 1 to derive normalized mathematical relationships between the electrophysiological properties listed in Table 1 and MN size SMN. First, the normalized size-dependent relationship ACV=4.0SMN0.7, reported in Table 3, was derived from the trendline fitting of the normalized ACV;SMN global dataset, which was obtained from three studies (Cullheim, 1978; Kernell and Zwaagstra, 1981; Burke et al., 1982) and is represented in the upper-left panel in Figure 3. This resulted in a statistically significant relationship (r2=0.58, p-value<10-5). A statistically significant (r2=0.59, p-value<10-5) inverse relationship SMN=kdACVd was also fitted to this dataset. Then, the AHP;ACV dataset, represented in the second panel upper row in Figure 3 and obtained from four studies (Eccles et al., 1958a; Zwaagstra and Kernell, 1980; Gustafsson and Pinter, 1984b; Foehring et al., 1987), was transformed into a new AHP;SMN1 dataset by applying the priory-derived SMN=kdACVd relationship to the list of ACV values. This transformed AHP;SMN1 dataset was then merged with the AHP;SMN2 dataset (third panel upper row in Figure 2) obtained from Zwaagstra and Kernell, 1980, yielding the final AHP;SMNf dataset of N=492 data points shown in Figure 4, second panel. A statistically significant (r2=0.58, p-value<10-5) power trendline AHP=2.5104SMN-1.5 was fitted to the AHP;SMNf dataset and is reported in Table 3. As before, a statistically significant (r2=0.38, p-value<10-5) inverse relationship SMN=kdAHPd was also fitted to this dataset for future use. A similar procedure was applied to derive the normalized final relationships between R, Ith, C,τ and SMN reported in the last column of Table 3. A statistically significant (p-value<10-5) correlation was obtained between each electrophysiological property ACV, AHP, R, Ith, C,τ and SMN as reported in Table 3. With r20.28;0.58, it was obtained that {Ith, C, ACV} and R, AHP, τ respectively increased and decreased with increasing MN sizes SMN, with slopes reported in Table 3 and plotted in Figure 5.

Normalized motoneuron (MN) size-related final datasets.

These were obtained from the 19 studies reporting cat data that concurrently measured at least two of the morphometric and electrophysiological properties listed in Table 1. For each {A;SMN} pair, the property A is read on the y-axis and SMN on the x-axis. The power trendlines A=kcSMNc (red dotted curves) are fitted to each dataset and are reported in Table 3. The 95% confidence interval of the regression is also displayed for each dataset (green dotted lines). For each {A;SMN} plot, the constitutive sub-datasets {A;SMN} that were obtained from different global {A;B} datasets are specified with the following symbols identifying the property B: SMN  Δ, ACV ×, AHP □, R  +, and Ith.

Normalized size-dependent behaviour of the motoneuron (MN) properties ACV, AHP, R, Ith,C, and τ.

For displaying purposes, the MN properties are plotted in arbitrary units as power functions (intercept k=1) of SMN : A=SMNc according to Table 3. The larger the MN size, the larger ACV, C, and Ith in the order of increasing slopes, and the lower AHP, τ and R in the order of increasing slopes.

Normalized mathematical relationships between electrophysiological properties

The final size-dependent relationships reported in Table 3 were mathematically combined to yield normalized relationships between all electrophysiological properties listed in Table 1. As an example, Ith=9.410-4SMN2.5 and R=9.6105SMN-2.4 were obtained in Table 3. When mathematically combined (with the latter mathematically inverted as SMN=2.9102R-0.4), these relationships yielded the normalized relationship between rheobase and input resistance Ith=1.5103R-1. All normalized relationships between MN electrophysiological properties hence obtained are provided in Appendix 1—table 2.

Validation of the normalized relationships

The normalized relationships reported in Appendix 1—table 2 were obtained from the complete global datasets represented in Figure 3. A fivefold cross-validation of these relationships was performed, as described in ‘Methods’, and assessed with calculations of the nME, nRMSE, and rpred2 validation metrics averaged across the five permutations. As displayed in Figure 6, the normalized relationships obtained from the 17 training datasets predicted, in average across the five permutations, the experimental data from the test datasets with average nME between 52 and 300%, nRMSE between 13 and 24%, and coefficients of determination rpred2 between 0.20 and 0.74 and of the same order of magnitude as the experimental rexp2 values.

Fivefold cross-validation of the normalized mathematical relationships.

Here are reported for each dataset the average values across the five permutations of (A) the normalized maximum error (nME), (B) the normalized root mean square error (nRMSE), and (C) coefficient of determination (rpred2, grey bars), which is compared with the coefficient of determination (rexp2, black bars) of the power trendline fitted to the log–log transformation of global experimental datasets directly.

Scaling of the normalized relationships

The normalized relationships provided in Appendix 1—table 2 were finally scaled to typical cat data obtained from the literature following the procedure described in ‘Methods’. Dsoma and Sneuron were found to vary in cats over an average qSE=2.4-fold range according to a review of the morphometric studies reported in the two first lines of Appendix 1—table 3. qSE may however be underestimated for reasons discussed in Appendix 1. Then, the empirical qAE and theoretical qAT ratios, defined in ‘Methods’, were calculated for each MN electrophysiological property A and are reported in Appendix 1—table 4. For example, MN resistance R was found to vary over an average qRE=10.9-fold range in an MN pool according to the literature, while the theoretical fold range qRT=qSEcR=2.42.4=8.4 was obtained from the R=9.6105SMN-2.4 normalized relationship previously derived when a 2.4-fold range is set for SMN. As shown in Appendix 1—table 4, qAEqAT0.8;1.3 for the MN properties ACV, AHP, C, and τ, while the theoretical ranges for R and Ith span over a narrower domain than expected from experimental results (qRT=8.410.9=qRE and qIthT=9.116.3=qIthE). This suggests that the range of variation of R and Ith is not entirely explained by the variation in MN size SMN, and that MN excitability is not only determined by SMN, as reviewed in Powers and Binder, 2001.

The normalized relationships were finally scaled using the theoretical ranges reported in the last column of Appendix 1—table 4. Taking the Ith;R pair as example, as Ith=1.4103R-1 (normalized relationship derived previously and provided in Appendix 1—table 2), R0.5;4.0106Ω and Ith3.9;35.010-9A (Appendix 1—table 4), it was directly obtained from ‘Methods’ that Ith=2.710-2R in SI base units. A similar approach yielded the final mathematical relationships reported in Table 4 between all MN electrophysiological and morphometric properties ACV, AHP, R, Ith,C,τ, and SMN. All constants and relationships are given in SI base units.

Table 4
Mathematical empirical cat relationships between the motoneuron (MN) properties Sneuron, Dsoma, R, Rm, C, τ, Ith,AHP, and ACV.

Each column provides the relationships between one and the eight other MN properties. If one property is known, the complete MN profile can be reconstructed by using the pertinent line in this table. All constants and properties are provided in SI base units (metres, seconds, ohms, farads, and amperes). The relationships involving Rm were obtained from theoretical relationships involved in Rall’s cable theory (see ‘Discussion’).

Sneuron[m2]Dsoma[m]R[Ω]Rm[Ωm2]C[F]τ[s]Ith[A]AHP[s]ACV[ms1]
Sneuron[m2]Dsoma=1.8102SneuronR=1.71010Sneuron2.43Rm=1.01010Sneuron1.43C=1.3102Sneuronτ=1.01012Sneuron1.48Ith=3.8108Sneuron2.52AHP=1.01011Sneuron1.51ACV=3.0106Sneuron0.69
Dsoma[m]Sneuron=5.5103DsomaR=5.1105Dsoma2.43Rm=1.7107Dsoma1.43C=7.9105Dsomaτ=2.3109Dsoma1.48Ith=7.8102Dsoma2.52AHP=2.7108Dsoma1.51ACV=8.1104Dsoma0.69
R[Ω]Sneuron=9.5105R0.41Dsoma=1.7102R0.41Rm=5.7105R0.59C=1.5106R0.40τ=9.6107R0.61Ith=2.7102R1.04AHP=1.3105R0.62ACV=4.9103R0.29
Rm[Ωm2]Sneuron=1.5107Rm0.70Dsoma=2.6105Rm0.70R=6.8106Rm1.70C=1.8109Rm0.69τ=1.4102Rm1.04Ith=2.2109Rm1.76AHP=2.2101Rm1.06ACV=5.4101Rm0.48
C[F]Sneuron=1.2102CDsoma=2.1104CR=1.51015C2.48Rm=1.71013C1.46τ=8.61016C1.51Ith=6.51013C2.57AHP=7.71015C1.54ACV=8.2107C0.71
τ[s]Sneuron=8.4109τ0.67Dsoma=1.5106τ0.67R=7.1109τ1.64Rm=3.6101τ0.96C=1.61010τ0.66Ith=1.61012τ1.70AHP=1.7101τ1.02ACV=7.5τ0.47
Ith[A]Sneuron=4.0104Ith0.40Dsoma=7.1102Ith0.40R=3.1102Ith0.96Rm=7.4106Ith0.57C=5.9106Ith0.39τ=1.2107Ith0.59AHP=1.5106Ith0.60ACV=1.3104Ith0.27
AHP[s]Sneuron=5.5108AHP0.66Dsoma=9.8106AHP0.66R=7.5107AHP1.61Rm=2.5AHP0.95C=9.71010AHP0.65τ=6.2102AHP0.98Ith=1.81010AHP1.67ACV=2.7101AHP0.46
ACV[ms1]Sneuron=4.61010ACV1.44Dsoma=8.3108ACV1.44R=8.31012ACV3.50Rm=2.3103ACV2.06C=9.01012ACV1.41τ=7.4101ACV2.14Ith=1.11015ACV3.64AHP=1.4103ACV2.18

Predicting correlations between MN properties

With this approach, the correlations between MN properties that were never concurrently measured in the literature are predicted, as displayed in Figure 2B. Such unknown relationships were indirectly extracted from the combination of known relationships (Table 3) and typical ranges of values obtained from the literature for these properties. Due to the prior fivefold cross-validation of the relationships in Table 4, these findings are reliable as indirectly consistent with the literature data processed in this study.

Extrapolation to rats and mice

It was assessed whether the scaled relationships obtained from cat data and reported in Table 4 could be successfully applied to rat and mouse data. A global Ith;R dataset was obtained from five studies focussing on rats reported in Appendix 1—table 5, while four mice studies reported in Appendix 1—table 5 provided processable data for the pairs Ith;R, τ;R, and C;R. Only data from control groups of wild-type animals were considered while data on mutated, operated, or trained animals were disregarded. Because the information on the age distribution of the retrieved datasets is not available, data from adults rats and mice of unknown age were used, despite well-known age-related variations in the MN electrophysiological properties in these animals (Highlander et al., 2020; Huh et al., 2021). Like the cat data, low variability was observed between the R and Ith distributions reported by different studies in the rat and mouse Ith;R datasets. Figure 7 reports the global datasets for all four pairs of properties and the predictions obtained with the cat relationships reported in Table 4. As displayed in Appendix 1—table 5, rpred2 and rexp2 values were close for all four datasets, with a maximum 12.5% difference for τ;R in mice. nRMSE values remained low in both rat and mouse Ith;R datasets, while being substantially higher in both τ;R and C;R datasets in mice because of an inter-species-specific offset in the relationships.

Global datasets for rat and mouse species and predictions of the motoneuron (MN) properties with the cat mathematical relationships (Table 4).

They were obtained from the five studies reporting data on rats and the four studies presenting data on mice reported in Appendix 1—table 5 that measured the {Ith;R}, {τ;R}, and {C;R} pairs of MN electrophysiological properties. Ith, R, τ, C are given in nA, MΩ, ms, and nF, respectively. The experimental data (○ symbol) is fitted with a power trendline (red dotted curve) and compared to the predicted quantities obtained with the scaled cat relationships in Table 4 (■ symbol).

Relationships between MN and mU properties

As shown in Figure 2A, eight pairs of one MN and one mU property were investigated in 14 studies in the literature in cats and rats. As remarked by Heckman and Enoka, 2012 and Manuel et al., 2019, most recent studies were performed on mice, while most dated were dominantly performed on cats. One study on the rat gastrocnemius muscle (Kanda and Hashizume, 1992) indicated no correlation between IR and ACV. However, after removing from the dataset two outliers that fell outside 2 standard deviations of the mean IR data, a statistically significant correlation (p<0.05) between IR and ACV was successfully fitted with a power trendline (r2=0.43). Eight studies, dominantly focussing on the cat soleus and medial gastrocnemius muscles, found a strong correlation between Ftet and ACV, while one study showed a significant correlation for the pair Ftet;R in both the cat tibialis anterior and extensor digitorum longus muscles. One study (Burke et al., 1982) on the cat soleus, medial and lateral gastrocnemius muscles inferred a statistically significant correlation for Ftet;SMN. In mice, both R and Ith were significantly related to muscle tetanic and twitch forces. As Ftet, Ftw, and IR are reliable indices of SmU as discussed in ‘Methods, the MN size-dependent relationships in Table 3 return eight SmU=kSMNc relationships between mU and MN properties (Table 5). A 2.8-fold range in positive c-values c[2.4;6.6] (mean 3.8±1.5 sd) was obtained between studies. This result infers that mU and MN size are positively related in cats, mice, and potentially in rats.

Table 5
Fitted experimental data of pairs of one muscle unit (mU) and one motoneuron (MN) property and subsequent SmUSMNc relationships.
SpeciesA=kBc(fitted relationships)(final relationships)
Relationshipcr2p-valueReference studiesc
RatIR=kACVc3.40.456.10-3Kanda and Hashizume, 19922.4
CatFtw=kACVc9.40.43<10-5Knott et al., 19716.6
Ftet=kACVc7.20.37<10-5Mcphedran et al., 1965; Wuerker et al., 1965; Appelberg and Emonet-Dénand, 1967; Proske and Waite, 1974; Bagust, 1974; Jami and Petit, 1975; Stephens and Stuart, 1975; Burke et al., 1982; Emonet-Dénand et al., 19885.0
Ftet=kRc-1.30.276.10-5Dum and Kennedy, 19803.2
Ftet=kSMNc2.00.212.10-2Burke et al., 19822.0
MouseFtw=kRc-2.10.42<10-5Manuel and Heckman, 2011; Martínez-Silva et al., 20185.1
Ftw=kIthc1.30.642.10-4Manuel and Heckman, 20113.3
Ftet=kIthc1.00.806.10-2Manuel and Heckman, 20112.5
Mean ± sd3.8 ± 1.5
Table 5—source data 1

Numerical data used to derive the relationships presented in Table 5.

https://cdn.elifesciences.org/articles/76489/elife-76489-table5-data1-v2.xlsx

Discussion

We processed the normalized data (Figure 3) from previous cat experimental studies to extract mathematical relationships (Table 4) between several MN electrophysiological and morphological properties (Table 1), providing a clear summary of previously known qualitative inter-relations between these MN properties. Table 4 is a new convenient tool for neuroscientists, experimenters, and modellers. Besides obtaining significant quantitative correlations between MN electrophysiological properties and MN size (Table 3, Figure 4, Figure 5), we established (Table 4) and validated (Figure 6) a comprehensive mathematical framework that models an extensive set of experimental datasets available in the literature, quantitatively links all the pairs of the MN properties in Table 1, and is consistent with Rall’s cable theory, as discussed in the following. This framework is directly applicable, although with limitations, to other mammalian species (Figure 7), and is used along with other measurements of mU properties (Table 5) to discuss the Henneman’s size principle.

Consistency of the relationships with previous empirical results and Rall’s theory

The mathematical relationships derived in Table 4 are first in strong agreement with the relative variations between MN properties that were reported in review papers and are listed in ‘Introduction’. Some combinations of the relationships in Table 4 are besides consistent with other cat experimental data that were not included in the global datasets and were not used for deriving the relationships. For example, the relationships Ith=kτ-1.7, 1C=kτ0.7, 1SMN =kACV-1.4, and 1R=kACV3.5 in Table 4 yield, when combined, IthC=kτ-1 and 1RSMN=kACV2.1, which are consistent with the relationship IthC=kτ-1 experimentally observed by Gustafsson and Pinter, 1985 and 1RSMN=kACV1.9 measured by Kernell and Zwaagstra, 1981. The stronger-than-linear inverse relationship R=kSMN-2.4 is also consistent with the phenomenological conclusions by Kernell and Zwaagstra, 1981. Moreover, some standard equations resulting from Rall’s cable theory can be derived with a combination of the relationships in Table 4, assuming the MNs as equivalent cylinders with spatial uniformity of both Cm and Rm. For example, taking tanhLL=0.5, which is in the typical 0.48;0.76 range reported in the literature (Powers and Binder, 2001), the relation C=τRLtanhL (eq. 1. in Gustafsson and Pinter, 1984a) is obtained when combining the R-Sneuron, C-Sneuron, and τ-Sneuron relationships. Also, taking tanhLL=0.6 (standard value provided in Powers and Binder, 2001), the R-Sneuron equation in Table 4 applied to R=RmSneuronLtanhL in Rall, 1977, yields Rm=1.010-10Sneuron1.4. This Rm-Sneuron relationship is consistent with the observations in Powers and Binder, 2001 that smaller MNs have larger Rm values. Also, considering that qSE=2.4 and Sneuron0.18;0.44 mm2 (Appendix 1—table 3), Rm=1.010-10SMN1.4 makes Rm vary over the 3.5-fold theoretical range 0.12;0.44 Ωm2 in the MN pool. This is highly consistent both with the qAE=3.2-fold range and the AminE; AmaxE=0.16;0.59 Ωm2 mean range previously reported in the literature (Albuquerque and Thesleff, 1968; Barrett and Crill, 1974; Burke et al., 1982; Gustafsson and Pinter, 1984b). Also, Rm=1.010-10Sneuron1.4 yields Rm=2.5AHP0.9 from the relationships in Table 4, which is consistent with the indirect conclusion of a positive correlation between Rm and AHP in Gustafsson and Pinter, 1984b. Finally, R=0.6RmSneuron and Rm=1.010-10Sneuron1.4 support the statement that the variations in Rm in the MN pool may be as important as cell size Sneuron in determining MN excitability (Gustafsson and Pinter, 1985), although the variation of R in an MN pool may be also related to cell geometry (Gustafsson and Pinter, 1984b), which was not investigated in this study. This is consistent with the results in Appendix 1—table 4, which demonstrate that the range of variation qSE of SMN does not entirely explain the ranges of variation qRE and qIthE of R and Ith, and that MN excitability is not only determined by MN size. Because of the aforementioned consistency with previous findings and with Rall’s theory, the remaining relationships between Rm and the other MN properties were calculated from Rm=1.010-10Sneuron1.4, following the same methods as before, and added to Table 4. Importantly, Rm might however not be uniform across the somatodendritic membrane, according to results from completely reconstructed MNs (Fleshman et al., 1988), and might be positively correlated to the somatofugal distance (Fleshman et al., 1983), with dendritic Rm being higher than somatic Rm (Barrett and Crill, 1974; Powers and Binder, 2001; Olinder et al., 2006) by a factor 100–300 (Clements and Redman, 1989). This contradicts the assumptions of membrane isopotentiality necessary to apply Rall’s cable theory. The relationships in Table 4, which were obtained from measurements in the soma where Rm is mainly constant, may thus not directly extrapolate to the dendritic regions of the MNs. Yet, the dimensionless variations of the average across MN surface of Rm with the other MN properties may still be valid (Olinder et al., 2006).

The equations in Table 4 support the notion that MNs behave like resistive-capacitive systems, as demonstrated by the combination of the R-τ and C-τ relationships in Table 4, which yields τ0.9=0.8RC, close to the theoretical equation τ=RCtanhLL derived from Rall’s theory. Also, Table 4 directly yields, from the combination of Rm=1.010-10Sneuron1.4 with the other relationships, τ=1.410-2Rm, which is exactly eq. 2.15 in Rall, 1977 (τ=RmCm) with a constant value for Cm=1.410-2Fm-2 that is consistent with the literature.

The relationship C=1.310-2Sneuron obtained in Figure 4 and reported in Table 4 is in agreement with the widely accepted proportional relationship between MN membrane capacitance and surface area, when spatial uniform values of Rm and Cm are assumed. It is worth noting that this proportional relationship was obtained in this study without simultaneous empirical measurements of SMN and C. The values of C were obtained from Rall’s C=τRLtanhL equation (Gustafsson and Pinter, 1984b) and the corresponding values of SMN were indirectly obtained from the final datasets ACV;SMN, AHP;SMN, R;SMN and Ith;SMN that involved the data from 17 selected studies as described in Figure 1. Although consistent with typical values reported in the literature (Lux and Pollen, 1966; Albuquerque and Thesleff, 1968; Barrett and Crill, 1974; Adrian and Hodgkin, 1975; Sukhorukov et al., 1993; Major et al., 1994; Solsona et al., 1998; Thurbon et al., 1998; Gentet et al., 2000), the Cm=1.310-2Fm-2 constant value reported in this study is slightly higher than the widely accepted Cm=1.010-2Fm-2 value (Ulfhake and Kellerth, 1984; Gustafsson and Pinter, 1984b). This high value is consistent with the conclusions from Gustafsson and Pinter, 1984a, who obtained slightly conservative values for Sneuron when approximating the MN surface area as Sneuron=CCm with Cm=1.010-2Fm-2.

Finally, the empirical relationship Ith=2.710-2R1.0 (Table 4) provides Vth=27mV, consistently with typical reported values (Brock et al., 1952; Eccles et al., 1958a), despite uncertainties in the value of the membrane resting potential (Heckman and Enoka, 2012) and lower values reported by Gustafsson and Pinter, 1984a. This supports the conclusions that the relative voltage threshold Vth may be assumed constant within the MN pool in first approximation (Coombs et al., 1955; Gustafsson and Pinter, 1984a; Powers and Binder, 2001), and that Ohm’s law is followed in MNs for weak subthreshold excitations (Glenn and Dement, 1981; Ulfhake and Kellerth, 1984; Spruston and Johnston, 1992; Powers and Binder, 2001; Olinder et al., 2006), as discussed in ‘Methods’. However, voltage thresholds Vth tend to be lower in in MNs of large R and long AHP (Gustafsson and Pinter, 1984a), suggesting that Vth might be inversely correlated to MN size SMN. Moreover, because of the voltage activation of persistent inward currents (PICs) near threshold which add to the external stimulating current and increases the MN input resistance (Fleshman et al., 1988), with a prominent effect on small MNs (Powers and Binder, 2001), the experimental values of Vth are generally greater than that directly predicted from the product IthR (Gustafsson and Pinter, 1984b), and a variant of Ohm’s law such as VthV=IthRV should be considered during MN discharging events. This voltage-dependent variation in the values of R might partly explain why the range of the values of R reported in the literature for weak subthreshold currents (Appendix 1—table 4) is generally lower than that of Ith. While Vth may be always size-dependent, and size-voltage-dependent close to threshold and during firing events, no correlation between VthRIth and C, AHP, or ACV was found and reported in Table 4 in subthreshold conditions, consistent with measurements performed by Gustafsson and Pinter, 1984b and Ulfhake and Kellerth, 1984, substantiating that the dynamics of MN recruitment dominantly rely on R, Ith, and Vth (Heckman and Enoka, 2012).

Relevance of the relationships

Table 4 provides the first quantitative description of an extensive set of cat experimental data available in the literature. These inter-related relationships were validated, successfully reconciliate the conclusions formulated by the selected experimental studies, and are in strong agreement with the theoretical equations describing Rall’s cable theory. This robust framework advances our general understanding of the MN neurophysiology by inferring quantitative relationships for pairs of MN properties that were never concurrently measured in experiments (Figure 2B) by crossmatching other relationships involving these properties. For example, Ith=3.8108Sneuron2.5 (Table 4) quantitatively supports the size dependency of Ith that was predicted by Powers and Binder, 2001 from Rall’s equations. Furthermore, modellers can directly and realistically tune with Table 4 the physiological parameters of phenomenological MN models, such as leaky fire-and-integrate models, and/or build pools of MNs displaying a continuous distribution of realistic MN profiles of MN-specific electrophysiological and morphometric properties, as it has been previously attempted (Dong et al., 2011; Negro et al., 2016). Models tuned with Table 4 can, for example, better replicate the discharging behaviour of MNs, obtained from decomposed high-density EMGs, than models scaled with generic property values (Caillet et al., 2022). In this view, such tuned models display MN-specific behaviours that should be easier to interpret. The mathematical relationships in Table 4 can also support future experimental investigations performed on spinal MNs in adult mammals in vivo by completing experimental datasets for which properties that are difficult to measure were not obtained, such as the MN membrane time constant. Experimenters can therefore choose to reduce the workload of measuring all the MN properties in Table 1 to favour the identification of a maximum number of MNs and obtain extensive datasets describing larger MN populations that provide an informative window on the continuous distribution of MN properties in an MN pool. The unknown electrophysiological MN properties, typically obtained in vivo, of cadaveric specimens can also be estimated with the size-related relationships in Table 4 from in vitro measurements of the somal diameter Dsoma in a slice preparation of the spinal cord. The mathematical relationships in Table 4 describe the experimental data published in the literature and provide a convenient metric for experimenters to compare their measurements to previous findings.

Extrapolation of the relationships to other mammals

It is demonstrated in Figure 7 and with the calculation of the r2 values in Appendix 1—table 5 that the mathematical equations in Table 4 obtained from cat data adequately predict the normalized associations between pairs of MN properties in rats and mice. However, the scaling factor K applied to the intercept ke to scale the normalized relationships A=KkeBe is species-specific, being, for example, around three times smaller in mice than in cat for the τ;R and C;R pairs (Figure 7) and explaining the large nRMSE values in these cases. Absolute values of the Ith;R pair were nevertheless accurately predicted in both rats and mice from the cat relationships, as explained by the larger values of R which counterbalances the respectively lower Ith values in mice and rats than in cats (Manuel et al., 2019; Highlander et al., 2020). Despite the age-related data variability (Highlander et al., 2020) in the rodent dataset displayed in Figure 7, these findings advance our understanding of the systematic inter-mammalian-species variations in MN properties. While the mathematical equations in Table 4 are specific to cats, the normalized equations in Appendix 1—table 2 can be scaled with species-specific data if investigating other mammals.

Henneman’s size principle of motor unit recruitment

Table 5 reports statistically significant power relationships of positive power values between MN and mU indices of size. These results substantiate the concept that SMN and SmU are positively correlated in a motor unit (MU) pool and that large MNs innervate large mUs (Henneman, 1981; Heckman and Enoka, 2012), a statement that has never been demonstrated from the concurrent direct measurement of SMN and SmU. Besides, considering the positive Ith-SMN correlation (Table 4), and that the mU force recruitment threshold Fth is positively correlated to Ftet (Heckman and Enoka, 2012) and thus to SmU (Table 2), larger MUs have both larger current and force recruitment thresholds Ith and Fth than relatively smaller MUs, which are thus recruited first, consistently with the Henneman’s size principle of MU recruitment (Henneman, 1957; Wuerker et al., 1965; Henneman et al., 1965a; Henneman et al., 1965b; Henneman et al., 1974; Henneman, 1981; Henneman, 1985). The terminologies ‘small MU’, ‘low-force MU’, and ‘low-threshold MU’ are thus equivalent. Henneman’s size principle thus relies on the amplitude of the MN membrane resistance (Binder et al., 1996; Powers and Binder, 2001; Heckman and Enoka, 2012). Finally, the relationships SMNCV1.4τ-0.7AHP-0.7 (Table 4) suggest that high-threshold MUs rely on relatively faster MN dynamics, which might partially explain why large MNs can attain relatively larger firing rates than low-thresholds MNs, for example, during ballistic contractions or during events close to maximum voluntary contractions (Powers and Binder, 2001; Heckman and Enoka, 2012).

It has been repeatedly attempted to extend Henneman’s size principle and the correlations between the MU properties in Table 1 to the concept of ‘MU type’ (Burke and ten Bruggencate, 1971; Burke, 1981; Bakels and Kernell, 1993; Powers and Binder, 2001). While a significant association between ‘MU type’ and indices of MU size has been observed in some animal (Fleshman et al., 1981; Burke et al., 1982; Zengel et al., 1985) and a few human (Milner-Brown et al., 1973; Stephens and Usherwood, 1977; Garnett et al., 1979; Andreassen and Arendt-Nielsen, 1987) studies, it has however not been observed in other animal studies (Bigland-Ritchie et al., 1998 for a review) and in the majority of human investigations (Sica and McComas, 1971; Goldberg and Derfler, 1977; Yemm, 1977; Young and Mayer, 1982; Thomas et al., 1990; Nordstrom and Miles, 1990; Elek et al., 1992; Macefield et al., 1996; Van Cutsem et al., 1997; Mateika et al., 1998; Fuglevand et al., 1999; Keen and Fuglevand, 2004). Moreover, the reliability of these results is weakened by the strong limitations of the typical MU-type identification protocols. Sag experiments are irrelevant in humans (Buchthal and Schmalbruch, 1970; Thomas et al., 1991; Bakels and Kernell, 1993; Macefield et al., 1996; Bigland-Ritchie et al., 1998; Fuglevand et al., 1999) and lack consistency with other identification methods (Nordstrom and Miles, 1990). MU-type identification by twitch contraction time measurements is limited by the strong sources of inaccuracy involved in the transcutaneous stimulation, intramuscular microstimulation, intraneural stimulation, and spike-triggered averaging techniques (Taylor et al., 2002; Keen and Fuglevand, 2004; McNulty and Macefield, 2005; Negro et al., 2014; Dideriksen and Negro, 2018). Finally, as muscle fibres show a continuous distribution of contractile properties among the MU pool, some MUs fail to be categorized in discrete MU types in some animal studies by histochemical approaches (Reinking et al., 1975; Tötösy de Zepetnek et al., 1992). Owing to these conflicting results and technical limitations, MU type may not be related to MN size and the basis for MU recruitment during voluntary contractions (McNulty and Macefield, 2005; Duchateau and Enoka, 2011).

Limitations

The mathematical relationships derived Table 4 and the conclusions drawn in ‘Discussion’ present some limitations, most of which were previously discussed and are summarised here. A first limitation relates to the experimental inaccuracies in measuring the MN properties in Table 1. As discussed, the true values of R and τ, and thus of Rm and C, may be underestimated because of a membrane leak conductance around the electrode and some voltage-activated membrane nonlinearities. Yet, the selected studies reduced the impact of these phenomena, which besides may not affect the association between MN properties (Gustafsson and Pinter, 1984a). Also, τ, C, and Rm were estimated assuming that the MN membrane is isopotential. This is valid in the soma, where τ, C, and Rm were measured, and the relationships in Table 4 are adequate for an MN simplified to one equivalent cable. However, as the membrane resistivity may increase with somatofugal distance, as discussed previously, the relationships in Table 4 may be offset towards larger R values in dendritic areas, which should be accounted for when modelling the MN dendrites in compartmental MN models or completely reconstructed dendritic trees. Finally, the selected cat studies are relatively dated, as observed in Manuel et al., 2019. Thus, recent computer-assisted techniques of MN morphometric measurements were notably not applicable in those studies, yielding important sources of inaccuracies, as discussed in Appendix 1. Therefore, the relationships involving Sneuron are subjected to a higher level of inaccuracy. A second limitation arises from the inter-study variability of the experimental datasets reported in the selected studies. As discussed, animal preparations slightly diverged between the selected studies. Also, because of the scarcity of available data in the literature, measurements obtained from MNs innervating different muscles from adult animals of different sex and age were merged into unique datasets. This inter-study variability, added to the ‘unexplained’ variance reported in Highlander et al., 2020, was however reduced by the normalization of the experimental datasets, as investigated in Appendix 1—figures 2 and 4. The normalization approach is however only valid if the experimental studies identified the same largest MN relatively to the MN populations under investigation. This cannot be verified but is likely as most studies returned similar normalized distributions of their measured parameters, as displayed in Appendix 1—figure 1, inferring that a similar portion of the MN pools was identified in the selected studies. Finally, the results proposed in this study may be affected by methodological bias from the three research groups that provided most of the data points processed in this study, even if all the selected studies reported similar methods to measure the properties, as discussed in ‘Methods’.

A third limitation is due to the limited processable experimental data available in the literature, especially in rats and mice, as also discussed by Heckman and Enoka, 2012. The extrapolation of the relationships in Table 4 to rat and mouse species could only be assessed against a small set of studies and for three associations only. The conclusions on the associations between SmU and SMN also relied on few studies providing few measurements due to the considerable amount of work required to measure both MN and mU properties in single MUs. In cats, some pairs of MN properties were investigated in only one study, preventing inter-study comparisons, and stressing the usefulness of Table 4 for reconciliating the data available in the literature. Within the experimental datasets, the scarcity of data points in some regions of the data distributions (e.g., see the skewed Ith and R distributions in the density histograms in Appendix 1—figure 1) transfers during processing to the final size-dependent datasets (Figure 4). This data heterogeneity may affect the reliability of the relationships in Table 4 for the extreme MN sizes (see, in Figure 4, the 95% confidence interval of the regressions widening for small SMN in the AHP;SMN and R;SMN datasets and for large SMN in the Ith;SMN dataset). In the validation procedure, the highest nME values reported in Figure 6 were obtained for these regions of limited data for most of the global datasets. This heterogeneous density of the data may be explained by the natural skewness of the MN property distributions in the MN pool towards many small MNs and few large MNs (Heckman and Enoka, 2012) and/or a systematic experimental bias towards identifying and investigating relatively large MUs.

A fourth limitation arises from the relatively low r2 values obtained from the power regressions in Table 3. Although these power trendlines are statistically significant (see p-values in Table 3) and globally provide a better description of the data than other fits (see Appendix 1—table 1), they cannot entirely describe the associations between MN property distributions, in line with the conclusions reported in the selected studies. According to these first four limitations, the mathematical relationships in Table 4 best reproduce published cat data but include the level of inaccuracies of the original experimental approaches.

A fifth limitation is the restriction of the study to the MN properties reported in Table 1. Other MN-specific properties, such as the electronic length, the relationship between discharge rate and excitatory current and its hysteresis, or the amplitude of the inhibitory and excitatory postsynaptic potentials, were omitted because they were rarely measured or reported along with SMN, ACV, AHP, R, Ith, C, and τ in the selected studies. Moreover, as the cell geometry was not investigated in this study, other MN properties (Clements and Redman, 1989) such as the effects of the distributed synaptic integration along the dendritic tree were overlooked. Finally, apart from Ith, ACV, and AHP, all the MN properties investigated in this study are passive properties, that is, baseline properties of the cell at rest (Powers and Binder, 2001), as they were measured using weak sub-threshold current pulses. Therefore, other MN active properties, such as voltage-dependent ionic channel-related properties including PICs, which activate close or above threshold, were not considered in this study. While this limits the applicability of Table 4 to a functional context or in comprehensive Hodgkin–Huxley-like models of MNs, some relationships in Table 4 remain pertinent in those conditions if they are linked with other observations in the literature on MN active properties, for example, MNs of small SMN having longer lasting total PICs of greater hyperpolarized activation than MNs of large SMN (Heckman and Enoka, 2012). Despite the aforementioned limitations, the relationships in Table 4 can be directly used in phenomenological RC approaches like LIF models (Izhikevich, 2004; Teeter et al., 2018), which rely on the properties reported in Table 1, to derive profiles of inter-consistent MN-specific properties and describe realistic continuous distributions of the MN properties R, Ith, C, Vth, and τ in the MN pool, as attempted in previous works (Dong et al., 2011; Negro et al., 2016). An example of a computational modelling application of Table 4 is provided in Caillet et al., 2022, where the relationships tune a cohort of LIF models which estimate, from the firing activity of a portion of the MN pool obtained from decomposed high-density EMGs, a realistic distribution of the MN properties in the MN pool and the firing behaviour of a complete MN pool in an isometrically contracting human muscle.

As a final limitation, the relationships in Table 4 were obtained from a regression analysis and therefore provide correlations between some MN properties in an MN pool but cannot be used to draw conclusions on the causality behind these associations.

Conclusion

This study provides in Table 4 a mathematical framework of quantitative associations between the MN properties SMN,ACV, AHP, R, Rm, Ith,C, and τ. This framework, which is consistent with most of the empirical and theoretical conclusions in the literature, clarifies our understanding of the association between these MN properties and constitutes a convenient tool for neuroscientists, experimenters, and modellers to generate hypotheses for experimental studies aiming at investigating currently unreported relationships, support experimentations, and build virtual MN profiles of inter-consistent MN-specific properties for MN modelling purposes.

Appendix 1

Inter-study data variability

Inter-study variability greater than the threshold defined in ‘Methods’ was observed in three, two, and four property distributions for the metrics range, mean, and CoV metrics, respectively. Indeed, sensibly different ranges of property values were observed to be reported between studies of the ACV;SMN, AHP;ACV, and R;AHP datasets for the SMN, ACV, and AHP properties, respectively, with sdg10;15% and/or sdgmeang0.15;0.33. Also, the studies of the R;AHP and τ;R datasets returned some variability in the mean of the reported data for properties AHP and R, respectively, with sdg10;11% and/or sdgmeang0.15;0.23. The ACV;SMN, R;ACV, and R;AHP global datasets showed variable data dispersion around the mean for properties SMN ,R and ACV, and AHP, respectively, with sdg10;11% and/or sdgmeang0.15;0.30.

Discussion on the potential underestimation of qSE

The minimum and maximum values of Sneuron and Dsoma reported in Appendix 1—table 1 may underestimate the true distribution of MN sizes in an MN pool for several reasons, and thus underestimate the true value of qSE. First, the studies reported in Appendix 1—table 1 investigated MN populations of small size, so that the largest and smallest MNs of the investigated MN pools may not have been systematically identified and measured. Also, dendritic labelling by staining may only identify dendrites up to the third dendritic branch (Issa et al., 2010; Brandenburg et al., 2018) and full dendritic trees may not be identified, thus underestimating Sneuron for the largest MNs that exhibit the highest dendritic tree complexity (Brandenburg et al., 2020). Finally, the distribution of gamma- and alpha-MN sizes is typically bimodal and the size distributions of the two populations are reported to overlap (Moschovakis et al., 1991; Vult von Steyern et al., 1999; Friese et al., 2009; Deardorff et al., 2013). However, most studies reported in Appendix 1—table 1 identify the alpha MNs from the gamma MNs by visually splitting the bimodal distributions of MN sizes in two independent domains (Donselaar et al., 1986; Ishihara et al., 2001), thus neglecting the overlap. Similarly, some studies identify as gamma MNs any MN showing an axonal conduction velocity typically less than ~60ms-1 ( Zwaagstra and Kernell, 1980; Ulfhake and Kellerth, 1984). These approaches lead to overestimating the lower limit for the reported Dsoma and Sneuron values. Consequently, true ranges in MN sizes may be larger than the theoretical ones reported in Appendix 1—table 1. A similar comment can be made for the electrophysiological properties reported in Appendix 1—table 2. However, the qSE=2.4-fold range obtained for cats is consistent with the fold ranges obtained from rat and mouse data for Dsoma, as reported in Appendix 1—table 1, while it underestimates the 4.4-fold range in Sneuron values observed in mice. Appendix 1—table 1 finally shows that typical cat MN sizes are larger than in rats and mice, as reviewed by Manuel et al., 2019.

Appendix 1—figure 1
Density histograms of the data distributions reported in the experimental studies included in the global datasets.

Depending on the typical range over which each property spans, the distributions are divided in steps of 10 or 20%. The frequency distribution is provided in percentage of the total number of reported data points in a study. Different studies are displayed with different colours in each graph. For each dataset {A;B}, the frequency distributions of property A and B are provided in the left and right plot, respectively.

Appendix 1—figure 2
Assessment of data variability between the experimental studies that constitute the eight global datasets that include at least two experimental studies.

For each experimental study included in the global dataset {A;B}, the range, mean, coefficient of variation (CoV=standard deviationMean), and the ratio MeMd=MeanMedian of the experimental A values measured in this study were computed. Then, the average (bars) and standard deviation (error bars) across experimental studies of these four metrics (range, mean, CoV, MeMd) were calculated for each global dataset independently. For example, the first bar in the ‘ SMN parameter’ plot represents the average across three studies (Kernell and Zwaagstra, 1981; Cullheim, 1978; Burke et al., 1982) of the range of SMN values reported in these studies, for the global dataset {ACV;SMN}. The computed standard deviations (error bars) express for each global dataset the inter-study variability of (1) the length of the identified bandwidth of the motoneuron (MN) pool (range metric), (2) the spread of values around the mean (CoV metric), (3) the skewness of the distributions, and (4) whether the distributions from different studies are centred (mean metric). A global dataset {A;B} reporting narrow error bars for parameter A for the four metrics infers that the experimental studies constituting this dataset measured similar distributions of property A.

Appendix 1—figure 3
Distribution of the size of the experimental datasets constituting the global datasets for assessment of the variability in the input data.

The histogram is divided between global datasets (half vertical lines), grouped as final size-dependent datasets (full vertical lines). For each global dataset, the total number of data points is reported (label: ‘N’), while the size of the constitutive experimental studies is reported in decreasing order.

Appendix 1—figure 4
Assessment of data variability between the global datasets.

This figure compares the distributions of the SMN, ACV, AHP, R, Ith, and τ properties between the normalized global datasets built in this study. For each property, the range, mean, coefficient of variation (oV=standard deviationMean), and the ratio MeMd=MeanMedian of its distribution in a global dataset is compared to the other global datasets it appears in. For each property, the standard deviation sdG between global datasets of these four metrics is computed. A low sdG value reflects the low variability of the property distribution between global datasets.

Appendix 1—table 1
r2 values obtained for each experimental dataset {A;B} when performing a linear regression analysis on {A;B} directly (‘Linear’) and on the ln(A)ln(B) (‘Power’) and ln(A)B (‘Exponential’) transformations of {A;B}.

The r2 values returned by the three types of regression cannot be directly compared to estimate the best model. However, the power fit returned r2>0.5 for relatively more experimental datasets than the linear and exponential fits.

StudiesLinearPowerExponential
{ACV;SMN}Cullheim, 19780.420.450.43
Kernell and Zwaagstra, 19810.620.630.6
Burke et al., 19820.230.260.24
{AHP;SMN}Zwaagstra and Kernell, 19800.330.370.35
{AHP;ACV}Eccles et al., 1958b0.530.540.54
Zwaagstra and Kernell, 19800.420.440.42
Gustafsson and Pinter, 1984a0.710.680.7
Foehring et al., 19870.150.140.14
{R;SMN}Kernell and Zwaagstra, 19810.680.710.7
Burke et al., 19820.410.50.46
{R;ACV}Kernell, 19660.700.650.64
Burke, 19680.280.310.31
Kernell and Zwaagstra, 19810.680.710.71
Fleshman et al., 19810.30.210.2
Gustafsson and Pinter, 1984a0.450.440.42
Sasaki, 19910.680.680.66
{R;AHP}Gustafsson, 19790.430.410.42
Gustafsson and Pinter, 1984a0.710.680.67
Foehring et al., 19870.330.320.33
Pinter and Vanden Noven, 19890.660.530.59
Sasaki, 19910.570.490.52
{Ith;R}Kernell, 19660.520.520.68
Fleshman et al., 19810.280.340.38
Gustafsson and Pinter, 1984b0.630.830.83
Zengel et al., 19850.490.60.62
Munson et al., 19860.340.420.47
Foehring et al., 19870.340.470.47
Krawitz et al., 20010.320.490.44
{Ith;ACV}Kernell and Zwaagstra, 19810.230.240.24
Gustafsson and Pinter, 1984b0.40.520.5
{Ith;AHP}Gustafsson and Pinter, 1984b0.550.720.69
{C;R}Gustafsson and Pinter, 1984a0.550.570.57
{C;Ith}Gustafsson and Pinter, 1984a0.230.340.24
{C;AHP}Gustafsson and Pinter, 1984a0.250.260.26
{C;ACV}Gustafsson and Pinter, 1984a0.170.220.19
{τ;R}Burke and ten Bruggencate, 19710.040.070.03
Gustafsson, 19790.60.630.61
Gustafsson and Pinter, 1984a0.630.690.59
Zengel et al., 19850.390.330.36
Pinter and Vanden Noven, 19890.650.690.62
Sasaki, 19910.680.720.64
{τ;Ith}Gustafsson and Pinter, 1984b0.630.590.57
{τ;ACV}Gustafsson and Pinter, 1984a0.690.760.76
{τ;AHP}Gustafsson and Pinter, 1984a0.340.320.31
Appendix 1—table 2
Mathematical empirical normalized relationships between the motoneuron (MN) properties Dsoma, R, C, τ,Ith,AHP, and ACV.

Each column provides the relationships between one and the six other MN properties. All constants and properties are normalized up to a theoretical 100% maximum value. To scale the normalized relationships A=kBc for a specific mammalian species, the normalized intercept k must be scaled with a pair of values for properties A and B obtained from the same motoneuron in that species.

DsomaRCτIthAHPACV
DsomaR=9.6105Dsoma2.4C=1.2Dsomaτ=2.6104Dsoma1.5Ith=9.0104Dsoma2.5AHP=2.5104Dsoma1.5ACV=4.0Dsoma0.7
RDsoma=2.9102R0.4C=2.7102R0.4τ=5.8R0.6Ith=1.5103RAHP=4.7R0.6ACV=2.0102R0.3
CDsoma=8.6101CR=1.4106C2.5τ=3.3104C1.5Ith=6.2104C2.6AHP=3.1104C1.6ACV=3.6C0.7
τDsoma=9.5102τ0.7R=5.6102τ1.6C=8.5102τ0.7Ith=2.9104τ1.7AHP=7.8101τACV=4.6102τ0.5
IthDsoma=1.6101Ith0.4R=1.1103IthC=1.7101Ith0.4τ=4.2102Ith0.6AHP=3.7102Ith0.6ACV=2.7101Ith0.3
AHPDsoma=8.1102AHP0.7R=8.4102AHP1.6C=7.3102AHP0.6τ=1.3AHPIth=1.9104AHP1.7ACV=4.1102AHP0.5
ACVDsoma=1.4101ACV1.4R=1.2108ACV3.5C=1.7101ACV1.4τ=5.0105ACV2.1Ith=5.9106ACV3.6AHP=5.0105ACV2.2
Appendix 1—table 3
Typical ranges of physiological values for Dsoma and Sneuron in cat rat and mouse species.

SMN is found to vary over an average qSE=2.4-fold range, which sets the amplitude of the theoretical ranges. Absolute {min; max} reports the minimum and maximum values retrieved in the reference studies for Dsoma and Sneuron, while average {min; max} is obtained as the average across reference studies of minimum and maximum values retrieved per study.

Appendix 1—table 4
Typical ranges of physiological values for the motoneuron (MN) properties R, Rm, C, τ, Ith, AHP, and ACV.

As described in ‘Methods’, qAE is the average among reference studies of the ratios of minimum and maximum values; the properties experimentally vary over a qAE -fold range. This ratio compares with the theoretical ratio qAT=qSEc (with c taken from Table 3), which sets the amplitude of the theoretical ranges. Absolute and average {min; max} are obtained as described in the main sections.

MN propertyUnitAbsolute exp {min;max}Average exp {min; max}qAEReference studiesqATqAEqATTheoretical range
R[MΩ]{0.2;8.1}{0.4;4.0}10.9Kernell, 1966; Burke, 1968; Burke and ten Bruggencate, 1971; Barrett and Crill, 1974; Gustafsson, 1979; Kernell and Zwaagstra, 1981; Fleshman et al., 1981; Burke et al., 1982; Ulfhake and Kellerth, 1984; Gustafsson and Pinter, 1984a; Zengel et al., 1985; Foehring et al., 1986; Munson et al., 1986; Foehring et al., 1987; Sasaki, 1991; Krawitz et al., 20018.41.3[0.5;4.0]
C[nF]{2.2;8.5}{2.2;8.5}3.9Gustafsson and Pinter, 1984a; Gustafsson and Pinter, 19852.31.6[3.2;7.5]
τ[ms]{2.0;14.2}{2.9;10.2}3.5Burke and ten Bruggencate, 1971; Barrett and Crill, 1974; Gustafsson, 1979; Ulfhake and Kellerth, 1984; Gustafsson and Pinter, 1984a; Gustafsson and Pinter, 1985; Sasaki, 19913.70.9[2.8;10.3]
Ith[nA]{1.7;52.7}{2.3;36.6}16.3Fleshman et al., 1981; Kernell and Monster, 1981; Ulfhake and Kellerth, 1984; Gustafsson and Pinter, 1984b; Zengel et al., 1985; Foehring et al., 1986; Munson et al., 1986; Foehring et al., 1987; Krawitz et al., 20019.11.8[3.9;35.0]
AHP[ms]{10.6;266.6}{44.2;158.7}4.1Eccles et al., 1957; Gustafsson, 1979; Dum and Kennedy, 1980; Zwaagstra and Kernell, 1980; Ulfhake and Kellerth, 1984; Gustafsson and Pinter, 1984a; Zengel et al., 1985; Foehring et al., 1987; Sasaki, 19913.81.1[42.6;160.2]
ACV[m·s-1]{51.1;127.2}{65.3;114.8}1.8Eccles et al., 1957; Mcphedran et al., 1965; Kernell, 1966; Appelberg and Emonet-Dénand, 1967; Burke, 1968; Barrett and Crill, 1974; Proske and Waite, 1974; Bagust, 1974; Stephens and Stuart, 1975; Cullheim, 1978; Dum and Kennedy, 1980; Zwaagstra and Kernell, 1980; Kernell and Zwaagstra, 1981; Fleshman et al., 1981; Glenn and Dement, 1981; Burke et al., 1982; Gustafsson and Pinter, 1984a; Zengel et al., 1985; Foehring et al., 1986; Foehring et al., 19872.30.8[63.5;116.6]
Appendix 1—table 5
Validation of the cat relationships (Table 4) against rat and mouse data.

nME, normalized maximum error; nRMSE, normalized root mean square error; rpred2, coefficient of determination between experimental and predicted quantities; rexp2, coefficient of determination of the power trendline directly fitted to the experimental data.

AnimalDatasetReference studiesnME (%)nRMSE (%)rpred2rexp2
RatIth;RGardiner, 1993; Bakels and Kernell, 1993; Lee and Heckman, 1998; Button et al., 2008; Turkin et al., 2010; Krutki et al., 2015352160.530.54
MouseIth;RDelestrée et al., 2014; Martínez-Silva et al., 2018; Huh et al., 2021385160.460.46
τ;RManuel et al., 200912191680.350.40
C;RManuel et al., 20091915980.380.38

Data availability

Figure 3 - Source Data 1 contains the numerical data used to generate Figure 3. Figure 7 - Source Data 1 and Figure 7 - Source Data 2 contain the numerical data used to generate Figure 7. Table 5 - Source Data 1 contains the numerical data used to compute the mathematical relationships presented in Table 5. Please note that our study used exclusively data from previous investigations, for which public datasets were not available. The data were manually digitised by the authors from published figures and are made available as supplementary materials.

References

    1. Adrian RH
    2. Hodgkin AL
    (1975) Conduction velocity and gating current in the squid giant axon
    Proceedings of the Royal Society of London. Series B, Biological Sciences 189:81–86.
    https://doi.org/10.1098/rspb.1975.0043
  1. Book
    1. Binder MD
    2. Heckman CJ
    3. Powers RK
    (1996)
    The physiological control of motoneuron activity
    In: Rowell LB, editors. Handbook of Physiology: Exercise, Regulation and Integration of Multiple Systems. New York: Oxford Univ Press. pp. 3–53.
  2. Book
    1. Burke RE
    (1981)
    Motor units: anatomy, physiology, and functional organization
    In: Brooks VB, editors. Handbook of Physiology, the Nervous System, Motor Control. American Physiological Society. pp. 345–422.
  3. Book
    1. Chollet F
    (2021)
    Deep Learning with Python
    Simon and Schuster.
  4. Conference
    1. Fleshman JW
    2. Segev I
    3. Cullheim S
    4. Burke RE
    (1983)
    Matching electrophysiological with morphological measurements in cat α-motoneurons
    In Soc. Neurosci. Abstr (Vol. 9, p. 431).
    1. Henneman E
    (1981)
    Recruitment of motoneurons: the size principle
    Motor Unit Types, Recruitment and Plasticity in Health and Disease pp. 26–60.
    1. Major G
    2. Larkman AU
    3. Jonas P
    4. Sakmann B
    5. Jack JJ
    (1994)
    Detailed passive cable models of whole-cell recorded CA3 pyramidal neurons in rat hippocampal slices
    The Journal of Neuroscience 14:4613–4638.
    1. Obregón G
    2. Ermilov LG
    3. Wen-Zhi Z
    4. Sieck GC
    5. Mantilla CB
    (2009)
    Modeling dendritic arborization based on 3D-reconstructions of adult rat phrenic motoneurons
    Revista Ingeniería Biomédica 3:47–54.
  5. Book
    1. Rall W
    (1977)
    Core conductor theory and cable properties of neurons
    In: Poeter R, editors. The Nervous System. American Physiological Society. pp. 39–97.

Decision letter

  1. Hugo Merchant
    Reviewing Editor; National Autonomous University of Mexico, Mexico
  2. Barbara G Shinn-Cunningham
    Senior Editor; Carnegie Mellon University, United States
  3. Randall K Powers
    Reviewer; University of Washington, United States
  4. Leonardo Abdala Elias
    Reviewer; University of Campinas - UNICAMP, Brazil

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Mathematical Relationships between Spinal Motoneuron Properties" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by Barbara Shinn-Cunningham as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Randall K Powers (Reviewer #2); Leonardo Abdala Elias (Reviewer #3).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

1. The format of the paper should follow the structure of eLife. The legend of Figure 1 should describe very explicitly the steps used to create the final size-related datasets, maybe giving an example. The circles and symbols in Figure 2 should be larger.

2. The Introduction should be more inclusive and respectful with previous studies. It is unfair to state that some inferences are "speculative" since some authors provide linear regression curves (Eccles et al., 1958, for instance) while others consider the theoretical analysis by Rall. Also, when mentioning computational models, since the present study is focused on motor neurons, consider citing the extensive literature on models of motor neurons (see p. 3, l. 51 and 52).

3. In both Intro and Discussion, the authors should be careful when declaring the relevance of the current study. Several statements cannot be supported by analysis, for example, the comments that the method "can accelerate future research in the behaviour of individual MNs" and that mathematical equations can be scaled "to estimate the MN property values that cannot be measured in humans in vivo". How can you scale the equations without prior information on human motor neuron electrophysiology and morphology?

4. The methods section could be expanded to explain some of the measured electrophysiological properties in more detail. Most of the direct measurements (ACV, AHP duration, RN, Ith) are fairly straightforward, with the possible exception of time constant (see comments in Public Review). However, the derived measurements need to be explained in the Methods section. Both estimates of specific membrane resistivity (Rm) and whole cell capacitance (C) can be estimated using an equivalent cable model of the motoneuron along with an estimate of electrotonic length, and the formulas for these should be stated (see Ulfhake and Kellerth, Brain Res, 1984 equations 2 and 3 and Gustafsson and Pinter, JPhysiol, 1984a, equation 1). It should be mentioned in Methods that whole cell capacitance can be used as an estimate of cell surface area by assuming a value for specific capacitance of 1 microF/cm2 and that these area estimates are roughly in line with direct measurements although they tend to be a bit high (see Figure 6 of Gustafsson and Pinter 1984a).

5. Data used in the analysis are both from in vivo and in vitro recordings, not exclusively from in vivo experiments as declared in several parts of the manuscript (e.g., Abstract). Morphometric data cannot be recorded in vivo.

6. The authors should revise their method to include some information on data variability. A simple way is to provide the confidence interval of the regression. In this case, all parameters of the mathematical equations will have a range (95% confidence interval). The values provided in the current manuscript are only average values.

7. It is not clear why the authors used a 70-30% scheme for crossvalidation and why only 5 validations were performed. A larger set of validations are probably needed and at least another training-testing proportion should be carried out.

8. Tables 4, 5 and 7 could be moved to the Supplementary material section.

9. In the merged data sets shown in Figure 4 it would be useful to use different symbols for each data set (i.e., different symbols for {AHP;Smn}1 and {AHP;Smn}2).

10. There is no oval box with fc(C) in Figure 1. Please, revise. Also, consider changing Figure 1 with a pseudo-code of the algorithm to improve clarity.

11. Ith was estimated using different approaches in the selected studies. Some studies have used triangular-shaped currents, while others used step currents. How would different methods for estimating similar parameters influence your analysis (inter-study variability)?

12. Some relations are not unexpected. The proportional relation between cell capacitance and cell size is obvious. Again, you should consider previous theoretical studies based on Rall's theory.

13. The authors should provide more clear evidence on how the proposed method should be translated to future studies on synaptic integration. Since the analysis did not consider any active property, I am not confident that the mathematical relations can help a more comprehensive computer simulation study.

14. Please use ACV or CV, not a mixture of both abbreviations.

15. The section on specific resistance (starting on Line 536) should be expanded to consider the fact that specific resistance may be different in different parts of the neuron (see Public review comments). This is important for modelers who want to explicitly represent dendrites in their motoneuron models (either as a separate compartment, an equivalent cable or completely reconstructed trees).

16. The Discussion section on relevance for modelers contains the following statement: "This supports the conclusions that the relative voltage threshold is constant within the MN pool and that Ohms law is followed in MNs". This is not strictly correct. Gustafsson and Pinter (1984b) showed that voltage threshold tended to be lower in low input conductance cells with long AHPs (see their Figure 4C and D). Also, the voltage threshold predicted by Ohms Law from the product of the measured rheobase and input resistance Vth = Ith*R tended to be lower than the measured voltage threshold, suggesting the activation of an inward conductance near threshold. This requires a variant of Ohm's Law in which the effective resistance is voltage-dependent: Vth=Ith*R(v).

17. Finally, although the vast majority of readers will know that correlation does not imply causation, it is worth stating explicitly in the Discussion section on limitations.

18. p. 33, l. 721: Actually, you should consider that in some experimental conditions, large (high-threshold) MNs attain lower firing rates than small (low-threshold) MNs. Maybe you have to explain in what conditions large MNs attain high firing rates.

https://doi.org/10.7554/eLife.76489.sa1

Author response

Essential revisions:

1. The format of the paper should follow the structure of eLife. The legend of Figure 1 should describe very explicitly the steps used to create the final size-related datasets, maybe giving an example. The circles and symbols in Figure 2 should be larger.

The authors modified the manuscript as suggested by the Reviewer.

To answer both this and comment #10, Figure 1 and its caption were changed to display a concrete example of how the two successive {R; S;MN} and {ith; sMN} final datasets were generated. This concrete example made it possible to go into the detail of the explicit steps used to create these datasets (L.308 to L.322). Figure 2 was adapted, and the circles and symbols were made larger.

2. The Introduction should be more inclusive and respectful with previous studies. It is unfair to state that some inferences are "speculative" since some authors provide linear regression curves (Eccles et al., 1958, for instance) while others consider the theoretical analysis by Rall. Also, when mentioning computational models, since the present study is focused on motor neurons, consider citing the extensive literature on models of motor neurons (see p. 3, l. 51 and 52).

We carefully considered this important comment, echoed by Reviewer 3. It was definitively not our intention to be disrespectful to previous work. We did not wish to mean that the correlations between MN properties were speculative, which is of course not true according to the multiple empirical and theoretical reports from previous studies of significant correlations between MN properties. Our message was rather that, although global conclusions on the relative variations of two MN properties (A is positively correlated with B) were made in previous review papers, these reviews relied on scattered experimental data from independent studies that used different mathematical models to describe the association between the same two MN properties (see L.42 to L.46 of the new Introduction). This weakened the possibility of previous papers to reach quantitative relationships between the MN properties. To avoid misinterpretations and the impression of not well acknowledging previous work, the Introduction has been entirely reshaped. The Introduction was also shortened and oriented to a broader audience, as advised by Reviewer 1.

In more details, the 1st and 3rd paragraphs of the first version of the Introduction were re-written (now L.20 to L.38). This new paragraph provides an extensive list of the strong associations between MN properties that were provided in previous review papers. Contrary to the original version, it is also now underlined that these empirical associations are consistent with the theoretical functions from Rall's cable theory, as extensively reviewed in Powers and Binder (2001), and detailed later in the Discussion section.

Then, a new section was added to the Introduction Section (L.39 to L.54). This section highlights the current limitations in the field that make it difficult to derive quantitative relationships between MN properties. Finally, the end of the introduction was condensed to the main outcomes and applications of this study (L.55 to L.74).

As detailed later in this document, the other sections of the manuscript that were not inclusive with previous studies, such as the 1st section of the Discussion section, were also adjusted.

To follow Reviewer 1’s advice of shortening the Introduction, less emphasis was given on computational MN models in the Introduction; therefore, the authors consider that an extensive citation of the literature on models of MNs is less suitable in this new version of the Introduction.

3. In both Intro and Discussion, the authors should be careful when declaring the relevance of the current study. Several statements cannot be supported by analysis, for example, the comments that the method "can accelerate future research in the behaviour of individual MNs" and that mathematical equations can be scaled "to estimate the MN property values that cannot be measured in humans in vivo". How can you scale the equations without prior information on human motor neuron electrophysiology and morphology?

Thank you for pointing this out. In both the Introduction (L.58 to L.68) and the Discussion sections (L.783 to L.810), the relevance of the study was narrowed for conciseness and clarity to the most pertinent applications of the mathematical relationships. We also provided an example of application from our own work. Debatable or less evident applications such as the ones cited above were removed from the manuscript.

4. The methods section could be expanded to explain some of the measured electrophysiological properties in more detail. Most of the direct measurements (ACV, AHP duration, RN, Ith) are fairly straightforward, with the possible exception of time constant (see comments in Public Review). However, the derived measurements need to be explained in the Methods section. Both estimates of specific membrane resistivity (Rm) and whole cell capacitance (C) can be estimated using an equivalent cable model of the motoneuron along with an estimate of electrotonic length, and the formulas for these should be stated (see Ulfhake and Kellerth, Brain Res, 1984 equations 2 and 3 and Gustafsson and Pinter, JPhysiol, 1984a, equation 1). It should be mentioned in Methods that whole cell capacitance can be used as an estimate of cell surface area by assuming a value for specific capacitance of 1 microF/cm2 and that these area estimates are roughly in line with direct measurements although they tend to be a bit high (see Figure 6 of Gustafsson and Pinter 1984a).

The authors thank the Reviewing Editor and Reviewer 2 for underlining this imbalance between the measurements of SMN and the electrophysiological properties in the description of the experimental approaches and their sources of error.

To extend this section, the selected experimental studies were reviewed to add to the manuscript a detailed description (L.148 to L.177) of the experimental approaches these studies took to measure ACV, AHP, R, Ith, τ and c. The main sources of error in measuring these properties, that were reported in the selected studies and in other studies cited in the manuscript, were identified to be a variable membrane leak conductance around the electrode, voltage-induced membrane nonlinearities, and manual curve fitting, and are now detailed from L.178 to L.217. The equation taken from Gustafsson and Pinter (1984a) to C(C=τRLtanh(L))

was added in L.177, in the context of an equivalent cable model.

As the values of Rm and the indirect predictions of sneuron from electrophysiological values (Gustafsson and Pinter (1984a)) were not included in the global datasets, the two other equations the reviewer refers to

Rm=RSneurontanh(L)L and Sneuron= CCm

were added in the Discussion section, where they were discussed as advised by the Reviewer. From L.712 to L.739, it is explained how Rm=RSneurontanh(L)L is consistent with the sMN − ACV − AHP − R − I − C − τ relationships in Table 4 considering other findings reported in the literature, after which the relationships between these properties and rm were calculated and added to Table 4, as done in the original version of the study. From L.757 to L.760, it is now discussed that the Sneuron\ =CCm relationship can be used to approximate the MN area if the value of c is known and an adequate constant value for cm is taken.

To maintain consistency in the revised manuscript and clarify the Results section, the first sub-section of the Results section of the original version of the manuscript was moved to the Methods section from L.125 to L.147. This block, which has been unchanged since V1, covers the protocols taken by the selected studies to prepare the animals for experimentations and to perform morphometric measurements. It now introduces the added block, described above, about the electrophysiological measurements (from L.148 to L.177).

Sections of this answer were reproduced to answer Reviewer 2’s comments, with the intent of being made public with the preprint.

5. Data used in the analysis are both from in vivo and in vitro recordings, not exclusively from in vivo experiments as declared in several parts of the manuscript (e.g., Abstract). Morphometric data cannot be recorded in vivo.

The authors thank the Reviewing Editor for underlining this important source of confusion for the reader. In the first version of the manuscript, in vivo implicitly referred exclusively to the electrophysiological properties, and not to the morphological parameters which are indeed only obtainable in vitro. in vivo was being repeated through the manuscript to underline that the experimental studies reporting in vitro electrophysiological data were not included in this study. Where this could have led to some confusion, i.e. referred to in vivo morphological data, the related sentences were clarified, or the term ‘in vivo’ was removed in the revised version of the manuscript.

6. The authors should revise their method to include some information on data variability. A simple way is to provide the confidence interval of the regression. In this case, all parameters of the mathematical equations will have a range (95% confidence interval). The values provided in the current manuscript are only average values.

This is an important comment, in line with the comments from Reviewer 3 about data variability, as the inter-study data variability might have been a crucial factor in this analysis. Therefore, we considered important to provide a detailed analysis of data variability, that is now described in the Methods section from L.255 to L.279 and applied in the Results section from L.496 to L.533. The results of this analysis, which include the 95% confidence intervals the reviewer mentions, show in the revised manuscript that the data variability was low between the input experimental datasets and between the global datasets. These results are provided in Appendix Figures 1 to 4 of Appendix 1, by adding the ranges of the parameters for the 95% confidence interval in Table 3, and by plotting the 95% confidence intervals in Figure 3. Below we provide a summary of the procedures we followed.

First, we assessed the variability of the normalized property distributions between the experimental studies that constituted the same global dataset, for each of the 8 global datasets that included at least 2 experimental studies, i.e. {ACV; SMN}, {AHP; ACV}, {R; SMN}, {R; ACV}, {R; AHP}, {Ith; R}, {Ith; ACV}. Therefore, the inter-study data variability of 16 property distributions was analysed. ANOVA and Leven’s tests were not retained to perform this analysis of data variability as most of the property distributions were strongly not normal, according to a Jarque-Bera test for data normality which rejected the null hypothesis with p-values as low as 10−9. The Jarque-Bera test was preferred to other common tests (Shapiro-Wilk for instance) as it was identified to perform well for slightly skewed distributions with long tails (Thadewald, T., and Büning, H. (2007) Journal of applied statistics, 34(1), 87-105), which corresponds well to the data distributions in Figure FSM1. It was rather assessed for each global dataset whether the constitutive property distributions showed high inter-dataset variability in (1) the dataset size, (2) the range of the measured properties, (3) the mean of the distributions, (4) the dispersion of the distributions around their mean, and (5) the ratio between mean and median values of the distributions. The methodological approach to assess such variability was added to the Methods section. With the calculation of the standard deviations across the constitutive datasets of four metrics that described points (2) to (5) above, it could be concluded that in all global datasets, the constitutive studies (1) provided enough data points to have a significant impact on the regressions fitted to the global datasets, and reported data distribution that (2) spanned over similar normalized property ranges, (3) were centred around similar means, (4) displayed similar data dispersion, and (5) displayed similar skewness, as visually displayed with frequency histograms in Appendix Figure 1 of Appendix 1. These results were added to the Results section and the related graphs were added to the Appendix 1 section in Appendix Figure 2 and Appendix Figure 3. The 95% confidence intervals of the regressions were then added to the plots in Figure 2, and the ranges of the parameters of the mathematical equations added to Table 3. The reported narrow ranges of these intervals and the low variability of the regression parameters demonstrated that the constitutive studies showed little interstudy variability in associating the MN property distributions. The low observed inter-study variability made it pertinent to merge the normalized data from the selected studies into global datasets.

Further, we considered important to assess the variability between the distributions of the properties that appeared in different global datasets, for example the variability of the Ith distributions between the global datasets {Ith; R}, {Ith; ACV}, {Ith; AHP} and {τ; Ith}. This was assessed with the standard deviation across global datasets of the same four metrics indicated before, as now described in the Methods Section. Low variability was observed, which made the global datasets fit to be processed and transformed into final size-dependent datasets.

7. It is not clear why the authors used a 70-30% scheme for cross-validation and why only 5 validations were performed. A larger set of validations are probably needed and at least another training-testing proportion should be carried out.

Thank you for this comment. In the previous version of the manuscript, 5 iterations of 70-30 random splits of the global datasets were performed to constitute complementary training and test sets. The normalized relationships were derived from the training sets and validated against the test sets with three metrics (nME, NRMSE, r2) which were averaged across the five trials.

Rather than pursuing with this non-standard approach of randomly 70-30 splitting the data between the 5 iterations and improving this approach by repeating the shuffle-and-split process with different train-test ratios or an increased number of random splits, we have preferred to change the validation method and have opted for a standard and widely-accepted K-fold cross-validation approach, which is now described in the Methods from L.329 to L.349 and applied in the Results section from L.591 to L.599. The validation plots in Figure 6 were updated accordingly.

In brief, we chose K=5 for the K-fold cross-validation advised in ‘Deep Learning in Python’ from François Chollet (p. 87 section 3.6.4 Validating your approach using K-fold validation) when validating models on rather small datasets, as in our study (16 of the 17 datasets have less than 500 points). This 5-fold procedure is also the default validation method in the scikit-learn toolbox in Python. In this new approach, the data of each global dataset was shuffled once. Then, 5 non-overlapping partitions p1to p5 were created in each dataset, each containing 20% of the data (see Author response table 1). These partitions were permutated to create five different pairs of one training set including four partitions and one test set including one partition. For each permutation, and as before, the training set was used to derive the final relationships as described in the manuscript in the Methods. The final relationships were then validated against the corresponding test set. The final metrics (nME, nRMSE, r2) were averaged across the five permutations.

Author response table 1
For one {A; B} global datasetTraining setTest set
Permutation 1p1, p2, p3, p4p5
Permutation 2p2, p3, p4, p5p1
Permutation 3p1, p3, p4, p5p2
Permutation 4p1, p2, p4, p5p3
Permutation 5p1, p2, p3, p5p4

The new validation results are strongly in agreement with the ones reported with the previous approach, as shown in Author response table 2.

Author response table 2
DatasetsPrevious validationNew validation
nMEnRMSER2nMEnRMSER2
AHP = ka ∙ SMNa70140,4067120,46
AHP = ka ∙ ACVa158200,43147190,43
R = ka ∙ SMNa255200,56312310,51
R = ka ∙ ACVa246230,43205180,46
R = ka ∙ AHPa153210,63143210,66
Ith = ka ∙ Ra410190,37300190,36
Ith = ka ∙ ACVa171260,34155200,42
Ith = ka ∙ AHPa84150,5978170,61
C = ka ∙ Ra61120,5855130,57
C = ka ∙ ITaℎ47170,5352190,49
C = ka ∙ AHPa75160,2662160,22
C = ka ∙ ACVa82210,1674210,20
τ = ka ∙ Ra86130,5473130,51
τ = ka ∙ AHPa88140,6277130,64
τ = ka ∙ ITaℎ83150,6865130,74
τ = ka ∙ ACVa81170,3884160,44

8. Tables 4, 5 and 7 could be moved to the Supplementary material section.

We thank the Reviewing Editor and Reviewer 1 for their advice for a clearer and simpler framing of the paper. We modified the manuscript as suggested by the Reviewer and moved these three tables to Appendix 1, now Tables 3 to 5. The references in the text to these tables were adapted accordingly.

9. In the merged data sets shown in Figure 4 it would be useful to use different symbols for each data set (i.e., different symbols for {AHP;Smn}1 and {AHP;Smn}2).

Thank you for this comment which made Figure 4 more informative. The figure was adjusted, as suggested, and the sub-datasets constituting the final {A; SMN} datasets were identified with symbols specified in the caption of Figure 4 (L.576-577).

10. There is no oval box with fc(C) in Figure 1. Please, revise. Also, consider changing Figure 1 with a pseudo-code of the algorithm to improve clarity.

In line with the first comment of the reviewer, Figure 1 was revised to be a more detailed box-and-arrow diagram with a detailed caption explicating the step-by-step process applied to a specific example. The oval box issue was corrected.

11. Ith was estimated using different approaches in the selected studies. Some studies have used triangular-shaped currents, while others used step currents. How would different methods for estimating similar parameters influence your analysis (inter-study variability)?

As advised in the fourth comment from the Reviewing Editor and in the reviewer 2’s review report, we added in the revised manuscript a review of the experimental approaches taken by the selected studies for the measurements of the electrophysiological properties investigated in the manuscript (L.148 to L.177). When screening the selected studies, all the selected studies which provided Ith data for the global datasets in Figure 3 used step currents to estimate Ith.

– Kernell (1966): ‘cells were stimulated to steady repetitive firing by long-lasting injected currents’

– Fleshman (1981): ‘Rheobase, defined as the minimal 50-ms current pulse necessary to produce an action potential, was measured by slowly increasing the current intensity until discharge occurred intermittently’.

– Kernell and Monster (1981): ‘The threshold current for maintained repetitive impulse firing was determined by aid of long-lasting injected currents’

– Gustafsson and Pinter (1984b): ‘The current magnitude necessary to depolarize from resting potential to the critical firing level using long (50-60 ms) current pulses injected through the recording micro-electrode’

– Zengel et al. (1985): ‘Rheobase was determined by slowly increasing the intensity of long (50 or 100 ms) current pulses until discharge occurred intermittently.’

– Munson et al. (1986) – Foehring et al. (1987): similar protocol as in Zengel et al. (1985)

– Krawitz et al. (2001): ‘Rheobase was defined as the minimum amplitude of a depolarizing (50 ms duration) current pulse that evoked an action potential.’

However, it must be noted that Krawitz et al., (2001) reported that ‘In two cells rheobase was estimated using a slowly rising current ramp’. However, the influence on the analysis of these two Ith data points being obtained with a different experimental approach, among the hundreds included in the global {Ith; R} dataset, is negligeable. Therefore, this different experimental approach was not discussed in the manuscript.

12. Some relations are not unexpected. The proportional relation between cell capacitance and cell size is obvious. Again, you should consider previous theoretical studies based on Rall's theory.

Thank you for this comment, which echoes with Comment #4 and other comments by Reviewer 3. We realise that the Discussion section of the first version of the manuscript was not framed correctly and could have misled the reader to think that we were stating that we obtained unexpected and ‘new’ relationships when actually referring to well-known associations, such as the proportional relation between cell capacitance and cell size. We now avoid this issue in two ways. First, by providing in the Introduction section a list of well-known empirical associations, as discussed in Comment #2. And second by assessing in detail in the Discussion section in a rearranged and newly-named sub-section ‘Consistency of the relationships with previous empirical results and Rall’s theory’ (L.698 to L.760), the consistency of various combinations of the mathematical relationships in Table 4 with some of the equations that describe Rall’s theory of a MN equivalent cable model:

CCmSMN

τ=RmCm

C=τRLtanh(L)

R=RmSneuronLtanh(L)

Yet, we would like to highlight that the proportional relation between cell capacitance and cell size obtained in Table 4, while not unexpected in the literature, was not necessarily evident in the procedure taken in this study (L.748 to L.753). In the global datasets, the values of C were calculated as C=τRLtanh(L) and empirically related to the four parameters ACV, AHP, R, Ith, i.e. never to direct measures of SMN. Then, the inverted relationships SMN = f(ACV, AHP, R, Ith) obtained from 17 of the selected studies were used to derive the final {C; SMN} dataset, which happened to yield a proportional relationship, consistently with Rall’s theory, and without relying on direct measurements of SMN.

13. The authors should provide more clear evidence on how the proposed method should be translated to future studies on synaptic integration. Since the analysis did not consider any active property, I am not confident that the mathematical relations can help a more comprehensive computer simulation study.

This is an important comment concerning one of the proposed applications of this study. Consistently with one of Reviewer 3’s comments, we added in the Limitations section between L.924 and L.948 that this study mainly focused on passive MN electrophysiological properties; the membrane properties R, τ, C were measured with weak sub-threshold current pulses. Other properties like ACV, AHP, Ith were however measured in a functional context of elicited action potentials. Yet, some active MN membrane properties, such as PICs, were disregarded in this work, as they were not quantified in the selected studies simultaneously with the 9 MN properties investigated in this work. While quantitative relationships like those in Table 4 are missing for such active properties, the results from Table 4 can still be used in more comprehensive computer simulation studies involving Hodgkin-Huxley-like models for example, if they are linked with other observations in the literature that involve active properties, such as MNs of small SMN having longer lasting total PICs of greater hyperpolarized activation than MNs of large SMN (Heckman and Enoka, 2012).

The relationships in Table 4 can however be directly used in standard phenomenological MN models that typically neglect these active properties in a functional context, like the phenomenological RC LIF-like models Reviewer 3 refers to, which are known for their capacity to predict spiking behaviour (Teeter et al., 2018). Like in Watanabe et al., (2011) or Negro and Farina (2011) which involve comprehensive MN models, some previous studies using RC phenomenological models have attempted to build profiles of interconsistent MN-specific properties using maximum likelihood methods (Dong et al., 2011) and to describe a realistic continuous distribution of the MN properties in the MN pool (Negro et al., 2016) to estimate the proportion of common input to a motor neuron pool in humans. Such studies, and others that experimentally measured the distribution in the MN pool of some of the MN properties from Table 1 for modelling purposes (Raikova et al., 2018), would gain from the quantitative relationships provided in Table 4. An example of a direct application of Table 4 for modelling purposes is provided in Caillet et al., (2022) (https://www.biorxiv.org/content/10.1101/2022.02.21.481337v1.abstract), where the relationships in Table 4 are used to scale, after SMN calibration, a cohort of LIF models to predict the firing behaviour of a theoretical MN pool from a subset of firing MNs identified with decomposed HDEMGs. The predicted spike trains then drive an N-MUs Hill-type model to predict the isometric muscle force trajectory.

14. Please use ACV or CV, not a mixture of both abbreviations.

Thank you for having spotted these inconsistencies. The manuscript was revised to maintain the ‘ACV’ abbreviation throughout the manuscript.

15. The section on specific resistance (starting on Line 536) should be expanded to consider the fact that specific resistance may be different in different parts of the neuron (see Public review comments). This is important for modelers who want to explicitly represent dendrites in their motoneuron models (either as a separate compartment, an equivalent cable or completely reconstructed trees).

Both this and Reviewer 2’s are important comments because a non-constant Rr value across the MN membrane is a limitation to using Rall’s simplest model of a MN (soma + dendrites) as a unique equivalent cylinder of constant diameter. In this view, the assumption of isopotentiality across the membrane surface was repeated in the manuscript whenever Rall’s approach was considered. We now discuss the variability of Rr across the somatodendritic surface and the subsequent limitations to the relationships in Table 4 in the Discussion section from L.731 to L.739 and in the Limitations Section from L.874 to L.880, with mention to MN modelling.

It should be noted that, after the comments from Reviewer 1 for a simpler framing of the paper, this section on Rr was moved from the Results section to the Discussion section (L.713 to L.739). In this updated paragraph, some sentences were left unchanged since the original version, others were rephrased or added.

16. The Discussion section on relevance for modelers contains the following statement: "This supports the conclusions that the relative voltage threshold is constant within the MN pool and that Ohms law is followed in MNs". This is not strictly correct. Gustafsson and Pinter (1984b) showed that voltage threshold tended to be lower in low input conductance cells with long AHPs (see their Figure 4C and D). Also, the voltage threshold predicted by Ohms Law from the product of the measured rheobase and input resistance Vth = Ith*R tended to be lower than the measured voltage threshold, suggesting the activation of an inward conductance near threshold. This requires a variant of Ohm's Law in which the effective resistance is voltage-dependent: Vth=Ith*R(v).

The authors thank the Reviewer for this comment that helped in enriching the discussion of the relationships in Table 4 (L.761 and L.781 of the Discussion section). It is now stated that, consistently with the relationships in Table 4 and in first approximation, Ohm’s law is followed in MNs when MNs are excited by weak subthreshold currents. This is how R was measured in the selected studies to avoid voltage-activated membrane nonlinearities, as discussed in the Methods section. It is then added in the manuscript that Vth may however be MN size-dependent in the MN pool, whatever the strength of the excitation current. It is finally mentioned that Vth is underestimated by the product RIth near and above threshold, due to voltage-activated inward conductance. It is concluded that Ohm’s law must be revised both in the passive and the functional context to be a function of MN voltage and size.

17. Finally, although the vast majority of readers will know that correlation does not imply causation, it is worth stating explicitly in the Discussion section on limitations.

As advised by the Reviewer, this was added as the concluding limitation of the Limitations Section. L.949951:

‘As a final limitation, the relationships in Table 4 were obtained from a regression analysis and therefore provide correlations between some MN properties in a MN pool, but cannot be used to draw conclusions on the causality behind these associations.’

18. p. 33, l. 721: Actually, you should consider that in some experimental conditions, large (high-threshold) MNs attain lower firing rates than small (low-threshold) MNs. Maybe you have to explain in what conditions large MNs attain high firing rates.

Thank you for this comment. To remove any confusion regarding the ‘onion skin’ effect, it was added L.841843 that large MNs could attain higher firing rates than small MNs in the events of ballistic contractions or force generations close to maximum voluntary contractions.

References

Caillet, A., Phillips, A. T., Farina, D., and Modenese, L. (2022). Estimation of the firing behaviour of a complete motoneuron pool by combining EMG signal decomposition and realistic motoneuron modelling. bioRxiv.

Chollet F (2021) Deep Learning with Python. Simon and Schuster

Dong, Y., Mihalas, S., Russell, A., Etienne-Cummings, R., and Niebur, E. (2011). Estimating parameters of generalized integrate-and-fire neurons from the maximum likelihood of spike trains. Neural computation, 23(11), 2833-2867.

Eccles JC, Eccles RM, Lundberg A (1958) The action potentials of the α motoneurones supplying fast and slow muscles. J.Physiol.(Lond.) 142:275

Fleshman JW, Segev I, Burke RB (1988) Electrotonic architecture of type-identified α-motoneurons in the cat spinal cord. J.Neurophysiol. 60:60-85

Fleshman JW, Munson JB, Sypert GW, Friedman WA (1981) Rheobase, input resistance, and motor-unit type in medial gastrocnemius motoneurons in the cat. J Neurophysiol 46:1326-1338

Foehring RC, Sypert GW, Munson JB (1987) Motor-unit properties following cross-reinnervation of cat lateral gastrocnemius and soleus muscles with medial gastrocnemius nerve. II. influence of muscle on motoneurons. J Neurophysiol 57:1227-1245.

Gustafsson, Pinter MJ (1984a) An investigation of threshold properties among cat spinal α‐motoneurones. J. Physiol. 357:453-483.

Gustafsson, Pinter MJ (1984b) Relations among passive electrical properties of lumbar α-motoneurones of the cat. J. Physiol. 356:401-431.

Heckman CJ, Enoka RM (2012) Motor unit. Compr. Physiol. 2:2629-2682

Huh, S.,Heckman, C.J.,Manuel, M.(2021) Time course of alterations in adult spinal motoneuron properties in the SOD1(G93A) mouse model of ALS. eNeuro, 8, ENEURO

Kernell D, Monster AW (1981) Threshold current for repetitive impulse firing in motoneurones innervating muscle fibres of different fatigue sensitivity in the cat. Brain Res. 229:193-196

Kernell D (1966) Input resistance, electrical excitability, and size of ventral horn cells in cat spinal cord. Science 152:1637-1640

Krawitz S, Fedirchuk B, Dai Y, Jordan LM, McCrea DA (2001) State-dependent hyperpolarization of voltage threshold enhances motoneurone excitability during fictive locomotion in the cat. The Journal of physiology 532:271-281

Krutki P, Hałuszka A, Mrówczyński W, Gardiner PF, Celichowski J (2015) Adaptations of motoneuron properties to chronic compensatory muscle overload. Journal of neurophysiology 113:2769-2777.

Manuel M, Iglesias C, Donnet M, Leroy F, Heckman CJ, Zytnicki D (2009) Fast kinetics, high-frequency oscillations, and subprimary firing range in adult mouse spinal motoneurons. The Journal of neuroscience 29:11246-11256

Munson JB, Foehring RC, Lofton SA, Zengel JE, Sypert GW (1986) Plasticity of medial gastrocnemius motor units following cordotomy in the cat. Journal of neurophysiology 55:1454

Negro F, Yavuz U, Farina D (2016) The human motor neuron pools receive a dominant slow‐varying common synaptic input. J. Physiol. 594:5491-5505.

Negro, F., and Farina, D. (2011). Linear transmission of cortical oscillations to the neural drive to muscles is mediated by common projections to populations of motoneurons in humans. The Journal of physiology, 589(3), 629-637.

Powers RK, Binder MD (2001) Input-output functions of mammalian motoneurons. Rev. Physiol. Biochem. Pharmacol. 143:137-263

Raikova, R., Celichowski, J., Angelova, S., and Krutki, P. (2018). A model of the rat medial gastrocnemius muscle based on inputs to motoneurons and on an algorithm for prediction of the motor unit force. Journal of Neurophysiology, 120(4), 1973-1987.

Sasaki M (1991) Membrane properties of external urethral and external anal sphincter motoneurones in the cat. The Journal of physiology 440:345-366

Teeter C, Iyer R, Menon V, Gouwens N, Feng D, Berg J, Szafer A, Cain N, Zeng H, Hawrylycz M, Koch C, Mihalas S

(2018) Generalized leaky integrate-and-fire models classify multiple neuron types. Nat Commun 9:1-15

Ulfhake B, Kellerth JO (1984) Electrophysiological and morphological measurements in cat gastrocnemius and soleus α-motoneurones. Brain Res. 307:167-179

Watanabe, R. N., Magalhães, F. H., Elias, L. A., Chaud, V. M., Mello, E. M., and Kohn, A. F. (2013). Influences of premotoneuronal command statistics on the scaling of motor output variability during isometric plantar flexion. Journal of neurophysiology, 110(11), 2592-2606.

Zengel JE, Reid SA, Sypert GW, Munson JB (1985) Membrane electrical properties and prediction of motor-unit type of medial gastrocnemius motoneurons in the cat. J Neurophysiol 53:1323-1344

https://doi.org/10.7554/eLife.76489.sa2

Article and author information

Author details

  1. Arnault H Caillet

    Department of Civil and Environmental Engineering, Imperial College London, London, United Kingdom
    Contribution
    Conceptualization, Data curation, Formal analysis, Validation, Methodology, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6146-1829
  2. Andrew TM Phillips

    Department of Civil and Environmental Engineering, Imperial College London, London, United Kingdom
    Contribution
    Supervision, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6618-0145
  3. Dario Farina

    Department of Bioengineering, Imperial College London, London, United Kingdom
    Contribution
    Conceptualization, Supervision, Writing – review and editing
    Contributed equally with
    Luca Modenese
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-7883-2697
  4. Luca Modenese

    1. Department of Civil and Environmental Engineering, Imperial College London, London, United Kingdom
    2. Graduate School of Biomedical Engineering, University of New South Wales, Sydney, Australia
    Contribution
    Supervision, Methodology, Writing – review and editing
    Contributed equally with
    Dario Farina
    For correspondence
    l.modenese@unsw.edu.au
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-1402-5359

Funding

Imperial College London (Skempton Scholarship)

  • Arnault H Caillet

Imperial College London (Imperial College Research Fellowship)

  • Luca Modenese

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

AC was supported by the Skempton Scholarship and LM by an Imperial College Research Fellowship granted by Imperial College London. The authors want to thank Dr Marin Manuel for making available the experimental data published in Huh et al., 2021, Time course of alterations in adult spinal motoneuron properties in the SOD1 (G93A) mouse model of ALS, Eneuro.

Senior Editor

  1. Barbara G Shinn-Cunningham, Carnegie Mellon University, United States

Reviewing Editor

  1. Hugo Merchant, National Autonomous University of Mexico, Mexico

Reviewers

  1. Randall K Powers, University of Washington, United States
  2. Leonardo Abdala Elias, University of Campinas - UNICAMP, Brazil

Publication history

  1. Preprint posted: August 5, 2021 (view preprint)
  2. Received: December 17, 2021
  3. Accepted: July 13, 2022
  4. Accepted Manuscript published: July 18, 2022 (version 1)
  5. Version of Record published: October 27, 2022 (version 2)

Copyright

© 2022, Caillet et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 542
    Page views
  • 219
    Downloads
  • 2
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Arnault H Caillet
  2. Andrew TM Phillips
  3. Dario Farina
  4. Luca Modenese
(2022)
Mathematical relationships between spinal motoneuron properties
eLife 11:e76489.
https://doi.org/10.7554/eLife.76489

Further reading

    1. Developmental Biology
    2. Neuroscience
    Emily L Heckman, Chris Q Doe
    Research Advance

    The organization of neural circuits determines nervous system function. Variability can arise during neural circuit development (e.g. neurite morphology, axon/dendrite position). To ensure robust nervous system function, mechanisms must exist to accommodate variation in neurite positioning during circuit formation. Previously we developed a model system in the Drosophila ventral nerve cord to conditionally induce positional variability of a proprioceptive sensory axon terminal, and used this model to show that when we altered the presynaptic position of the sensory neuron, its major postsynaptic interneuron partner modified its dendritic arbor to match the presynaptic contact, resulting in functional synaptic input (Sales et al., 2019). Here we investigate the cellular mechanisms by which the interneuron dendrites detect and match variation in presynaptic partner location and input strength. We manipulate the presynaptic sensory neuron by (a) ablation; (b) silencing or activation; or (c) altering its location in the neuropil. From these experiments we conclude that there are two opposing mechanisms used to establish functional connectivity in the face of presynaptic variability: presynaptic contact stimulates dendrite outgrowth locally, whereas presynaptic activity inhibits postsynaptic dendrite outgrowth globally. These mechanisms are only active during an early larval critical period for structural plasticity. Collectively, our data provide new insights into dendrite development, identifying mechanisms that allow dendrites to flexibly respond to developmental variability in presynaptic location and input strength.

    1. Epidemiology and Global Health
    2. Neuroscience
    Lorenza Dall'Aglio, Hannah H Kim ... Henning Tiemeier
    Research Article Updated

    Background:

    Associations between attention-deficit/hyperactivity disorder (ADHD) and brain morphology have been reported, although with several inconsistencies. These may partly stem from confounding bias, which could distort associations and limit generalizability. We examined how associations between brain morphology and ADHD symptoms change with adjustments for potential confounders typically overlooked in the literature (aim 1), and for the intelligence quotient (IQ) and head motion, which are generally corrected for but play ambiguous roles (aim 2).

    Methods:

    Participants were 10-year-old children from the Adolescent Brain Cognitive Development (N = 7722) and Generation R (N = 2531) Studies. Cortical area, volume, and thickness were measured with MRI and ADHD symptoms with the Child Behavior Checklist. Surface-based cross-sectional analyses were run.

    Results:

    ADHD symptoms related to widespread cortical regions when solely adjusting for demographic factors. Additional adjustments for socioeconomic and maternal behavioral confounders (aim 1) generally attenuated associations, as cluster sizes halved and effect sizes substantially reduced. Cluster sizes further changed when including IQ and head motion (aim 2), however, we argue that adjustments might have introduced bias.

    Conclusions:

    Careful confounder selection and control can help identify more robust and specific regions of associations for ADHD symptoms, across two cohorts. We provided guidance to minimizing confounding bias in psychiatric neuroimaging.

    Funding:

    Authors are supported by an NWO-VICI grant (NWO-ZonMW: 016.VICI.170.200 to HT) for HT, LDA, SL, and the Sophia Foundation S18-20, and Erasmus University and Erasmus MC Fellowship for RLM.