TY - JOUR AU - Nishimura,, Yukio AB - Abstract After spinal cord injury (SCI), the motor-related cortical areas can be a potential substrate for functional recovery in addition to the spinal cord. However, a dynamic description of how motor cortical circuits reorganize after SCI is lacking. Here, we captured the comprehensive dynamics of motor networks across SCI in a nonhuman primate model. Using electrocorticography over the sensorimotor areas in monkeys, we collected broadband neuronal signals during a reaching-and-grasping task at different stages of recovery of dexterous finger movements after a partial SCI at the cervical levels. We identified two distinct network dynamics: grasping-related intrahemispheric interactions from the contralesional premotor cortex (PM) to the contralesional primary motor cortex (M1) in the high-γ band (>70 Hz), and motor-preparation-related interhemispheric interactions from the contralesional to ipsilesional PM in the α and low-β bands (10–15 Hz). The strengths of these networks correlated to the time course of behavioral recovery. The grasping-related network showed enhanced activation immediately after the injury, but gradually returned to normal while the strength of the motor-preparation-related network gradually increased. Our findings suggest a cortical compensatory mechanism after SCI, where two interdependent motor networks redirect activity from the contralesional hemisphere to the other hemisphere to facilitate functional recovery. cortical reorganization, electrocorticography (ECoG), motor cortex, network connectivity, spinal cord injury (SCI) Introduction Spinal cord injury (SCI) induces neuronal plasticity in the central nervous system not only in the spinal cord but also in the brain. This widespread plasticity is the basis for the nervous system to repair or compensate for the damaged sensory and motor pathways. Therefore, understanding of the pathways involved and their responses to the injury is crucial for developing therapeutic strategies to guide regeneration or reorganization of neural circuits for functional restoration (Nishimura and Isa 2012; Isa and Nishimura 2014; Isa 2017). Moreover, exploration of how plasticity after SCI influences sensorimotor behavior could be vital for the development of noninvasive biomarkers for use in prognosis and clinical evaluation of recovery (Ellaway et al. 2011). Extensive and long-term changes in brain activity after SCI have been reported in functional imaging studies in humans (Roelcke et al. 1997; Bruehlmeier et al. 1998; Curt et al. 2002; Jurkiewicz et al. 2007) and in nonhuman primates (Nishimura et al. 2007, 2011). Moreover, the degree of these activity changes, particularly in the cerebral cortex, has been correlated with spinal cord atrophy and also recovery of functionality (Freund et al. 2011; Lundell et al. 2011). However, crucial information on how reorganization of brain networks influences the fast neural dynamics underlying movement execution cannot be accessed by functional imaging techniques, simply due to their low temporal resolution. Recording modalities with higher temporal resolution have also been used in SCI studies, such as the local field potential (LFP) and multiunit activity (MUA) (Aguilar et al. 2010) and electroencephalography (EEG) (Green et al. 1998; Fallani et al. 2007), but only provide a snapshot view of SCI-induced changes in the brain. To acquire a comprehensive picture of how brain networks change immediately after the impact of SCI and how they reorganize toward functional recovery, longitudinal recordings of brain activity both before and after the injury are essential. To achieve this, we used a nonhuman primate model of partial SCI, which assesses functional recovery of hand dexterity in macaques after a unilateral transection of the lateral corticospinal tract at the C4/C5 spinal segments. Compared with human studies, our model distinctively provides a total access to neural activity across the injury. Moreover, this model targets dexterous hand movement, one the most affected functions in tetraplegic patients (Anderson 2004), and provides a controlled experimental platform, where precise lesions with similar extent and location can be induced in different subjects to robustly reproduce nearly-complete functional recovery several weeks after a complete impairment of dexterous finger movements (Sasaki et al. 2004; Nishimura et al. 2007; Higo et al. 2009; Alstermark et al. 2011; Sawada et al. 2015; Tohyama et al. 2017). To acquire the fast neural dynamics underlying long-term and large-scale reorganization in motor networks, we used electrocorticography (ECoG) arrays over the sensorimotor areas of both hemispheres in macaques. Compared with other recording modalities, ECoG offers wider spatial coverage than does LFP or MUA, broader bandwidth than EEG, and superior temporal resolution to fMRI and PET. In addition, ECoG enables long-term stable recording (several months or more) (Chao et al. 2010), which is ideal for monitoring cortical activity throughout recovery. The balanced spatiotemporal resolution and long-term stability enable access to large-scale cortical dynamics during movement over the recovery period, which can last for months. After recording, we quantified the dynamics of local cortical activation and directional corticocortical interactions within the motor network and employed a data-driven analytical approach (Chao and Fujii 2013; Chao et al. 2015) to monitor how network dynamics underlying dexterous finger movements changed across SCI and evolved during its recovery. Materials and Methods Subjects and Task Design Two macaque monkeys (Monkeys M and R) were used in the present study. The data presented in this study include those shown in a previous study (Sawada et al. 2015), but with a completely different focus. The previous study investigated subcortical modulation, while the present study focused on cortical reorganization. The experimental protocols followed the guidelines set forth by the Ministry of Education, Culture, Sports, Science, and Technology (MEXT) of Japan, and were approved by the Institutional Animal Care and Use Committee of the National Institutes of Natural Sciences (Approval No. 14A126). To assess finger dexterity, the monkeys were trained to perform a reaching-and-grasping task (Fig. 1A). The subjects were seated in a monkey chair with their heads fixed and were required to put their left hand on a board positioned at the height of their abdomen. A nontransparent gate in front of the subjects was opened after the monkey held their left hand on the board for at least 2 s, and the subjects could retrieve a piece of food that was presented in the slit. Thirty trials were recorded on each experimental day. A successful precision grip was identified as successful removal of a cube of sweet potato using the pads of the index finger and thumb. A digital video camera (33 frames/s) was used to record the reaching-and-grasping sequence, and infrared sensors placed on the board and window allowed us to determine when the monkey’s hand was in the rest position and when the food piece was being retrieved. The reaching onset was measured as the time when the subject’s left hand was lifted from the board; the grasping onset was measured as the time when the subject picked up the piece of food from the slit. Figure 1. View largeDownload slide Experimental setup and functional recovery of finger dexterity. (A) A reach-and-grasp task to assess finger dexterity. Two macaque monkeys (Monkeys M and R) were trained to perform a reach-and-grasp task. A successful trial was defined as successful removal of a cube of sweet potato by precision grip with the pads of the index finger and thumb. (B) Spinal cord lesion. The left lateral corticospinal tract was transected at the C4/C5 level. The black area represents the lesion site, and the proportion of the lesioned tissue is indicated. The lesion extent was calculated as 100 × (1--α/β), where α and β are areas shown in the inset. (C) Arrangement of the electrocorticography (ECoG) electrodes. Two 15-channel ECoG electrode arrays were implanted, each over a hemisphere. Electrode numbers are labeled. Channel 25 showed strong artifacts and thus was excluded. Different symbols in the contralesional hemisphere indicate the movement evoked by threshold electrical stimulation. The arcuate sulcus and central sulcus are labeled as as, and cs, respectively. Panels (A) and (B) are adapted from (Sawada et al. 2015). Figure 1. View largeDownload slide Experimental setup and functional recovery of finger dexterity. (A) A reach-and-grasp task to assess finger dexterity. Two macaque monkeys (Monkeys M and R) were trained to perform a reach-and-grasp task. A successful trial was defined as successful removal of a cube of sweet potato by precision grip with the pads of the index finger and thumb. (B) Spinal cord lesion. The left lateral corticospinal tract was transected at the C4/C5 level. The black area represents the lesion site, and the proportion of the lesioned tissue is indicated. The lesion extent was calculated as 100 × (1--α/β), where α and β are areas shown in the inset. (C) Arrangement of the electrocorticography (ECoG) electrodes. Two 15-channel ECoG electrode arrays were implanted, each over a hemisphere. Electrode numbers are labeled. Channel 25 showed strong artifacts and thus was excluded. Different symbols in the contralesional hemisphere indicate the movement evoked by threshold electrical stimulation. The arcuate sulcus and central sulcus are labeled as as, and cs, respectively. Panels (A) and (B) are adapted from (Sawada et al. 2015). Spinal Cord Lesion and Histological Analysis All surgeries were performed using sterile techniques under general anesthesia, induced with a combination of ketamine (10 mg/kg, intramuscular injection [i.m.]) and xylazine (1 mg/kg, i.m.) and supplemented with intubation, isoflurane (1–1.5%) inhalation to maintain stable, deep anesthesia throughout the surgery. Heart rate, peripheral capillary oxygen saturation, and end-expiratory carbon dioxide pressure were monitored during the surgery. Ringer’s solution was continually administered through an intravenous (i.v.) drip. Dexamethasone (0.825 mg/body weight) and ampicillin (40 mg/kg) were administered after anesthesia. Lateral corticospinal tract (l-CST) lesions were established as described previously (Sasaki et al. 2004; Nishimura et al. 2007, 2009; Sawada et al. 2015) (Fig. 1B). First, the border between the C4 and C5 segments (C4/C5) was exposed by laminectomy of the C3 and C4 vertebrae, and a transverse opening was made in the dural membrane. A small opening was then made in the pia mater at the lateral convexity of the spinal cord. A horizontal strip was made in a mediolateral direction relative to the lateral funiculus by inserting a minute L-shaped hook that could not be inserted >5-mm deep, which corresponds to the distance from the lateral convexity of the spinal cord to the midline. The dorsal part of the lateral funiculus was transected from the dorsal root entry zone ventrally to the level of the horizontal strip lesion. The lesion was extended ventrally at the most lateral part of the lateral funiculus. The skin and back muscles were sutured with nylon or silk. At the end of all experiments, the monkeys were deeply anaesthetized with an overdose of sodium pentobarbital (50–75 mg/kg, i.v.) and perfused transcardially with 0.1 M phosphate buffered saline (PBS, pH 7.3), followed by 4% paraformaldehyde in 0.1 M PBS (pH 7.3). The spinal cord was immediately removed and immersed in 30% sucrose solution of 0.1 M PBS (pH 7.3), and serial sections (50μm) of the spinal cord (axial) were cut on a freezing microtome. All sections were processed for Nissl-staining with 1% cresyl violet. The photomicrographs of the spinal cord lesion were captured and then reconstructed using Adobe Illustrator. The extent of the SCI lesion in the C4/C5 spinal cord was calculated from spinal sections using the following equation: R = 100 × (1 − α/β), where R is the percentage of the lesion extent, α is the area of white matter remaining in the lateral and ventral funiculi on the lesioned side, and β is the area of total white matter in the lateral and ventral funiculi on the intact side (Sugiyama et al. 2013) (see inset in Fig. 1B). As results, 54% lesion extent was found in both Monkeys M and R (Fig. 1B). ECoG Implant and Recording Cortical oscillatory activities were recorded from the sensorimotor cortex, including the premotor cortex, primary motor cortex, and somatosensory cortex. First, a median linear skin incision was performed on the head, and the skull was exposed over the bilateral frontal cortices. Craniotomies were located around the central sulcus, and the cortex around the central sulcus was exposed bilaterally in the monkeys. Two platinum ECoG arrays, each comprised of 15-channel (5 × 3 grid) electrodes, were placed on the digit, hand, and arm areas of the primary motor cortex (M1), the primary somatosensory cortex (S1), and the premotor cortex (PM) on both hemispheres (Fig. 1C). The electrodes had a diameter of 1 mm and an interelectrode distance of 4 mm center-to-center. Small titanium–steel screws were attached to the skull as anchors. The skull and screws were completely covered with acrylic resin. Two titanium–steel tubes were mounted in parallel over the frontal and occipital lobes for fixation of the head. To estimate the motor areas covered by ECoG electrodes, microstimulation current through a custom made stimulator was delivered to each ECoG electrode on the contralesional hemisphere, and movements in body parts were examined. Each electrode was labeled based on the body parts of which movements could be evoked (Fig. 1C). To verify the position of the ECoG arrays, the brain was removed at the end of all experiments after the transcardial perfusion of 4% paraformaldehyde, photographed, and the electrode locations were reconstructed. The electrode locations were further estimated by comparing the MRI images and the physical dimension of the arrays. ECoG recordings were done on different days before and after the spinal cord lesion. In Monkey M, data were recorded on Days −47, −31, −27, −19, −7, 3, 4, 8, 9, 10, 15, 16, 21, 22, 23, 25, 29, 32, 35, 38, 78, 80, 85, 86, 91, 93, 94, 95, 100, 101, 105, 107, 109, 113, 151, 260, and 274. In Monkey R, data were recorded on Days −27, −23, −19, −13, −3, 10, 11, 16, 18, 22, 23, 25, 28, 31, 37, 38, 42, 44, 53, 63, 67, 71, 73, 84, 85, 118, 126, 128, 156, and 158. Signals were recorded with a Cerebus™ data acquisition system at a sampling rate of 2 000 Hz. ECoG signals were extracted using multichannel amplifiers with a 0.3-Hz high-pass analog filter and with band-pass digital filters in the Neural Signal Processor (0-Hz high-pass and 250-Hz low-pass). Data Analysis Preprocessing The 50-Hz line noise was removed from raw ECoG data using the Matlab Chronux toolbox (Bokil et al. 2010). The data were then downsampled eight times, resulting in a sampling rate of 250 Hz. Channels and trials with abnormal spectra were rejected using an automated algorithm from the EEGLAB library (Delorme et al. 2011), which has been suggested as the most effective method for artifact rejection (Delorme et al. 2007). For both subjects, signals from Channel 25 were artifactual and were excluded from analyses. The ECoG signals from each channel were then aligned using the timing of reaching onset (Time 0; see Fig. 2A), where grasping onset varied but typically occurred 0.3–1 s after the reaching onset (Supplementary Fig. S1). The further analyses of cortical activation and corticocortical interactions were done using the Matlab FieldTrip toolbox (Oostenveld et al. 2010). Figure 2. View largeDownload slide Quantification of the dynamics of cortical activation and corticocortical interactions. (A) A schematic showing raw ECoG signals from a channel (Channel 12), where the trials are presented as a stack. Each trial was aligned by the corresponding reaching onset (Time 0; indicated by the vertical line), while the grasping onset varied across trials (indicated by the vertical dashed line). Actual data from Monkey M are used. (B) A time–frequency representation (TFR) was calculated for each trial. (C) Event-related activity (ERA) was calculated for each trial by dividing each TFR value by the baseline value (baseline period indicated as the horizontal bar), and shown as ratios. (D) ERAs from every three consecutive experimental days (n = 90 trials) were averaged to represent the in-trial activity dynamic at the time of recovery. The vertical dashed line represents the averaged timing of grasping onset (also see Supplementary Fig. S1). (E) A schematic showing raw ECoG signals from two channels (Channels 12 and 7), where trials are presented as a stack. (F) Spectral Granger causality (GC) was measured from every three consecutive experimental days (n = 90 trials); event-related causality (ERC) was calculated by dividing each GC value by the baseline value (baseline period indicated as the horizontal bar). ERCs from two opposite connections, Channel 12–7 and Channel 7–12, are shown. Figure 2. View largeDownload slide Quantification of the dynamics of cortical activation and corticocortical interactions. (A) A schematic showing raw ECoG signals from a channel (Channel 12), where the trials are presented as a stack. Each trial was aligned by the corresponding reaching onset (Time 0; indicated by the vertical line), while the grasping onset varied across trials (indicated by the vertical dashed line). Actual data from Monkey M are used. (B) A time–frequency representation (TFR) was calculated for each trial. (C) Event-related activity (ERA) was calculated for each trial by dividing each TFR value by the baseline value (baseline period indicated as the horizontal bar), and shown as ratios. (D) ERAs from every three consecutive experimental days (n = 90 trials) were averaged to represent the in-trial activity dynamic at the time of recovery. The vertical dashed line represents the averaged timing of grasping onset (also see Supplementary Fig. S1). (E) A schematic showing raw ECoG signals from two channels (Channels 12 and 7), where trials are presented as a stack. (F) Spectral Granger causality (GC) was measured from every three consecutive experimental days (n = 90 trials); event-related causality (ERC) was calculated by dividing each GC value by the baseline value (baseline period indicated as the horizontal bar). ERCs from two opposite connections, Channel 12–7 and Channel 7–12, are shown. Dynamics of Cortical Activation The dynamics of cortical activation in each channel were quantified by the time–frequency representation (TFR) generated by Morlet wavelet transformation, at 111 different center frequencies (10–120 Hz) with the half-length of the Morlet analyzing wavelet set at the coarsest scale of five samples (Fig. 2B). Each TFR represented the in-trial dynamics from a channel, during the time between −1.5 s and 1 s (125 time points), and across the frequency between 10 Hz and 120 Hz (111 frequency bins). To quantify event-related activity (ERA) during reaching-and-grasping, we further divided each TFR value by the baseline value (the mean TFR value at the corresponding frequency during the resting period from −1.5 s to −1 s) (Fig. 2C). To provide a stable representation of in-trial dynamics during the recovery, we then averaged ERAs across trials collected from every three consecutive experimental days (n = 90 trials) (Fig. 2D). The trial averaging was done separately for data before and after the lesion, which resulted in 33 and 26 averaged ERAs for Monkeys M and R, respectively. Dynamics of Corticocortical Interactions The dynamics of corticocortical interactions were quantified by spectral Granger causality (GC), which can represent phase differences between signals from two cortical areas to provide their asymmetric causal dependence (Kamiński et al. 2001; Brovelli et al. 2004). Between signals from two channels (Fig. 2E), in-trial spectral GCs were calculated from every three consecutive experiments (n = 90 trials), where each GC represents a unidirectional connectivity from one channel to another, during times between −1.5 s and 1 s (125 time points), and across frequencies between 10 Hz and 120 Hz (111 frequency bins). Event-related causality (ERC) was further quantified by dividing each GC value by the baseline value (the mean GC value at the corresponding frequency during the resting period from −1.5 s to −1 s) (Fig. 2F). Three preparation steps were performed for the spectral GC calculation (1) Preprocessing: detrending, temporal normalization, and ensemble normalization were performed to achieve local stationarity of the data (Ding et al. 2000). Detrending, which is the subtraction of the best-fitting line from each time series, removes the linear drift in the data. Temporal normalization, which is the subtraction of the mean of each time series and division by the standard deviation, ensures that all variables have equal weights across the trial. These processes were performed on each trial for each channel. Ensemble normalization, which is the pointwise subtraction of the ensemble mean and division by the ensemble standard deviation, dramatically improves the local stationarity of the data (Ding et al. 2000; Bressler and Seth 2011). (2) Window length selection: the length and the step size of the sliding-window for segmentation were set as 250 ms and 50 ms, respectively. (3) Model order selection: the model order, which is related to the length of the signal in the past that is relevant to the current observation, was determined by the Akaike information criterion (AIC) (Akaike 1974). In both subjects, a model order of 10 samples (equivalent to 10 × 4 = 40 ms of history) resulted in minimal AIC and was selected. The selected model order also passed the Kwiatkowski–Phillips–Schmidt–Shin (KPSS) test (Lee and Schmidt 1996), thus maintaining local stationarity. Furthermore, the vector autoregression model was validated by the consistency test (Ding et al. 2000). Parallel Factor Analysis To extract structured information from high-volume datasets of cortical activation or corticocortical interactions, we used PARAFAC, a generalization of principal component analysis for higher order arrays (Harshman and Lundy 1994), to enable the computational extraction of latent structures in functional network dynamics (Chao et al. 2015). PARAFAC was performed using the N-way toolbox (Andersson and Bro 2000), with the non-negativity constraint on all dimensions. The convergence criterion (i.e., the relative change in fit for which the algorithm stops) was set to 1e−6. The initialization method was set to be direct trilinear decomposition, which was considered the most accurate method (Cichocki et al. 2009). To determine the number of structures hidden in the dataset, we performed the core consistency diagnostic (CORCONDIA) to identify the appropriate latent structures where adding other latent structures does not considerably improve the model fit (Bro and Kiers 2003). Results Recovery of Finger Dexterity After Partial Spinal Cord Lesion For both Monkeys M and R, the precision grip was completely impaired immediately after the SCI, but gradually recovered (Fig. 3), as described previously (Sasaki et al. 2004; Nishimura et al. 2007, 2009; Sawada et al. 2015). The success rate reached nearly 100% after 32 days in Monkey M and 44 days in Monkey R (defined as “full recovery” time). Figure 3. View largeDownload slide Recovery of finger dexterity. The time course of precision grip success rate for each subject is shown. Full recovery occurred ~1 month after the lesion. The time scale is used for comparison. Adapted from Sawada et al. (2015). Figure 3. View largeDownload slide Recovery of finger dexterity. The time course of precision grip success rate for each subject is shown. Full recovery occurred ~1 month after the lesion. The time scale is used for comparison. Adapted from Sawada et al. (2015). In-trial Dynamics of Cortical Activation During Recovery To evaluate whether and how cortical activity reorganized during the recovery, we first examined activity in motor networks during reaching-and-grasping, and how this in-trial activity dynamics changed across the SCI and evolved over the course of recovery. In-trial event-related cortical activation at different stages of recovery was quantified by averaged ERA (see Fig. 2D). Examples of averaged ERAs before and after lesion in Monkey M are shown in Figure 4, each representing unique activity dynamics. The activity in Channel 12 (in M1) arose during grasping (grasping-related), primarily in the high-γ frequency band (>70 Hz) (Fig. 4A). In contrast, the activity in Channel 7 (in PM) arose before reaching (motor-preparation-related), primarily in the α and low-β bands (10–15 Hz) and the high-γ band (>70 Hz) (Fig. 4B). The magnitudes of grasping- and preparation-related activity changed over the course of the recovery. For example, the grasping-related activity in Channel 12 was enhanced immediately after the lesion and gradually weakened during the recovery (Fig. 4A). In Channel 7, the preparation-related activity in the α and low β bands was considerably reduced after the lesion except during days 14–29 (Fig. 4B). Significant changes in cortical activity of different frequency bands across different recovery stages are shown in Supplementary Figure S2. Figure 4. View largeDownload slide Example of grasping- and preparation-related activity before and after lesion. (A) In-trial grasping-related activity dynamics in Channel 12 (electrode location marked on the brain map) in Monkey M from 35 days before the lesion (Day −35) to 228 days after the lesion (Day 228). Each panel represents an averaged ERA (ratios of TFR values relative to the corresponding baseline values, as in Fig. 2D), whose title is either shaded in green (before lesion) or red (after lesion). The day of lesion is indicated. (B) In-trial preparation-related activity dynamics in Channel 7. Note that both Channels 12 and 7 were on top of the finger movement areas (see Fig. 1C). Figure 4. View largeDownload slide Example of grasping- and preparation-related activity before and after lesion. (A) In-trial grasping-related activity dynamics in Channel 12 (electrode location marked on the brain map) in Monkey M from 35 days before the lesion (Day −35) to 228 days after the lesion (Day 228). Each panel represents an averaged ERA (ratios of TFR values relative to the corresponding baseline values, as in Fig. 2D), whose title is either shaded in green (before lesion) or red (after lesion). The day of lesion is indicated. (B) In-trial preparation-related activity dynamics in Channel 7. Note that both Channels 12 and 7 were on top of the finger movement areas (see Fig. 1C). To obtain a comprehensive view of the total activity data from all the channels, we deconvolved the dataset to identify its latent structures by performing PARAFAC. We first pooled ERAs from all channels and all experiments to create a broadband library of activity during the recovery. To organize and visualize this dataset, we created a tensor with three dimensions (3D): “Channel” (cortical area), “Time–Frequency” (in-trial dynamics), and “Day” (recovery time course), for the anatomical, dynamic, and functional aspects of the data, respectively. For Monkey M, the dimensionality of the tensor was 29 (channels) by 13 875 (125 time windows and 111 frequency bins) by 33 (days of experiment); for Monkey R, the dimensionality was 29 by 13 875 by 26. We then extracted latent structures from the pooled activity by using PARAFAC. The core consistency diagnostic of the PARAFAC models indicated that the pooled activity could be robustly decomposed into two latent structures, but not more (see the results of the core consistency diagnostic in Supplementary Fig. S3). The two latent activity structures are shown in Figure 5. Each structure represented a unique fingerprint of network anatomy, dynamics, and function, described by its compositions in the three tensor dimensions. For Monkey M, the first structure was grasping-related, due to its occurrence during grasping, which appeared around the contralesional M1 and S1 (Fig. 5A, top), during the grasping and in the high-γ frequency band (>70 Hz) (Fig. 5B, top). Furthermore, this grasping-related activity pattern increased immediately after the lesion, and gradually returned to its pre-lesional state around the time when the finger dexterity fully recovered (Fig. 5C, top). The second structure was preparation-related, due to its occurrence right before movement, which appeared in both hemispheres (Fig. 5A, bottom), before reaching and in both the α and low-β bands (10–15 Hz) and the high-γ band (>70 Hz) (Fig. 5B, bottom), and showed no significant changes across the lesion. Remarkably, similar grasping- and preparation-related activity structures were also found in Monkey R (Fig. 5D–F). Figure 5. View largeDownload slide Two principal activity structures. (A) The anatomical dimension of the two principal activity structures in Monkey M. The size and color of each circle represent the activation level (arbitrary unit) at the corresponding electrode. (B) The dynamic dimension of the two principal structures in Monkey M. The solid vertical line represents the reaching onset (Time 0), and the dashed vertical line represents the average grasping onset. (C) The functional dimension of the two principal structures in Monkey M. The solid vertical line represents the lesion (Day 0), and the dashed vertical line represents the full recovery day. The horizontal line represents the average activation level before the lesion. (D–F) Results from Monkey R. Figure 5. View largeDownload slide Two principal activity structures. (A) The anatomical dimension of the two principal activity structures in Monkey M. The size and color of each circle represent the activation level (arbitrary unit) at the corresponding electrode. (B) The dynamic dimension of the two principal structures in Monkey M. The solid vertical line represents the reaching onset (Time 0), and the dashed vertical line represents the average grasping onset. (C) The functional dimension of the two principal structures in Monkey M. The solid vertical line represents the lesion (Day 0), and the dashed vertical line represents the full recovery day. The horizontal line represents the average activation level before the lesion. (D–F) Results from Monkey R. In both subjects, the activation time course of the grasping-related activity structure (Fig. 5C and F, top) was negatively correlated with the time course of behavioral recovery (the success rate of precision grip in Fig. 3) (correlation coefficient in Monkey M: r = −0.78, P < 0.05, n = 33; in Monkey R: r = −0.89, P < 0.05, n = 26), while the activation time course of the preparation-related activity structure (Fig. 5C and F, bottom) showed no significant correlation (Monkey M: P = 0.15; Monkey R: P = 0.30). These results suggest that the sudden increase of grasping-related activity in the contralesional M1 was directly induced by the lesion, and the progressive reduction of this enhanced activity after the lesion was associated with the recovery of finger dexterity. In-trial Dynamics of Corticocortical Interactions During Recovery We then evaluated whether neuronal interactions between cortical areas, or corticocortical interactions, were also reorganized during the recovery. In-trial event-related corticocortical interactions at different stages of recovery were quantified by ERC (see Fig. 2F). Examples of ERCs before and after lesion are shown in Figure 6. Each represents a unique connectivity dynamic: the intrahemispheric interaction from Channels 7 to 12 arose during grasping (grasping-related), primarily in the high-γ frequency band (>70 Hz), with an additional enhancement of a low frequency component which appeared during the later stage of the recovery (Fig. 6A). In contrast, the interhemispheric interaction from Channels 2 to 16 arose before reaching (motor-preparation-related), primarily in the α and low-β bands (10–15 Hz) (Fig. 6B). The magnitudes of grasping- and preparation-related connectivity changed over the course of the recovery. For example, the high-γ grasping-related connectivity from Channels 7 to 12 was enhanced immediately after the lesion and gradually weakened during the recovery (Fig. 6A), and the preparation-related connectivity from Channels 2 to 16 was enhanced after the lesion (Fig. 6B). Significant changes in corticocortical interactions of different frequency bands across different recovery stages are shown in Supplementary Figure S4. Figure 6. View largeDownload slide Examples of grasping- and preparation-related connectivity before and after lesion. (A) In-trial grasping-related connectivity dynamics from Channel 7 to 12 (an intrahemispheric connection marked on the brain map) in Monkey M from 35 days before the lesion (Day −35) to 228 days after the lesion (Day 228). Each panel represents an ERC (ratios of GC values relative to the corresponding baseline values as in Fig. 2F), whose title is shaded either in green (before lesion) or red (after lesion). The day of lesion is indicated. The day of lesion is indicated. (B) In-trial preparation-related connectivity dynamics from Channel 2 to 16 (an interhemispheric connection). Figure 6. View largeDownload slide Examples of grasping- and preparation-related connectivity before and after lesion. (A) In-trial grasping-related connectivity dynamics from Channel 7 to 12 (an intrahemispheric connection marked on the brain map) in Monkey M from 35 days before the lesion (Day −35) to 228 days after the lesion (Day 228). Each panel represents an ERC (ratios of GC values relative to the corresponding baseline values as in Fig. 2F), whose title is shaded either in green (before lesion) or red (after lesion). The day of lesion is indicated. The day of lesion is indicated. (B) In-trial preparation-related connectivity dynamics from Channel 2 to 16 (an interhemispheric connection). To obtain a more comprehensive view of the total dynamics, we pooled ERCs from all connections and all experiments to create a broadband library of connectivity during the recovery, and a tensor was organized in three dimensions: “Channel-Channel” (corticocortical connection), “Time–Frequency” (in-trial dynamics), and “Day” (recovery time course). For Monkey M, the dimensionality of the tensor was 812 (unidirectional connections) by 13 875 (125 time windows and 111 frequency bins) by 33 (days of experiment); and for Monkey R, the dimensionality was 812 by 13 875 by 26. We then deconvolved the 3D connectivity tensor into multiple components by performing PARAFAC, and two latent structures from the pooled connectivity were observed (see the results of the core consistency diagnostic in Supplementary Fig. S5). The two connectivity structures are shown in Figure 7. For Monkey M, the first structure was grasping-related, due to its occurrence during grasping, which represented intrahemispheric interactions, primarily within the contralesional motor cortex (Fig. 7A, top), during grasping and in the high-γ frequency band (Fig. 7B, top). Furthermore, this grasping-related connectivity pattern increased immediately after the lesion, and gradually returned to its pre-lesional state around the time when the finger dexterity fully recovered (Fig. 7C, top). The second structure was preparation-related, due to its occurrence right before movement, which represented interhemispheric interactions, primarily from the contralesional to ipsilesional motor cortices (Fig. 7A, bottom), before reaching and in the α and low-β frequency bands (Fig. 7B, bottom). Furthermore, these interhemispheric interactions were gradually enhanced after the lesion, and peaked when the finger dexterity fully recovered (Fig. 7C, bottom). Similar grasping- and preparation-related connectivity structures were also found in Monkey R (Fig. 7D–F), although the α- and low-β-band connectivity in the preparation-related connectivity structure continued during the movement in this monkey. Figure 7. View largeDownload slide Two principal connectivity structures. (A) The anatomical dimension of the two principal connectivity structures in Monkey M. Each arrow indicates the directionality of the interaction; its color and width indicate the activation level (arbitrary unit). (B) The dynamic dimension of the two principal connectivity structures in Monkey M. The solid vertical line represents the reaching onset (Time 0), and the dashed vertical line represents the average grasping onset. (C) The functional dimension of the two principal connectivity structures in Monkey M. The solid vertical line represents the lesion (Day 0), and the dashed vertical line represents the full recovery day. The horizontal line represents the average activation level before the lesion. (D–F) Results from Monkey R. Figure 7. View largeDownload slide Two principal connectivity structures. (A) The anatomical dimension of the two principal connectivity structures in Monkey M. Each arrow indicates the directionality of the interaction; its color and width indicate the activation level (arbitrary unit). (B) The dynamic dimension of the two principal connectivity structures in Monkey M. The solid vertical line represents the reaching onset (Time 0), and the dashed vertical line represents the average grasping onset. (C) The functional dimension of the two principal connectivity structures in Monkey M. The solid vertical line represents the lesion (Day 0), and the dashed vertical line represents the full recovery day. The horizontal line represents the average activation level before the lesion. (D–F) Results from Monkey R. To further visualize the connectivity patterns in the structures (Fig. 7A and D), we quantified the causal outflow and causal inflow for each channel, which are the sums of all outgoing interactions from the channel and all incoming interactions to the channel, respectively. Therefore, causal outflows and inflows revealed the source and sink areas of the interactions, respectively. The source and sink areas of the grasping-related connectivity structure indicated that intrahemispheric interactions were primarily from the contralesional PM to the contralesional M1 (Fig. 8A). In both subjects, the sink areas were located on the areas for finger and wrist controls (see Fig. 1C). On the other hand, the source and sink areas of the preparation-related connectivity structure indicated that interhemispheric interactions were primarily from the contralesional PM to the ipsilesional PM (Fig. 8B). Figure 8. View largeDownload slide Source and sink areas of the connectivity structures. (A) The grasping-related connectivity structure represents intrahemispheric PM–M1 networks within the contralesional hemisphere. The source (left) and sink (right) areas for Monkeys M (top) and R (bottom) are shown. The source and sink areas were measured from the top 10% connections shown in Figure 7. Darker colors indicate stronger contribution. (B) The preparation-related connectivity structure represents interhemispheric PM–PM networks from the contralesional to ipsilesional hemisphere. Figure 8. View largeDownload slide Source and sink areas of the connectivity structures. (A) The grasping-related connectivity structure represents intrahemispheric PM–M1 networks within the contralesional hemisphere. The source (left) and sink (right) areas for Monkeys M (top) and R (bottom) are shown. The source and sink areas were measured from the top 10% connections shown in Figure 7. Darker colors indicate stronger contribution. (B) The preparation-related connectivity structure represents interhemispheric PM–PM networks from the contralesional to ipsilesional hemisphere. In both subjects, the activation time course of the grasping-related connectivity structure (Fig. 7C and F, top) was negatively correlated with the time course of the precision grip success rate (Fig. 3) (Monkey M: Pearson’s r = −0.91, P < 0.05, n = 33; Monkey R: r = −0.72, P < 0.05, n = 26, see Fig. 9A), while the activation time course of the preparation-related connectivity structure (Fig. 7C and F, bottom) showed a positive correlation (Monkey M: r = 0.69, P < 0.05; Monkey R: r = 0.49, P < 0.05, see Fig. 9B). This showed that the transient enhancement of grasping-related neuronal interactions in the contralesional motor cortex could be directly linked to the decline in finger dexterity, and interhemispheric interactions could be the driving force for the recovery. This functional correlation between the grasping- and preparation-related connectivity structures was further shown by a negative correlation between their activations after the lesion Monkey M: r = −0.84, P < 0.05; Monkey R: r = −0.52, P < 0.05 (see Fig. 9C). Figure 9. View largeDownload slide Correlation between connectivity structures and behavioral recovery. (A) Correlations between the activation time course of the grasping-related connectivity structure and the time course of the success rate of precision grip in Monkeys M (top) and R (bottom). Each circle represents an experimental day, and the red line indicates the linear fit, where the corresponding Pearson’s correlation coefficient (r) and P-value (P) are shown. (B) Correlations between the activation time course of the preparation-related connectivity structure and the time course of the precision grip success rate. (C) Correlations between the activation time courses of the grasping- and preparation-related connectivity structures. Figure 9. View largeDownload slide Correlation between connectivity structures and behavioral recovery. (A) Correlations between the activation time course of the grasping-related connectivity structure and the time course of the success rate of precision grip in Monkeys M (top) and R (bottom). Each circle represents an experimental day, and the red line indicates the linear fit, where the corresponding Pearson’s correlation coefficient (r) and P-value (P) are shown. (B) Correlations between the activation time course of the preparation-related connectivity structure and the time course of the precision grip success rate. (C) Correlations between the activation time courses of the grasping- and preparation-related connectivity structures. The results from the activity and connectivity analyses both demonstrated that the contralesional motor cortex became hyperactive after SCI and then gradually returned to the pre-lesional state during recovery. On the other hand, a significant increase in interhemispheric connectivity was found after SCI, while no significant change in ipsilesional activity was found (see Fig. 5). This indicated that ipsilesional activity remained at the same level but became more in synchrony with contralesional activity during the recovery. Nonetheless, activity and connectivity analyses were complementary with each other, and collectively demonstrate the existence of two distinct mesoscopic network dynamics in the motor cortex during the functional recovery. The intrahemispheric network in the contralesional motor cortex could reflect enhanced neuronal activation and interactions directly induced by the lesion, and the interhemispheric network could underlie a compensatory mechanism to reshape or counterbalance the hyperactive intrahemispheric network. Discussion In this study, we monitored cortical activation and corticocortical interactions across SCI and during long-term SCI recovery, and revealed a rich landscape of how neuronal interactions in the motor cortex evolve during month-long recovery from SCI in a nonhuman primate model. Our results demonstrated that neuronal plasticity during SCI recovery took place in two distinct motor networks: an intrahemispheric grasping-related network from the contralesional PM to the contralesional M1, and an interhemispheric preparation-related network from the contralesional PM to the ipsilesional PM. From the neuroanatomical perspective, the intrahemispheric network could be the PM−M1 corticocortical pathway, while the interhemispheric network could be the PM−PM commissural pathway, or some other PM−PM pathways mediated by subcortical networks. We theorize that the interhemispheric interactions modulate the ipsilesional motor cortex to facilitate functional recoupling between the contralesional motor cortex and muscles after the injury, and subserve a compensatory mechanism to shape enhanced neuronal activation and interactions in the contralesional motor cortex for functional recovery. Nonhuman Primate Model of Partial SCI Studies on partial SCI model in macaques have shown that the lesion of the lateral corticospinal tracts affects the direct cortico-motoneuronal connection, but nearly-complete functional recovery can be achieved by indirect pathways from the corticospinal tracts to motoneurons (Sasaki et al. 2004; Nishimura et al. 2007, 2009; Higo et al. 2009; Sawada et al. 2015; Tohyama et al. 2017). The advantage of this lesion-recovery model is that it enables a controllable recovery time course and the direct examination of how motor networks continuously evolve across SCI and toward its recovery, which are both difficult to obtain in human studies. With the superior controllability of the recovery time course and accessibility of brain activity, we could identify the high-resolution network dynamics underlying SCI recovery in this study as a stepping stone to exploration of the lesion-dependent variability of recovery dynamics which has been a major issue in human studies. Interplay Between Intrahemispheric and Interhemispheric Networks During SCI recovery, the strengths of the grasping-related intrahemispheric and preparation-related interhemispheric networks changed in an anticorrelated fashion (Fig. 9), suggesting their functional interdependence. In Figure 10, we theorize a step-by-step reorganization of cortical networks during functional recovery based on our unbiased data-driven results. Before the injury (Fig. 10A), the contralesional PM interacts with M1 (in high-γ frequency band), which then sends motor commands to the finger muscle. The functional coupling between M1 and the muscle is primarily through a direct cortico-motoneuronal connection (Mima et al. 2001; Nishimura et al. 2009) and partially via the spinal interneuronal pathways (Kinoshita et al. 2012). Immediately after lesion of the direct pathway in the spinal cord (Fig. 10B), the contralesional motor cortex and the muscles are decoupled. This induces enhanced activation in the contralesional M1 (see the recovery time course of the grasping-related activity structure in Fig. 5) and enhanced interactions from the contralesional PM to M1 (see the grasping-related connectivity structure in Fig. 7). During early recovery (Fig. 10C), the interhemispheric interactions modulate the ipsilateral PM (in α and low-β frequency bands) (see the preparation-related connectivity structure in Fig. 7), which further enhances the connection between the ipsilesional M1 and the muscles and gradually recouples the contralesional motor cortex to the muscles. This transient pathway further shapes the enhanced activation and interactions in the contralesional motor cortex (see the grasping-related activity and connectivity structures in Figs 5 and 7, respectively). Toward full recovery (Fig. 10D), the intensifying interhemispheric network continues to recouple the contralesional motor cortex with the finger muscles and facilitate the indirect pathway in the spinal cord, and reaches its maximal activation when the indirect pathway becomes fully functional (see the preparation-related connectivity structure in Fig. 7). After finger dexterity fully recovers (Fig. 10E), the interhemispheric network gradually fades away (see the preparation-related connectivity structure in Fig. 7). Figure 10. View largeDownload slide Coordination of intrahemispheric and interhemispheric networks during SCI recovery. The theorized step-by-step reorganization of cortical networks before the injury (A), immediately after lesion of the direct pathway in the spinal cord (B), during early recovery (C), toward full recovery (D), and after full recovery (E). The blue pathway represents the grasping-related fast-oscillatory intrahemispheric network; the red pathway represents the preparatory slow-oscillatory interhemispheric network; the green pathway represents the hypothesized network from the ipsilesional M1 to the finger muscles that was modulated by the interhemispheric network. The thickness of the arrows represents the strength of the interactions, and the thickness of the shapes (round rectangle: neurons; oval: muscles) represents the level of the activation. Figure 10. View largeDownload slide Coordination of intrahemispheric and interhemispheric networks during SCI recovery. The theorized step-by-step reorganization of cortical networks before the injury (A), immediately after lesion of the direct pathway in the spinal cord (B), during early recovery (C), toward full recovery (D), and after full recovery (E). The blue pathway represents the grasping-related fast-oscillatory intrahemispheric network; the red pathway represents the preparatory slow-oscillatory interhemispheric network; the green pathway represents the hypothesized network from the ipsilesional M1 to the finger muscles that was modulated by the interhemispheric network. The thickness of the arrows represents the strength of the interactions, and the thickness of the shapes (round rectangle: neurons; oval: muscles) represents the level of the activation. Enhanced Activation in Intrahemispheric Network as a Substrate for Plasticity Transient enhancements of activation and intrahemispheric interactions found in the contralesional motor cortex could reflect hyperexcitable neuronal circuitry (see the grasping-related activity in Figs 4A and 5, also the high-γ band activity in Supplementary Fig. S2A), which has been observed after neural injuries. In a rat model of focal brain damage, hyperexcitability appeared immediately after cerebral infarction and gradually decreased over subsequent months (Buchkremer-Ratzmann et al. 1996). Similar hyperexcitability is also found in patients with stroke, which usually subsides a few months after stroke (Shimizu et al. 2002; Ward 2004). Hyperexcitable neuronal tissue is thought to be a potential environment for cortical reorganization after neural injury (Nudo 2003; Palva et al. 2005; Wieloch and Nikolich 2006), where it could embed various substrates for neural plasticity, such as enhanced long-term potentiation (Hagemann et al. 1998), downregulated type A gamma-aminobutyric acid (GABAA) receptors (Witte et al. 1997), an increased number of N-methyl-d-aspartate (NMDA) receptor binding sites (Mittmann et al. 1998), and promoted axonal sprouting (Witte et al. 2000; Higo et al. 2009). The impact of SCI on cortical neurons could be subtler and less direct than stroke pathology, nonetheless, enhanced activation found in the intrahemispheric network could represent a constructive rewiring of the dense horizontal interconnections that facilitates the re-establishment of lost corticomuscular coupling, rather than a pathological condition that needs to be alleviated for recovery. Information Relay via Interhemispheric Interactions The causal role of the ipsilesional motor cortex in SCI recovery has been demonstrated by a reduced success rate of precision grip during inactivation of ipsilesional M1 or PM (Nishimura et al. 2007). The next question was whether the involvement of the ipsilesional motor cortex is mediated by emergent interhemispheric neuronal interactions from the contralesional motor cortex, and what kind of information is transferred to the ipsilesional hemisphere. Interhemispheric interactions, via the corpus callosum, have been demonstrated to regulate the brain’s ability to establish new neuronal connections, acquire new functions, and compensate for injuries (Dayan and Cohen 2011; Hosp and Luft 2011; Takeuchi and Izumi 2012). To control neural plasticity for functional recovery, reorganization of the interhemispheric interactions found in the present study (Fig. 7 and Supplementary Fig. S4) could play a role in relaying motor commands from the contralesional motor cortex to the finger muscles after the original descending connection is damaged. Consequently, the interhemispheric network could mediate the recoupling of the contralesional motor cortex and the finger muscles to strengthen residual descending pathways, likely through the propriospinal and/or reticulospinal neurons (Alstermark and Isa 2012; Isa and Nishimura 2014; Isa 2017; Tohyama et al. 2017). In nonpathological settings, transferring explicit motor commands through the corpus callosum has been linked to interlimb coordination during bimanual actions (Gooijers et al. 2013), where an efference copy of the planned motor program is thought to be sent to the opposite hemisphere to enable the optimal timing of movements (Geffen et al. 1994; Carson 2005). More abstract forms of motor signals could also be transferred through the corpus callosum; for example, callosal transfers of motor knowledge are found to contribute to enhanced performance after motor learning (Perez et al. 2007; Lee et al. 2010). Our results suggest the possibility that the signals critically involved in the online control of hand movements are conveyed through the callosal fibers. Information Modulation via Interhemispheric Interactions Instead of directly relaying motor signals to the ipsilesional hemisphere, the interhemispheric interactions could subserve a gating function that modulates motor signals within the ipsilesional hemisphere. This view is supported by the observation that activation of the interhemispheric network occurred during motor preparation but not during movement. Furthermore, slow oscillations, particularly in the α band, have been thought to subserve an inhibitory control that gates the information flow in the task-relevant brain networks (Klimesch et al. 2007; Jensen and Mazaheri 2010; Palva and Palva 2011). Based on this view, our finding of enhanced α oscillations after the lesion could imply that the interhemispheric interactions reduced the inhibition that existed before the lesion, and therefore activated neuronal resources in the ipsilesional hemisphere to process motor-related information. The emergence of interhemispheric networks was consistently observed in both subjects, which suggested that a common mechanism was involved for functional recovery. However, it’s important to note that the functions of the ipsilesional motor cortex might be determined by multiple variables, such as the stage of recovery, degree of injury, task complexity, and type of rehabilitation therapy, as in the case of stroke (Hoyer and Celnik 2011; Dodd et al. 2017). Therefore, data from subjects with a variety of spinal cord lesions will be valuable to examine the full spectrum of possible neuronal solutions for functional recovery. Frequency Channels of Intrahemispheric and Interhemispheric Networks Our results showed that the intrahemispheric network transmitted signals predominantly in the high-γ frequency band, while the interhemispheric network used the α and low-β frequency bands (see Fig. 7B and E). High-γ activity has been commonly found during motor controls (Rickert et al. 2005; Ball et al. 2008; Chao et al. 2010), and has been directly linked to local neuronal firing and synchrony (Ray et al. 2008; Siegel et al. 2012). This is consistent with our view that the high-γ interactions localized within the contralesional motor cortex underlies the motor control of grasping. On the other hand, slow oscillations have been associated with the regulation of long-distance inter-areal communication (Von Stein et al. 2000; Palva et al. 2005; Palva and Palva 2007), which agrees with our view that the slow-oscillatory interactions relay or regulate information across hemispheres. Slow oscillations have also been consistently found to correlate with cortical plasticity after neuronal injury and during development. It has been reported that in patients with stroke, motor recovery is related to large-scale changes in α-band functional connectivity, both in the perilesional and distant cortical areas (Dubovik et al. 2012; Westlake et al. 2012). Changes in long-range coherence in the α-band have also been found to accompany developmental trends in children (Srinivasan 1999; Uhlhaas et al. 2010). Moreover, the functional correlation observed between the slow- and fast-oscillatory networks suggested the significance of cross-frequency interactions in functional recovery. Further understandings of the neuronal mechanisms underlying this coordination will be crucial for developing strategies for its control. Time-Dependent Activation in the Motor Cortex Changes in activations in the motor cortex during SCI recovery were demonstrated previously in a PET study (Nishimura et al. 2007), where activations in bilateral M1 were found during the early stage of recovery (1 month after the lesion); activation in contralesional M1 was further enhanced during the late stage of recovery (3–4 months); and activations in bilateral PM were also observed during the late stage. In the present study, we confirmed the activation of the bilateral M1 during the early recovery stage (see Fig. 5 and the high-γ band activity in Supplementary Fig. S2A). However, we also found different sequence of events, despite a similar experimental protocol being used. In particular, we did not find enhanced contralesional M1 activation during late recovery; instead we found that activity and connectivity in contralesional M1 gradually returned to its pre-lesional state (see the grasping-related activity in Fig. 5C and F, and the grasping-related connectivity in Fig. 7C and F). Since the recording modalities differed in these two studies, one explanation for the inconsistent observation is that the motor cortex could adapt its pre-lesional control strategy and return to the pre-lesional level of activation after recovery (as shown in Figs 5 and 7) but a higher level of cerebral blood flow was required to maintain the same level of activation. Dissecting Complex Multidimensional Brain Activity The main challenge to identify concurrent and interdependent cortical processes during SCI recovery is to extract latent structures from multidimensional high-volume data and examine their interdependence without unnecessary restrictions. Here, we show how a comprehensive description of complex brain responses can be achieved automatically using the PARAFAC method of tensor decomposition. PARAFAC is one of several methods to decompose a multidimensional data into a set of modular structures that can describe the data in a more condensed form. PARAFAC has proven to be a powerful analytical tool for EEG (Miwakeichi et al. 2004; Mørup et al. 2006), ECoG (Yanagawa et al. 2013; Chao et al. 2015), and fMRI (Beckmann and Smith 2005), and is ideal for dissecting brain responses that consist of multiple superimposed network dynamics (Chao et al. 2015). A Data-driven Approach for Testable Hypotheses With a data-driven approach (Yanagawa et al. 2013; Chao et al. 2015), we successfully visualized the comprehensive and robust cortical dynamics from the high-volume and high-dimensional data. This unbiased exploration enabled us to identify novel network structures, which led us to a two-part hypothesis: (1) interhemispheric interactions are critical for SCI recovery, and (2) interhemispheric interactions redirect motor signals to facilitate functional recoupling between the brain and muscles. To experimentally test the first part of the hypothesis, interhemispheric interactions should be manipulated either by a physical lesion to irreversibly block the interactions, or by electrical, optogenetic, or chemical genetic techniques to selectively and reversibly control the pathways by using viral vectors (Kinoshita et al. 2012; Tohyama et al. 2017). To test the second part of the hypothesis, we need to simultaneously record brain activity from sensorimotor areas in both hemispheres and the muscle activity from the fingers, and use decoding techniques to examine whether the activity flowing from the contralesional to ipsilesional motor cortices indeed carries motor information about finger muscle control. These hypothesis-driven investigations could be crucial for understanding the circuit mechanism of functional recovery and furthermore, for development of future therapeutic strategies to guide plasticity in neural circuits toward functional restoration. Supplementary Material Supplementary material is available at Cerebral Cortex online. Notes We thank N. Takahashi, M. Togawa, and K. Isa for technical help. The monkeys used in this study were provided by the National Institute of Natural Sciences through the National Bio-resource Project of the Ministry of Education, Culture, Sports, Science, and Technology of Japan (MEXT). This work was supported by grants from the Brain Machine Interface Development and performed under the Strategic Research Program for Brain Sciences from MEXT and the Japan Agency for Medical Research and Development, a Ministry of Education, Culture, Sports, Science, and Technology Grant-in-Aid for Scientific Research (KAKENHI 23680061) to Y.N., a Grant-in-Aid for Scientific Research on Innovative Areas Adaptive Circuit Shift to T.I. (Project no. 26112008). All data are stored at the Department of Developmental Physiology at the National Institutes for Physiological Sciences, Okazaki, Japan. Conflict of Interest: None declared. References Aguilar J , Humanes-Valera D , Alonso-Calviño E , Yague JG , Moxon KA , Oliviero A , Foffani G . 2010 . Spinal cord injury immediately changes the state of the brain . J Neurosci . 30 : 7528 – 7537 . Google Scholar Crossref Search ADS PubMed WorldCat Akaike H . 1974 . A new look at the statistical model identification . IEEE Trans Automat Contr . 19 : 716 – 723 . Google Scholar Crossref Search ADS WorldCat Alstermark B , Isa T . 2012 . Circuits for skilled reaching and grasping . Annu Rev Neurosci . 35 : 559 – 578 . Google Scholar Crossref Search ADS PubMed WorldCat Alstermark B , Pettersson L , Nishimura Y , Yoshino-Saito K , Tsuboi F , Takahashi M , Isa T . 2011 . Motor command for precision grip in the macaque monkey can be mediated by spinal interneurons . J Neurophysiol . 106 : 122 – 126 . Google Scholar Crossref Search ADS PubMed WorldCat Anderson KD . 2004 . Targeting recovery: priorities of the spinal cord-injured population . J Neurotrauma . 21 : 1371 – 1383 . Google Scholar Crossref Search ADS PubMed WorldCat Andersson CA , Bro R . 2000 . The N-way Toolbox for MATLAB . Chemometr Intell Lab Syst . 52 : 1 – 4 . Google Scholar Crossref Search ADS WorldCat Ball T , Demandt E , Mutschler I , Neitzel E , Mehring C , Vogt K , Aertsen A , Schulze-Bonhage A . 2008 . Movement related activity in the high gamma range of the human EEG . Neuroimage . 41 : 302 – 310 . Google Scholar Crossref Search ADS PubMed WorldCat Beckmann CF , Smith SM . 2005 . Tensorial extensions of independent component analysis for multisubject FMRI analysis . Neuroimage . 25 : 294 – 311 . Google Scholar Crossref Search ADS PubMed WorldCat Bokil H , Andrews P , Kulkarni JE , Mehta S , Mitra PP . 2010 . Chronux: a platform for analyzing neural signals . J Neurosci Methods . 192 : 146 – 151 . Google Scholar Crossref Search ADS PubMed WorldCat Bressler SL , Seth AK . 2011 . Wiener–Granger causality: a well established methodology . Neuroimage . 58 : 323 – 329 . Google Scholar Crossref Search ADS PubMed WorldCat Bro R , Kiers HA . 2003 . A new efficient method for determining the number of components in PARAFAC models . J Chemom . 17 : 274 – 286 . Google Scholar Crossref Search ADS WorldCat Brovelli A , Ding M , Ledberg A , Chen Y , Nakamura R , Bressler SL . 2004 . Beta oscillations in a large-scale sensorimotor cortical network: directional influences revealed by Granger causality . Proc Natl Acad Sci USA . 101 : 9849 – 9854 . Google Scholar Crossref Search ADS PubMed WorldCat Bruehlmeier M , Dietz V , Leenders K , Roelcke U , Missimer J , Curt A . 1998 . How does the human brain deal with a spinal cord injury? Eur J Neurosci . 10 : 3918 – 3922 . Google Scholar Crossref Search ADS PubMed WorldCat Buchkremer-Ratzmann I , August M , Hagemann G , Witte OW . 1996 . Electrophysiological transcortical diaschisis after cortical photothrombosis in rat brain . Stroke . 27 : 1105 – 1111 . Google Scholar Crossref Search ADS PubMed WorldCat Carson R . 2005 . Neural pathways mediating bilateral interactions between the upper limbs . Brain Res Rev . 49 : 641 – 662 . Google Scholar Crossref Search ADS PubMed WorldCat Chao ZC , Fujii N . 2013 . Mining spatio-spectro-temporal cortical dynamics: a guideline for offline and online electrocorticographic analyses. In: Ogawa H , Oka K , editors . Advanced methods in neuroethological research. Tokyo: Springer . p. 39 – 55 . Google Preview WorldCat Chao ZC , Nagasaka Y , Fujii N . 2010 . Long-term asynchronous decoding of arm motion using electrocorticographic signals in monkeys . Front Neuroeng . 3 . doi:10.3389/fneng.2010.00003 . WorldCat Chao ZC , Nagasaka Y , Fujii N . 2015 . Cortical network architecture for context processing in primate brain . eLife . 4 : e06121 . Google Scholar Crossref Search ADS WorldCat Cichocki A , Zdunek R , Phan AH , Amari S-i . 2009 . Nonnegative matrix and tensor factorizations: applications to exploratory multi-way data analysis and blind source separation . John Wiley & Sons . Google Preview WorldCat Curt A , Bruehlmeier M , Leenders K , Roelcke U , Dietz V . 2002 . Differential effect of spinal cord injury and functional impairment on human brain activation . J Neurotrauma . 19 : 43 – 51 . Google Scholar Crossref Search ADS PubMed WorldCat Dayan E , Cohen LG . 2011 . Neuroplasticity subserving motor skill learning . Neuron . 72 : 443 – 454 . Google Scholar Crossref Search ADS PubMed WorldCat Delorme A , Mullen T , Kothe C , Acar ZA , Bigdely-Shamlo N , Vankov A , Makeig S . 2011 . EEGLAB, SIFT, NFT, BCILAB, and ERICA: new tools for advanced EEG processing . Comput Intell Neurosci . 2011 . doi: 10.1155/2011/130714 WorldCat Delorme A , Sejnowski T , Makeig S . 2007 . Enhanced detection of artifacts in EEG data using higher-order statistics and independent component analysis . Neuroimage . 34 : 1443 – 1449 . Google Scholar Crossref Search ADS PubMed WorldCat Ding M , Bressler SL , Yang W , Liang H . 2000 . Short-window spectral analysis of cortical event-related potentials by adaptive multivariate autoregressive modeling: data preprocessing, model validation, and variability assessment . Biol Cybern . 83 : 35 – 45 . Google Scholar Crossref Search ADS PubMed WorldCat Dodd KC , Nair VA , Prabhakaran V . 2017 . Role of the contralesional vs. ipsilesional hemisphere in stroke recovery . Front Hum Neurosci . 11 : 469 . Google Scholar Crossref Search ADS PubMed WorldCat Dubovik S , Pignat J-M , Ptak R , Aboulafia T , Allet L , Gillabert N , Magnin C , Albert F , Momjian-Mayor I , Nahum L . 2012 . The behavioral significance of coherent resting-state oscillations after stroke . Neuroimage . 61 : 249 – 257 . Google Scholar Crossref Search ADS PubMed WorldCat Ellaway P , Kuppuswamy A , Balasubramaniam A , Maksimovic R , Gall A , Craggs M , Mathias C , Bacon M , Prochazka A , Kowalczewski J . 2011 . Development of quantitative and sensitive assessments of physiological and functional outcome during recovery from spinal cord injury: a clinical initiative . Brain Res Bull . 84 : 343 – 357 . Google Scholar Crossref Search ADS PubMed WorldCat Fallani FDV , Astolfi L , Cincotti F , Mattia D , Marciani MG , Salinari S , Kurths J , Gao S , Cichocki A , Colosimo A . 2007 . Cortical functional connectivity networks in normal and spinal cord injured patients: evaluation by graph analysis . Hum Brain Mapp . 28 : 1334 – 1346 . Google Scholar Crossref Search ADS PubMed WorldCat Freund P , Weiskopf N , Ward NS , Hutton C , Gall A , Ciccarelli O , Craggs M , Friston K , Thompson AJ . 2011 . Disability, atrophy and cortical reorganization following spinal cord injury . Brain . 134 : 1610 – 1622 . Google Scholar Crossref Search ADS PubMed WorldCat Geffen GM , Jones DL , Geffen LB . 1994 . Interhemispheric control of manual motor activity . Behav Brain Res . 64 : 131 – 140 . Google Scholar Crossref Search ADS PubMed WorldCat Gooijers J , Caeyenberghs K , Sisti HM , Geurts M , Heitger MH , Leemans A , Swinnen SP . 2013 . Diffusion tensor imaging metrics of the corpus callosum in relation to bimanual coordination: effect of task complexity and sensory feedback . Hum Brain Mapp . 34 : 241 – 252 . Google Scholar Crossref Search ADS PubMed WorldCat Green J , Sora E , Bialy Y , Ricamato A , Thatcher R . 1998 . Cortical sensorimotor reorganization after spinal cord injury an electroencephalographic study . Neurology . 50 : 1115 – 1121 . Google Scholar Crossref Search ADS PubMed WorldCat Hagemann G , Redecker C , Neumann-Haefelin T , Freund HJ , Witte OW . 1998 . Increased long‐term potentiation in the surround of experimentally induced focal cortical infarction . Ann Neurol . 44 : 255 – 258 . Google Scholar Crossref Search ADS PubMed WorldCat Harshman RA , Lundy ME . 1994 . PARAFAC: parallel factor analysis . Comput Stat Data Anal . 18 : 39 – 72 . Google Scholar Crossref Search ADS WorldCat Higo N , Nishimura Y , Murata Y , Oishi T , Yoshino-Saito K , Takahashi M , Tsuboi F , Isa T . 2009 . Increased expression of the growth‐associated protein 43 gene in the sensorimotor cortex of the macaque monkey after lesioning the lateral corticospinal tract . J Comp Neurol . 516 : 493 – 506 . Google Scholar Crossref Search ADS PubMed WorldCat Hosp JA , Luft AR . 2011 . Cortical plasticity during motor learning and recovery after ischemic stroke . Neural Plast . 2011 : 871296 . Google Scholar Crossref Search ADS PubMed WorldCat Hoyer EH , Celnik PA . 2011 . Understanding and enhancing motor recovery after stroke using transcranial magnetic stimulation . Restor Neurol Neurosci . 29 : 395 – 409 . Google Scholar PubMed WorldCat Isa T . 2017 . The brain is needed to cure spinal cord injury . Trends Neurosci . 40 : 625 – 636 . Google Scholar Crossref Search ADS PubMed WorldCat Isa T , Nishimura Y . 2014 . Plasticity for recovery after partial spinal cord injury–hierarchical organization . Neurosci Res . 78 : 3 – 8 . Google Scholar Crossref Search ADS PubMed WorldCat Jensen O , Mazaheri A . 2010 . Shaping functional architecture by oscillatory alpha activity: gating by inhibition . Front Hum Neurosci . 4 : 186 . Google Scholar Crossref Search ADS PubMed WorldCat Jurkiewicz MT , Mikulis DJ , McIlroy WE , Fehlings MG , Verrier MC . 2007 . Sensorimotor cortical plasticity during recovery following spinal cord injury: a longitudinal fMRI study . Neurorehabil Neural Repair . 21 : 527 – 538 . Google Scholar Crossref Search ADS PubMed WorldCat Kamiński M , Ding M , Truccolo WA , Bressler SL . 2001 . Evaluating causal relations in neural systems: Granger causality, directed transfer function and statistical assessment of significance . Biol Cybern . 85 : 145 – 157 . Google Scholar Crossref Search ADS PubMed WorldCat Kinoshita M , Matsui R , Kato S , Hasegawa T , Kasahara H , Isa K , Watakabe A , Yamamori T , Nishimura Y , Alstermark B . 2012 . Genetic dissection of the circuit for hand dexterity in primates . Nature . 487 : 235 – 238 . Google Scholar Crossref Search ADS PubMed WorldCat Klimesch W , Sauseng P , Hanslmayr S . 2007 . EEG alpha oscillations: the inhibition–timing hypothesis . Brain Res Rev . 53 : 63 – 88 . Google Scholar Crossref Search ADS PubMed WorldCat Lee M , Hinder MR , Gandevia SC , Carroll TJ . 2010 . The ipsilateral motor cortex contributes to cross‐limb transfer of performance gains after ballistic motor practice . J Physiol . 588 : 201 – 212 . Google Scholar Crossref Search ADS PubMed WorldCat Lee D , Schmidt P . 1996 . On the power of the KPSS test of stationarity against fractionally-integrated alternatives . J Econom . 73 : 285 – 302 . Google Scholar Crossref Search ADS WorldCat Lundell H , Christensen MS , Barthélemy D , Willerslev-Olsen M , Biering-Sørensen F , Nielsen JB . 2011 . Cerebral activation is correlated to regional atrophy of the spinal cord and functional motor disability in spinal cord injured individuals . Neuroimage . 54 : 1254 – 1261 . Google Scholar Crossref Search ADS PubMed WorldCat Mima T , Toma K , Koshy B , Hallett M . 2001 . Coherence between cortical and muscular activities after subcortical stroke . Stroke . 32 : 2597 – 2601 . Google Scholar Crossref Search ADS PubMed WorldCat Mittmann T , Qü M , Zilles K , Luhmann H . 1998 . Long-term cellular dysfunction after focal cerebral ischemia: in vitro analyses . Neuroscience . 85 : 15 – 27 . Google Scholar Crossref Search ADS PubMed WorldCat Miwakeichi F , Martınez-Montes E , Valdés-Sosa PA , Nishiyama N , Mizuhara H , Yamaguchi Y . 2004 . Decomposing EEG data into space–time–frequency components using parallel factor analysis . Neuroimage . 22 : 1035 – 1045 . Google Scholar Crossref Search ADS PubMed WorldCat Mørup M , Hansen LK , Herrmann CS , Parnas J , Arnfred SM . 2006 . Parallel factor analysis as an exploratory tool for wavelet transformed event-related EEG . Neuroimage . 29 : 938 – 947 . Google Scholar Crossref Search ADS PubMed WorldCat Nishimura Y , Isa T . 2012 . Cortical and subcortical compensatory mechanisms after spinal cord injury in monkeys . Exp Neurol . 235 : 152 – 161 . Google Scholar Crossref Search ADS PubMed WorldCat Nishimura Y , Morichika Y , Isa T . 2009 . A subcortical oscillatory network contributes to recovery of hand dexterity after spinal cord injury . Brain . 132 : 709 – 721 . Google Scholar Crossref Search ADS PubMed WorldCat Nishimura Y , Onoe H , Morichika Y , Perfiliev S , Tsukada H , Isa T . 2007 . Time-dependent central compensatory mechanisms of finger dexterity after spinal cord injury . Science . 318 : 1150 – 1155 . Google Scholar Crossref Search ADS PubMed WorldCat Nishimura Y , Onoe H , Onoe K , Morichika Y , Tsukada H , Isa T . 2011 . Neural substrates for the motivational regulation of motor recovery after spinal-cord injury . PLoS One . 6 : e24854 . Google Scholar Crossref Search ADS PubMed WorldCat Nudo R . 2003 . Adaptive plasticity in motor cortex: implications for rehabilitation after brain injury . J Rehabil Med-Suppl . 41 : 7 – 10 . Google Scholar Crossref Search ADS WorldCat Oostenveld R , Fries P , Maris E , Schoffelen J-M . 2010 . FieldTrip: open source software for advanced analysis of MEG, EEG, and invasive electrophysiological data . Comput Intell Neurosci . 2011 : 156869 . Google Scholar PubMed WorldCat Palva S , Palva JM . 2007 . New vistas for α-frequency band oscillations . Trends Neurosci . 30 : 150 – 158 . Google Scholar Crossref Search ADS PubMed WorldCat Palva S , Palva JM . 2011 . Functional roles of alpha-band phase synchronization in local and large-scale cortical networks . Front Psychol . 2 : 204 . Google Scholar Crossref Search ADS PubMed WorldCat Palva JM , Palva S , Kaila K . 2005 . Phase synchrony among neuronal oscillations in the human cortex . J Neurosci . 25 : 3962 – 3972 . Google Scholar Crossref Search ADS PubMed WorldCat Perez MA , Wise SP , Willingham DT , Cohen LG . 2007 . Neurophysiological mechanisms involved in transfer of procedural knowledge . J Neurosci . 27 : 1045 – 1053 . Google Scholar Crossref Search ADS PubMed WorldCat Ray S , Crone NE , Niebur E , Franaszczuk PJ , Hsiao SS . 2008 . Neural correlates of high-gamma oscillations (60–200 Hz) in macaque local field potentials and their potential implications in electrocorticography . J Neurosci . 28 : 11526 – 11536 . Google Scholar Crossref Search ADS PubMed WorldCat Rickert J , de Oliveira SC , Vaadia E , Aertsen A , Rotter S , Mehring C . 2005 . Encoding of movement direction in different frequency ranges of motor cortical local field potentials . J Neurosci . 25 : 8815 – 8824 . Google Scholar Crossref Search ADS PubMed WorldCat Roelcke U , Curt A , Otte A , Missimer J , Maguire R , Dietz V , Leenders K . 1997 . Influence of spinal cord injury on cerebral sensorimotor systems: a PET study . J Neurol Neurosurg Psychiatry . 62 : 61 – 65 . Google Scholar Crossref Search ADS PubMed WorldCat Sasaki S , Isa T , Pettersson L-G , Alstermark B , Naito K , Yoshimura K , Seki K , Ohki Y . 2004 . Dexterous finger movements in primate without monosynaptic corticomotoneuronal excitation . J Neurophysiol . 92 : 3142 – 3147 . Google Scholar Crossref Search ADS PubMed WorldCat Sawada M , Kato K , Kunieda T , Mikuni N , Miyamoto S , Onoe H , Isa T , Nishimura Y . 2015 . Function of the nucleus accumbens in motor control during recovery after spinal cord injury . Science . 350 : 98 – 101 . Google Scholar Crossref Search ADS PubMed WorldCat Shimizu T , Hosaki A , Hino T , Sato M , Komori T , Hirai S , Rossini PM . 2002 . Motor cortical disinhibition in the unaffected hemisphere after unilateral cortical stroke . Brain . 125 : 1896 – 1907 . Google Scholar Crossref Search ADS PubMed WorldCat Siegel M , Donner TH , Engel AK . 2012 . Spectral fingerprints of large-scale neuronal interactions . Nat Rev Neurosci . 13 : 121 – 134 . Google Scholar Crossref Search ADS PubMed WorldCat Srinivasan R . 1999 . Spatial structure of the human alpha rhythm: global correlation in adults and local correlation in children . Clin Neurophysiol . 110 : 1351 – 1362 . Google Scholar Crossref Search ADS PubMed WorldCat Sugiyama Y , Higo N , Yoshino-Saito K , Murata Y , Nishimura Y , Oishi T , Isa T . 2013 . Effects of early versus late rehabilitative training on manual dexterity after corticospinal tract lesion in macaque monkeys . J Neurophysiol . 109 : 2853 – 2865 . Google Scholar Crossref Search ADS PubMed WorldCat Takeuchi N , Izumi S-I . 2012 . Maladaptive plasticity for motor recovery after stroke: mechanisms and approaches . Neural Plast . 2012 : 359728 . Google Scholar PubMed WorldCat Tohyama T , Kinoshita M , Kobayashi K , Isa K , Watanabe D , Kobayashi K , Liu M , Isa T . 2017 . Contribution of propriospinal neurons to recovery of hand dexterity after corticospinal tract lesions in monkeys . Proc Natl Acad Sci . 114 : 604 – 609 . Google Scholar Crossref Search ADS PubMed WorldCat Uhlhaas PJ , Roux F , Rodriguez E , Rotarska-Jagiela A , Singer W . 2010 . Neural synchrony and the development of cortical networks . Trends Cogn Sci . 14 : 72 – 80 . Google Scholar Crossref Search ADS PubMed WorldCat Von Stein A , Chiang C , König P . 2000 . Top-down processing mediated by interareal synchronization . Proc Natl Acad Sci . 97 : 14748 – 14753 . Google Scholar Crossref Search ADS PubMed WorldCat Ward NS . 2004 . Functional reorganization of the cerebral motor system after stroke . Curr Opin Neurol . 17 : 725 – 730 . Google Scholar Crossref Search ADS PubMed WorldCat Westlake KP , Hinkley LB , Bucci M , Guggisberg AG , Findlay AM , Henry RG , Nagarajan SS , Byl N . 2012 . Resting state alpha-band functional connectivity and recovery after stroke . Exp Neurol . 237 : 160 – 169 . Google Scholar Crossref Search ADS PubMed WorldCat Wieloch T , Nikolich K . 2006 . Mechanisms of neural plasticity following brain injury . Curr Opin Neurobiol . 16 : 258 – 264 . Google Scholar Crossref Search ADS PubMed WorldCat Witte OW , Bidmon H-J , Schiene K , Redecker C , Hagemann G . 2000 . Functional differentiation of multiple perilesional zones after focal cerebral ischemia . J Cereb Blood Flow Metab . 20 : 1149 – 1165 . Google Scholar Crossref Search ADS PubMed WorldCat Witte OW , Buchkremer-Ratzmann I , Schiene K , Neumann-Haefelin T , Hagemann G , Kraemer M , Zilles K , Freund H-J , Payne B , Lomber S . 1997 . Lesion-induced network plasticity in remote brain areas. Authors’reply . Trends Neurosci . 20 : 348 – 349 . Google Scholar Crossref Search ADS PubMed WorldCat Yanagawa T , Chao ZC , Hasegawa N , Fujii N . 2013 . Large-scale information flow in conscious and unconscious states: an ECoG study in monkeys . PLoS One . 8 : e80845 . Google Scholar Crossref Search ADS PubMed WorldCat © The Author(s) 2018. Published by Oxford University Press. All rights reserved. For Permissions, please e-mail: journals.permissions@oup.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model) TI - Dynamic Reorganization of Motor Networks During Recovery from Partial Spinal Cord Injury in Monkeys JF - Cerebral Cortex DO - 10.1093/cercor/bhy172 DA - 2019-07-05 UR - https://www.deepdyve.com/lp/oxford-university-press/dynamic-reorganization-of-motor-networks-during-recovery-from-partial-KLfacuePDh SP - 3059 VL - 29 IS - 7 DP - DeepDyve ER -