Time-Varying Functional Connectivity of Rat Brain during Bipedal Walking on Unexpected Terrain

The cerebral cortex plays an important role in human and other animal adaptation to unpredictable terrain changes, but little was known about the functional network among the cortical areas during this process. To address the question, we trained 6 rats with blocked vision to walk bipedally on a treadmill with a random uneven area. Whole-brain electroencephalography signals were recorded by 32-channel implanted electrodes. Afterward, we scan the signals from all rats using time windows and quantify the functional connectivity within each window using the phase-lag index. Finally, machine learning algorithms were used to verify the possibility of dynamic network analysis in detecting the locomotion state of rats. We found that the functional connectivity level was higher in the preparation phase compared to the walking phase. In addition, the cortex pays more attention to the control of hind limbs with higher requirements for muscle activity. The level of functional connectivity was lower where the terrain ahead can be predicted. Functional connectivity bursts after the rat accidentally made contact with uneven terrain, while in subsequent movement, it was significantly lower than normal walking. In addition, the classification results show that using the phase-lag index of multiple gait phases as a feature can effectively detect the locomotion states of rat during walking. These results highlight the role of the cortex in the adaptation of animals to unexpected terrain and may help advance motor control studies and the design of neuroprostheses.


Introduction
In the natural environment, animals often encounter various complex terrains when walking. Affected by environmental factors, such as a low-light environment, animals may not be able to visually detect changes in terrain when walking, making animals need to adjust their movements to adapt to unexpected terrain changes. Compared with normal walking, this un predictable movement may cause considerable differences in the motor control mechanism. A previous study showed that humans take more cautious strategies when walking on unpredictable uneven terrain, such as increasing step length and reducing step width [1]. In addition, when subjects accidentally touched an obstacle during walking, they would choose different coping strategies according to the movement state of their limbs. For example, when a leg touches an obstacle in the process of swinging, the swinging leg will be immediately lifted up and over the obstacle in the early stage and be placed in front of the obstacle in the later stage [2]. It is worth mentioning that adapting to unexpected terrain seems to be not only an instinctive reaction, and its process was also influenced by previous experience [3,4]. With the development of neuroscience, more and more studies have shown that the brain plays an im portant role in animals' responses to unexpected disturbances. A recent study showed that the spectral power of the premotor cortex and supplementary motor area increased significantly when humans encountered unexpected obstacles while walking or running on a treadmill [5]. When subjects unexpectedly lost their balance during walking, the theta-band energy in multiple brain regions, such as the anterior cingulate and parietal, increased significantly [6]. When the sudden acceleration of the support plane disturbs the posture stability of the subjects, the activities of some cerebral cortex (such as the anterior prefrontal cortex) change differently according to the coping strategies taken by the subjects [7]. Some animal experiments have shown that when the standing platform suddenly tilts, the motor cortex become more active [8,9]. However, these studies mainly focused on analyzing the activation of each brain region and did not analyze the interaction between brain regions. The study of the interaction between cortices can help us better understand the role of the brain in the process of animal adaptation to unexpected terrain.
Brain connectivity estimation aims to describe the interaction between the cortex, including structural connectivity, functional connectivity, and effective connectivity [10]. Structural connectivity refers to the anatomical connection between brain regions [11]; functional connectivity was defined as the temporal correlation between spatially remote neurophysiologic events [12]; and effective connectivity describes the direct or indirect influence of one nervous system on another, which can be derived from functional connectivity [13]. Functional connectivity was considered the most central and challenging of the 3 concepts of brain connectivity and has been widely used in neuroscience research [14]. Currently, most network-based studies use the data of an entire scan session to build a static brain function network, which can obtain the brain connectivity characteristics of the whole process. However, walking was a cycle process consisting of stance and swing phases, and the activity of the cerebral cortex couples with gait phases [15,16]. Therefore, the static function network may not be able to capture changes in brain connectivity with gait phases. Some studies analyze the time-varying properties of networks during cognitive tasks and motor imagery by scanning the data using sliding time windows and estimating the functional connectivity within the window [17,18]. However, these studies did not apply this dynamic analysis method to the actual walking process. Considering the continuity of the walking process, analyzing the network properties on a time scale may be more helpful in understanding the cortical information interaction process during walking.
As mentioned above, our current understanding of the information interaction among brain regions in animal adaptation to unexpected terrain was limited. This study aimed to investigate animals' dynamic brain network properties when encountering unexpected terrain while walking. To address this question, we used rats as the model and built a posturerestrained platform based on a treadmill to help the rat walk in a bipedal posture. The rats' vision was blocked to remove all visual cues on the upcoming change in the terrain, and an uneven area was randomly placed on the treadmill belt to simulate the unexpected complex terrain. An implantable multimodal recording platform equipped with high-speed cameras was used to simultaneously record electroencephalography (EEG) signals and behavioral data in rats. Whole-brain EEG signals from 6 rats were scanned by time windows, and functional connectivity within the window was quantified by phase-lag index (PLI). Furthermore, the global and local metrics of the functional network were calculated on the basis of PLI. Finally, support vector machine (SVM), naive Bayes (NB) classifier, and K-nearest neighbor (KNN) classifier were used to classify rat motion states based on the PLI matrix, global metrics, and local metrics. We hypothesize that dynamic functional connectivity of the rat brain varied with gait phase and terrain conditions. Based on the differences in the time-varying function connectivity, the locomotion state of the rat can be detected.

