Reverse engineering of motor unit discharge in multiple sclerosis reveals heterogeneity of voluntary motor commands

Reviewer #1 (Public review):

This study applies a recently developed reverse-engineering framework for motor unit discharge to a large dataset from individuals with multiple sclerosis (MS) and neurologically intact controls. The authors aim to determine whether abnormalities in voluntary motor control in MS can be attributed to different patterns of excitatory, inhibitory, and neuromodulatory input to spinal motoneurons. A major conceptual emphasis of the study is on heterogeneity: rather than asking only whether people with MS differ from controls on average, the authors examine whether the distributions of motor unit discharge features and derived physiological variables are broader and more diverse across affected individuals.

A major strength of the work is the size and richness of the dataset. The study includes 89 participants with MS and 34 controls, with high-density surface electromyography (HDsEMG) used to obtain large populations of motor unit discharge patterns from the tibialis anterior (TA) and soleus (SOL) muscles. The resulting dataset contains many thousands of motor unit recordings and allows the authors to examine both group means and the shapes and variances of participant-level distributions. The finding that several motor unit discharge characteristics are more broadly distributed in MS than in controls is convincing and potentially important. In particular, the observation that affected individuals can occupy both high and low extremes of these distributions provides a useful empirical description of the diversity of motor unit behavior in this disease.

The principal limitation concerns the physiological interpretation assigned to these discharge patterns. The excitation, inhibition, and neuromodulation variables are not directly measured physiological inputs. They are composite variables derived from several features of motor unit discharge, with each feature weighted according to relationships identified in simulations reported previously by Chardon and colleagues (reference 1). Thus, the step from observed discharge behavior to specific underlying synaptic mechanisms is necessarily model-dependent. The distinction between these two levels of inference is important for interpreting the main conclusions of the study.

This issue is especially relevant because several different physiological processes can plausibly influence the same discharge features. Motor unit firing patterns reflect not only excitatory and inhibitory synaptic inputs and neuromodulation, but also intrinsic motoneuron properties, persistent inward currents (PICs), after-hyperpolarization (AHP) properties, tonic inhibition, the time course of synaptic excitation, and afferent input. Some of these factors are acknowledged as limitations of the modeling framework in the manuscript. The current data therefore provide strong evidence for heterogeneous motor unit discharge phenotypes in MS, but more indirect evidence that this heterogeneity can be uniquely attributed to distinct patterns of excitatory, inhibitory, and neuromodulatory input.

The interpretation of inhibition is a particularly clear example of this general inverse problem. In the underlying modeling framework, inhibition is represented along a continuum from proportional or balanced inhibition to reciprocal or push-pull inhibition. These different patterns influence PICs and consequently alter firing-rate nonlinearity and rate modulation. This provides a plausible forward-model relationship between inhibitory organization and motor unit discharge. However, observing a particular firing pattern in vivo does not necessarily identify the organization of inhibitory input uniquely, because similar changes in firing-rate modulation or hysteresis could arise from altered neuromodulation, intrinsic motoneuron properties, or other changes in synaptic drive. The inhibition composite is therefore best interpreted as a discharge phenotype that is consistent with a particular inhibitory organization under the assumptions of the model, rather than as a direct measure of inhibitory synaptic input.

A related methodological issue concerns the construction of the composite variables. The authors use mutual-information (MI) values from the previous simulation study as weights in signed linear combinations of normalized discharge features. Mutual information quantifies how informative a feature is about a modeled input parameter, but it is not itself a regression coefficient or a measure of the magnitude of a physiological effect. In addition, different discharge features may contain overlapping information about the same underlying process. The resulting composites are therefore useful summary measures of patterns associated with the modeled physiological variables, but their quantitative interpretation as direct estimates of those variables is less certain. This distinction is particularly relevant because the manuscript sometimes moves from describing the composite variables to describing the corresponding physiological inputs themselves.

The study's emphasis on variability also raises an important measurement issue. Because increased between-participant variance is itself one of the central biological findings, differences in measurement precision between the MS and control groups are more consequential here than in a conventional comparison of group means. Participants with MS sometimes had greater difficulty producing smooth triangular contractions, and motor unit yield and decomposition quality may plausibly vary more across affected participants. If measurement or decomposition uncertainty were more heterogeneous in the MS group, this could broaden participant-level distributions and thereby amplify the appearance of biological heterogeneity. The manuscript uses established decomposition and quality-control procedures, so this is not a general challenge to the validity of HDsEMG. Rather, it is a consideration that is particularly important when increased distributional spread is itself the primary result.

