Introduction

Array electrodes that allow simultaneous recording of firing patterns of populations of neurons in mammals have transformed understanding of the computations implemented by neural net-works during natural behaviors and the utilization of these computations for brain machine interfaces (Buonomano and Maass, 2009; Chandrasekaran et al., 2021; Stevenson et al., 2012; Stevenson and Kording, 2011; Vyas et al., 2020). However, to understand how a given neural network generates the firing patterns that implement these computations, it is necessary to understand the transformation of inputs to outputs by the constituent neurons. The net effect of ionotropic excitatory inputs in activating a neural circuit is likely to be strongly impacted by the pattern of inhibitory inputs arising from local interneurons and the effects of neuromodulatory inputs acting on G-protein coupled receptors (Binder et al., 2020; Goaillard and Marder, 2021; McCormick and Nusbaum, 2014). In-depth analyses of the resulting complex interactions require intracellular electrode techniques, but these are not yet feasible simultaneously in many neurons during a natural behavior. An alternative is reverse engineering (RE) the firing patterns recorded by arrays to identify the neurons’ inputs and properties.

Thus far, RE has been used to achieve several different goals, including estimation of circuit structures within simulations with simple neuron models (Lim et al., 2011; Pisokas, 2021; Rostro-Gonzalez et al., 2012) and estimation local field potentials from measurements of spiking patterns (Telenczuk et al. (2020); see also Yochum et al. (2019)). In this study, our goal was to harness the computational power provided by implementation of realistic neuron models on a supercomputer to investigate the feasibility of reverse engineering these model firing patterns to identify the underlying organization of their simulated excitatory, inhibitory and neuromodulatory inputs. Successful RE of firing patterns into input organization would greatly advance the effectiveness of array recordings in identifying how neural circuits implement their computations.

RE of neural outputs into inputs faces a formidable barrier. Multiple neural systems have been shown to exhibit non-uniqueness, which is manifest as closely similar outputs being produced by many combinations of parameters specifying neuron properties and inputs (Edelman and Gally, 2001; Marder, 2012; Prinz et al., 2004). This set of parameters forms the “solution space” for that neural system (Prinz, 2010). Although likely advantageous for stability and resilience of neurons and circuits (Goaillard and Marder, 2021), a large solution space that contain many combinations of excitation, inhibition and neuromodulation may limit the effectiveness of RE. In addition, the parameters within the solution space may interact in complex, nonlinear ways (Mukunda and Narayanan, 2017). Neuromodulatory inputs, for example, may have highly nonlinear effects on how neurons process their ionotropic inputs (Binder et al., 2020; Marder, 2012; McCormick and Nusbaum, 2014). On the other hand, when neuronal output behaviors are complex, non-uniqueness may be reduced by the constraints of matching multiple outputs (Mukunda and Narayanan, 2017; Yang et al., 2022).

To reduce the impact of non-uniqueness, we chose a well-understood neural system for evaluation of RE, the motoneurons that transform all motor commands into signals for control of muscles (Heckman and Enoka, 2012). Motoneuron properties and inputs have been extensively studied using in situ voltage clamp methods in animal preparations, with a focus on neuromodulatory actions (Heckman and Enoka, 2012; Henneman and Mendell, 1981). This data set has enabled us to generate highly realistic motoneuron models (Powers et al., 2012; Powers and Heckman, 2017), which we have successfully implemented on a supercomputer for the present study. Furthermore, recent advances in flexible array electrode methods now allow measurement of firing patterns of populations of motoneurons in humans, due to the one-to-one relation between motoneuron action potentials and those of its innervated muscle fibers (Holobar et al., 2010; Hug et al., 2021; Negro et al., 2016). Consequently, RE of these population firing patterns may generate deep insights into the synaptic organization of motor commands in humans as well as other animals.

We used an ensemble modeling approach (Prinz, 2010) to determine how non-uniqueness and nonlinearity affected RE identification of the pattern of excitatory, inhibitory and neuromodulatory inputs that produces motor output. All studies were carried out in silico, so that both inputs and outputs were known and thus the effectiveness of RE could be accurately quantified. The su-percomputer implementation of motoneuron models allowed us to thoroughly explore this input solution space, via more than 300,000 simulations of the effects of many input combinations on firing patterns. Our previous intracellular studies of motoneuron input-output functions suggest that motoneuron firing patterns contain substantial information about their inputs (Hyngstrom et al., 2008; Kuo et al., 2003; Lee and Heckman, 1999, 1998a,b). We therefore hypothesized that RE of motoneuron firing patterns would effectively deal with the problems of non-uniqueness and would therefore identify a reasonably small solution space of input combinations. To assess the importance of information contained within neuron firing patterns for achieving a small solution space, we compared RE of motoneuron firing patterns to RE of the cumulative spike train (CST), an overall measure of motor output generated by summing all neuron firing patterns together, which has been shown to closely replicate the muscle electromyogram (EMG), which is directly related to muscle force. Our simulation results supported the hypothesis and demonstrated a much smaller solution space for RE of firing of neuron populations than for the CST.

Results

Performance of the motoneuron models for integration of inputs

