Authors: A. Del Vecchio, F. Sylos-Labini, V. Mondì, P. Paolillo, Y. Ivanenko, F. Lacquaniti, D. Farina
Categories: Research Articles, Neurophysiology, Neuroscience, Research Article, SciAdv r-articles
Source: Science Advances
Human neonatal motoneurons show high synchronized discharge timings that modulate the generation of fast movements.
Mammalian neonates show a repertoire of coordinated movements important for survival, produced by central pattern generator (CPG) networks (1) under sensory and brainstem/spinal cord control. Alpha motoneurons are the final common pathway, interfacing CPGs with the external environment. In neonatal rodents, motoneurons show a millisecond synchronization (2, 3). Highly synchronized neural activations might be crucial to perform fast actions and to reorganize neural circuits for complex behaviors (2, 3). Synchronized activity of motoneurons has been associated to motoneuron path-finding (2), synapse maturation [e.g., establishing the initial polyneuronal innervation and subsequent synapse elimination (4)], refinement of pattern-generating circuits, and modulation of cortical somatosensory maps (5). The activation mode of motor units (i.e., synchronization) rather than their intensity (e.g., discharge rate) seems to control synaptic integrations and connections at different levels of the neuromuscular system in a Hebbian fashion (4, 6).
Although humans have limited motor abilities at birth, CPGs allow newborns to eat, breathe, and step or kick automatically (7, 8). Moreover, neonates are endowed with a set of primitive reflexes involving fast movements, which are important for protection, nutrition, and survival or are simply phylogenetic rudiments (9). However, very little is known about the behavior of individual motor units in human neonates due to the lack of adequate recording methods (10). It is unknown whether the human neonatal spinal cord behaves in a way similar to other mammals. An interface with the output circuitries of the human neonatal spinal cord would allow us to understand the behavior and development of the neural cells that ultimately control force generation. Here, we show and validate a technique that identifies individual motoneuron discharges from leg muscles of human neonates (Fig. 1). With this technique, we investigate the discharge rate and synchronization of motoneurons in neonates, as associated to the generation of fast leg movements. Because of the relatively slow twitch forces exerted by neonates’ muscle units, we hypothesized that fast leg movements in neonates are mediated by high motoneuron synchronization.
Fig. 1 Interfacing the neonatal spinal cord.(A) Alpha motoneurons in the ventral horn of the spinal cord. (B) A representative neonatal movement that was recorded via 128 high-density EMG electrodes, covering most of the lower leg of the neonates. The spatial information in the high-density EMG recording allowed decomposition via blind source separation techniques into the constituent motoneuron action potential discharge timings (C). The box in (C) shows an example of high synchronization of the neonatal motoneurons.
We developed a noninvasive, high-density neuromuscular interface to decode the motor unit discharge timings in freely moving human neonates via gradient convolution kernel compensation [Figs. 1 and 2; for review, see (11)]. We applied this technique to analyze the electromyographic (EMG) activity of tibialis anterior muscle in four babies (2 to 14 days old) laying supine on a table, and we were able to identify the activity of 31 motor units (Figs. 1 to 3, figs. S1 and S2, and movie S1).
Fig. 2 High-density EMG activity and motor unit discharge times in neonates.(A) Two grids of 64 electrodes (128 electrodes in total) covered most of the anterior lower leg of the neonates. The EMG grids were aligned with the direction of the tibia bone. The interelectrode distance was 4 mm. Photo Alessandro Del Vecchio, Imperial College London. (B) Ten EMG channels showing common myoelectric activity during spontaneous movements. Vertical red lines indicate the time intervals when the RMS values were computed for each EMG channel. RMS was projected in 2D maps. (C) Four interpolated EMG-amplitude maps. The time windows correspond to the red lines in (B). These maps show distinct patterns of activation, resulting from activation of different muscles or distinct motor unit innervation zones. We then applied blind source separation on the EMG signals and identified the constituent motor unit action potentials. (D) Average EMG activity across all 128 channels with a 2-Hz low-pass filter (redline). (E) Raster plot of the discharges of six identified motor units. (F) The instantaneous discharge rate of these neurons displayed in a 10-s interval.
Fig. 3 Validity of blind source separation on high-density EMG signals.(A) Raster plot of four identified motor units during neonatal spontaneous movements. (B) The average discharge rate for the neurons. (C) The output of the gradient convolution kernel compensation algorithm applied to the high-density EMG signals. We used the time instants identified by the algorithm to spike-trigger average the high-density EMG signals in 25-ms time intervals. (D) The spike-triggered average of the motor unit number #3 (green). The EMG channels in this example are displayed in monopolar derivation with each individual spike superimposed (black line denotes the average action potential). (E) The average motor unit action potential across the full duration of the trial for the 128 electrodes. The colormaps in (E) and (F) represent the RMS amplitude of the action potentials. Note that each motor unit has its unique territory. AU, arbitrary units.
Figure 2 shows an example of 10 EMG signals recorded during a 60-s trial. After projecting the amplitude of the EMG [root mean square (RMS) value over 250 ms; see Materials and Methods] in the two-dimensional (2D) maps of the high-density grids (Fig. 2, B and C), we observed distinct patterns of muscle activation. To identify the discharge timings of individual motor units, we decomposed the full high-density EMG signal (128 channels) via blind source separation. These procedures have been detailed previously (11) and validated under a vast range of conditions involving healthy adult individuals (12), patients with Parkinson’s disease (13), and following nerve transfers (14). Here, we tailored these methods for neonates, and we succeeded in reliably identifying motor unit activity (see Fig. 3 and fig. S1).
We tested the accuracy and validity of decomposition by different metrics. First, we designed a custom software that would plot the motor unit spike rasters (Fig. 3, A to C), innervation pulse trains (output of the convolutive blind source decomposition; Fig. 3C), and the 2D waveforms and shimmer plots of the motor unit action potentials (Fig. 3D) extracted by spike-triggered average. The shimmer plots represented action potentials of the same motor units superimposed in a 25-ms time window centered at the motor unit spike instants extracted by decomposition (see Fig. 3D and Materials and Methods). This same time window was applied to all EMG channels. Figure 3D shows shimmer plots for four electrodes. Moreover, for validation purposes, the motor unit spike trains (binary code) were altered with random perturbations (noise), with the assumption that even small perturbations (5 to 20% of interpulse intervals) would significantly corrupt the motor unit template identified by decomposition (fig. S1). We also calculated the pulse-to-noise ratio (PNR) of the motor unit innervation spike trains, as proposed in (15), and only motor units showing >26 dB were considered for further analysis. The motor unit action potential waveforms identified in each neonate showed a 2D cross-correlation >70% when compared to the random identified motor units (<30%; see fig. S1).
Figure 3 (D to F) shows the shimmer plot and the location of representative motor units (RMS amplitude of each electrode) plotted in 2D grids. The location of motor units was consistent across the whole duration of the task (fig. S1), with distinct hotspots (Fig. 3, E and F). Validity of the decomposition was therefore verified by the PNR, the identification of the same motor unit action potentials across the task duration, and by the analysis of sensitivity to added noise.
After testing the EMG decomposition accuracy, we evaluated the properties of the discharge patterns of the motoneurons and compared them to those of 12 adults performing a wide range of contractions. Figure 4 shows representative data for an adult and a neonate and the motor unit characteristics across all subjects. The motor unit discharge rates observed in the neonates were similar to those of the adults performing contractions at 10 to 30% of the maximal voluntary force (MVC; Fig. 4H). Average motor unit instantaneous discharge rates were 16.0 ± 4.8 pps in neonates and 10.0 ± 1.4 pps (low forces) or 12.7 ± 1.7 pps (medium forces) in the adults (adult versus neonate motor unit discharge rate, Kruskal-Wallis P = 0.3). The average cumulative discharge rate for the motor units during the slow and fast movements showed similar values in all neonates (Fig. 4H).
Fig. 4 Estimates of motoneuron output for the adult and neonatal human spinal cord.(A) Twelve adults performed isometric ankle-dorsiflexion contractions at 10 and 30% of maximal force for 50 s. (B) The experimental setup and a representative isometric contraction with the corresponding motor unit activity. (C) Motor unit spike raster and after convolution with a 1-s Hann window [smoothed spikes per second (D)] for a representative newborn during spontaneous movements. The black line in (C) represents the average discharge frequency calculated after convolving the sum of the spike timings. Note the highly synchronized discharge timings in the scaled version (E and F). The cross-correlation value from the smoothed discharge timings was 0.52. (G) The cross-correlation function for 50 permutations of two equally sized groups of motoneurons (gray lines) and the corresponding average (black line). (H) The average motor unit discharge rate of the neonate’s individual movements and for the adults during the two tasks. The discharge rate was averaged in the intervals corresponding to slow or fast neonatal movements. Moreover, the average value of the instantaneous discharge rate across the full duration of the tasks (neonate inst.) is also reported. (I) The average motoneuron synchronization for each neonate and the adults. The neonates showed a higher synchronization value compared to the adults during both slow and fast movements (*P < 0.001). Moreover, the neonates exhibited higher synchronous firing during fast movements. (J) The average coherence across all neonates. The dashed red line indicates the significance threshold.
We then computed time-domain cross-correlations between the discharge timings of the pools of identified motoneurons (Fig. 4, F to H). It is known that, in adults, the estimated degree of synchronization between motor unit action potentials is positively correlated to the average discharge rate of the motoneurons population (16). Critically, however, we found that, although the neonate motoneurons showed similar discharge rates as those of the adults (Fig. 4H), neonate motoneurons exhibited a higher level of synchronization during both low- and high-speed movements (Fig. 4I). The average cross-correlation value was 0.19 ± 0.09 for the neonates and 0.044 ± 0.023 or 0.040 ± 0.021 for the adults at 10 or 30% of MVC, respectively (Kruskal-Wallis P < 0.001). In previous studies (2, 3), neonatal rat motoneurons showed synchronizations similar to our human newborns [values between 0.1 and 0.3; bin width of cross-correlogram of 100 ms (2, 3)].
We further computed the value of synchronization in two frequency bandwidths. Motoneuronal output shows common oscillations within the 0- to 40-Hz bandwidth (17). However, the bandwidth in which muscles produce force is much smaller (0 to 5 Hz) (18). Thus, to smooth motoneuron spike times, we used either a short Hann window (40 ms; Fig. 4, C and D) or a longer window (1 s; Fig. 4, C to F). The Hann window is a smoothing function with shape proportional to the square of the cosine. The smoothing depends on the window duration. The short duration of 40 ms corresponds to retaining most of the data of the neural oscillations (0 to 40 Hz) while the longer duration of 1 s retained mainly the low frequencies (0 to 5 Hz) that correspond to the electromechanical filtering obtained by the muscular apparatus. We hypothesized that the correlation of motoneuron firings in the low-frequency range would show a higher synchronization than in other bands because it directly translates to force due to the filtering by the muscle dynamics. Accordingly, we observed a higher correlation at the low-frequency bandwidth with a value of 0.56 ± 0.14 (P < 0.0001; Fig. 4C) with respect to the overall correlation for the entire 0- to 40-Hz band (shown above). For all correlations, including neonate and adult data, the time lag was ~0 ms (maximum of <1 ms). The width of the cross-correlogram was also limited to short time frames [average range of the correlation peak (−26.30 ± 1.60 to 27.42 ± 2.52 ms)], and the values were not statistically different for the adults (Kruskal-Wallis P = 0.4). These small ranges indicate high synchronization within short time scales of the neuronal spinal cord output.
The values of synchronization computed above were averaged across the full duration of the movements of the neonates, and silent periods greater than 500 ms (19) were removed from the analysis (see Materials and Methods). Since there is evidence that synchronization modulates the development of the neuromuscular system in a Hebbian fashion, we further hypothesized that high synchronization may be needed for neonates to produce fast rates of force development. While in healthy young adults the main determinant of rate of force development is the fast recruitment of motor units and high discharge rates of the motoneurons (12, 20, 21), motor unit discharge rates in neonates are relatively low during fast movements (see Fig. 4H). We therefore checked whether high values of synchronization corresponded to specific movements of the neonates, with the hypothesis that higher synchronization generated faster rates of force development.
Thus, we analyzed the neonates’ movements recorded with a digital camera synchronized with EMG recordings. The recordings involved either spontaneous leg movements or leg movements elicited by light mechanical stimulation of the plantar surface of the foot or other body parts. Figure 5A shows the estimates of neonate movements in the horizontal and vertical planes along with the concurrently recorded EMG signals. Notice that the strength of EMG and motor unit activity was highly correlated with the speed of foot movements (Fig. 5A and movie S1).
Fig. 5 2D pose estimation of the neonates with EMG and motoneuron activity.(A) The change in position of the big toe across the x (red) and y (blue) positions [as shown in (B)] is plotted as a function of time concurrently with the EMG signal (green line) and spike raster of the motor unit discharge timings. Note the synchronized activity between the neural stimuli and movements. The movements of the neonates are shown for two representative images corresponding to the dotted magenta and black lines, respectively, in (B) and (C). Photo Alessandro Del Vecchio, Imperial College London. We tracked the positions of nine markers corresponding to the five toes, two EMG connectors, and two locations on the neonate belly as shown in (D). (E) Bivariate correlations between the finger positions (R > 0.9 for all cases, P < 0.0001). (F) The maximum speed, mean speed, and average distance covered by the fingers of the neonates are shown. *P < 0.0001.
Although several movements were relatively slow (1 to 10 cm/s), some movements reached maximal speeds of up to 400 cm/s (see Fig. 5, A and F, and movie S1). Even considering the low inertia of the neonatal limbs, these movements are very fast. For comparison, kicking by adult taekwondo athletes reaches speeds of 700 to 900 cm/s (22). For all neonates, the synchronization values were significantly higher (up to ~50% times as high) during faster with respect to slower foot movements (P < 0.001; Fig. 4I). These values were significantly higher compared to the adults for both slow and fast neonatal movements (Fig. 4I). It must be noted that the higher values of synchronization in neonates were not due to a difference in discharge rate between contraction speeds. We did not find significant differences in discharge rate during fast and slow movements of the neonates (Fig. 4H). Thus, contrary to previous observations in adults (12), an increase in contraction speed in the neonates was not mediated by an increase in discharge rate, but rather, it was exclusively associated to an increase in motoneuron synchronization.
We further assessed correlated input to motoneurons in the frequency domain by coherence analysis across the unfiltered motor unit data (spike trains) (17, 23, 24). Significant peaks in the coherence reflect common fluctuations of the motoneuronal discharge timings. We found that all neonates exhibited significant coherent oscillations at low frequencies (<3 Hz; Fig. 4J), which are known to directly translate into muscle force during steady (25) and oscillatory contractions (26), and are in agreement with the time-domain analysis reported above. We did not find any significant coherence, except for one neonate, in the high-frequency band (beta band, 15 to 30 Hz), which is thought to be associated to cortical dynamics (24).
We presented a noninvasive neural interface that can decode individual neuronal activity in freely behaving human neonates. Our interface was able to reliably decode the activity of more than 30 motoneurons and to associate the features of motoneuron discharge to kinematic variables obtained from a digital camera and machine learning processing. We found that human neonates, in a similar way as neonatal rodents, exhibit strong synchronization of motoneuron output, greater than that found in adults during constant-force contractions. Moreover, contrary to adults, the performance of fast contractions in neonates is not mediated by a modulation of motoneuron discharge rate. These results are crucially important for our understanding of the development of spinal neural networks. Although it has been largely shown that synchronization of neurons at different levels of the nervous system modulates plasticity, here, we show that in newborns, synchronization is also critical for generating fast contractions, as those involved in primitive reflexes important for protection and survival.
In adults, fast contractions are generated by high motoneuronal discharges (12). The adult muscular apparatus produces fast movements by recruiting the pool of spinal motoneurons in a very short time frame and at discharge rates that are much greater than during controlled low-medium force contractions [200 versus 20 spikes/s; (12, 27, 28)]. The neonatal motoneurons typically showed low discharge rates (~10 to 25 spikes/s) that did not vary with the speed of contraction, despite the very large range of movement speeds we observed (~10 to 400 cm/s). These low discharge rates may depend on immature projections to motoneurons from supraspinal pathways, such as the corticomotoneuronal projections. Therefore, the main strategy that the nervous system uses to increase the rate of force development in neonates is the increase in synchronization of the motoneuron outputs to the muscular apparatus.
The observed synchronization in the neonatal motor circuitries is presumably due to a high common synaptic input from premotor interneurons of CPGs and from afferent feedback. In addition, electrical gap-junction coupling between motor neurons has been suggested to account for synchronization at least in neonatal rodents (2, 3). Here, even during slow neonatal movements (1 to 10 cm/s; see the initial part of movie S1), the synchronization was significantly greater when compared to adults exerting low or moderate isometric forces (30% of maximal). Although it is very challenging to directly compare force levels between adults and neonates, our results suggest that the neonatal spinal cord shows higher values of synchronization than in adults. Motoneuron synchronization was significant at the delta band (0 to 5 Hz), which is the bandwidth that reflects the generation of movement. This bandwidth is also dominant in adult motoneurons, as observed in a large number of studies in different muscles (29, 30). Instead, in neonates, there was little coherence in the beta band (15 to 30 Hz), which is thought to be associated to cortical dynamics (24). This result suggests either motoneuronal and interneuron filtering of CPG inputs or that the CPG dynamics are transmitted to motoneurons mainly within the low-frequency range.
Comparison of the results of synchronization in neonates with previous work on very fast (ballistic) contractions in adults requires careful considerations. We have previously reported the synchronization and discharge rate levels for motoneurons of adults during rapid, maximal isometric contractions (12, 20). The adult spinal cord, during rapid “ballistic” contractions, shows motoneuron discharge rates as high as 200 Hz. Moreover, the human ability to produce force rapidly in adulthood is explained at the individual subject level by the discharge rate of the motoneurons (12). When we estimated the synchronization of motoneurons during ballistic contractions in adults in previous work (20), we observed a value close to the saturation level because of the confounding factor of the extremely high discharge rates. At these levels of discharge rates, computational studies have previously shown that the estimates of motoneuron synchronization by correlation analysis in the time or frequency domain are insensitive to the strength of common synaptic inputs (20). Therefore, the appropriate comparison of neonates’ motoneuron synchronization is with moderate force contractions in adults, where the discharge rates are matched, as we have done in this study. Although the comparison is still challenging because of lack of measures of relative levels of muscle forces, our results on the comparison between adults and neonates clearly indicate diverging strategies to produce fast movements. Neonates do not change discharge rates with changing contraction speed while they change synchronization level (Fig. 4). Conversely, adults rely on modulation of discharge rate to increase the speed of contractions. Therefore, in addition to promoting plasticity, as classically observed, synchronization in neonates may also have the functional role of facilitating fast movements in the absence of a strong corticospinal drive—a mechanism that is important for some neonatal behaviors.
The discharge characteristics of neonates’ motoneurons may also be partly related to intrinsic muscle fiber properties that are different in neonates from that in adults due to the absence of fast fibers at birth (31, 32). Motoneurons innervating slow motor units typically have much longer after-hyperpolarization following the discharge of an impulse, which limits their high discharge rates. Also, mean discharge rates of 8 spikes/s in neonates (Fig. 4H) presumably correspond to fused tetanic contractions in slow neonatal muscles and to unfused tetanus in adult muscles, implying different regimes of muscle force production. Overall, the results point toward a simpler and less flexible control of neonatal muscles with a lack of the “asynchronous” mechanism of motor unit recruitment [which represents a more efficient regime of muscle force production than synchronous activation; (33, 34)], a lack of modulation of motoneuron discharge rate with increasing muscle contraction speed, and a lack of heterogeneity in the muscle fiber composition (only slow fibers in neonates, while an adult muscle contains various types of fast and slow motor units allowing their flexible task-dependent involvement).
In sum, we showed that it is possible to decipher the spinal cord output in human newborns by decoding noninvasive recordings from high-density grids of electrodes placed on limb muscles. We found a high level of motoneuron synchronization, likely the result of strong common excitatory inputs from CPGs and sensory afferents to the motoneuron pools. Moreover, we found that high synchronization in neonates is associated with fast limb movements while discharge rates in neonate do not vary with movement speed. Our results with human newborns are closely reminiscent of previous observations in neonatal rats (2–4). From a translation perspective, with our approach, one might tackle the difficult problem of early diagnosis of motoneuron dysfunction in developmental motor disorders.
The neonates were 2, 3, 10, and 14 days old (two males and two females; gestational age, 37 to 40 weeks; length at birth, 49.25 ± 0.96 cm; weight at birth, 3330 ± 422 g; 5-min Apgar score of ≥ 8). The experiments were carried out at the Casilino Hospital in Rome (Italy), approved by the Research Ethics Committee of Azienda Sanitaria Locale (Local Health Centre) Roma C (protocol 27593, study n. 38.15), and conformed to the Declaration of Helsinki. The consent form was signed by the parents of the children. The adult data were collected at the University of Rome “Foro Italico” (Ethical Committee approval no. 44 680) and conformed to the requirements of the Declaration of Helsinki. After being informed of the purpose and experimental procedures of the study, a written informed consent was signed by all adult participants (12 men; age, 25.1 ± 2.9 years; height, 1.78 ± 0.06 m; weight, 73.3 ± 8.0 kg; International physical activity 2161 ± 1128 metabolic equivalent of task min week^−1^).
The neonates were placed on a massage table in a supine position. After gently rinsing the lower leg with an alcoholic solution for pediatric use, we applied two high-density EMG electrodes to cover most of the tibialis anterior muscle. The electrodes of the high-density grids (Fig. 2) were separated by a 4-mm interelectrode distance, fully gold coated, and of 1 mm in diameter. The EMG signals were acquired in monopolar derivation at 2048 samples/s and filtered at the hardware level with a 10 high-pass and 900-Hz low-pass filter by the hardware system, Quattrocento, OT Bioelettronica (Turin, Italy).
Because the aim of the experiment was to record motor unit activities from the ankle dorsiflexor muscles, we either recorded spontaneous ankle-dorsiflexion (when present) or triggered leg movements by means of light mechanical stimulation of the plantar surface of the foot or other body parts. In the latter case, after the stimulation, the neonate typically showed several bursts of spontaneous activity (movie S1). Neonate movements were filmed with a camera (Realtek integrated camera) at 30 frames/s. We observed that most neonates were developing spontaneous EMG activity with bursts of activations (Figs. 2 to 4 and movie S1).
The experimental apparatus recording ankle dorsiflexion for the adults is shown in Fig. 4A. After a warm-up, the adults performed three maximal voluntary isometric contractions separated by 30 s of rest. The highest peak in force was used to define the threshold for the 50-s long isometric steady-state force contractions (Fig. 4B). The force level for these contractions was either 10 or 30% of maximal force. The participants were instructed to follow with the force biofeedback a trajectory on a screen. The force signal was synchronized with the EMG signals at source by the OT Bioelettronica system Quattrocento. The electrode grids used for the adults had 8-mm interelectrode spacing and 64 electrodes. The signals were also acquired in monopolar mode, and the decomposition procedure is described below and common to both adult and neonate high-density EMG data.
The 128 monopolar EMG signals were acquired at 2048 samples/s and digitally filtered between 10 and 500 Hz and resampled at 10,240 Hz. EMG channels showing poor signal-to-noise ratios were removed by a semiautomatic function that displays the EMG time series for visual inspection after highlighting the channels with an integral of the power spectrum above 3 SDs from the mean (power spectrum averaged across all signals in the band 10 to 500 Hz). For all neonates and under all conditions, 95% of the total channels showed high signal-to-noise ratios.
After manual and visual inspection of the signals, we decomposed the high-density EMG signals via a gradient convolution kernel compensation algorithm (11). This procedure takes advantage of the high dimensionality and uniqueness in the shape of the motor unit action potentials to converge in their unique time-series representation of the discharge timings of the alpha motoneurons (the binary code of movement). The mathematical description of the algorithm has been described in detail previously (11). We also briefly describe here the general steps of decomposition.
Each channel of the EMG signal represents the convolution of the motoneuron discharge timings (sources) by the muscle tissue (muscle unit action potentials). Therefore, the EMG can be mathematically described in a matrix x form (e.g., when recorded with high-density EMG grids) as
where s¯¯(k) = [s1(k), s2(k),…, sn(k)]^T^ represents the n motor unit discharge times that generate the EMG signal (x¯¯), and n¯¯ is the noise for each electrode. The sources (s) are the axonal action potentials reaching the muscle units and can be written as Dirac delta function
where φjr represents the spike times of the jth motor unit. The matrix H¯¯ (l) in Eq. 1 holds the information of the motor unit action potential and has size m × l with lth sample of motor unit action potentials for the n motor units and m channels.
The high spatial sampling given by the 128 electrodes further enhanced by extending the observation numbers (11) allows to recover the sources in a blind way with a gradient descent update rule that maximizes the sparsity (the motor unit waveform; Fig. 3 and fig. S1) to converge in the unique time-space representation of each motor unit discharge timings (the sources, s). Because this process is obtained in a completely blind way, the output of decomposition can be then inspected by spike-triggered averaging.
We designed a custom software for testing the accuracy of decomposition in neonates. The general purpose of this software was to display the motor unit action potential for visual inspection, 2D cross-correlation of the motor unit action potentials, and random perturbations of the motor unit trains (Eq. 2) with Gaussian distributed noise. Figure 3 and fig. S1 show the general overview of this process. After the identification of a portion of the signal that was successfully decomposed on the basis of a high PNR of the extracted motor unit waveform (11), we performed cross-correlation of the action potential to converge in the motor unit waveform. When a motor unit was detected, the shape of the motor unit was cross-correlated with a 2D cross-correlation function. The 2D cross-correlation ranges from 0 to 1, where 1 indicates maximal similarity. We then perturbed the motor discharge timings, with noise ranging from 5 to 20% of the interspike interval variability (5% = ~2 ms; see fig. S1 for the full distribution of values). This approach has been used both in healthy human (35) and feline motor units (36) and is based on the assumption that each motor unit has a unique waveform that changes minimally throughout the contraction. This assumption has been largely confirmed in both simulated and experimental EMG data by comparison of gold-standard intramuscular EMG signals and surface EMG [see (11) for review]. After testing the validity and reliability of each motor unit waveform, we extracted several physiological parameters of the motor unit discharge timings, as discussed below.
Figure 3 and fig. S1 show the steps for obtaining the spatiotemporal representation of the motor unit action potential. After obtaining high reproducibility in the discharge times obtained by blind source separation (motor unit innervation trains with PNRs above 30 dB), we performed spike-triggered averaging on the discharge times across the movements performed by the neonates. This step guarantees the uniqueness of each identified motor unit firing across the tasks.
The uniqueness of each motor unit was assessed with 2D correlation analysis of the motor unit waveforms that were extracted via spike-triggered averaging the high-density EMG in 25-ms windows, centered at the spike instants. Figure 3 and fig. S1 show the unique spatiotemporal representations of the action potentials of the four identified motor units across the task in one neonate. Each motor unit displayed a specific territory and amplitude of the action potential. There was high similarity in the discharge times obtained by decomposition. The 2D cross-correlation value across the same motor unit during the task was 0.8 ± 0.2 while between different motor unit was 0.3 ± 0.2 (see fig. S1). We then injected noise in the pulse trains for increasing time windows and calculated the 2D cross-correlation of the motor unit action potentials across the identified motor unit pool (fig. S1). We consistently observed a decrease in the 2D cross-correlation value with small perturbations of the action potential. The possibility to find motor units by chance, similarly to what has been observed in adult healthy humans (35) and feline motor units (36), is negligible.
The instantaneous discharge rate was computed as the inverse of the interspike interval. The motor unit spikes were stored in binary format with each row corresponding to a motor unit and in the columns either zeros or ones corresponding to the firing instance of a unit. We then extracted the cross-correlogram (estimate of motoneuron synchronization) from this binary signal.
Because the motor units in the pool have distinct synaptic connections and therefore unique inputs, some motor units inevitably show a higher or smaller degree of synchrony that is not representative of the full common synaptic input of the population of motoneurons. It has been previously shown that permutating the cross-correlation function across many motor units improves the estimates of the common input that the population of motoneuron receives. For this purpose, two separate vectors were generated with two different smoothing functions. The smoothing was applied with a convolution with a 25-ms and 1-s Hann window to the motor unit discharge times (Eq. 2, binary signal). The 40-Hz signal represents the full bandwidth of the neural drive to the muscle (17). Conversely, the smoothing of the motor unit spike times with a 1-s Hann window reflects the passive filtering of the musculotendinous unit (18). Previous motor unit correlation studies in neonates adopted smoothing functions of 100 ms (2).
We then divided and summed the number of identified motor units in two groups, generating two independent vectors. This process was iterated with 100 random permutations or corresponding to the product of all integers of the total number of motor units (when less than four motor units were identified; see Fig. 2F for an example). The cross-correlation function was then applied to the two independent cumulative discharge timings. The motoneuron synchronization value corresponded to the peak of the averaged cross-correlogram (maximal value consistently at lag 0). The width of the cross-correlogram was defined as the onset and offset of the peak of the cross-correlogram (Fig. 4G).
Once the n motor unit discharge times were identified, we extracted the synchronization of the motor unit spike times in the signal s as follows
for k = r, if n > 4, r = 100, else r = n!
[N1,N2] = randperm (n).
s (n
(1,2),j) = conv (sum (s (Nx)),window_length) – s¯
[R (k),lag (k)] = cross_correlation (s (n
(1),j), s (n
(2),j))
end
where s is the binary signal corresponding to the motor unit action potentials at time j and window length is Hann window either of 40 or 1000 ms, and s¯ is the average of that signal. The cross-correlation was also extracted in time frames of 500 ms. Note that when the motoneuron discharge timings showed gaps in the cumulative sum of their firing, above 500 ms were removed from the analysis, as described above and in (19). For this reason, the values obtained in the full-time frames may yield a higher value of cross-correlation power. Within the same 500-ms time windows, we also calculated the average movement speed of the neonate, as described below. First, we calculated the average of the movement speed in time windows of 500 ms as the first derivative divided by time for the movement trajectories of the foot. This was computed as M(k)=|x¯+y¯|, for slow movements corresponding to M(i)<M¯ and fast for M(k)≥M¯, where M¯ is the average of the movement speed, x¯ and y¯ are the average change in position of the foot along the x and y coordinates at time k calculated in 500-ms time frames. The synchronization value R¯ (k) was computed for the motoneurons in the same 500-ms time windows and the average R¯ (m) values below the average speed (M¯). In these time frames, we also calculated the average discharge rate of the motor units. This value, which is different from the instantaneous discharge rate, corresponds to the average number of discharges per motor unit per second.
We also estimated the coherent oscillations among the spike trains of the motoneurons for both neonates and adults. For this purpose, we used coherence analysis within the groups of motor units (as described above for the cross-correlation). The coherence analysis is a cross-correlation in the frequency domain. The estimated peaks in the coherence function represent the common motoneuronal fluctuations at the given peak frequency. This type of analysis has been previously used to indirectly infer the common synaptic input received by the motoneurons (17). There are three commonly observed dominant peaks in the adult spinal cord corresponding to low-frequency range (0 to 5 Hz), ~5 to 12 Hz or tremor frequency, and cortical beta band (15 to 30 Hz). The coherence function was computed as follows
The coherence function (C) was computed in 1-s nonoverlapping windows and permutated across two equally sized groups of motor units, as described for the time-domain cross-correlation. The coherence bias was estimated as the maximum value of coherence above the 100-Hz band. The cross-spectrum between x and y [Gxy (f)] and the autospectral densities [Gxx(f) and Gyy(f)] correspond to the unfiltered motor unit data (spike trains). Therefore, the coherence can be significant at any bandwidth within the sampling frequency resolution. Moreover, because the coherence is dependent on the resolution of the signal, the larger the number of spike instants, the greater the accuracy of the measure. For this reason, some neonates did not show high coherent values but high time-domain cross-correlation values, mainly because the frequency content (average of ~1-s time windows) may not show prominent peaks. It is usually advised to use at least 30 s of steady contractions with constant discharge of the motor units (37). Despite these constraints, we consistently observed a peak in the low-frequency band for the neonates (see Results), which is associated to force generation.
For the adult data, we computed the cross-correlation and coherence in a similar way as described above. The cross-correlation function was applied only to the steady-state portion of the contractions (50 s, after removing the ramp-up period) since the motor unit in this time interval is not firing concurrently due to recruitment, and when recruited together, they would inevitably show a high degree of correlation that is due to a shared supraspinal common input sent by the central nervous system. Because the correlation between motoneuron spike instants depends on the number of identified motoneurons, we randomly permutated the number of motor unit identified across the neonates for the cross-correlation function in adults. It is challenging to match the neonate and adult spinal cord output. The cross-correlation and therefore the motoneuron synchronization depend on the discharge rate of motoneurons (16), and there are no methods to compensate for this association. Therefore, we chose steady isometric contractions at 10 and 30% of MVC, since from pilot data, these would relatively match the instantaneous discharge rates of motor units for both adults and neonates (Fig. 4H). We previously described in details the methods and experimental protocols used to acquire adult motor unit data (35).
The video was analyzed with DeepLabCut (38), a Python implementation of a residual neural network that is trained on user-defined labeled data (see Fig. 5). We trained the neural network on the entire duration of the video of the neonates, extracting up to 100 labels per video. The labels were extracted with k-means, and then a residual neural network (resnet50) was trained on the labels. The training was stopped when the loss of neural network plateaued in the learning rate. The network was evaluated by computing the RMS error function (Euclidean distance) and then visually tested on all videos (four videos in total). The error in the estimated position of the markers was very low (1.6 ± 1.2 cm, means ± SD).
Because the number of neonates was low, we performed nonparametric Kruskal-Wallis tests between the motor unit characteristics in adults and neonates, and Bonferroni correction was applied on multiple comparisons. The shape similarity across motor unit was tested via 2D cross-correlation function (MATLAB). All the statistical and software implementation for the analysis of the motor units was performed in MATLAB 2018a (MathWorks Inc.).
Funding: This work was supported by the Italian Ministry of Health (Ricerca corrente, IRCCS Fondazione Santa 486 Lucia), Italian Space Agency (grants I/006/06/0 and ASI-MARS-PRE DC-VUM 2017-006), and Italian University Ministry (PRIN grant 2017CBF8NJ_005) and has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (project NaturalBionicS; grant agreement no. 810346). **Author ** All authors conceived the project and designed the experiments. P.P. provided key populations and comments and edits. A.D.V., F.S.-L., Y.I., and V.M. performed the experiments; A.D.V. analyzed the data; A.D.V. and D.F. drafted the paper; A.D.V., F.L., P.P., and D.F. revised the manuscript. All authors discussed and interpreted the results. **Competing ** The authors declare that they have no competing interests. **Data and materials ** All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
View/request a protocol for this paper from Bio-protocol.