Mathematical relationships between spinal motoneuron properties

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.


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 ( I th ) of an MN are both positively correlated to the axonal conduction velocity ( ACV ) and afterhyperpolarization duration ( AHP ). Also, I th 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 ).
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 , 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 andKernell, 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 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.
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 Table 1. The motoneuron (MN) and muscle unit (mU) properties investigated in this study with their notations and SI base units. S MN is the size of the MN. As reproduced in Table 2, the MN size S MN is adequately described by measures of the MN surface area Sneuron and the soma diameter Dsoma . R and Rm define the MNspecific electrical resistance properties of the MN and set the value of the MN-specific current threshold I th (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). ∆V th is the amplitude of the membrane voltage depolarization threshold relative to resting state required to elicit an action potential. I th is the corresponding electrical current causing a membrane depolarization of ∆V th . 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. S mU is the size of the mU. As indicated in Table 2, the mU size S mU is adequately described by measures of (1) the sum of the cross-sectional areas (CSAs) of the fibres composing the mU CSA tot , (2) the mean fibre CSA CSA mean , (3) the innervation ratio IR, that is, the number of innervated fibres constituting the mU, and (4) the mU tetanic force F tet . F tw is the MU twitch force. Twitch force F tw [N] 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, I th was obtained as the minimal intensity of a 50-100-ms-long depolarizing rectangular current pulse required to produce regular discharges. I th 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 I th . 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 ( I in ) 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 − I in plots. In all studies, R was obtained by calculating the slope of the linear part of the Vm − I in 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 dVm dt 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): 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 Binder, 2001 andOlinder 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 − I in 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 ln ( 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 dV dt -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 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, I th , 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 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 ln transformation of each global dataset yielding relationships of the type ln ( A ) = a · ln ( B ) + k , which were then converted into power relationships of the type A = k · B a (eq. 1). The adequacy of these global power trendlines and the statistical significance of the correlations were assessed with the coefficient of determination r 2 (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 r 2 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 crossstudy analysis. For this preliminary analysis and the rest of the study, power regressions (eq. 1) were preferred to linear ( A = k + a · B ) or exponential ( A = k · e a·B ) 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 = k · I −1 th . Also, a least-squares linear regression analysis was preliminary performed for each normalized experimental dataset and for its ln and ln ( A ) − B transformations, yielding linear, power, and exponential fits to the data. The power regressions overall returned r 2 > 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 Me Md 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 Me Md metrics, that is, when sdg < 10% and sdg meang < 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 = k · B a . 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, Me Md ) for each global dataset, and then computing the mean G ± sd G 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 sd G < 10% and sdG meanG < 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 S MN 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 S MN 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, I th , 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 S MN . For each electrophysiological property A and each { A; B } dataset, we considered two cases. If B = S MN , the global dataset was not further processed as the electrophysiological property A was already related to MN size S MN from experimental measurements.
+ kc , which were converted into the power relationships: 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 r 2 (squared value of the Pearson's correlation coefficient) and a threshold of 0.01 on the p-value of the regression analysis, respectively. If r 2 > 0.3 (i.e., r > 0.55) and p<0.01, another power trendline S MN = k d · A d was fitted to the final dataset to describe the inverse relationship for { S MN ; A } and used in the processing of the next-in-line property.

Normalized mathematical relationships between electrophysiological properties
When a final relationship with S MN was obtained for two MN properties A and B by the procedure described above, that is, A = k cA · S cA MN and B = k cB · S cB MN , one of the expressions was mathematically inverted and a third empirical relationship was derived for the property pair 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 Figure 1. Detailed example for the process adopted to successively create the two final {R; SMN} and {I th ; 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) 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 = kc · S c MN and relationships between electrophysiological properties A = ke · B e obtained with the training set were applied to the S MN 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 r 2 pred 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 r 2 pred values were finally averaged across the five permutations. For each global dataset, the average r 2 pred was compared to the r 2 exp 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 = kc · S c MN and relationships between electrophysiological properties A = ke · B e were computed for the complete global datasets.

Scaling of the normalized relationships
The A = kc · S c MN and A = ke · B e 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. S MN varied over a q E S -fold range, taken as the average across cat studies of the ratios of minimum and maximum values measured for S MN : An experimental ratio q E A was similarly obtained for each electrophysiological property A . As was similarly built over the q E S -fold range previously derived. Finally, the intercept kc in the size-dependent relationships A = kc · S c MN was scaled as A similar approach was used to mathematically derive the intercept ke and scale the relationships A = ke · B e between electrophysiological properties.

Extrapolation to rats and mice
It was finally assessed whether the A = kc · S c MN and A = ke · B e scaled relationships derived from cat data accurately predicted rat and mouse quantities reported in other experimental studies. These mathematical relationships were applied to the S MN 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 r 2 between predicted and experimental values of A.

Definitions of MN size S MN and mU size S mU MN size S MN
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 S MN = Dsoma in the described methodology. In the literature on spinal MNs in adult mammals, the size of an MN S MN 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 D dendrite (Amendola and Durand, 2008;Carrascal et al., 2009;Mantilla et al., 2018), individual dendritic surface area S dendrite (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  defined as the summed soma and dendritic surface areas.
In the references Ulfhake and Kellerth, 1981), a linear correlation between Dsoma and the average diameter of the stem dendrites D dendritemean 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 D dendrite and the membrane surface area of the corresponding dendritic tree S dendrite (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 S dendritemean . Moreover, the number of dendritic trees N dendrites per cell increases with increasing soma surface Ssoma (Brandenburg et al., 2018), and a linear correlation between Dsoma or S dendrites with N dendrites has been observed (Zwaagstra and 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 S tot dendrite are linearly correlated. This assumption/approximation is consistent other conclusions of a linear correlations between Dsoma and S tot dendrite (Amendola and Durand, 2008;Brandenburg et al., 2020). Then, according to typical values of Dsoma , Ssoma , and S tot dendrites obtained from the studies previously cited, yielding S tot dendrites ≫ Ssoma , it was also assumed S tot dendrites ≈ 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 = kc · S c MN derived previously, with S MN = Dsoma , were extrapolated with a gain to relationships between A and S MN = Sneuron , notably permitting the definition of surface specific resistance Rm and capacitance Cm . Following the same method as described above, theoretical ranges were derived from additional morphometric studies on adult cat spinal alpha-MNs. The new relationships were extrapo- min · S c neuron . It must however be highlighted that the conclusion of a linear correlation between Dsoma and S neuron , 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 S dendrite . 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 N dendrites 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 S mU
To enable future comparisons between MN and mU properties, we here assess potential indices of the mU size S mU suggested in the literature. The size of an mU ( S mU ) can be defined as the sum CSA tot of the CSAs of the innervated fibres composing the mU. CSA tot depends on the mU innervation ratio ( IR ) and on the mean CSA ( CSA mean ) of the innervated fibres: S mU = CSA tot = IR · CSA mean . CSA tot 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 Rafuse et al., 1997). The mU tetanic force F tet 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 F tet = σ · IR · CSA mean (Burke, 1981;Enoka, 1995) returns a linear correlation between F tet and IR · CSA mean = S mU 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 F tet and IR and between F tet and CSA mean . Consequently, F tet , IR , and CSA mean are measurable, consistent, and linearly related measures of S mU in animals, as summarized in Table 2.

Relationships between MN and mU properties
To assess whether the size-dependent relationships A = kc · S c MN 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 B MN (Table 1) and an mU property A mU for the same MU. The normalized global datasets obtained for the pairs { A mU ; B MN } were fitted with power trendlines, as previously described for MN properties, yielding A mU = k · B b MN relationships. Using both the definition of S mU ( Table 2) and the A = kc · S c MN relationships derived previously, the 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.

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 included the data from two or more experimental studies.
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 Me Md 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 sdg meang < 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 sdg meang < 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 = ka · B a (p-value<10 -5 and r 2 ∈  Figure 3, the confidence interval remained narrow for all datasets with the value of power a varying less than 0.4, except for the { I th ; 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 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.
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 S MN , ACV, AHP, R, I th , 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, sd G < 10% and sdG meanG < 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.

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, I th = 9.4 · 10 −4 · S 2.5 MN and R = 9.6 · 10 5 · S −2.4 MN were obtained in Table 3. When mathematically combined (with the latter mathematically inverted as S MN = 2.9 · 10 2 · R −0.4 ), these relationships yielded the normalized relationship between rheobase and input resistance I th = 1.5 · 10 3 · R −1 . All  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 = kc · S c MN (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 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 r 2 pred 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 r 2 pred between 0.20 and 0.74 and of the same order of magnitude as the experimental r 2 exp values.

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 q E S = 2.4 -fold range according to a review of the morphometric studies reported in the two first lines of Appendix 1-table 3. q E S may however be underestimated for reasons discussed in Appendix 1. Then, the empirical q E A and theoretical q T A 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 q E R = 10.9 -fold range in an MN pool according to the literature, while the theoretical fold range q T R = ( q E S ) cR = 2.4 2.4 = 8.4 was obtained from the R = 9.6 · 10 5 · S −2.4 MN normalized relationship previously derived when a 2.4-fold range is set  for the MN properties ACV, AHP, C , and τ , while the theoretical ranges for R and I th span over a narrower domain than expected from experimental results ( q T R = 8.4 ≪ 10.9 = q E R and q T Ith = 9.1 ≪ 16.3 = q E I th ). This suggests that the range of variation of R and I th is not entirely explained by the variation in MN size S MN , and that MN excitability is not only determined by S MN , 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-  (Appendix 1-table 4), it was directly obtained from 'Methods' that 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, I th , C, τ , and S MN . All constants and relationships are given in SI base units.

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 . 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 I th distributions reported by different studies in the rat and mouse { I th ; 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 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 {I th ; 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).
The online version of this article includes the following source data for figure 7: Source data 1. Datasets used in the rat plot.
Source data 2. Datasets used in the mouse plots.

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 Enoka, 2012 andManuel 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 ( r 2 = 0.43 ). Eight studies, dominantly focussing on the cat soleus and medial gastrocnemius muscles, found a strong correlation between F tet and ACV , while one study showed a significant correlation for the pair { F tet ; R } in both the cat tibialis anterior and extensor digitorum longus muscles. One study  on the cat soleus, medial and lateral gastrocnemius muscles inferred a statistically significant correlation for { F tet ; S MN } . In mice, both R and I th were significantly related to muscle tetanic and twitch forces. As F tet , F tw , and IR are reliable indices of S mU as discussed in 'Methods, the MN size-dependent relationships in Table 3 return eight S mU = k · S c MN 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.

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.  Mcphedran 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., 1988 5.0 The online version of this article includes the following source data for table 5: Source data 1. Numerical data used to derive the relationships presented in Table 5.

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 I th = k · τ −1.7 , 1 C = k · τ 0.7 , 1 SMN = k · ACV −1.4 , and 1 R = k · ACV 3.5 in Table 4 yield, when combined, I th C = k · τ −1 and 1 R·SMN = k · ACV 2.1 , which are consistent with the relationship I th C = k · τ −1 experimentally observed by Gustafsson and Pinter, 1985 and 1 R·SMN = k · ACV 1.9 measured by Kernell and Zwaagstra, 1981. The stronger-than-linear inverse relationship R = k · S −2.4 MN 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 range reported in the literature (Powers and Binder, 2001), (eq. 1. in Gustafsson and Pinter, 1984a) is obtained when combining the R − Sneuron, C − Sneuron , and τ − Sneuron relationships. Also, taking tanh ( L ) L = 0.6 (standard value provided in Powers and Binder, 2001), the R − Sneuron equation in Table 4 Rall, 1977,  ] Ω · m 2 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.0·10 −10 S 1.4 neuron yields Rm = 2.5 · AHP 0.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.6Rm Sneuron and Rm = 1.0·10 −10 S 1.4 neuron 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 q E S of S MN does not entirely explain the ranges of variation q E R and q E Ith of R and I th , 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.0·10 −10 S 1.4 neuron , 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.8 · RC , close to the theoretical equation τ = RC · tanh ( L ) L derived from Rall's theory. Also, Table 4 directly yields, from the combination of Rm = 1.0·10 −10 S 1.4 neuron with the other relationships, τ = 1.4 · 10 −2 · Rm , which is exactly eq. 2.15 in Rall, 1977 ( τ = RmCm ) with a constant value for Cm = 1.4 · 10 −2 F · m −2 that is consistent with the literature. The relationship C = 1.3 · 10 −2 · Sneuron 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 S MN and C . The values of C were obtained from Rall's C = τ R · L tanh ( L ) equation (Gustafsson and Pinter, 1984b)  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.3 · 10 −2 F · m −2 constant value reported in this study is slightly higher than the widely accepted Cm = 1.0 · 10 −2 F · m −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 = C Cm with Cm = 1.0 · 10 −2 F · m −2 . Finally, the empirical relationship I th = 2.7·10 −2 R 1.0 ( Table 4) provides ∆V th = 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 ∆V th 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 ∆V th tend to be lower in in MNs of large R and long AHP (Gustafsson and Pinter, 1984a), suggesting that ∆V th might be inversely correlated to MN size S MN . 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 ∆V th are generally greater than that directly predicted from the product I th · R (Gustafsson and Pinter, 1984b), and a variant of Ohm's law such as ∆V th 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 I th . While ∆V th may be always size-dependent, and size-voltage-dependent close to threshold and during firing events, no correlation between ∆Vth R·Ith and C , AHP , or ACV was found and reported in Table 4 in subthreshold conditions, consistent with measurements performed by Gustafsson andPinter, 1984b andKellerth, 1984, substantiating that the dynamics of MN recruitment dominantly rely on R, I th , and ∆V th (Heckman and Enoka, 2012). 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, I th = 3.8 · 10 8 · S 2.5 neuron ( Table 4) quantitatively supports the size dependency of I th 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 sizerelated 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 r 2 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 = K · ke · B e 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 { I th ; 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 I th 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 S MN and S mU 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 S MN and S mU . Besides, considering the positive I th − S MN correlation ( Table 4), and that the mU force recruitment threshold F th is positively correlated to F tet (Heckman and Enoka, 2012) and thus to S mU (Table 2), larger MUs have both larger current and force recruitment thresholds I th and F th 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 S MN ∝ CV 1.4 ∝ τ −0.7 ∝ AHP −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 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 S mU and S MN 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 I th 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 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 r 2 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 S MN , ACV, AHP, R, I th , 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 I th , 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 S MN having longer lasting total PICs of greater hyperpolarized activation than MNs of large S MN (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, I th , C, ∆V th , 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 S MN , ACV, AHP, R, Rm, I th , 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. 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. The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.   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.

Author contributions
values reported in these studies, for the global dataset {ACV; S MN } . 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. The r 2 values returned by the three types of regression cannot be directly compared to estimate the best model. However, the power fit returned r 2 > 0.5 for relatively more experimental datasets than the linear and exponential fits.