The manuscript also makes a stronger interpretive step from broad group distributions to patient-specific pathophysiology. The data convincingly show that motor unit discharge-derived measures are more heterogeneous among people with MS. They do not yet establish whether this variation represents distinct mechanisms in individual patients, continuous variation in a common mechanism, identifiable pathophysiological subgroups, differences in disease severity or lesion distribution, or some combination of these factors. The manuscript itself recognizes this distinction when it identifies the separation of individual-level variation from potential subgroups as an important goal for future work.

Overall, this is a valuable study with an unusually large motor unit dataset and a compelling demonstration that motor unit discharge behavior is markedly heterogeneous in MS. The work also provides an informative application of a model-based reverse-engineering framework to a clinically diverse human population. The evidence is strongest for the descriptive conclusion that motor unit discharge phenotypes are heterogeneous and altered in MS. The more specific attribution of these phenotypes to excitatory, inhibitory, and monoaminergic inputs is plausible and potentially useful, but remains contingent on the assumptions and identifiability of the underlying model. With this distinction in mind, the dataset and analytical approach should be useful to researchers interested in motor unit physiology, disease-related variability in motor control, and the possibilities and limitations of inferring latent physiological mechanisms from human motor unit discharge.

Reference 1: Chardon, M. K. et al. Supercomputer framework for reverse engineering firing patterns of neuron populations to identify their synaptic inputs. eLife 12, RP90624 (2024).

Reviewer #2 (Public review):

The Beauchamp and Chardon model, like any computational model, is only as informative as the physiological parameters it includes. By varying only three synaptic related dimensions, namely neuromodulatory (i.e., PIC strength from 5HT inputs, etc.), the pattern of inhibition relative to excitation, and the distribution of excitatory input across smaller versus larger motoneurons, the model necessarily holds constant many other properties that may differ substantially between individuals and between controls and patients with MS. These properties include intrinsic membrane conductance properties (with some conductances like NaV that may differ in MS not included in the Chardon model), afterhyperpolarization properties, tonic inhibition, dendritic and axonal structure, synaptic kinetics, glial changes, etc., etc. The present paper then moves one step further away from direct physiological estimation because it does not fit these three model parameters to each participant's data, but instead combines control-normalized firing pattern features using weights derived from the original simulations. As a result, the resulting "excitation," "inhibition," and "neuromodulation" scores should be regarded as indirect similarity scores within a restricted model space, with a substantial risk that changes caused by unmodeled physiology are misattributed to one of the three modeled components. An additional limitation is that the underlying motoneuron models were originally tuned to intracellular recordings from medial gastrocnemius motoneurons in decerebrate cats and then manually modified to generate more human-like firing rates and hysteresis, rather than being formally fitted to human motor unit recordings.

The authors should justify their conclusions using the Chardon reverse engineering method. In addition, more of the raw firing rate profiles should be presented to give a better sense of the data and the quality of the motor unit identification.

Reviewer #3 (Public review):

Summary:

The authors set out to describe how the three basic ingredients of a voluntary motor command, that is excitation, inhibition and neuromodulation, are altered in people with multiple sclerosis. They recorded high-density surface electromyograms from tibialis anterior and soleus during slow triangular contractions in 89 patients and 34 control participants, decomposed the signals into the discharge times of individual motor units, and extracted seven features of the resulting discharge patterns. These features were then combined into three composite scores intended to represent the three ingredients, and the authors asked whether the scores in patients differ from those in controls in their average value, in their spread, and in the shape of their distribution. The working hypothesis was that the motor command is disturbed in different ways in different patients rather than in one consistent direction, and the authors report findings that they consider consistent with that hypothesis.

Strengths:

The dataset is the main strength of this work, and it is quite a considerable one. Studies of motor unit behaviour in neurological populations are usually built on ten or twenty participants, whereas here there are 89 patients spanning the full range of disability together with 34 controls, two muscles per person, and more than 12,000 unique motor units. The clinical characterisation is thorough and includes disease subtype, symptom duration, disability score, walking tests and lesion locations from clinical imaging. The experimental protocol is appropriate and was clearly demanding to deliver in a population of this kind. The treatment of the motor unit data in the statistical models is also more careful than is common in this literature, in particular the recognition that motor unit labels are arbitrary and the correct nesting of units within muscle and within participant. The question itself is well worth asking, and the central observation, that average discharge rates and rate modulation are reduced while a minority of patients sit far above anything seen in controls, is interesting.

Weaknesses:

Three issues limit how far the conclusions can be taken as the paper stands.

First, the framework that gives the paper its title is not the framework that was used. The published approach searches a large database of simulations to find the combination of synaptic inputs that reproduces an observed discharge pattern. What is applied here is a weighted average of standardised discharge features, in which the weights are mutual information values taken from that earlier work and the signs of the weights are supplied by the authors on the basis of the literature. Mutual information measures how much a feature tells you about a parameter. It does not carry the direction of that association, and it does not become a regression coefficient by having a sign attached to it. Since the direction of every conclusion in the paper depends on those assigned signs, the reader has no way of judging how faithfully the composite scores track the physiological quantities they are named after. I do not understand why this was not done with simulations. The authors should expand or clarify this point.

The simulations have known inputs, so the composites could be computed on the simulated discharge patterns and their accuracy reported directly. Two of the features were also modified relative to those simulations, being computed against joint torque rather than against the synaptic drive, and using a normalised version of the hysteresis measure, while the weights derived for the original features were retained. The authors should justify this or specify why this was not done and how it affects the underlying physiology.

A related difficulty is that the three composite scores are not independent of one another. The hysteresis measure contributes to all three of them and several other features contribute to two. Finding abnormality in all three components of the motor command may therefore reflect a single underlying signal expressed three times over. The correlations between the composites are not reported, and without them the reader cannot tell which of these two readings is correct.

Second, the comparison with stroke and with spinal cord injury, which carries much of the novelty of the paper, is asserted rather than demonstrated. A good deal of recent motor unit work in spinal cord injury and in stroke is also omitted, which is an important weakness, and I would encourage the authors to engage with it directly rather than treat those populations as a settled contrast. The abstract and the discussion state that the variability seen here is fundamentally different from the consistency seen in those populations, but no stroke or spinal cord injury data are presented, and no quantitative comparison with published values is offered. This matters because the individual patterns illustrated in the paper, that is, reduced peak discharge rate, compressed rate modulation, a narrowed recruitment range, synchronisation between units, and continued firing after the end of the task, are all well-described features of spastic paresis of other causes. This is commonly observed in people with spinal cord injury as well. One of the three illustrated patients is described as having mild hemiparesis and spasticity. Most of this cohort also carries cervical and thoracic cord lesions in addition to lesions above the spinal cord, so a simpler reading of the spread in these data is a mixture of the mechanisms already known from stroke and from spinal cord injury, present in varying proportions in different patients. That would be a different and considerably less novel conclusion, and it deserves to be considered explicitly. There is a statistical asymmetry in the comparison as well, since the populations described as consistent have been studied in samples of ten to twenty people, and small samples cannot reveal the tails of a distribution.

Third, the central finding of greater variability between patients has plausible alternative explanations that have not been excluded. The value for each participant is the median across the motor units identified in that person, and the number of units varies widely between participants and between the two groups. The precision of a median depends on how many units contribute to it, so precision that differs systematically between groups will inflate the spread of the participant values on its own. In addition, the leg studied was the more affected leg in the patients but the dominant leg in the controls. Choosing the more affected of two limbs is a selection on an extreme value, and it will both shift the patient average and widen the patient distribution for reasons that are purely statistical. Since the Methods state that both legs were recorded, perhaps the authors could consider adding this as an additional control. Finally, the hypothesis is framed in terms of multiple peaks and possible subgroups of patients, yet no test of multimodality is performed, and the density estimates shown do not obviously support it.

The feature carrying by far the largest weight in the neuromodulation score shows no group difference at all, and depending on which version of the hysteresis measure entered the composite, either one or none of its six constituent features differs between groups. The group difference in neuromodulation is also absent in the primary model and emerges only when maximal strength is added as a covariate, a change that the text reports but describes as leaving the results generally unchanged.

Moreover, medication may be assessed more clearly considering the large cohort. Thirty-four of the 89 patients take antispastic drugs, and both tizanidine and baclofen act on the system that the neuromodulation score measures, while others take reuptake inhibitors that act in the opposite direction. A quantitative comparison, with the obvious caveat of confounding by indication, would be needed before the neuromodulation findings can be read as disease-related. The authors should consider this to strengthen the manuscript or justify why this was not done.

Appraisal:

The authors achieve their descriptive aim. They show convincingly that motor unit discharge is altered in multiple sclerosis, that the average change is towards lower discharge rates and reduced rate modulation, and that the patient group is more dispersed than the control group on most measures. I do not think they establish the two claims that give the work its stated significance, that is, that the abnormalities can be attributed specifically to excitatory, inhibitory, and neuromodulatory drive, and that the resulting pattern distinguishes multiple sclerosis from other conditions affecting the same pathways.

