MRI data were collected from 20 healthy right-handed subjects (10 males/females, 24.5 $\pm$ 2.16 years).

The fMRI experiment contained six successive sessions for every volunteer. The first session had a duration of 480 s, and resting-state fMRI data were collected. The second session tested a self-paced, free-speed walking task to simulate a usual walk without pressure following an audio alert “alternating left and right sides in a free speed and keeping a constant speed” repeated five times for 30 s. Subjects were instructed to perform ankle dorsiflexion and plantarflexion alternately and repeatedly to simulate rhythmic, bilaterally alternating muscle activities in the lower limbs during walking. The following four sessions tested ankle movements at 0.5, 1, and 2 Hz frequencies or a changing speed in a counterbalanced order for every four subjects. Subjects paced the movements following the audio alert “left–right–left–right…” at preset frequencies. The 0.5, 1, and 2 Hz conditions consisted of 15, 30, and 60 audio alerts, respectively, during a task block. The audio alert had a fixed order in every task block of changing-speed conditions: 4 audio alerts with 1 Hz, 6 alerts with 1.5 Hz, 8 alerts with 2 Hz, 6 alerts with 1.5 Hz, 8 alerts with 2 Hz, 8 alerts with 1.5 Hz, and 4 alerts with 1 Hz with an average frequency of 1.47 Hz.

The following four sessions were counterbalanced for every four subjects in the order of arrival to the MRI room.

During fMRI scanning, each subject was instructed to close his/her eyes but remain awake while lying in a supine position with a pad under the knees to support the thighs to relax the spine and prevent motion from transferring to the head during ankle movements [10].

A monitoring system displayed the real-time amplitude of translational head motions of the subject during MR scanning. Data scanning was restarted if head movements surpassed 1 mm in the x-, y-, or z-axis.

T2* images were collected using the Siemens Prisma 3T MRI (Siemens Healthineers Ltd., Germany) system (time of repetition [TR], 2000 ms; time of echo [TE], 30 ms; flip angle, 90${}^{\circ }$; and voxel size, 2.0 $\times$ 2.0 $\times$ 2.0 mm) consisting of 64 interleaved axial slices with 0.2 mm gaps between slices. The 3D MPRAGE sequence was used to acquire T1-weighted images (TR, 2530 ms; TE, 2.98 ms; flip angle, 7${}^{\circ }$; inversion time, 1100 ms; and voxel size, 0.5 $\times$ 0.5 $\times$ 1.0 mm). Field mapping images were also collected (TR, 635 ms; TE, 4.92 and 7.38 ms; flip angle, 60${}^{\circ }$; and voxel size 2.0 $\times$ 2.0 $\times$ 2.0 mm).

MRI data were obtained from a previous study [26], which focused on brainstem and brain activations by using the model-based approach and investigated the region of interest-based functional connectivity measured with Pearson correlation coefficients.

Data and the informed consent of all participants were obtained. The institutional review board of the National Research Center for Rehabilitation Technical Aids approved the protocol (code: 2018YFC2001401).

