Experimental setup.

Participants were seated with the tested foot secured to the dynamometer. High-density surface EMG arrays were affixed to the TA and soleus muscles. Real-time visual feedback of torque generation during experimental trials (inset) was displayed on a TV monitor.

Calculation of reverse engineering features and interpretation of composite variables.

A. Calculation of reverse engineering features from the smoothing firing rate vs. time function: torque at recruitment and de-recruitment, duration, recruitment range, delta-F, maximal firing rate, firing rate at recruitment and de-recruitment, and firing rate range. The delta-F calculation shown is for the higher threshold blue motor unit, whereas the other features are shown for the lower threshold purple motor unit. B. Calculation of reverse engineering features from the smoothed firing rate vs. torque function: braceheight and rate attenuation slope. C. Neuromodulatory input facilitates persistent inward currents PICs) that drastically affect motor unit firing patterns. Top: during a triangular contraction, motor units with large PICs produce firing patterns with a characteristic non-linear shape. Three distinct phases can be seen: rapid initial acceleration, rate attenuation marked by a more gradual rate increase, and recruitment hysteresis where the motor unit de-recruits at a lower level of excitatory synaptic input than when it was recruited. The example shown is a motor unit from the TA of a control participant. Bottom: Motor units with small (or no) PICs produce more linear firing patterns. D. Physiologic interpretation of excitatory, inhibitory, and neuromodulatory composite variables.

Weights for each reverse engineering feature in calculating excitation, inhibition, and neuromodulation composite variables.

Motor unit firing patterns from single trials.

Smoothed motor unit discharge rates from one control participant (left) and three MS participants. Each colored line shows the firing pattern of a different motor unit, and torque over time is shown in black. Data from the TA during dorsiflexion are shown on the top row, and data from the soleus during plantarflexion are shown on the bottom row.

Group distribution shapes for the excitation, inhibition, neuromodulation, and firing rate composite variables.

Probability density functions are shown for the control (blue) and MS (pink) groups for the TA (top) and soleus (bottom). The individual participant median values that make up the density functions are shown at the bottom of each plot, colored by group. Estimated marginal means (i.e., adjusted for covariates) are reported in the text.

Summary of reverse engineering feature data.

Median values for control participants are shown as blue circles and those for MS participants are shown as pink circles. Group mean values for each Group and Muscle combination are indicated by horizontal black lines. For parameters with a significant effect of Group or Group x Muscle interaction, the associated p-value is displayed on the plot.

The mean ± SD (total) number of motor unit spike trains decomposed per trial (top) and the mean ± SD (total) number of unique motor units per participant (bottom) for each group, muscle, and sex combination.

Maximum voluntary torque (MVT) in dorsiflexion and plantarflexion for females (left) and males (right).

Each circle represents the MVT value for one participant, and horizontal black lines indicate the group mean. Circle size corresponds with Age, and circle color in the MS group corresponds to the participant EDSS score. We compared maximal strength between groups for dorsiflexion and plantarflexion separately using linear regression with a fixed effect of Group and covariates of Age_c (grand mean-centered Age) and Sex. Maximal strength was significantly lower for the MS group compared with the control group for both dorsiflexion (30.1 Nm vs. 36.7 Nm, p = 0.0008) and plantarflexion (86.2 Nm vs. 108.8 Nm, p = 0.0004). For the five participants whose TA data could not be included, their maximal strength values, reasons for exclusion, and sex were as follows: 23.9 Nm (insufficient yield, F), 38.5 Nm (insufficient yield, F), 0.75 Nm (weakness, M), 0.00 Nm (weakness, F), 0.78 Nm (weakness, F).

Schematic of the mixed-effect model structure for neurophysiological variables collected at the motor unit level.

For any given outcome (e.g., delta F), each subject has multiple observations per motor unit and trial. Thus, there are three sources of variance we generally want to account for in the random effects (participant, trial, and motor unit). However, motor unit labels are arbitrary: i.e., motor unit “Φ” is same for one participant from trial to trial, but Φ in the soleus is not the same as Φ in the TA, and Φ for Participant 1 is not the same as Φ for Participant 2. Thus, we concatenate participant, muscle, and motor unit labels (e.g., “1TAΦ”), to obtain a random intercept that appropriately nests motor units within a given participant and muscle. Similar to the previous model, fixed effects then allow us to test the study’s main hypotheses with respect to Group and Group x Muscle. Continuous fixed effects were mean centered and categorical fixed effects were contrast coded using orthogonal polynomials. Notably, torque at recruitment was grand-mean centered (TQ_recrt_c), as opposed to being centered within a participant, and serves a motor unit-level covariate (except for the model with TQ_recrt as the outcome measure).

Linear mixed effect model results for the composite variables and reverse engineering features: Effects of Group, Muscle, Group x Muscle

Linear mixed effect model results for the composite variables and reverse engineering features: Effects of Age, Sex, and Torque at Recruitment

Linear mixed effect model results for the composite variables when controlling for maximal strength

Tests for equality of variance and differences in distribution shape for the composite variables.

Tests for equality of variance for the reverse engineering features.