Impact and utility:

The findings and dataset are very novel and, once shared, will be a resource for the field, and I would encourage the authors to release the analysis code and the exact weights alongside it. Moreover, it would have been valuable to see more relationships reported between the clinical measures and motor unit behaviour, in particular disability, walking speed, lesion location and medication.

Author response:

We sincerely thank the reviewers and editors for their thorough and constructive comments. Our goal was to study characteristics of excitatory, inhibitory, and neuromodulatory components of voluntary motor commands in participants with multiple sclerosis. To do so, we drew on a novel reverse-engineering framework designed to estimate characteristics of these components from a selection of features extracted from motoneuron discharge patterns (Chardon et al., 2024).

We appreciate the reviewers’ recognition of the considerable size and richness of our dataset of motor unit firing patterns recorded from the tibialis anterior and soleus of 89 participants with MS and 34 neurologically intact controls, the rigor of our statistical analyses, the novelty and importance of our descriptive findings, and the potential value of our dataset and analytical approach to other researchers.

We address reviewer comments and discuss our planned revisions to the paper below.

Concerns about physiological inference of findings from the reverse engineering framework. Reviewer 1, for example:

The principal limitation concerns the physiological interpretation assigned to these discharge patterns. The excitation, inhibition, and neuromodulation variables are not directly measured physiological inputs. They are composite variables derived from several features of motor unit discharge, with each feature weighted according to relationships identified in simulations reported previously by Chardon and colleagues. Thus, the step from observed discharge behavior to specific underlying synaptic mechanisms is necessarily model-dependent. The distinction between these two levels of inference is important for interpreting the main conclusions of the study.

The issue is especially relevant because several different physiological processes can plausibly influence the same discharge features. … The interpretation of inhibition is a particularly clear example of this general inverse problem.

We acknowledge all three reviewers’ concerns that, as with any computational model, assumptions of the Chardon reverse engineering framework place limits on the physiological interpretation of our findings.

However, two key strengths of the reverse engineering framework warrant further emphasis when considering the extent to which model assumptions limit physiological inference: (1) the extensive experimental foundation and physiological realism of the motoneuron model, and (2) the demonstrated ability to address the problem of non-uniqueness when predicting characteristics of synaptic and neuromodulatory inputs from discharge patterns of model motoneurons. We will discuss these strengths and clarify the physiological properties that are modeled.

(1) Physiological foundation of the motoneuron model

We first highlight the extensive experimental foundation underlying the motoneuron model. The Heckman laboratory has devoted considerable effort to the development of a physiologically realistic motoneuron model. The motoneuron model is grounded in over 20 years of in situ voltage-clamp experiments that characterized motoneuron intrinsic properties and their responses to excitatory, inhibitory, and neuromodulatory inputs (e.g., Lee & Heckman, 1998a, 1998b, 1999; Kuo et al., 2003; Hyngstrom et al., 2008). The resulting model motoneurons used for simulations were described by the authors as “designed to closely recreate behaviors documented in our extensive database of current and voltage clamp studies in motoneurons within animal preparations” (Chardon et al., 2024).

Moreover, the motoneuron model did not originate with the reverse-engineering framework used in the present study. Earlier versions were refined and evaluated across multiple prior publications (Kim et al., 2009; Kim & Jones, 2012; Powers et al., 2012; Powers & Heckman, 2015, 2017; Beauchamp et al., 2023) and optimized to balance realism with computational speed.

While the motoneuron model’s extensive experimental foundation does not eliminate the inherent limitations of model-based inference, it constrains the model substantially and reduces the extent to which model predictions depend on arbitrary or purely theoretical assumptions.

(2) Addressing non-uniqueness in the inverse problem

We next discuss whether characteristics of motoneuronal inputs can be inferred with reasonable accuracy from motoneuron discharge patterns. This is the inverse problem for which non-uniqueness is a fundamental challenge.

Non-uniqueness must be considered when making physiological inferences from motor unit discharge. Many motoneuron discharge features, when examined on their own, are sensitive to changes in more than one type of motoneuronal input. For example, discharge hysteresis is affected by both the level of neuromodulation and the pattern of inhibition relative to excitation. Reviewers were concerned that the reverse engineering framework also suffered from non-unique solutions when predicting motoneuronal inputs from discharge patterns, suggesting as an example that changes to motoneuron firing patterns resulting from a change in the pattern of inhibition could also result from changes to the level of neuromodulation.