The SPM software (SPM12, Wellcome Centre for Human Neuroimaging, University College London, https://www.fil.ion.ucl.ac.uk/spm/software/spm12) was used for temporal and spatial preprocessing. Differences in slice acquisition times of functional images were corrected first. The FieldMap Toolbox was applied while performing realignment and unwarped simultaneously to correct head motion and geometric distortion. T2* and T1 images were coregistered and normalized to the standard MNI space as introduced in [27].

ICA was performed using the GIFT4.0b (http://mialab.mrn.org/software/gift/index.html). The number of ICs from each condition was estimated with the minimum description length criteria [28]. The Infomax algorithm [29] was repeated 10 times with the ICASSO [30] to separate ICs from the fMRI data set of each condition. The one-sample t test was used to acquire the contrast images and t map of each component (against the null of zero task-related change either positive or negative, i.e., a two-sided test) [31]. The spatial map of each component was then acquired using the SPM software at a threshold of p 0.001 with a cluster extent of more than 10 contiguous voxels.

Motion artifacts were identified in accordance with previous studies [32, 33, 34, 35] if the main part of a component was located on large vessels (e.g., superior sagittal or transverse sinus), margins of the brain (displaying ring-like shapes), or regions full of cerebrospinal fluid, such as that surrounding the brainstem.

Task-evoked ICs were identified if the Pearson correlation coefficient $\ge$0.3 when performing “temporal sorting” function for task fMRI data sets with the GIFT software. ICs with correlation coefficient 0.3 might present a spatial pattern similar to that of the default mode network, which is found to be uncoupled from task waveform.

Spatial correlation coefficients between each IC map of each condition and all other IC maps of the remaining conditions were calculated to identify task-evoked components that commonly appeared in all conditions. For a component, strong and positive correlation values indicated that the spatial map of this component was stable across all conditions, and the brain regions in the spatial map of this component constituted a distributed brain network.

The estimated numbers of ICs in each condition were 25 (0.5 Hz), 26 (1 Hz), 26 (2 Hz), 25 (changing-speed condition), 26 (free-speed condition), and 27 (resting-state condition). Three task-related components that commonly appeared in all conditions presented positive correlation $>$0.3 for each task condition. On the basis of their spatial maps, these ICs were named as cerebellum, PAC-S2-IPC-preSMA, and medial M1S1-SMA-PreM-SPL networks. In addition, one component presented negative correlations with the task time course in 0.5 and 1 Hz conditions but positive correlations with the task time course in the 2 Hz, changing-speed (average frequency = 1.47 Hz), and free-speed (average frequency = 1.46 $\pm$ 0.67 Hz) conditions. This component was also detected in the resting state and named the lateral S1-SPL-PMd network on the basis of its spatial map.

The common spatial map of each of the four components across all conditions (Fig. 1) was determined using the inclusive masking function in the SPM software at primary thresholds of p 0.001 and p 0.05 family-wise error corrected cluster extent provided by the SPM toolbox for the whole search volume.

Fig. 1.

Spatial maps of four independent components in all conditions. The spatial map of each commonly appearing independent component was projected on the average T1 images by using the MRIcron software (v1.0.20190902, https://www.nitrc.org/projects/mricron). The numbers indicate MNI coordinates in the z-axis. Cerebellum network: red; PAC-S2-IPC-preSMA network: yellow; medial M1S1-SMA-PreM-SPL network: blue; S1-SPL-PMd network: purple.

The cerebellum network was distributed in the cerebellum (from z = –55 to z = –10 in Fig. 1). The PAC-S2-IPC-preSMA network (Supplementary Figs. 1,2) was mainly distributed in the PAC and S2. Some parts of the network were distributed in the IPC and preSMA. The PAC is located at Heschl’s gyrus, including cytoarchitectonic subdivisions TE 1.0, 1.1, and 1.2 [36] and at the medial part of the temporal pole lateral to the posterior insula [37].

The lateral S1-SPL-PMd network (Supplementary Figs. 3,4) was mainly distributed in the parietal cortex (lateral parts of S1 and SPL). Some parts of the network were distributed in the PMd. The extensions of preSMA, SMA, and PMd were determined using a meta-analysis [38].

The medial M1S1-SMA-PreM-SPL network was mainly distributed in the medial parts of the primary sensorimotor cortex (M1 and S1). This component also covered bilateral SMA, PreM, and SPL (SupplementaryFigs. 5,6). The scopes of these brain regions except that of SMA were determined using the SPM Anatomy toolbox.

For each of the four components (networks), the pairwise spatial correlation values among different conditions are shown in Table 1. Under most conditions, the medial M1S1-SMA-PreM-SPL network presented the lowest correlation value, ranging from 0.30 to 0.38 between the resting-state and other conditions. However, the pairwise correlation values of the components between the task conditions were not less than 0.54. These results indicated that this network exhibited much more variations between resting-state and task states than between different task conditions. As shown in Supplementary Fig. 7, the network covered lateral parts of M1 and S1 during the resting state and medial parts of these brain regions during task conditions.

Correlation values were ranked in descending order in each row. The cerebellum network always presented the highest correlation value, that is, the most consistent spatial maps across all conditions. Under most conditions, the lateral S1-SPL-PMd network presented the second highest correlation value, whereas the medial M1S1-SMA-PreM-SPL network presented the lowest correlation value.

The cerebellum network covers most parts of the cerebellum and exhibits the highest level of consistency across all conditions (Table 1). These findings suggest that this component plays a fundamental role in walking control. Compared with simple limb movements, the cerebellum is more likely to participate in modulating complicated movements [39]. Additionally, the cerebellum is a place to modulate movements with somatosensory inputs [40]. The motor tasks applied in this study are bilaterally involved alternating ankle dorsiflexion and plantarflexion, which are more complex than unilateral ankle movements and require more somatosensory information for feedforward and feedback controls. Therefore, this component covers extensive areas of the cerebellum.

The cerebellum is responsible for the coordination of movements, especially skilled voluntary motor tasks. Cerebellar lesions lead to incoordination (ataxia) of volitional movement and gait disorders [41]. Therefore, the bilaterally alternating ankle dorsiflexion–plantarflexion movements depend on cerebellar functions.

The PAC-S2-IPC-preSMA network mainly covers the PAC, indicating a primary role in processing audio information. This feature is consistent with the audio cues used to prompt ankle movements for all tasks in this study. This network also covers some parts of the preSMA. When performing sequential movements, the preSMA participates in chunking and related processes, such as task switching, response selection, and response inhibition [42]. Thus, for sequential movements, such as simulated walking movements, preSMA can play a role in separating the motor task into minor fragments and in switching these fragments. The preSMA uses audio information from the PAC to achieve these functions. Other brain regions in these networks, including S2 and IPC, can provide sensory information for this function. The inferior temporal cortex preferentially participates in the cognitive aspects of discriminating shape and object, whereas the IPC plays roles in object-oriented action, object recognition, and motor planning [43]. This network controls the rhythm of the simulated walking tasks by integrating audio inputs (PAC), chunking/switching function (preSMA), and sensory processing in high-order sensory cortices (S2 and IPC).

The lateral S1-SPL-PMd network is mainly distributed in the parietal lobe, including the lateral parts of S1 and SPL. This network also covers a small part of PMd. This network presents a negative temporal correlation value with the time course of motor tasks with slow speeds (0.5 and 1 Hz) but a positive temporal correlation value with the time course of high-speed motor tasks (2 Hz, free- and changing-speed conditions). Thus, this network plays a role in differentiating walking speed.

The SPL participates in sensorimotor integration, such as tactile exploration [43] and processing kinesthetic cues during action-related somatosensory inputs [44]. The PMd controls proximal limb muscles for positioning the limb to perform a motor task during the preparation stage [45]. According to a previous study [46], the lateral part of S1 in this network covers the representation area of body parts above the knees but not the lower leg. We hypothesize that the lateral S1-SPL-PMd network integrates somatosensory inputs from the rest of the body (other than the moving part) and the functions of the PMd (movement preparation) and SPL (sensorimotor integration) to facilitate the involvement of the whole body during fast walking. This network does not directly control movement because M1 is not included in this network.

The medial M1S1-SMA-PreM-SPL network is mainly distributed over the medial parts of M1 and S1 and covers the representation areas of the lower leg, including ankles and feet. These regions are commonly reported in previous fMRI studies on foot movements and are fundamental brain regions for motor execution and related sensory processing. Hence, this network is a basic actor network to control the walking movement.

This basic actor network can be detected in all motor task conditions and in the resting state. The map during the resting state covers more lateral parts of M1 and S1 than the map during task conditions (Supplementary Fig. 7). Such spatial inconsistency leads to low correlation values (Table 1). The spatial map in the task conditions is only the medial part of the map in the resting state because only the lower leg is involved in motor tasks.