Animals
The experiments were conducted on 6 male Sprague Dawley rats (aged 8 to 10 weeks and weighing 200 to 250 g). Each rat was individually maintained in a ventilated animal facility with 12-h light/12-h dark cycles. The temperature (25 ± 2 °C), humidity, and ventilation conditions in the cage were appropriate, and food and water were provided ad lib. Rats were trained to walk in bipedal posture after 1 week of dietary control, and food rewards will be given on the basis of the performance of the rats during training. Each rat was trained for about half an hour per day. After 2 weeks of training, a home-designed 32-channel flexible electrode array (Fig. 1A) was surgically fixed on the surface of the rat skull (see [19] for surgical details). After 1 week of recovery, the rats continued to train for 2 weeks and began the formal experiment. All animal procedures were approved by the Institutional Animal Care and Use Committee of the Beijing Institute of Technology. Figure 1D shows the platform used in this study. Rats walked on a moving treadmill belt (6.5 to 7.5 cm/s) in a bipedal posture with the help of a suspension device, and a black tube was used to block all frontal visual cues on the upcoming terrain. An 80-fps camera was placed on each side of the treadmill to record images of the rats. First, the rat stands still on the treadmill for 10 s, then the treadmill starts, and the rat walks on a flat terrain for 20 s. Next, the rat walks for another 30 s, during which we randomly choose to place or not to place an uneven area (the size of the area was 11 × 11 cm, and hard hemispherical protrusions with a diameter of 10 mm were randomly distributed in it) on the treadmill belt in front of the rat. After that, it was regarded as the end of a trial, and the uneven area was removed. In this process, the placement time of uneven area was also random, and rats were able to pass through the uneven area and completely transition to the flat terrain. It was worth noting that uneven areas did not support continuous walking of rats (rats can leave the uneven terrain in the next step after a complete transition to the uneven terrain). There was a 30-s rest time between each trail, and each rat can complete 50 to 60 trails (about 25 to 30 trials with uneven area) on average every day.

Experimental procedures
Three types of locomotion tasks (Fig. 1E) were performed by the rats in the experiment: (a) the rat walked from the flat terrain to the flat terrain (FF); (b) the rat walked from the flat terrain to the uneven terrain (FU); and (c) the rat walked from the uneven terrain to the flat terrain (UF). In addition, 4 events were defined on the basis of the rat's gait: right hind paw coming off the ground (RO), right hind paw contacting the ground (RC), left hind paw coming off the ground (LO), and left hind paw contacting the ground (LC). Benefiting from the slow speed of the treadmill belt, the rats can stand on both hind paws for a period after completing a continuous stride (e.g., RO→RC→LO→LC) to prepare for the next continuous stride. Therefore, we defined the period before the rat stride continuously as the preparation phase and the process of striding continuously as the walking phase. In this study, the right hind paw of the rat was selected as the leading paw, so the preparation phase was from event LC to RO, the walking phase was from event RO to event LC, and the complete gait cycle was defined as LC→RO→RC→LO→LC. We used Simi Motion v9.2.2 (Simi Inc., Unterschlessheim, Germany) kinematics analysis software to extract gait events from rat walking videos to segment different gait phases and manually mark them in EEG signals. By averaging the time interval of each gait phase for all rats under all task conditions (1,846 complete gait cycles in total, each rat contributed a similar number of gait cycles under each condition), the duration of the preparation stage was 747 ± 38 ms (mean ± SD, LC→RO), and the time interval of each gait event The putative electrode positions on the rat's brain surface, as determined by Brainstorm3. The colors of the electrodes correspond to the different brain regions (blue, somatomotor area; red, somatosensory area; yellow, retrosplenial area; green, visual area). (D) Platform used in this study. Rats walked on a treadmill in a bipedal posture with the help of a suspension device, and a black tube was used to block all frontal visual cues on the upcoming terrain. A self-made shell was installed on the rat skull to protect and fix the interface of the electrode array. A removable uneven area was randomly placed on the treadmill belt. (E) Task conditions in this study.
in the walking stage was: 238 ± 15 ms (mean ± SD, RO→RC), 245 ± 18 ms (mean ± SD, RC→LO), 258 ± 16 ms (mean ± SD, LO→LC). For the convenience of comparative analysis, we define the duration of the preparation phase as 750 ms, and the interval between gait events in the walking phase was uniformly defined as 250 ms.