However, we must distinguish the general non-uniqueness problem from the performance of the reverse engineering framework specifically. Non-uniqueness is not an unaddressed weakness of the reverse-engineering framework. Rather, it is the central problem that the framework was explicitly designed to address.

The ability of the reverse-engineering framework to address and substantially reduce non-uniqueness was explicitly evaluated in Chardon et al., (2024) using simulated, known inputs to a population of model motoneurons. Characteristics of motoneuronal inputs (i.e., the distribution of excitatory inputs across motoneurons of different size, the pattern of inhibition relative to excitation, and the level of neuromodulation) were systematically varied in all combinations, generating a “library” of simulated spike trains from each model motoneuron. Seven reverse engineering features were calculated on the simulated discharge patterns and used to train linear and non-linear machine learning models. When applied to test data, the machine learning models predicted the known inputs with substantial accuracy. 

The key finding of Chardon et al. is that the discriminatory power of the reverse engineering framework comes from the collection of discharge features used in the inverse problem. Because individual features can be sensitive to more than one type of motoneuronal input, considering the features as an ensemble drastically reduces the inverse problem’s solution space. This allows the known inputs to the simulated model motoneurons to be predicted with substantial accuracy.

(3) Modeled and unmodeled parameters

Having discussed the unique strengths of the reverse engineering framework, we also acknowledge its limitations, as do the authors of the approach. There are aspects of motoneuronal inputs and intrinsic properties that were held constant or not explicitly represented in the Chardon et al. (2024) simulations. Therefore, application of the current framework to human motor unit data introduces uncertainty regarding the extent to which variation in unmodeled parameters may influence discharge patterns and, by extension, the inferred characteristics of excitatory, inhibitory, and neuromodulatory inputs.

However, we consider this uncertainty in the context of the physiological information gained by applying the framework to human data. Specifically, the framework allows us to interpret motor unit discharge patterns recorded in vivo using relationships established through systematic manipulation of physiologically grounded model parameters.

Thus, despite its limitations, the framework provides substantially greater physiological insight than can be obtained by interpreting individual motor unit discharge features in isolation and represents an important step toward identifying the combinations of motoneuronal inputs that contribute to typical and pathological motor unit discharge in humans. Further variation in these parameters could be incorporated in future iterations of the framework, as is possible with more complex models like that of Mousa and Elbasiouny (2026).

With regard to the inclusion of specific parameters, we discuss two points:

(A) It is necessary to clarify that some physiological properties identified in the reviewer comments as unmodeled were, in fact, represented in the motoneuron model and simulations used in Chardon et al. (2024).

Perhaps the most crucial is the presence of persistent inward currents (PICs), whose realistic representation is a defining feature of the motoneuron model due to its four dendritic compartments with varied Ca²⁺ channel densities (Powers & Heckman, 2017). In fact, in Chardon et al. (2024), the PIC-induced non-linearities in simulated firing patterns were not merely present in the model, but were a key factor contributing to the success of the reverse engineering framework.

Contributions from afferent inputs were also considered within the synaptic inputs to the model. The simulated distributions of excitatory input and patterns of inhibitory input both encompass potential contributions from afferent sources in addition to descending and spinal sources (Binder et al., 2002; Johnson et al., 2017; Chardon et al., 2024). Further, the opposing influence of inhibitory input from any source on facilitation of PICs is incorporated in the simulations.

(B) We agree with the reviewers that it is important to consider potential between-group differences in intrinsic membrane properties that affect PIC behavior in addition to neuromodulatory input. This is a point that we did not adequately discuss in the initial version of the paper. However, we see this primarily as a lack of precision in our discussion of the physiology rather than a deficit related to unmodeled parameters in the Chardon et al. model.

In the Chardon et al. (2024) simulations, the level of neuromodulatory input was varied with the intent of scaling PIC amplitude accordingly. This was done by changing the density of dendritic PIC channels to simulate how the number of channels activated scales with the amount of neuromodulatory input (i.e., serotonin and norepinephrine) present. Thus, the reverse-engineering features that predicted the simulated level of neuromodulatory input in the Chardon study more directly predicted the resulting variation in simulated PIC amplitude. Because other intrinsic membrane properties (e.g., NaV conductance) were held constant, changes in PIC amplitude could be attributed specifically to the simulated level of neuromodulatory input.