All simulations are run using a set of 20 model motoneurons designed to closely recreate behaviors documented in our extensive database of current and voltage clamp studies in motoneurons within animal preparations. These models were first developed in NEURON (Hines and Carnevale, 1997) for a previous study (Powers and Heckman, 2017) and have been implemented on a supercomputer for the present work (see Methods). Our goal of investigating the feasibility of RE of motoneuron firing patterns depends on accuracy of these motoneuron models. Figure 1 illustrates the performance of our motoneuron models in realistically capturing the integration of excitatory (Figure 1A), neuromodulatory (Figure 1B) and inhibitory (Figure 1C) inputs. Each of these simulated behaviors are considered in the following subsections.

Figure 1.

Common input structure, differences in motoneuron properties

A fundamental point in achieving accurate modeling of the outputs of a pool of motoneurons innervating a muscle is that they are subject to “common drive” (De Luca and Erim, 1994; Farina and Negro, 2015; Heckman and Enoka, 2012). As a result, motoneuron firing patterns for a specific task exhibit an overall similarity, which is likely essential for efficient force control (Farina and Negro, 2015). This point is important for RE, as common drive precludes independent inputs to each motoneuron. Firing patterns, however, are not identical. As in any biological system, noise is present, though this is relatively low in motoneurons (coefficient of variation during steady firing at about 10-20% (Heckman and Enoka, 2012)). We add noise to our excitatory synaptic conductances to match this level, as illustrated in Figure 1A1, which illustrates firing patterns of 5 of the 20 motoneuron models in response to a linearly rising and falling excitatory conductance. This figure also shows systematic differences in activation (i.e. recruitment) thresholds. These differences are due to systematic differences in intrinsic electrical thresholds for spiking of motoneurons, which generates an orderly activation (i.e. recruitment) sequence known as Henneman’s size principle (Henneman and Mendell, 1981). The size principle is essential for effective utilization of the slow vs fast (F) muscle fibers, with motoneurons innervating slow fibers having the lowest thresh-olds. S units tend to reach higher firing rates than F units, which require more input to reach their higher thresholds so less is available for suprathreshold rate modulation (Beauchamp et al., 2022; De Luca and Contessa, 2012; Hu et al., 2013).

Distribution of excitatory input to S vs F motoneurons