Data acquisition
EEG signals were recorded by a home-designed 32-channel flexible electrode array and Plexon OmniPlex system (Plexon Inc., Dallas, TX, USA). Figure 1A and B shows the actual electrode image and the position of the electrode on the rat skull. The electrode array contains 34 channels, including 32 working channels, a reference channel (REF), and a ground channel (GND). It was worth noting that the signals from the REF and GND channels were excluded, and only the signals from the 32 working channels were used for subsequent analysis. The electrode array was connected to the Digital Headstage processor (Plexon Inc., Dallas, TX, USA) through a cable, and the EEG signal was sampled at 2,000 Hz for storage.

EEG data analysis
EEG data preprocessing Figure 2 shows the processing procedure for EEG data. First, the raw EEG signals from the 32-channel electrode array were notch filtered at 50 Hz to reduce line noise. Our previous study showed that cortical activity in rats was concentrated in the 3-to 50-Hz frequency band when encountering unexpected terrain [19], so the signal was filtered at 3 to 50 Hz. Next, The EEG data were re-referenced to the common average and down-sampled to 500 Hz. Subsequently, the EEG data were decomposed into independent component (IC) sources by an independent component analysis algorithm. We removed those non-brain IC sources by visually inspecting each IC scalp projection and power spectrum to obtain clean EEG data. After the above processing, we extracted epoch from the EEG data of each rat according to the mean hind limb kinematics data of all rats (each epoch started 1 s before the event RC and the total duration was 2.5 s, including 2 complete preparation phases and 1 complete walking phase). After that, epochs from all rats were divided into 3 groups according to different task conditions (FF: 685 epochs, FU: 586 epochs, and UF: 575 epochs). Each rat contributed a similar number of epochs under each condition). The above EEG data analysis was performed using custom scripts written in MATLAB v2018b (The MathWorks, Inc., Nedick, MA, USA) containing functions from EEGLAB [20].

Functional connectivity computation
Each electrode node's EEG signal reflects the neurons' activity in a specific range of brain areas below the node [21]. The functional connectivity can be used to assess the functional relationship between brain regions, quantified with measures of statistical dependencies [22]. Phase synchronization was an effective measurement of functional connectivity. PLI (a phase synchronization-based functional connectivity index) was widely used in EEG signal correlation analysis due to its robustness to the presence of common sources [23]. Therefore, the PLI between the signals of each electrode node was calculated to evaluate the functional relationship among different brain areas. For 2 given EEG signals x, y, PLI values can be calculated by the following equations: The PLI value between 0 and 1: 0≤ ≤1, where 0 indicates no coupling or coupling with a phase difference cantered around 0 mod, and 1 indicates completed phase locking at a value of Δϕ xy (t n ) different from 0 mod. In Eq. 1, | X | represents the absolute value of X, where sign( ) was a signum function. If X > 0, sign(X) = 1; if X < 0, sign(X) = −1; if X = 0, sign(X) = 0. In Eq. 1, Δϕ xy (t n ) denotes the instantaneous phase difference between signals x and y, which can be computed by Eq. 2. where x H (t n ) and y H (t n ) were the Hilbert transforms of x(t n ) and y(t n ), respectively. In this study, we based on the duration of the epoch and gait phase to select a 250-ms long time window to scan the EEG signal. At the same time, there was a 50% overlap between the windows for a more continuous change in connectivity. The HERMES toolbox [23] was used to calculate the signals' PLI [24] matrix within each time window. Next, we got where N was the number of epochs, T was the number of time windows, and CH was the number of electrode nodes. It is worth noting that the preparation phase was a gait phase with a long duration (from −1,000 to −250 ms, i.e, 750 ms before event RO), and using a 250-ms time window (50% overlap) scan resulted in 5 PLI matrices. To facilitate comparative analysis with other gait phases, we averaged the first, third, and fifth PLI matrices to obtain the mean PLI matrix (1,846 × 1 × 32 × 32) of the whole preparation phase.