In contrast, variation in PIC amplitude estimated from our experimental data cannot be attributed to differences in neuromodulatory input alone. We must also consider other factors that affect PIC amplitude that are unknown in our participants, including differences in intrinsic membrane properties. Importantly, several discharge features incorporated into the reverse engineering framework, including braceheight and delta-F, are commonly used to estimate PIC amplitude in human motor unit recordings, independently of the Chardon et al. (2024) framework (e.g., Mesquita et al., 2024).

Thus, to the extent that brace height, delta-F, and the “neuromodulation” composite variable reflect PIC amplitude in our human data, our findings reflect variation in the physiological factors that determine PIC amplitude. These include neuromodulatory input and intrinsic membrane properties, as well as, for delta-F in isolation, the pattern of inhibition. In our revised version of the manuscript, we will incorporate discussion of these additional factors in addition to our current discussion of neuromodulatory input and the evidence for its alteration in MS.

Finally, while identifying the specific motoneuron inputs and properties that contribute to altered PIC amplitude among participants with MS is an important goal for future work, uncertainty regarding those contributors does not preclude the potential functional significance of our finding that estimated PIC amplitude can be abnormally high or low among participants with MS. Because PIC amplitude directly influences motoneuron excitability, abnormally high or low PIC amplitudes could contribute substantially to motor deficits that emerge in MS and their variation across patients.

Calculation of composite variables using mutual information scores. Reviewer 3, for example:

Mutual information measures how much a feature tells you about a parameter. It does not carry the direction of that association, and it does not become a regression coefficient by having a sign attached to it. Since the direction of every conclusion in the paper depends on those associated signs, the reader has no way of judging how faithfully the composite scores track the physiological quantities they are named after. I do not understand why this was not done with simulations. … The simulations have known inputs, so the composites could be computed on the simulated discharge patterns and their accuracy reported directly.

A related difficulty is that the three composite scores are not independent of one another. The hysteresis measure contributes to all three of them and several other features contribute to two. Finding abnormality in all three components of the motor command may therefore reflect a single underlying signal expressed three times over. The correlations between the composites are not reported, and without them the reader cannot tell which of these two readings is correct.

We agree that mutual information scores are non-directional and that they are not regression coefficients. We did not use them as regression coefficients; rather, we used them as an informed way to create a weighted average of the reverse engineering features that, individually, were most informative about each type of input. Signs were assigned to each feature to reflect the direction of the relationship between the feature and each type of input. We took great care to base the assigned signs on published literature, and we will update the paper to more explicitly link each assigned sign to the supporting literature.

Thank you for the useful suggestion to calculate our composite scores from the simulated features from Chardon et al. and compare them with the known input values. We obtained these data from the authors and found that two of the signs in the excitation composite (for torque at recruitment and duration) needed to be updated based on the Chardon data. After correcting these signs, the neuromodulation and excitation composite variables were both highly correlated with their respective known inputs (r > 0.85). Correction of the two signs in the excitation composite did not change the general pattern of our results. The inhibition composite variable was moderately correlated with the known input (r = 0.50), consistent with the original Chardon et al. study showing that machine learning predictions of inhibition using non-linear regression greatly outperformed those using linear regression. We will include these validation results in the revised manuscript, as they provide a direct assessment of how well the composite variables reflect the known inputs.

We agree that it is important to demonstrate the independence of the composite variables, and it was an oversight on our part not to include that information. The correlations among the composite variables were very low, ranging from r = 0.006 to r = 0.18, with none reaching statistical significance. As discussed, many of the discharge features are sensitive to more than one type of input, which is why it is difficult to interpret them physiologically in isolation. However, the different types of input affect the features in different ways. For example, push-pull/reciprocal inhibition increases both delta-F and rate attenuation slope compared with uniform inhibition. In contrast, an increase in neuromodulation increases delta-F but decreases rate attenuation slope. The signs within the composite variable calculations reflect these differences in the relationships between each feature and the inputs. Therefore, the shared constituent features do not necessarily result in strongly correlated composite variables, as confirmed by the very low correlations observed in our data.

We are now collaborating with the authors of Chardon et al. (2024) to explore application of their reverse engineering framework (and/or its subsequent updates currently under development) to our MS and control data to supplement our mutual information-based composite variable analyses.

Comparison with stroke and spinal cord injury populations, Reviewer 3.