In addition to the pattern and amplitude of excitatory input, there is another important aspect of excitatory motor commands, which is that all excitatory inputs studied thus far project differentially to S vs. F motoneurons (Powers and Binder, 2001) (in contrast, inhibition tends to be uniform ((Powers and Binder, 2001; Lindsay and Binder, 1991) as does neuromodulation (Lee and Heckman, 1998a, 1999). Descending inputs tend to generate larger synaptic currents in the F’s, especially the corticospinal input (Powers and Binder, 2001) (see Discussion). As yet only one input favors S motoneurons, but it is an essential one, monosynaptic Ia excitation from muscles spindles (Heckman and Binder, 1988). These differences, though, do not provide independent control of F motoneu-rons, much less of individual motoneurons, as common drive still dominates and recruitment still follows the size principle (Heckman and Binder, 1993). Figure 1A2 however shows that greater relative excitatory synaptic current in high threshold F motoneurons does tend to compress the range of recruitment thresholds and increase the steepness of rate modulation as input levels rise. This compression/steepening may thus provide RE sufficient information to identify the distribution of excitatory input to S vs F motoneurons and provide information on the relative roles of descending versus sensory input.

Neuromodulation and excitatory input

We chose to focus on just one type of neuromodulation (see Methods), the monoaminergic input originating in the brainstem, because this system has proved to have particularly powerful effects on motoneuron excitability (Heckman et al., 2005; Heckman and Enoka, 2012). These monoamin-ergic axons form a dense and monosynaptic projection onto motoneurons throughout the spinal cord, with the synapses releasing either 5HT (axons originating in the caudal raphe nuclei) or NE (axons from the locus coeruleus) (Alvarez et al., 1998; Holstege and Kuypers, 1987; Maratta et al., 2015). Both 5HT and NE have multiple effects on motoneuron ion channels, but their strongest action is via facilitation of persistent inward currents (PICs) mediated by dendritic Ca channels (likely CaV 1.3) (Binder et al., 2020). Our motoneuron models have been fine tuned to accurately replicate the highly nonlinear interactions between PICs and excitatory inputs (Powers and Heckman, 2017). Figure 1B shows that, when brainstem neuromodulatory are moderately high and dendritic PICs become strong, three marked nonlinearities become apparent: an initial acceleration (due to PIC activation and consequent amplification of input), followed by attenuation (due to PIC depolarization of dendritic regions, which reduces excitatory driving force) and then hysteresis (offset at a lower level than onset, due to the prolongation of input). In this example, where a pure excitatory drive is applied with no background of inhibition, low threshold motoneurons continue firing long after the input returns to baseline (this is known as self-sustained firing). These nonlinearities are striking in comparison to the nearly linear behavior in the low neuromodulatory state (Figure 1A1).

Neuromodulation and inhibition

Our models also accurately capture the very strong interaction of the PIC with inhibition. Even a small background of inhibition reduces PIC amplitudes (Hyngstrom et al., 2007; Kuo et al., 2003). Figure 1C shows that modulation of inhibitory input that is either reciprocally or proportionally organized with respect to excitation also dramatically alters the pattern of rate modulation. Reciprocal inhibition (Figure 1C1) tends to allow strong PIC expression with acceleration, attenuation and hysteresis, whereas proportional inhibition (Figure 1C2) tends to flatten rate modulation and reduce hysteresis (in fact, of all the patterns shown in Figure 1, C1 is closest to that seen in many human muscles (Beauchamp et al., 2022; Heckman and Enoka, 2012). Thus the nonlinearities in firing patterns are not a straightforward reflection of neuromodulatory input but instead result from a complex interaction between neuromodulation and inhibition.

Given the above effects of the organization of excitation, inhibition and neuromodulation on firing patterns, our RE effort targeted identification of the following aspects of excitation, inhibition and neuromodulation: amplitude and time course of excitation, distribution of excitation to S versus F motoneurons, amplitude and time course of inhibition (including the baseline level and then whether inhibition was proportional, constant or balanced) and level of neuromodulation from low to high.

Ensemble modeling for Reverse Engineering

The concept for the RE approach we used is illustrated in Figure 2. The left portion shows the flow of motor commands from the brain, brainstem and spinal cord to motoneurons, which then activates a muscle. The muscle produces an output force. This force is due to the summed actions of recruited motor units, whose spike trains sum to produce an overall electrical signal (EMG). Our focus is on the spike trains themselves, which can be measured in both animals and humans. The pattern of average firing rate is obtained by summing the rates of all active units (see Methods) and is thus referred to as the cumulative spike train (CST) (righthand side, bottom) We then apply RE methods to both the CST and to individual motor unit firing patterns (right, middle), which, if successful, will allow us to estimate the temporal patterns and amplitudes of the three components of motor commands. Two hypothetical outcomes are illustrated. In both, a triangular form of excitation is applied. The push-pull scheme is a reciprocal organization in which a tonic level of inhibition decreases as excitation increases. The balanced scheme has the opposite organization, inhibition increasing in proportion to excitation. Neuromodulation is assumed to be constant but can varied widely in its level. Although the balanced scheme often is used to drive cortical neurons (Berg et al., 2007), it is unknown how these three components interact in motor commands to motoneurons.

Figure 2.

The ensemble modeling approach we utilize for RE begins with estimation of the trajectory of the excitatory motor command. The procedure is illustrated in Figure 3. We required all simulations to achieve the same overall output, a triangular waveform for CST (upper left; Methods for the match criterion). For each simulation, we specified the following properties of the input organization:

Figure 3.

1. Distribution of excitatory input on to S vs F motoneurons (range of 0.5 to 2.5; at 0.5 the highest threshold type F motoneuron received half the average input of lowest threshold type S; at 2.5 the F received 2.5 times the input of the S).

2. Level of neuromodulation: range 0.8 to 1.2, specified in terms of the maximum conductance of the PIC.

3. Baseline level of inhibition: This was set to the lowest level that prevented the PIC from generating self-sustained firing (Figure 1B shows an example of motoneurons with strong self-sustained firing due to lack of sufficient inhibitory baseline)

4. Pattern of inhibition relative to excitation: This pattern was set to be linearly proportional to excitation via multipliers ranging from +0.7 to −0.7. Figure 1C1 and Figure 1C2 illustrate firing patterns resulting from maximum proportional inhibition (multiplier of +0.7) to maximal push pull (−0.7). At 0.0, the baseline of inhibition remained constant.

5. A physiological level of synaptic noise (see Methods). The RMS value for this noise was the same in all simulations, but was regenerated 30 times to allow estimates of variance for each combination of the above 4 input properties.

Then, for each simulation, we set the input and noise properties as above and then applied a triangular pattern of excitation to the motoneuron models, starting with a purely triangular form and modifying it in an iterative process until the CST of the motoneurons reproduced the linear target CST (Figure 3, flow diagram). We repeated these simulations with different sets of input parameters 315,000 times to form the ensemble model of the motor pool (see section Data preparation for parameter combination details). We used the resulting large data set for successful matches to the target CST to assess the feasibility of RE for both the CST and individual motoneuron firing patterns.

The ensemble modeling approach we utilize for RE begins with estimation of the trajectory of the excitatory motor command. The procedure is illustrated in Figure 3. We required all simulations to achieve the same overall output, a triangular waveform for CST (upper left; Methods for the match criterion). For each simulation, we specified the following properties of the input organization:

1. Distribution of excitatory input on to S vs F motoneurons (range of 0.5 to 2.5; at 0.5 the highest threshold type F motoneuron received half the average input of lowest threshold type S; at .5 the F received 2.5 times the input of the S).

2. Level of neuromodulation: range 0.8 to 1.2, specified in terms of the maximum conductance of the PIC.

3. Baseline level of inhibition: This was set to the lowest level that prevented the PIC from generating self-sustained firing (Figure 1B shows an example of motoneurons with strong selfsustained firing due to lack of sufficient inhibitory baseline)

4. Pattern of inhibition relative to excitation: This pattern was set to be linearly proportional to excitation via multipliers ranging from +0.7 to −0.7. Figure 1C1 and Figure 1C2 illustrate firing patterns resulting from maximum proportional inhibition (multiplier of +0.7) to maximal push pull (−0.7). At 0.0, the baseline of inhibition remained constant.

5. A physiological level of synaptic noise (see Methods). The RMS value for this noise was the same in all simulations, but was regenerated 30 times to allow estimates of variance for each combination of the above 4 input properties.

Then, for each simulation, we set the input and noise properties as above and then applied a triangular pattern of excitation to the motoneuron models, starting with a purely triangular form and modifying it in an iterative process until the CST of the motoneurons reproduced the linear target CST (Figure 3, flow diagram). We repeated these simulations with different sets of input parameters 315,000 times to form the ensemble model of the motor pool (see section Data preparation for parameter combination details). We used the resulting large data set for successful matches to the target CST to assess the feasibility of RE for both the CST and individual motoneuron firing patterns.

Solution space for overall output

We hypothesized that the CST, which contains information only about the overall output trajectory, would have a large solution space. It was nonetheless surprising that our process for modifying the shape of the excitatory input command generated good matches to the target CST for every motor command combination we tried (all 315,000 of them (see section Data preparation)). However, as illustrated in by the progression of input forms in Figure 2 and further elaborated in Figure 4, the potent amplification and large nonlinearities in firing patterns induced by PICs usually resulted in highly nonlinear excitatory commands being required to generate the linear CST. The excitatory conductance input that leads to a good match to the target firing rate typically shows a rapid increase to get recruitment started, followed by a drop off during the acceleration phase due to wide-spread PIC activation across the active motor units, followed by a more rapid increase to keep the average firing rate rising in the face of PIC rate saturation. The falling phase of the command was roughly the mirror image of the rising phase. This general trend is evident in all the examples in Figure 4.

Figure 4.

As expected, as the level of neuromodulation increased, these nonlinear trends increased, but less excitation is required to achieve the same global output. The examples in upper half of Figure 4 assumed inhibition was constant. The lower left example in Figure 4 shows examples of the effects of different organizations of inhibition. Much higher levels of excitation are needed to counteract inhibition when it varies in proportion to the excitatory command (balanced inhibition, red traces), then when it varies inversely (push-pull inhibition, green traces). To quantify the effects of neuromodulation and inhibition on the overall amplitude of excitation, we calculated the area under the curve for the excitation patterns as neuromodulation and inhibition varied. Figure 5 shows that the amount of excitation required to generate the standard CST target decreased as neuromodulation increased or as inhibition progressively changed from proportional to push pull. Thus, although the solution space for reproducing the target CST is very wide, encompassing the full range of neuromodulation and inhibition investigated, achieving this wide range required large changes in the pattern and amplitude of the excitatory drive.

Figure 5.

The solution space for firing patterns of motor units

If each of patterns of excitation, neuromodulation and inhibition that generated the the standard triangular CST generated highly similar motoneuron firing patterns, then the solution spaces for each set of firing patterns would each be large. This firing pattern similarity however did not occur. The heat map in Figure 6 provides a quantitative analysis of one of the main differences in motor unit firing patterns, onset-offset hysteresis (this hysteresis was quantified by the ΔF method developed for studies of real motor unit firing patterns in humans (Gorassini et al., 2002) (see Methods and Figure 7)). The differences in neuromodulation and inhibition produced a huge range of hysteretic firing behaviors, ranging from < 1.0 at low neuromodulation levels coupled to strong proportional inhibition (lower left in Figure 6) to > 7.0 for high neuromodulation and strong reciprocal inhibition. Fig. 6 also shows the simulated firing patterns generated by a subset of these combinations (arrows link these patterns to the hysteresis values). As neuromodulation increases, PIC effects (acceleration, attenuation, hysteresis) increase, with push-pull inhibition emphasizing these effects and balanced inhibition suppressing them. These differences for each combination of neuromodulation and inhibition shown in Figure 6 demonstrate that the solution space for matching a particular set of motor unit firing is much more restricted than for matching the overall CST (see Figure 4 and Figure 5).

Figure 6.

Figure 7.

Note however there is some degree of trade-off for effects of inhibition and neuromodulation on hysteresis along the diagonals in Figure 6. Fortunately, hysteresis is not the only index of motoneuron firing behaviors. Figure 7 illustrates additional features of motoneuron firing patterns that are easily measured and likely to reflect other effects of excitation, inhibition and neuromodulation. These features include firing rate acceleration and attenuation (which are also driven by actions of neuromodulation and inhibition on PICs – see Figure 1C1 and Figure 1C2) as well as recruitment and de-recruitment behaviors (sensitive to the distribution of excitation and S vs F motoneurons – see Figure 1A2). Each of these other characteristics had patterns of variation with respect to inhibition and neuromodulation that varied in comparison to each other and to hysteresis. Therefore, to quantify the effectiveness of RE of motor unit firing patterns to predict the organization of synaptic inputs, we next used multiple regression methods to determine how well combinations of these characteristics of simulated firing properties predicted the level of neuromodulation, the patterns of inhibition and the distribution of excitation to low verse high threshold motoneurons.

Regression using motoneuron outputs to predict input organization

We utilized several different types of linear regression techniques (see Methods) and one nonlinear technique (Support Vector Regression using a radial basis function kernel (RBF)). In each case, the regression parameters were calculated from a subset of the data (∼70% of the total) and then used to predict the remaining subset (∼30%). For neuromodulation, all regression methods were similar in their effectiveness.

Figure 8A shows the results for standard linear regression (the dashed red line has slope = 1.0 to indicate perfect predictive performance). For each level of neuromodulation (0.8 to 1.2), the vertical scatter is due to various combinations of the pattern of inhibition (which ranged from 0.7 to −0.7, as is in Figure 6) and of the distribution of excitation to low versus high threshold motoneurons (which ranged 0.5 to 2.5). The scatter is considerable, but there is also a tendency for points to cluster near the line. As a result, the regression relation accounts for 76% of the variance (i.e. R² = 0.76; Table 1 shows the R² values for all regression methods for prediction of neuromodulation, inhibition and the distribution of excitation). Thus, prediction of the level of neuromodulation is good and indicates that this regression approach can identify the level of neuromodulation with a resolution of about 0.1 units of neuromodulation within the range 0.8 to 1.2.

Table 1.

Figure 8.

For inhibition, the nonlinear regression was most effective, giving the results shown in Figure 8B. Variance accounted for is 0.85, implying a prediction resolution of about 0.23 units of a range of 1.4 units of inhibition. Finally, Figure 8C shows the regression results for prediction of the relative weighting ratio for the motoneuron pool, implying that characteristics of motor unit firing patterns can predict input structure with variances accounted for ranging from 76-91%.

As a final step, we analyzed which motor unit firing characteristics had the most predictive power in identifying each of the three input organizations. Figure 9 provides two perspectives on this question. The left column shows results of step-wise regression, starting with the characteristic with the strongest predictive power and then progressively adding the others. For neuromodulation, all regression methods behaved similarly, in that the first two characteristics (brace-height and ΔF see Figure 7) was highly effective in reducing error in the regression prediction, with little addition benefit from adding the others. For inhibition prediction via the linear regression models, ΔF was most effective, followed by α_sat, with all others being ineffective. This is shown in the mutual information ranking (right side) in Figure 9, which ranks the characteristics (features) by contribution of new information to the regression model. The left side of Figure 9 ranks the features by the F-statistic, which does not necessarily imply that new information is introduced into the regression model. This is shown by the very subtle improvements in the error (mean squared error) as features are cumulatively added to the regression model, until t_drac and α_sat are added to the model.

Figure 9.

However, the nonlinear regression method (RBF SVR) was more effective than any of the linear methods. For this method, each of the firing characteristics had a notable effect in reducing error. For linear regression for prediction of excitation, t_drac (the average time when units were derecruited) was the most impactful and adding the other characteristics brought little improvement. The right column in Figure 9 shows mutual information for each characteristic, which reflects how well each one predicts the input behaviors on its own. As expected, the order of effectiveness of prediction in these plots is identical to that for the step-wise regression.

Overall, neuromodulation is best predicted by brace-height and ΔF. Inhibition is also predicted by ΔF while brace–height is ineffective. Excitation (excitatory weight ratio) is best predicted by parameters that assess recruitment and derecruitment timings and less informed in prediction by ΔF and brace–height. These different roles for RE prediction of input make sense in terms of the physiological effects of the three types of inputs on motoneuron firing patterns (see Discussion).

Discussion

Our simulation shows that the solution space for a motor output variable that reflects the net action of all active motoneurons is indeed very broad. The variable investigated here, the cumulative spike train (CST), has been shown to be closely proportional to muscle EMG and muscle force (Thompson et al., 2018), which are widely assessed in experimental paradigms in both animals and humans. Our simulations however indicate that RE analyses of CSTs, EMGs and forces are likely to fail due to the problem of non-uniqueness. In contrast, motoneuron firing patterns provide a much greater level of information about motoneuron inputs. Our simulations showed that the solution space for a set of motoneuron firing patterns is reasonably small, allowing predictions that account for 75-90% of the variance in the level of neuromodulation, the pattern of inhibition and the distribution of excitation to S vs F motoneurons.

Generalization

Array electrodes are extensively used in many parts of the CNS to record firing patterns of populations of neurons in awake behaving animals (see Introduction). Moreover, highly realistic modelsof neurons and their interconnections in many systems are already available and undergoing rapid development, including supercomputer implementation (Einevoll et al., 2019). Since models and firing patterns are the fundamental elements for success of RE in the simulations presented here, this same type of approach may be effective in many other neural systems. In each neural system, however, success of RE of firing patterns into input structures as simulated here is likely to depend on the particular interactions between neuromodulation, inhibition and excitation. For motoneurons the nonlinearities induced by PICs in motoneuron firing patterns have proven to be a major advantage for RE. These nonlinearities impart distinctive characteristics to firing patterns (acceleration, attenuation, hysteresis – see Figure 1C1, the firing pattern at the upper right in Figure 6), which are easy to quantify (see Figure 7). Moreover, because the PIC is highly sensitive to inhibition, these same firing characteristics provided effective indices of the temporal pattern of inhibition. Furthermore, because these relations between neuromodulation and inhibition were different for each firing characteristic (see Supporting Information Figure S1), multiple regression using these characteristics was effective. It seems likely that if motoneurons were nearly linear input output processors, non-uniqueness would be a much greater problem for RE.

Muscle EMG versus single motor unit firing patterns

Muscle EMG is the electrical signal generated by the summed action potentials of the muscle fibers of all active motor units. EMG is easy to measure in many types of motor behaviors in humans and other animals and provides an effective measurement of the timing and overall pattern of the neural output from the spinal cord to muscle. However, the simulations here and in our previous work (Powers et al., 2012; Powers and Heckman, 2017) show that EMG does not provide information about the organization of motor commands that produce spinal motor output. On the other hand, because motoneurons action potentials are one-to-one with those of their muscle fibers in healthy humans and other animals, motor unit firing patterns are equivalent to motoneuron firing patterns. While motor unit firing patterns are more technically difficult to measure than whole muscle EMG, the accessibility of muscle makes these measurements highly feasible. In fact, motor unit firing patterns in humans were recorded at the very onset of single neuron studies, over 90 years ago (Adrian and Bronk, 1929). The likelihood that this single neuron information provide deep insights into the structure of motor commands has thus long been appreciated (Duchateau and Enoka, 2011). Many such studies were essential for establishing the prevalence of Henneman’s size principle of recruitment (reviewed in Heckman and Enoka (2012)). Many studies also focused on inferring the time course of EPSPs and IPSPs driving motoneuron responses to transient ionotropic inputs during steady contractions (reviewed in Türker and Powers (2005)). The distinctive effects of PICs on firing have also been a focus, with the ΔF method developed by Gorassini and colleagues for quantifying PIC-induced hysteresis having become a standard in the field Gorassini et al. (2002). However, intracellular current and voltage clamp studies have demonstrated that the PIC is highly sensitive both to neuromodulatory input (Hounsgaard et al., 1988; Lee and Heckman, 2000, 1999) and inhibition (Hyngstrom et al., 2008; Kuo et al., 2003). Consistent with this, ΔF does not provide good discrimination of neuromodulation versus inhibition in our simulations (Figure 9). In contrast, and consistent with our previous study (Beauchamp et al. (2023)), the brace–height parameter was sensitive to neuromodulation but insensitive to inhibition, likely because it is normalized for peak firing. In general, our simulations show that no one motor unit firing characteristic is sufficient for estimation of input structure and thus a form of RE is necessary.

Limitations

Ensemble model approaches for RE are obviously reliant on the biological accuracy of the constituent models. Our models have been carefully tuned to replicate experimental data from years of current and voltage clamp studies of the integration of neuromodulatory, inhibitory and excitatory inputs motoneurons in feline preparations. We made a few adjustments to these models to more closely replicate the slow firing rates in human motor units but will likely need to further adjust parameters like the spike AHP and the PIC voltage threshold. Our RE was not applied for identification of the amplitude or time course of the excitation. The peak excitatory conductance for each simulation strongly depends on the level of neuromodulation and the pattern of inhibition (Figure 4). This correlation was built into our simulations to match the CST for each combination of neurmodulation/inhibition and so excitation was not an independent parameter. Further simulations in which we relax this constraint or allow inhibition and/or neuromodulation to vary for a given pattern of excitation may prove interesting. It will also be interesting to expand our range of neuromodulation. At the high end, this expansion will have to be limited, as controlling such large PICs is difficult. But exploring whether RE works as at lower neuromodulation levels where input-output processing is nearly linear would reveal whether this simpler behavior limits RE effectiveness.

Summary & Future Application

Although our models likely need modification before we begin RE of human motor unit firing patterns, they were well situated to the goal of this study, which was to investigate the feasibility of RE of neuron firing patterns despite the strong tendency for neural systems to be non-unique. In addition, it would seem that the success of our RE approach suggests that this technique can be used to gain insight into the motor commands underlying different motor unit firing patterns associated with different motor tasks with or without disease states (e.g., Mottram et al. (2014)).

Methods and Materials

Motoneuron pool model

The motoneuron pool consists of 20 model motoneurons with a range of intrinsic properties. This motonouron pool size was chosen to reflect the typical sample size of motor units discriminated based on surface EMG array recordings. Each model motoneuron consisted of a soma compartment and four dendritic compartments, each coupled to the soma. The passive properties of the compartments, such as, size, capacitance and passive conductance were derived from works by described by Kim and colleagues (Kim et al. (2009); Kim and Jones (2012)). Spike conductances (g_Na and g_K) and conductances mediating the medium AHP were inserted into the soma compartment and a calcium conductance mediating the slowly-activating persistent inward current (PIC) was inserted into each of the dendritic compartments. In addition, a hyperpolarization-activated mixed-cation (HCN) conductance was inserted into all compartments. Conductance densities, kinetics and steady-state activation curves were originally tuned to recreate the range of input-output behavior recorded in medial gastrocnemius (MG) motoneurons in decerebrate cats, as described in Powers and Heckman (2017). We made four modifications to model parameters for the cat MG pool in order to produce the lower firing rates and less than full recruitment typically observed during moderate voluntary contractions in human subjects. First, the range of values governing excitability were restricted to the first 75% of the original recruitment range (i.e., the parameters of motoneuron 20 in the new model corresponded to motoneuron 15 in the original parameter set). Second, AHP durations and amplitudes were increased by increasing the values of the time constant of calcium removal for the calcium-activated potassium conductance (from 60 to 10 ms in the original model to 90 to 57 ms). The larger and longer AHPs acted to oppose early PIC activation during increasing excitatory synaptic drive, so we hyperpolarized the PIC half-activation threshold (from −40 to −37 mV in the original model to −42 to −40.4 mV). Finally, in the original model the slow decay of PICs observed in high threshold MG motoneurons (Lee and Heckman, 1996) was replicated by including a voltage-dependent inactivation process. We found that this process limited the firing rate hysteresis to values below those typically seen in human motor unit recordings, so we eliminated PIC inactivation in the present model pool. All simulations were run using the Python interface to the NEURON simulator (Hines and Carnevale, 1997). NEURON files specifying motoneuron pool parameters, conductance mechanisms, and protocols on on ModelDB with accession number 239582.

Simulation protocol

Inputs to the motoneuron pool

The target motoneuron pool output, representing the average firing rate across the entire pool of 20 motoneurons was 22 seconds in length and consisted of 1 sec delay followed by a linear rise to a peak value of 16 imp/s over the next 10 sec, and a linear decline back to zero in 10 sec (Figure 3, Reference Ref). A scaled version of this command, I_ex = 0.6Ref, (Figure 3, light salmon colored with the value 0 to its bottom left) was used as the first guess of the time course of the excitatory input needed to produce the target output. For each iteration an inhibitory conductance input (I_in) was also applied according to the following equation:

where the excitatory input I_ex is multiplied by a gain G_in and biased by b_in as defined by:

where R_mod is the neuromodulation value, set in this work, from 0.8 to 1.2 at increments of 0.1 (i.e. [0.8:0.1:1.2]). The inhibition is therefore coupled to the neuromodulation through its bias, b_in, and was set to value just sufficient to ensure PIC deactivation at the end of the command. The inhibitory gain, G_in, is varied from −0.7 to 0.7 in increments of 0.1 thus traversing the possible control schemes push-pull inhibition to balanced inhibition. Example of both control schemes can be seen in Figure 2 (green for excitatory and red for inhibitory).

Changes in neuromodulation were simulated by changing the density of dendritic PIC channels: from 80% to 120% of the standard density in 10% increments. Mathematically, the MN model was altered by changing the current of the L-type calcium channel according to these equations:

where is the current for L-type calcium channel, is the conductance of the L-type calcium channel, is the activation gate, V is the compartment voltage, and is the reversal potential for the L-type calcium channel. R_mod is the gain factor that changes the maximum conductance of the L-type Calcium channel with the following values [0.8, 0.9, 1.0, 1.1, 1.2].

Finally, filtered Gaussian noise was added to both the excitatory and inhibitory commands with a standard deviation that varied in proportion to the square root of the mean level and a decay time constant of 20 ms to reflect the fact that most of the power in the synaptic input to motoneurons is thought to be in the low frequency (< 10 Hz) range (Farina et al., 2014).

Feedback and optimization

We used an iterative optimization procedure to find the motoneuron pool input needed to produce the “Reference” output as shown in Figure 3. We opted for a feedback type of optimization method to converge onto the solution.

For the feedback to work, an error term is calculated between the output of the MN pool and the “Reference” (see the “Controller” Figure 3). The output of the MN pool is the average of the 20 MN pool calculated by convolving each model spike train with a 2-s Hanning window with unit area, adding up the result across all 20 spike trains and dividing by 20:

where is the new excitatory input, n is the iteration in the feedback loop, F (t)_i is the firing rate of a MN i = [0 : 20] in the 20 MN pool, H(t) is the 2-s Hanning window, e(t) is the error and K is the proportional gain (usually 20%). The resulting (the average firing rate response) plus the error e(t) is then used as the new excitatory command to the MN Pool. This is shown by the series of traces leaving the “Plant System - Motoneuron (MN) pool” in Figure 3. Using a new inhibitory command is also calculated using equation 1 and is subsequently used in the new computational iteration n.

This process was repeated until the mean squared error between the output and the target was less than 1 Hz² (generally less than 0.5 Hz²). In the case shown in Figure 3, good convergence was reached by the fifth iteration.

This iterative matching procedure was repeated for 225 different input combinations: five levels of neuromodulation, three different distributions of excitatory input, and fifteen different levels of inhibition.

Excitatory distribution or MN weights

The density of the excitatory synaptic conductance was either uniform across all cells, weighted toward low threshold units (synaptic weight onto lowest threshold unit 2 X that of highest threshold unit), or weighted toward high threshold units (high weight 2X low). The parameters governing cell excitability were skewed to produce a greater proportion of low threshold units, so synaptic weighting showed a similar skew according to the following equation:

where w_MN# is the weight of the given motoneuron, #MN is the number of the MN in the pool (e.g. 20), w_end and w_start are the weight range. For weights favoring low threshold units, w_start = 2.5 and w_end = 1.0 and for the opposite w_start = 1.0 and w_end = 2.5. In total we used 7 different types of weight configurations: [w_start, w_end] ⟺ [1, 1], [1.5, 1], [2, 1], [2.5, 1], [1, 1.5], [1, 2], [2.5, 1].

Noise Seed

We chose to calculate the model using 30 different noise seeds. We assumed that each seed would represent a “unique” individual and allow for the use of averages and variance in our interpretation of the results. Note: The noise seed was not reset for each of motoneurons in the pool. The synaptic noise is thus coupled for each of the motoneurons which fits with the “common drive” argument (see section Common input structure, differences in motoneuron properties).

Massively-parallel Data Generation Methods

Since the majority of the simulation is able to operate independently as the simulations for each of the parameter sets (noise seeds, inhibitory gains, and neuromodulation) did not need to feed into each other, we were able to run this simulation in a perfectly parallel way. Each unique parameter set was simulated using one node (computer) in a supercomputing cluster, which produced a dataset file of the simulation results (see Algorithm 1 in Supporting Information). Each node simulated one motor pool, with the computation of each motoneuron run on each processor on the compute node (see Algorithm 2 in Supporting Information). To implement this setup, we utilized the MPI interface in NEURON 7.7 Hines and Carnevale (1997) and Python 3.6.7 running on the Bebop supercomputer (Intel Broadwell nodes) at Argonne National Laboratory, (Supplementary Figure 2). Each combination tuple tuple_comb contains a particular ginmult, nmmult, ginadd, and random_seed (if noise is used).

Feature Extraction

For each simulation input combination, we analyzed the individual motoneuron discharge patterns to extract 7 features: Recruitment Time (t_rec), De-Recruitment Time (t_drac), Activation Duration (t_dur), Recruitment Range (t_range), Firing Rate Saturation (α_sat), Firing Rate Hysteresis (ΔF) and brace–height (brace–height). The calculation of these quantitative measures is illustrated in Figure 7, which shows the instantaneous discharge rates (dots) of a higher (upper panel) and lower threshold (lower threshold) motoneuron, along with the smoothed firing rates calculated by convolving the spike times with a 2-s Hanning window (continuous lines).

Recruitment Time (t_rec)

Recruitment time is the time of the first instantaneous discharge rate of a given motorneuron.

De-Recruitment Time (t_drac)

De-Recruitment time is the time of the last instantaneous discharge of a given motoneuron.

Activation Duration (t_dur)

Activation duration is the difference between de-recruitment time and recruitment time:

Recruitment Range (t_range)

The recruitment range is the difference between the recruitment time of the last and the first motorneuron in the motorpool.

Firing Rate Saturation (α_sat)

The motorneuron firing rate is characterized by an initial acceleration phase that lasts about 1s followed by partial saturation of firing rate. The degree of this saturation is quantified by calculating the rate of firing rate increase from 1 sec after recruitment to the target peak (11 s):

where f₃ is the firing rate at 11 s, f₂ is the firing rate 1 s post recruitment (t_rec), t₃ is time at target peak 11 s and t₂ is 1 s post recruitment.

Firing Rate Hysteresis (ΔF)

A motorneuron rate hysteresis can be quantified using the ΔF measure introduced by Gorassini and colleagues (Gorassini et al., 2002). ΔF is defined as the difference in the smoothed firing rate of a lower threshold reporter motorneuron (Figure 7 - bottom) at the onset (t_rec, f₀) and offset (t_drac, f₁) of the higher threshold motorneuron (Figure 7 - top):

brace–height (brace − height)

Brace height is a measure of the non-linear portion of the rising phase of the firing rate of the motoneuron (see Figure 7 - top trace). It is the maximum perpendicular distance from the line linking the firing rates at recruitment time (t_rec) to t₃, the time at which the triangular “Reference” reaches its maximum (11 s). This distance is calculate by solving the geometric problem of finding the “Distance from a point to a line”. This analysis is further detailed by Beauchamp et al. (2023).

We restricted the calculation to motorneuron pairs for which the recruitment time difference between the two was more than 1 s. This restriction was imposed to ensure that PIC activation in the reporter unit had reached near steady state, so that further changes in reporter discharge rate primarily reflected changes in synaptic drive (cf. Powers et al. (2008)).

Machine learning inference of motor pool characteristics

Data preparation

We ran the model for a range of neuromodulation (5), inhibition (15), excitatory weights (7) and noise seeds (30) for a total of 15,750 combinations. Given 20 motoneuron per pool and 20 optimization iterations (see Figure 3 or section Feedback and optimization we ran a grand total of 6,300,000 simulations.

Post simulations, we extracted the features from the firing rates of each motoneuron (see section Feature Extraction and organized them according to the neuromodulation, inhibition, excitatory weight and noise seed. We then first calculated the average of these features accordingly and then took an average of these averages with respect to the noise seed. In final, the was an average of the features grouped by neurmodulation, inhibition and excitatory weight. Supplementary

Figure 1 shows the average of each feature in a heat map configuration for the excitatory weight distribution [w_start, w_end] = [1, 1]. This data was the used subsequently (see section Model fits and Ranking features).

Model fits

One of the goals of our RE approach is to infer the characteristics of the input to the motoneuron pool based on our measures of motor unit output. The numeric features extracted from the discharge patterns are used to create five-dimensional input data for the models. The dimensions were the mean ΔF value, times of recruitment (t_rec) and derecruitment (t_drac), the firing rate saturation degree (α_sat), the recruitment range (t_range) and brace–height (brace − height). The target values to predict were the simulation input values, namely level of neuromodulation, inhibition, and excitatory weight. Since all target values to be predicted in this study were continuous numeric values, the RE process was posed as a regression task.

Since the dataset is relatively small (<1000 instances) and has only 5 dimensions, linear methods and kernel methods were chosen for predictive analysis. More complex models, such as tree-based methods and neural networks, were omitted due to inadequate training data instances and interpretability. Linear regression models using standard least-squares, Lasso, Ridge, and ElasticNet methods, and support vector models using linear and radial basis function kernels were created for each simulation input value to be predicted using the mean squared error as the minimization cost function.

Ranking features

To measure the performance of linear regression against the firing features (see section Feature Extraction), we use a greedy approach towards incorporating more and more features into the linear regression model. For this approach, the F-statistic (linear regression test) is used as the scoring function to order-rank the features for incorporation into a linear regression model.

We used mutual information (MI) (Ross, 2014; Kraskov et al., 2004; Kozachenko and Leonenko, 1987) to rank the dependence of the firing rate features (see section Feature Extraction) to the inhibition, neuromodulation and excitatory weight ratio levels of the motoneurons in the pool model.

Acknowledgements

This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357. We gratefully acknowledge the computing resources provided on Bebop, a high-performance computing cluster operated by the Laboratory Computing Resource Center at Argonne National Laboratory. This research was supported by the National Institute of Neurological Disorders and Stroke of the National Institutes of Health under awards R01NS062200 and R01NS125863.

Appendix 1

Algorithm 1

Python pseudocode for reverse engineering proportional controller on each compute node

Algorithm 2

Python pseudocode for simulating each motoneuron independently

Appendix 2

Supplementary Figures

Appendix 2—figure 1.

Appendix 2—figure 2.