Functional network construction
We constructed the functional network of the rat brain on the spatial and temporal scales by using the spatial position of each electrode point on the rat skull and the time-varying PLI matrix. In order to build networks of different gait phases, we extracted the matrix belonging to the walking phase from the PLI matrix calculated before and merged them to get a final size of 1,846 × 4 × 32 × 32 PLI matrix. Afterward, we averaged the PLI matrices within each gait phase from the same conditions, resulting in the averaged PLI matrices under the 3 experimental conditions (matrix size was 3 × 4 × 32 × 32, where 3 represents three task conditions, 4 represents four gait phases, and 32 represents different electrode channels). After obtaining the average adjacency matrix for each gait phase under different task conditions, we arranged all electrode points by spatial location (Fig. 1C). Then, by querying the PLI matrix, it was determined whether each electrode point was connected to other electrode points and the weight of the connection. We connect the connected electrode points with lines and use the thickness and color of the lines to indicate the strength of the connection. Additionally, the size of each electrode point was used to indicate the number of connections the electrode point has. It was worth mentioning that the edge threshold of the network was set to 50% of the strongest connections (i.e., only keep edges that are more than 50% of the strongest connections) [25]. Finally, the BrainNet Viewer toolbox [26] was used to visualize the brain functional network.

Network metrics analysis
In this study, the Brain Connectivity Toolbox (http://www.nitrc. org/projects/bct/) was used to calculate network metrics. First, the mean global and local metrics of the functional network of all rats were calculated to analyze the changes in the functional connectivity of the rat brain over time under different task conditions. For the global metrics, we calculate the global efficiency, transitivity, and characteristic path length of the whole brain function network in each time window and average the global metrics belonging to the same condition. In addition, we also calculate the global metrics belonging to each gait phase according to the PLI matrix of each gait phase. In particular, the global metrics of the preparation phase were calculated on the basis of the average PLI matrix of the preparation phase (LC→RO). These global metrics measure the brain's ability to combine information from different brain regions rapidly. Global efficiency was an index to measure the efficiency of distant information transmission in the network. It was defined as the inverse of the average characteristic path length between all nodes in the network [27]. The characteristic path length was the average shortest path length between all node pairs in the network, indicating how easily information can be transferred across the network [28]. The transitivity of the network reflects the prevalence of clustered connectivity around individual nodes [29]. For local metrics, the local efficiency [30], node strength [31], clustering coefficient [32], and eigenvector centrality [33] of each node in each gait phase were calculated. Similar to global metrics, the local metrics of the preparation phase were also calculated on the basis of the average PLI matrix of the preparation phase (LC→RO). The local efficiency was an index to measure the information transfer efficiency of each node in the network. The strength of the relationship between the individual node and other nodes in the network was called node strength, which can be computed as the sum of the connectivity weights of the edges attached to each node. The clustering coefficient was the fraction of a node's neighbors that were neighbors of each other, and eigenvector centrality was a measure of the node's influence on a network. Then, we averaged the local metrics of all nodes belonging to the same brain area and the same condition to obtained the average local x H t n x t n − arctan y H t n y t n metrics of each brain area. It is worth mentioning that the brain region to which each electrode point belongs (Fig. 1C) was determined according to the stereotaxic coordinates of the rat brain [34], the size of the electrode (Fig. 1A), and the installation position of the electrode on the rat skull (Fig. 1B). Finally, IBM SPSS Statistics for Windows, version 27.0 (IBMCorp., Armonk, NY, USA) was used to perform a 1-way analysis of variance for the significance (P < 0.05) of the network metrics across different task conditions.

Machine learning algorithms for the classifying locomotion states
This study used 3 classifiers, SVM, NB, and KNN, for locomotion state detection. These classification algorithms were all implemented by the Scikit-Learn toolbox based on python v3.8 [35]. The basic principle of SVM was to find a hyperplane in the sample space to separate the samples of different classes and to maximize the minimum distance between the classified points and this plane [36]. The NB classifier was a statistical classifier based on Bayes' theorem, which can predict class membership probabilities [37]. The basic idea of KNN was to judge which category the sample belongs to according to the category of the K points closest to the sample [38].
In order to compare the difference between time-varying network analysis method and other analysis methods in detecting rat locomotion state, we used raw EEG signal, whole-time PLI matrix, and time-varying PLI matrix in a complete gait cycle (i.e., LC→RO→RC→LO→LC) as features to classify rat locomotion state. The time length of the raw EEG signal was 1.5 s (i.e., −1,000 to 500 ms). Since the sampling rate was 500 Hz, the matrix of the EEG signal is N × 32 × 750 (FF: N = 685, FU: N = 586, UF: N = 575, input feature: 32 × 750 = 24,000). The whole-time PLI matrix was calculated from the EEG signal within −1,000 to 500 ms, and the size was N × 32 × 32 (input feature dimension: 32 × 32 = 1,024). For the time-varying PLI, a matrix of size N × 4 × 32 × 32 (4 represents four gait phases, i.e., preparation phase: LC→RO; walking phase: RO→RC, RC→LO, and LO→LC) was obtained, and the input feature dimension was 4,096 (4 × 32 × 32 = 4,096). Then, we use t test to analyze the significance of the classification results of each classifier. Furthermore, we use the global and local metrics of the time-varying network for all gait phases (4 phases in total) as features to classify the locomotion state of rats. For the global metrics (global efficiency, transitivity, and characteristic path length), a matrix of size N × 4 × 1 × 1 was obtained, and the input feature dimension was 4. Moreover, for the local metrics (local efficiency, node strength, clustering coefficient, and eigenvector centrality), we get a matrix of size N × 4 × 32 × 1, and the input feature dimension was 128. Finally, the effect of the gait phase on the classification performance was also evaluated. First, we extract the PLI matrix of each gait phase (i.e., LC→RO, RO→RC, RC→LO, and LO→LC). The size of the PLI matrix for a single gait phase was N × 1 × 32 × 32, and the input dimension was 1,024. Next, we merged PLI matrices of different gait phases according to the sequence of gait events (i.e., LC→RO→RC, LC→RO→RC→LO, and LC→RO→RC→LO→LC). The merged PLI matrix sizes were N × 2 × 32 × 32, N × 3 × 32 × 32, and N × 4 × 32 × 32, and the input feature dimensions were: 2,048, 3,072, and 4,096, respectively. It is worth mentioning that we did not train classifiers to separate all these gait phases but manually extracted data belonging to different gait phases. Finally, we sequentially use the PLI matrix of a single gait phase and the combined PLI matrices of multiple gait phases as features to detect the locomotion state of rat.
In this study, we obtained the average of classification results through 5-fold cross-validation. Specifically, all the datasets were divided into 5 parts, 1 of which was taken as the test set without repetition, and the other 4 were used as the training set to train the model. After that, the classification results of the model on the test set were calculated each time, and all classification results were averaged to obtain the final classification results. In order to evaluate the performance of the classifier more comprehensively, we use the accuracy and F1 score as evaluation metrics. Accuracy was defined as the proportion of correctly classified samples to the total samples, while F1 score was the harmonic mean of precision and sensitivity. Precision was the proportion of true positives (TP) to the sum of TP and false positives (FP), and sensitivity was defined as the proportion of the TP to the sum of the TP and false negatives (FN) [39]. Figure 3 shows the network connections under 3 different task conditions for each locomotion phase. We noticed that connections with big weights (red line) in the preparation phase appeared in the left somatomotor area (LMO), right somatomotor area (RMO), and retrosplenial area (RSP), while the strongest connection occurred between RMO and RSP. It was worth noting that the edges of the right somatosensory area were more than the left. In addition, compared with other phases under the same task condition, the connectivity of the preparation phase was stronger (larger edge weight and greater degrees). Interestingly, during the right hind limb swing (RO→RC) after the rats finished the preparation phase, under all task conditions, connections between RMO and RSP also remained high level. However, the connection weights between LMO and RSP were weakened compared to the preparation phase and much smaller than the connection between RMO and RSP. In addition, compared with other conditions, the left somatosensory area (LSS) under the UF condition has more edges, and the weight of the edges was relatively larger. Interestingly, the degree and weight of LSS under FU and UF conditions were larger than those under FF condition in the RC→LO phase and more evident in the FU condition. At the end of the walking phase, during the swing of the rat's left hind limb (LO→LC), most connections appeared between the RSP and the LSS. It was worth noting that the largest connection weight appears between RMO and RSP in the FU condition, while the opposite was confirmed in the FF and UF conditions. Furthermore, the degree and connection weight of LSS under FF and FU conditions were larger than those under UF condition. The connection weights under the UF condition were smaller than those in the other conditions, which occurs in all gait phases. Figure 4 shows the time-varying global metrics of the functional network. We note that in the first preparation phase (purple window), global efficiency and transitivity ( Fig. 4A and B) keep increasing until the end of the preparation phase. In contrast, during the walking phase (yellow window), both decrease until the walking phase ends, and they gradually recovered to the same level as the first preparation phase in the second preparation stage. The characteristic path length shows the opposite trend to other global metrics (Fig. 4C). Interestingly, in a period of time (RC→LO) after the right hind paw contacted the uneven terrain (FU condition), both the global efficiency (Fig. 4D) and the transitivity (Fig. 4E) were significantly improved compared with the FF and UF conditions, while the characteristic path length (Fig. 4F) Figure 5 shows the local metrics of each brain region in all gait phases. Similar to the global metrics, the local metrics of multiple brain regions in the RC→LO and LO→LC phases showed obvious difference under 3 different task conditions. However, there was no significant difference in the local metrics of each brain region under the 3 task conditions of LC→RO phase (preparation stage) and RO→RC phase. Specifically, in the RC→LO phase, the local metrics of almost all brain regions under the FU condition were stronger than those under other conditions, while in the LO→LC phase, the local metrics under the FF condition were stronger than those under FU and FU conditions. It was worth noting that in the RC→LO phase, the values of local efficiency (Fig. 5A), node strength (Fig. 5B), and clustering coefficient (Fig. 5C) of LMO under UF condition were smaller than those under other conditions. In addition, in the LO→LC phase, the value of the eigenvector centrality of the RMO under the FU condition was larger than that in the other conditions. Furthermore, in the RC→LO phase, the value of eigenvector centrality (Fig. 5D) of LSS under FU and UF conditions was larger than that under FF condition. In the LO→LC phase, the values of all local metrics of LSS under UF condition were smaller than those under FF and UF conditions. Moreover, all local metrics of RSP maintained the highest level in all task conditions, higher than other brain regions under the same condition. All local metrics of the visual area were very close under all conditions. Figure 6 shows the classification accuracy and F1 score of the 3 classifiers when using the raw EEG signal, whole-time PLI, and time-varying PLI as features to classify rat locomotion states. Clearly, all classifiers perform significantly better (P < 0.05) when using time-varying PLI as input features compared to using raw EEG signals and whole-time PLI. In addition, when using time-varying PLI as a feature for classification, the accuracy and F1 score of the NB classifier exceeded 95%, Fig. 3. Time-varying brain function networks. Each row represents a time-varying brain functional network under a specific task condition (FF, FU, and UF). Each column represents a locomotion phase. The node size denotes the degree, and the color and size of the edges represent the connection weight. Connections below 50% of the strongest connection were removed.

Analysis of classification results
followed by the SVM classifier (about 79%), and the KNN classifier performs the worst (about 58%). Figure 7A and B shows the accuracy and F1 score of the classifier when using time-varying network global and local metrics as features. In general, among all the metrics, the classification result of the Eigenvector centrality was the best. Among the 3 classifiers, NB has the best performance, followed by SVM, and the worst performance was KNN (accuracy and F1 score for all metrics were below 50%). The largest accuracy and F1 score come from the classification of the Eigenvector centrality by NB, which were 65.01% and 65.42%, respectively. Figure 7C and D shows the classification accuracy and F1 score of the 3 classifiers for the PLI matrices of individual and combined gait phases. Similarly, NB still performs best among the 3 classifiers, and KNN still performs worst. We note that for the same classifier, there was little difference in classification performance (the accuracy and F1 score of NB were around 75%, the accuracy and F1 score of SVM were about 60%, and the accuracy and F1 score of KNN were about 45%) when the PLI matrix of a single gait phase (preparation phase, RO→RC, RC→LO and LO→LC) was used as the input feature. However, as the number of PLI matrices from different gait phases increases (LC→RO→RC, LC→RO→RC→LO, and LC→RO→RC→LO→LC), the accuracy and F1 score of all classifiers show an upward trend (the accuracy and F1 score of NB and SVM increased significantly, while the accuracy and F1 score of KNN increased slightly) and reach a maximum after merging the PLI matrices of all gait phases.

Discussion
In this study, we analyzed the time-varying functional connectivity of rats when encountering unexpected terrain and proposed a method for classifying rat locomotion states based on time-varying network properties. Our results show that the brain functional connectivity of the rat varied with gait phases and terrain conditions. Furthermore, the classification results show that combining the PLI matrices of all gait phases as classification features achieves better classification performance, and the NB classifier exhibited the best performance. Our findings have important implications for understanding animals' motor control mechanisms in adapting to unexpected terrain. Our results show that the brain functional connectivity varied with gait phases. Overall, the connectivity in the preparation phase was stronger than in the walking phase. Specifically, the degree and connection weights in the preparation phase were larger than those in the walking phase, and the global efficiency and transitivity reach the highest at the end of the preparation stage and gradually decrease in the walking stage, while the characteristic path length shows the opposite trend (Figs. 3 and 4). A previous study of human walking showed that brain functional connectivity during walking phase was weaker than during standing phase [40]. From our results, this rule also applies to rat bipedal walking on a treadmill and does not change with different terrain. It is worth noting that the degree of participation of the cortex was related not only to the complexity of the task but also to the degree of volitional engagement [41,42]. During the preparation phase, the rat stands on the treadmill in a bipedal position, and the rat's body begins to tilt as the treadmill belt turns. The gradual increase of the tilt angle seriously disrupted the body balance of the rats, which gives the rat a strong desire to adjust the body. In this process, high-level motor function was recruited, which makes the information interaction between different cortical areas stronger. A study of rat walking suggests that the details of hind limb motor control in rats were completed by the spinal cord [43]. We speculate that the decrease of cortical functional connectivity in rats during walking phase was probably due to the fact that the cortex enters the idling state after releasing the mo tor task in the preparation phase, while the spinal controls the movement according to the task during walking. Interestingly, we noticed that the cortex seemed to be more concerned with the regulation of the hind limb that support the body. We found that in the double-support phase (preparation phase and RC→LO), high-level connectivity appeared in the bilateral somatomotor areas. In the single-support phase (RO→RC and LO→LC), the connectivity of the somatomotor area on the opposite side of the supporting hind limb was much stronger than that on the other side. Similar to our results, a previous study showed that the motor cortex was more interested in stance control than the hind limb swing [43]. This may be related to the higher demand for muscle activity of the hind limb supporting the body, and the somatomotor area contributes directly to changes in muscle activity [44]. In the doublesupport phase, the weight of the body was shared by both hind limb, and the muscle strength required by both hind limb was similar, resulting in the similar connectivity level of the bilateral somatomotor areas. In the single-support phase, the weight of the body was completely borne by 1 hind limb, which requires not only stronger muscle strength but also different muscles to cooperate with each other to maintain stable posture. It is worth noting that the somatosensory cortex also shows a similar pattern, and the change of the somatosensory cortex was probably caused by the change of gait. In general, bipedal animals gradually transfer their body weight to the opposite side of the swinging hind limb in the double-support phase, and in the swing phase, the body weight was supported completely by 1 hind limb [45][46][47]. This leads to increased muscle strength from supporting hind limb, which leads to an increase in sensory information coming from the muscles. The transition from 2-foot support to single-foot support increases the pressure on the tactile receptors of the paw, which also increases the incoming sensory information. The increase in incoming sensory information makes the somatosensory area more active in interacting with other cortical areas so that the brain can generate appropriate control instructions. Notably, RSP was recruited at all gait phases and connected to all cortical areas, and eigenvector centrality (Fig. 5D) also indicates its importance in adapting rats to unexpected terrain. This may be related to the functions of RSP with multisensory information integration, terrain feature information recognition, and motion planning [48][49][50].
Previous studies had shown that the cerebral cortex was more active when humans encountered unexpected external disturbances [5][6][7]. We also found that the brain information transmission ability increased in a short time after the rat's accidental contact with uneven terrain (RC→LO), manifested as a significant increase in global efficiency and transitivity and a significant decrease in characteristic path length (Fig.  4). In addition, (Fig. 5A to C) shows that local metrics of multiple brain regions were significantly increased under FU condition during this process. Obviously, this difference comes from uneven terrain, which makes the incoming sensory information different. In addition, the input of sensory information was crucial for gait adjustment in unpredictable environments, and the sensory cortex undoubtedly plays an important role in the adaptation of animals to unfamiliar environments [51,52]. In this study, the functional connectivity of the somatosensory area changed obviously. Specifically, after the rats were accidentally in contact with uneven terrain (FU), the degree and edge weights of the LSS were far greater than in contact with flat terrain (Fig. 3). Changes in eigenvector centrality also indicated the importance of the somatosensory cortex in rat adaptation to terrain changes (Fig. 5D). We noticed that the eigenvector centrality of LSS was much greater than that of walking on flat terrain (FF) after rats accidentally contact with uneven terrain (FU). The possible reason for this phenomenon was that, compared with the flat surface, the pressure on the tactile receptor of rat paws after contacting the uneven surface increases, thus more tactile information was transmitted. In addition, due to the uneven support surface, the state of the hind limb muscles of rats has changed obviously compared with that before, which changes the sensory information transmitted by muscles. These sensory information from muscles and skin play an important role in the motor recovery and gait update after sudden disturbance [53][54][55][56]. In this process, the somatosensory cortex not only needs to process a large amount of sensory information in a short time but also needs to interaction with other cortical regions as soon as possible to help complete motion planning in a timely manner. Obviously, this puts forward higher requirements for the computing ability of the cortex and consumes more energy. Interestingly, compared to normal walking (FF), there was a obvious reduction in brain connectivity during the transition between different terrains (FU) by swinging the hind limb (LO→LC). This was reflected in significant decreases in global efficiency and transitivity and increases in characteristic path length, as well as significant reductions in local metrics across multiple brain regions (Fig. 5A to C). Previous studies have shown that the spinal play an essential role in coping with terrain change [53,56]. Therefore, in order to enable the brain to recover the consumption caused by the previous fast calculation as soon as possible, the motion control during terrain transition may rely more on the spinal than normal walking. To our surprise, the connectivity under UF conditions was much weaker than that under FF and FU conditions, as seen in Fig. 3. In addition, the local metrics of LMO (Fig. 5A to C) and LSS (Fig. 5) were much smaller in the UF condition than in the other 2 conditions. In the swing phase, the global metrics were also significantly lower than that of normal walking (Fig. 4). As the uneven area used in this study was only 11 × 11 cm in size, unable to support continuous walking of rats. Therefore, after the rat completely transits to the uneven terrain, the next step must leave the uneven terrain and enter the flat terrain, which enables the rat to predict the terrain ahead when leaving the uneven terrain. Compared with the unpredictable situation, rats rely more on the spinal to regulate movement when they can predict environmental changes [51,57]. We speculate that in the face of a predictable environment, the cortex was less involved to reduce consumption and relies on spinal regulation to complete the transition of the terrain.
Some studies have shown that using raw EEG signals or whole-time functional network metrics as features can effectively detect the locomotion state of subjects [25,58]. Our results (Fig. 6) show that using the PLI of the time-varying network as a feature was more effective in detecting the locomotion state of rats during adaptation to unexpected terrain than using raw EEG signals or whole-time functional network metrics. This shows that compared with the other 2 methods, using time-varying functional network to analyze the process of rats adapting to unexpected terrain can obtain more useful information. The classification results (Fig. 7) show that using the time-varying adjacency matrix of the network as a feature classification works best among all network metrics. In addition, the classification results of PLI matrices for different gait phases ( Fig. 7C and D) show that classification using PLI matrices of multiple gait phases combined in chronological order leads to higher accuracy than classification using PLI matrices of a single gait phase. A previous study showed that the classification accuracy tends to increase as the number of available features increases and remains constant after the number of features reaches a certain number [59]. Obviously, the improvement of classifier performance was directly related to the number of features, and the PLI matrix synthesizes the global and local features of the functional network, so that the classifier can extract more useful features for classification. In addition, our analysis shows that functional connectivity varies with gait phases, which makes each gait phase possess features that were different from other gait phases. Therefore, combining network metrics for multiple gait phases yields more useful features than using a single gait phase network metric.
There were also some limitations of this study. In this study, we did not consider the effect of rat learning on functional connectivity. However, considering the random occurrence of uneven areas, and the change of position and angle of contact between rat paws and uneven areas each time, we believe that the learning of rats was not enough to affect the experimental results. In the future work, we can further discuss the relationship between functional connectivity and the number and the order of trials to clarify the impact of rat learning on functional connectivity. In addition, we only analyzed functional connectivity within a wide frequency range, while analyzing functional connectivity in different frequency bands may get more information. Finally, this study only focuses on the influence of terrain changes caused by the accidental appearance of uneven terrain on the functional connectivity but does not analyze the impact of other kinds of terrain. This may hamper the possibility of generalizing the current results, so it will be necessary to add more types of terrain to verify the current results in the future work. It is worth mentioning that discussing the individual differences in adapting to unexpected terrain was also an important future work for us.
In the present study, we performed a time-varying analysis of functional connectivity in bipedal walking rats that encountered unexpected terrain. Dynamic functional connectivity maps, time-varying functional network global metrics, and local network metrics for different gait phases suggest that information interactions between cortical areas vary with gait phase and terrain conditions. In addition, the classification results emphasize the potential of EEG-based dynamic functional network analysis in animal locomotion detection. The findings of this study deepen our understanding of the role of the cortex in the adaptation of animals to unexpected terrain from the network viewpoint.