Second, the comparison with stroke and with spinal cord injury, which carries much of the novelty of the paper, is asserted rather than demonstrated. A good deal of recent motor unit work in spinal cord injury and in stroke is also omitted, which is an important weakness, and I would encourage the authors to engage with it directly rather than treat those populations as a settled contrast. The abstract and the discussion state that the variability seen here is fundamentally different from the consistency seen in those populations, but no stroke or spinal cord injury data are presented, and no quantitative comparison with published values is offered. This matters because the individual patterns illustrated in the paper, that is, reduced peak discharge rate, compressed rate modulation, a narrowed recruitment range, synchronisation between units, and continued firing after the end of the task, are all well-described features of spastic paresis of other causes.

Although our discussion of the literature closely reflects the quantitative results from those populations, we agree that directly presenting those values would better support the comparison. We therefore will include quantitative comparisons with published values in these populations, to the extent that they are available, in the revised paper.

We are uncertain which additional studies the reviewer has in mind when stating that a “good deal of recent motor unit work in stroke and spinal cord injury is omitted, which is an important weakness.” We cited the published papers most directly relevant to our specific line of inquiry. Nonetheless, we will conduct an additional review of recent literature and update our references and discussion to address any relevant omissions in the paper revision.

Finally, we agree that the motor unit firing characteristics mentioned above can be found in spinal cord injury and stroke populations, although their prevalence and expression differ between the populations. However, the presence of those characteristics in some of our participants with MS does not undermine the primary distinction we intended to make between our findings and those reported in stroke and spinal cord injury. One novel aspect of our findings is the marked between-participant heterogeneity in motor unit firing patterns within the MS group, as shown in Figure 3. In particular, among participants with MS, alterations in multiple discharge features were not just more variable than controls, but in some cases, they deviated from controls in opposite directions. We will revise the manuscript to make clear that it is this heterogeneity, rather than the presence of any individual discharge characteristic, that distinguishes the patterns observed in our MS sample from the predominantly group-mean difference patterns typically reported in stroke and spinal cord injury. We will support this comparison quantitatively where published data permit.

Further, the motor unit characteristics listed by the reviewer were not our only finding. Another novel aspect of our findings is that a substantial number of participants with MS demonstrated decreases in estimated PIC amplitude. In contrast, recent work in spinal cord injury and stroke reported increases in estimated PIC amplitude during voluntary drive in both groups (Hassan, 2021; Benedetto et al., 2026).

Alternative explanations to increased variability in MS. Reviewers 3, for example:

… the central finding of greater variability between patients has plausible alternative explanations that have not been excluded.

In the revised manuscript, we will address the comments from Reviewers 1 and 3 regarding potential alternative contributors to the increased between-participant variability observed in MS, including motor unit yield, decomposition quality, task performance, anti-spastic medications, and other factors raised by reviewers.

Group differences in the constituent features of the composite variables. Reviewer 3:

The feature carrying by far the largest weight in the neuromodulation score shows no group difference at all, and depending on which version of the hysteresis measure entered the composite, either one or none of its six constituent features differs between groups.

There are two main points to clarify. First, a statistically significant group-mean difference in a composite variable does not require the individual constituent features to also have statistically significant group mean differences. As discussed, many of the features are sensitive to more than one type of motoneuronal input. Further, Chardon et al. (2024) demonstrated that prediction of motoneuronal inputs improved as more features were included.

Second, the absence of a group-mean difference in an individual discharge feature or composite variable does not imply that the feature or composite variable is unchanged among participants with MS. This is particularly important given that many of our MS distributions included values that deviated in opposite directions from controls, which can result in little or no difference in the group mean despite substantial differences at the individual level. For this reason, our analyses compared MS and control distributions not only in terms of central tendency but also their spread and shape.

Relationships between motor unit discharge features, composite variables, and clinical measures. Reviewer 3, for example:

Moreover, it would have been valuable to see more relationships reported between the clinical measures and motor unit behaviour, in particular disability, walking speed, lesion location and medication.

We agree that exploring these relationships is a logical next step and will be especially valuable for understanding the heterogeneity of discharge features and composite variables among participants with MS. The primary goal of the present study was to characterize this heterogeneity, whereas determining how specific clinical characteristics relate to that heterogeneity represents a substantial additional question. We are currently preparing a follow-up paper focused specifically on these relationships so that we have sufficient space to present and discuss them thoroughly. As mentioned above, however, in the paper revision, we will include more information about data from participants who were and were not taking anti-spastic medications.

In our revision, we will also use more cautious terminology when discussing physiological inference vs. motor unit phenotypes, where appropriate, revise the denominator used in calculating the weighted averages, add a formal test of multimodality, add more raw motor unit firing traces, and address the reviewers’ other outstanding minor suggestions.

Beauchamp JA, Pearcey GEP, Khurram OU, Chardon M, Wang YC, Powers RK, Dewald JPA & Heckman C (2023). A geometric approach to quantifying the neuromodulatory effects of persistent inward currents on individual motor unit discharge patterns. J Neural Eng 20, 016034.

Benedetto A, Jenz S, Farley M, Heit B, Sangari S, Beauchamp JA, McPherson L, Heckman C, Perez M & Pearcey G (2026). Muscle-specific motor unit firing characteristics in elbow flexors and extensors after cervical spinal cord injury. ; DOI: 10.64898/2026.06.03.729825. Available at: http://biorxiv.org/lookup/doi/10.64898/2026.06.03.729825 [Accessed June 9, 2026].

Binder MD, Heckman CJ & Powers RK (2002). Relative Strengths and Distributions of Different Sources of Synaptic Input to the Motoneurone Pool. In Sensorimotor Control of Movement and Posture, ed. Gandevia SC, Proske U & Stuart DG, pp. 207–212. Springer US, Boston, MA. Available at: https://doi.org/10.1007/978-1-4615-0713-0_25.

Chardon MK, Wang YC, Garcia M, Besler E, Beauchamp JA, D’Mello M, Powers RK & Heckman CJ (2024). Supercomputer framework for reverse engineering firing patterns of neuron populations to identify their synaptic inputs. eLife. Available at: https://elifesciences.org/articles/90624 [Accessed November 2, 2024].

Hassan AS (2021). Changes in Motor Unit Firing Patterns as a Function of Age, Muscle, and Following a Unilateral Brain Injury: Ionotropic and Metabotropic Effects (Ph.D. thesis). Northwestern University, United States -- Illinois. Available at: http://libproxy.wustl.edu/login?url=https://www.proquest.com/dissertations-theses/changes-motor-unit-firing-patterns-as-function/docview/2572564921/se-2?accountid=15159.

Hyngstrom AS, Johnson MD & Heckman CJ (2008). Summation of Excitatory and Inhibitory Synaptic Inputs by Motoneurons With Highly Active Dendrites. J Neurophysiol 99, 1643–1652.

Johnson MD, Thompson CK, Tysseling VM, Powers RK & Heckman CJ (2017). The potential for understanding the synaptic organization of human motor commands via the firing patterns of motoneurons. J Neurophysiol 118, 520–531.

Kim H & Jones KE (2012). The retrograde frequency response of passive dendritic trees constrains the nonlinear firing behaviour of a reduced neuron model. PloS One 7, e43654.

Kim H, Major LA & Jones KE (2009). Derivation of cable parameters for a reduced model that retains asymmetric voltage attenuation of reconstructed spinal motor neuron dendrites. J Comput Neurosci 27, 321–336.

Kuo JJ, Lee RH, Johnson MD, Heckman HM & Heckman CJ (2003). Active Dendritic Integration of Inhibitory Synaptic Inputs In Vivo. J Neurophysiol 90, 3617–3624.

Lee RH & Heckman CJ (1998a). Bistability in Spinal Motoneurons In Vivo: Systematic Variations in Persistent Inward Currents. J Neurophysiol 80, 583–593.

Lee RH & Heckman CJ (1998b). Bistability in Spinal Motoneurons In Vivo: Systematic Variations in Rhythmic Firing Patterns. J Neurophysiol 80, 572–582.

Lee RH & Heckman CJ (1999). Enhancement of Bistability in Spinal Motoneurons In Vivo by the Noradrenergic α₁ Agonist Methoxamine. J Neurophysiol 81, 2164–2174.

Mesquita RNO, Taylor JL, Heckman CJ, Trajano GS & Blazevich AJ (2024). Persistent inward currents in human motoneurons: emerging evidence and future directions. J Neurophysiol 132, 1278–1301.

Mousa MH & Elbasiouny SM (2026). A Multi-Scale, High-Fidelity Computational Model of the Mouse Triceps Surae Motor Pool.

Powers RK, ElBasiouny SM, Rymer WZ & Heckman CJ (2012). Contribution of intrinsic properties and synaptic inputs to motoneuron discharge patterns: a simulation study. J Neurophysiol 107, 808–823.

Powers RK & Heckman CJ (2015). Contribution of intrinsic motoneuron properties to discharge hysteresis and its estimation based on paired motor unit recordings: a simulation study. J Neurophysiol 114, 184–198.

Powers RK & Heckman CJ (2017). Synaptic control of the shape of the motoneuron pool input-output function. J Neurophysiol 117, 1171–1184.