A somato-cognitive action network alternates with effector regions in motor cortex

Washington University adult participants

Data were collected from three healthy, right-handed, adult participants (aged 35, 25 and 27 years; 1 female) as part of a study investigating effects of arm immobilization on brain plasticity (previously published data^29,30,69). Written informed consent was obtained from all participants. The study was approved by the Washington University School of Medicine Human Studies Committee and Institutional Review Board. The primary data employed here were collected either prior to the immobilization intervention (participants 1 and 3) or two years afterwards (participant 2). Data collected immediately after the intervention are presented for within-participant replication in Extended Data Fig. 1b. For details concerning data acquisition and processing, see ref. ²⁹.

For two participants (participants 1 and 2), we collected additional fMRI data using the same sequence during performance of two motor tasks: a somatotopic mapping task and a motor control task.

Movement task battery

A block design was adapted from the motor task in ref. ³¹. In each run, the participant was presented with visual cues that directed them to perform one of five specific movements. Each block started with a 2.2-s cue indicating which movement was to be made. After this cue, a centrally presented caret replaced the instruction and flickered once every 1.1 s (without temporal jittering). Each time the caret flickered, participants executed the proper movement. Twelve movements were made per block. Each block lasted 15.4 s, and each task run consisted of 2 blocks of each type of movement as well as 3 blocks of resting fixation. Movements conducted within each run were as follows:

Run type 1: Close left (L) hand/Close right (R) hand/Flex L foot/Move tongue L and R (participant 1: 24 runs; participant 2: 20 runs).

Run type 2: Flex L elbow/Flex R elbow/Flex L wrist/Flex R wrist/Lift bilateral shoulders (participant 1: 10 runs; participant 2: 11 runs).

Run type 3: Flex L gluteus/Flex R gluteus/Tense abdomen/Open and close mouth/Swallow (participant 1: 10 runs; participant 2: 11 runs).

Run type 4: Flex L ankle/Flex R ankle/Bend L knee/Bend R knee/Flex bilateral toes (participant 1: 10 runs; participant 2: 11 runs).

Run type 5: Lift L eyebrow/Lift R eyebrow/Wink L eyelid/Wink R eyelid/Flare nostrils (participant 1: 10 runs; participant 2: 11 runs).

Action control and coordination task

An event-related design implemented using JSpsych toolbox v6.3 was used to discriminate planning and execution of limb movement. See Supplementary Fig. 6 for an illustration of this task. Within the run, the participant is prompted to move either a single limb or to simultaneously move two limbs. There are four possible motions—open–close of fingers or toes, left–right flexion of the wrist or ankle, clockwise rotation of the wrist or ankle, and anticlockwise rotation of the wrist or ankle—each of which may be executed by any of the four extremities (left or right upper or lower extremity). Each motion–extremity combination may be required in isolation, or in combination with a second simultaneous motion. The participant is cued to prepare the movement(s) when they see one or two movement symbols placed on a body shape in a grey colour (planning phase), and is then cued to execute the movement(s) when the grey symbol or symbols turn green (execution phase). Using a pseudorandom jitter, the planning phase can last from 2 to 6.5 s followed by 4 to 8.5 s of movement execution. Each movement trial (planning and execution) is followed by a jittered fixation of up to 5 s. A rest block of 8.6 s is implemented every 12 movements. Two possible movements are requested during the task run and practiced before the task. The movement pair is changed for each task run. 48 trials were collected in each run. Twelve total runs were acquired per participant.

Laryngeal mapping task

For the same two participants, as well as for one additional participant (participant 8, 40 years of age, from whom written informed consent was obtained) additional fMRI data were collected during performance of a laryngeal mapping task using a multiband five-echo blood oxygen level-dependent (BOLD) contrast sensitive gradient echo-planar sequence (flip angle = 68°, resolution = 2.0 mm isotropic, TR = 1,761 ms, multiband 6 acceleration, TE₁: 14.20 ms, TE₂: 38.93 ms, TE₃: 63.66 ms, TE₄: 88.39 ms, and TE₅: 113.12 ms), with each run lasting 3 min 52 s. A pair of spin-echo echo-planar images (EPI) with opposite phase-encoding directions (anterior→posterior (AP) and posterior→anterior (PA)) but identical geometrical parameters to the BOLD sequence were acquired. In participant 8, 15 min of resting-state fMRI were also acquired using the sequence above.

An additional set of ten movement task runs adapted from³¹ were collected to localize laryngeal phonation. In each run, the participant was presented with visual cues that directed them to perform one of six specific movements: Left hand, Right hand, Left foot, Right foot, Tongue or Voice. In the Voice condition, participants were required to briefly make the noise “eeee” without moving their jaw. Each block started with a 3.0 s cue indicating which movement was to be made. After this cue, a centrally presented caret replaced the instruction and flickered once every 1.0 s (without temporal jittering). Each time the caret flickered, participants executed the proper movement. Ten movements were made per block. Each block lasted 15.0 s, and each task run consisted of 2 blocks of each type of movement as well as 2 blocks of resting fixation. Each participant completed 10 runs.

Cornell adult participants

Data were collected from four healthy adult participants (ages 29, 38, 24 and 31; all male) as part of a previously published study⁷⁰. The study was approved by the Weill Cornell Medicine Institutional Review Board. Written informed consent was provided by each participant.

For details concerning data acquisition and processing, see ref. ⁷⁰.

Neonatal participant

Data were collected from one sleeping, healthy full-term neonatal participant beginning 13 days after birth, corresponding to 42 weeks post-menstrual age. The study was approved by the Washington University School of Medicine Human Studies Committee and Institutional Review Board. Written informed consent was provided by a parent.

MRI acquisition

The participant was scanned while asleep over the course of 4 consecutive days using a Siemens Prisma 3T scanner on the Washington University Medical Campus. Every session included collection of a high-resolution T2-weighted spin-echo image (TE = 563 ms, TR = 3,200 ms, flip angle = 120°, 208 slices with 0.8 × 0.8 × 0.8 mm voxels). In each session, a number of 6 min 45 s multi-echo resting-state fMRI runs were collected as a five-echo BOLD contrast sensitive gradient echo-planar sequence (flip angle = 68°, resolution = 2.0 mm isotropic, TR = 1,761 ms, multiband 6 acceleration, TE₁: 14.20 ms, TE₂: 38.93 ms, TE₃: 63.66 ms, TE₄: 88.39 ms, and TE₅: 113.12 ms). The number of BOLD runs collected in each session depended on the ability of the neonate to stay asleep during that scan; across the 4 days, 23 runs were collected in total. A pair of spin-echo EPI images with opposite phase-encoding directions (AP and PA) but identical geometrical parameters and echo spacing were acquired between every three BOLD runs or any time the participant was removed from the scanner.

MRI processing

Structural and functional processing followed the pipeline used for the Washington University dataset, with two exceptions. First, segmentation, surface delineation and atlas registration were conducted using a T2-weighted image (the single highest quality T2 image, as assessed via visual inspection) rather than a T1-weighted image, due to the inverted image contrast observed in neonates. Second, after the multi-echo BOLD data were unwarped and normalized to atlas space, it was optimally combined before nuisance regression and mapping to cifti space. All fMRI scans from the second day of scanning were excluded due to registration abnormalities.

Infant participant

Data were collected from one healthy sleeping infant aged 11 months. The study was approved by the Washington University School of Medicine Human Studies Committee and Institutional Review Board. Written informed consent was provided by a parent.

MRI acquisition

The participant was scanned while asleep over the course of three sessions using a Siemens Prisma 3T scanner on the Washington University Medical Campus. The first session included collection of a high-resolution T1-weighted MP-RAGE (TE = 2.24 ms, TR = 2,400 ms, flip angle = 8°, 208 slices with 0.8 × 0.8 × 0.8 mm voxels) and a T2-weighted spin-echo image (TE = 564 ms, TR = 3200 ms, flip angle = 120°, 208 slices with 0.8 × 0.8 × 0.8 mm voxels). The second and third sessions included collection of 26 total runs of resting-state fMRI, each collected as a 6 min 49 s-long BOLD contrast sensitive gradient echo-planar sequence (flip angle = 52°, resolution = 3.0 mm isotropic, TE = 30 ms, TR = 861 ms, multiband 4 acceleration). For each run, a pair of spin-echo EPI images with opposite phase-encoding directions (AP and PA) but identical geometrical parameters and echo spacing were acquired to correct spatial distortions.

MRI processing

Structural processing followed the DCAN Labs processing pipeline found in the ABCD BIDS Community Collection (ABCC; NDA Collection 3165) (https://github.com/DCAN-Labs/abcd-hcp-pipeline)⁷¹, which we found performed the best surface segmentation at this age. Functional processing followed the pipeline used for the Washington University adult dataset.

Child participant

Data were collected from one healthy awake male child age 9 years. The study was approved by the Washington University School of Medicine Human Studies Committee and Institutional Review Board. Written informed consent was provided by a parent and assent was given by the participant.

MRI acquisition

The participant was scanned repeatedly over the course of 12 sessions using a Siemens Prisma 3T scanner on the Washington University Medical Campus. These sessions included collection of 14 high-resolution T1-weighted MP-RAGE images (TE = 2.90 ms, TR = 2,500 ms, flip angle = 8°, 176 slices with 1 mm isotropic voxels), 14 T2-weighted spin-echo images (TE = 564 ms, TR = 3200 ms, flip angle = 120°, 176 slices with 1 mm isotropic voxels), and 26 total runs of resting-state fMRI, each collected as a 10 min-long BOLD contrast sensitive gradient echo-planar sequence (flip angle = 84°, resolution = 2.6mm isotropic, 56 slices, TE = 33 ms, TR = 1,100 ms, multiband 4 acceleration). In each session, a pair of spin-echo EPI images with opposite phase-encoding directions (AP and PA) but identical geometrical parameters and echo spacing were acquired to correct spatial distortions in the BOLD data.

MRI processing

Structural and functional processing followed the DCAN Labs processing pipeline found in the ABCD BIDS Community Collection (ABCC; NDA Collection 3165)⁷¹ (https://github.com/DCAN-Labs/abcd-hcp-pipeline).

Participant with perinatal stroke

PS1, a left-handed, 13-year-old male who played for a competitive youth baseball team, was referred to an orthopaedic physician because of difficulty using his right arm effectively. Ulnar neuropathy was considered and he was referred for physical therapy. However, PS1 was first seen by a child neurologist (N.U.F.D.) for further evaluation. Structural brain MRI revealed unexpectedly extensive bilateral cystic lesions consistent with perinatal infarcts. Review of PS1’s medical history revealed that the injury occurred in the perinatal period.

Data acquisition from PS1 were performed with the approval of the Washington University Institutional Review Board. Written informed consent was provided by PS1’s mother and assent was given by PS1 at the time of data acquisition.

For additional details regarding clinical history, neuropsychological evaluations, motor assessments, or MR image acquisition or processing, see ref. ³².

UMN macaque

Data were collected from a sedated adult female macaque monkey (Macaca fascicularis) aged 6 years. Experimental procedures were carried out in accordance with the University of Minnesota Institutional Animal Care and Use Committee and the National Institute of Health standards for the care and use of non-human primates. The subject was fed ad libitum and pair-housed within a light- and temperature-controlled colony room. The animal was not water restricted. The subject did not have any prior implant or cranial surgery. The animal was fasted for 14–16 h prior to imaging. On scanning days, anaesthesia was first induced by intramuscular injection of atropine (0.5 mg kg⁻¹), ketamine hydrochloride (7.5 mg kg⁻¹), and dexmedetomidine (13 μg kg⁻¹). The subject was transported to the scanner anteroom and intubated using an endotracheal tube. Initial anaesthesia was maintained using 1.0%–2% isoflurane mixed with oxygen (1 l min⁻¹ during intubation and 2 l m⁻¹ during scanning to compensate for the 12-m length of the tubing used). For functional imaging, the isoflurane level was lowered to 1%. The subject was placed onto a custom-built coil bed with integrated head fixation by placing stereotactic ear bars into the ear canals. The position of the animal corresponds to the sphinx position. Experiments were performed with the animal freely breathing. Continuous administration of 4.5 μg kg⁻¹ h⁻¹ dexmedetomidine using a syringe pump was administered during the procedure. Rectal temperature (~37.6 °C), respiration (10–15 breaths per min), end-tidal CO₂ (25–40), electro-cardiogram (70–150 bpm), and peripheral capillary oxygen saturation (SpO₂) (>90%) were monitored using an MRI compatible monitor (IRAD-IMED 3880 MRI Monitor). Temperature was maintained using a circulating water bath as well as chemical heating pads and padding for thermal insulation.

MRI acquisition

Data were acquired on a Siemens Magnetom 10.5 T Plus. A custom in-house built and designed RF coil was used with an 8-channel transmit/receive end-loaded dipole array of 18-cm length combined with a close-fitting 16- channel loop receive array head cap, and an 8-channel loop receive array of 50 × 100 mm² size located under the chin⁷². A B1+ (transmit B1) field map was acquired using a vendor provided flip angle mapping sequence and then power calibrated for each individual. Following B1+ transmit calibration, 3–5 averages (23 min) of a T1-weighted MP-RAGE were acquired for anatomical processing (TR = 3300 ms, TE = 3.56 ms, TI = 1,140, flip angle = 5°, slices = 256, matrix = 320×260, acquisition voxel size = 0.5 × 0.5 × 0.5 mm³, in-plane acceleration GRAPPA = 2). A resolution and field of view-matched T2-weighted 3D turbo spin-echo sequence was run to facilitate B1 inhomogeneity correction. Five images were acquired in both phase-encoding directions (R→L and L→R) for offline EPI distortion correction. Six runs of fMRI time series, each consisting of 700 continuous 2D multiband EPI^73,74,75 functional volumes (TR = 1,110ms; TE = 17.6 ms; flip angle = 60°, slices = 58, matrix = 108 × 154; field of view = 81 × 115.5 mm ; acquisition voxel size = 0.75 × 0.75 × 0.75 mm) were acquired with a left–right phase-encoding direction using in-plane acceleration factor GRAPPA = 3, partial Fourier = 7/8, and MB factor = 2. Since the macaque was scanned in sphinx position, the orientations noted here are what is consistent with a (head first supine) typical human brain study (in terms of gradients) but translate differently to the actual macaque orientation.

MRI processing

Processing followed the DCAN Labs non-human primate processing pipeline (http://github.com/DCAN-Labs/nhp-abcd-bids-pipeline), with minor modifications. Specifically, we observed that distortion from the 10T scanner was so extensive that the field maps did not fully correct it. Therefore, instead of field map-based unwarping, we used the computed field map-based warp as an initial starting point for Synth, a field map-less distortion correction algorithm that creates synthetic undistorted BOLD images for registration to anatomical images⁷⁶. Synth substantially reduced residual BOLD image distortion.

PRIME-DE macaque

Raw structural and functional data were provided from the Oxford dataset of the PRIMatE Data Exchange (PRIME-DE) consortium (https://fcon_1000.projects.nitrc.org/indi/PRIME/oxford.html)^77,78. The full dataset consisted of 19 (age 4.1 ± 0.98 years, weight 6.61 ± 2.94 kg) rhesus macaques (Macaca mulatta). The animal care, anaesthesia and MRI protocols were carried out in accordance with the UK Animals (Scientific Procedures) Act of 1986. Animals in the study were group-housed prior to scanning. Ketamine (10 mg kg⁻¹) was administered via intramuscular injection for induction, along with either xylazine (0.125–0.25 mg kg⁻¹), midazolam (0.1 mg kg⁻¹), or buprenorphine (0.01 mg kg⁻¹). Additionally, injections of atropine (0.05 mg kg⁻¹, intramuscular injection), meloxicam (0.2 mg kg⁻¹, intravenous injection), and ranitidine (0.05 mg kg⁻¹) were administered. A minimum of 15 min prior to being placed in the stereotaxic frame, animals also received local anaesthetics (5% lidocaine/prilocaine cream and 2.5% bupivacaine injected subcutaneously around ears). Finally, anaesthesia was maintained with isoflurane, and scanning began 1.5–2 h after the initial ketamine induction.

MRI acquisition

Anaesthetized animals were placed in the sphinx position into a stereotactic frame (Crist Instrument) and scanned in a horizontal 3T MRI scanner using a four-channel phased-array coil (Windmiller Kolster Scientific, Fresno, CA). Each animal received 53.33 min (1,600 volumes) of resting-state data, which was acquired at a 2.0 mm isotropic voxel resolution (TR = 2,000 ms, TE = 19 ms, Flip angle = 90°). A T1-weighted MP-RAGE sequence was used to acquire anatomical data (TR = 2,500 ms, TE = 4.01 ms, TI = 1,100, flip angle = 8°, acquisition voxel size = 0.5 × 0.5 × 0.5 mm, 128 slices).

MRI processing

Processing for structural data followed the DCAN Labs non-human primate processing pipeline (https://github.com/DCAN-Labs/nhp-abcd-bids-pipeline). Smoothing was applied with FWHM = 1.5 mm in both volume and surface space. The surface data were then down-sampled to a 10k surface to create the preprocessed cifti data. Finally, each animal’s data was closely visually inspected for quality. Following these inspections, data from 11 animals were excluded due to the presence of artefact in or near the central sulcus, leaving eight animals in the final data. This sample size of eight was chosen to include all available artefact-free data. No randomization or blinding was performed.

Group-averaged datasets

Resting-state fMRI data was averaged across participants within each of five large datasets.

UK Biobank

A group-averaged weighted eigenvectors file from an initial batch of 4,100 UKB participants aged 40–69 years (53% female) scanned using resting-state fMRI for 6 min was downloaded from https://www.fmrib.ox.ac.uk/ukbiobank/. This file consisted of the top 1,200 weighted spatial eigenvectors from a group-averaged principal component analysis. See ref. ⁷⁹ and documentation at https://biobank.ctsu.ox.ac.uk/crystal/ukb/docs/brain_mri.pdf for details of the acquisition and processing pipeline. This eigenvectors file was mapped to the Conte69 surface template atlas⁸⁰ using the ribbon-constrained method in Connectome Workbench⁸¹, and the eigenvector time courses of all surface vertices were cross-correlated.

Adolescent Brain Cognitive Development Study

Twenty minutes (4 × 5-min runs) of resting-state fMRI data, as well as high-resolution T1-weighted and T2-weighted images, were collected from 3,928 9- to 10-year-old participants (51% female), who were selected as the participants with at least 8 min of low-motion data from a larger scanning sample. Data collection was performed across 21 sites within the USA, harmonized across Siemens, Philips and GE 3T MRI scanners. See ref. ⁸² for details of the acquisition parameters. Data processing was conducted using the ABCD-BIDS pipeline found in the ABCD BIDS Community Collection (ABCC; NDA Collection 3165) (https://github.com/DCAN-Labs/abcd-hcp-pipeline)⁷¹; see ref. ⁸³ for details.

Human Connectome Project

A vertexwise group-averaged functional connectivity matrix from the HCP 1200 participants release was downloaded from https://db.humanconnectome.org. This matrix consisted of the average strength of functional connectivity across all 812 participants aged 22–35 years (410 female) who completed 4× 15-min resting-state fMRI runs and who had their raw data reconstructed using the newer recon 2 software. See refs. ^81,84,85,86 for details of the acquisition and processing pipeline.

Washington University 120

Data were collected from 120 healthy young adult participants recruited from the Washington University community during relaxed eyes-open fixation (60 females, ages 19–32). Scanning was conducted using a Siemens TRIO 3.0T scanner and included collection of high-resolution T1-weighted and T2-weighted images, as well as an average of 14 min of resting-state fMRI. See ref. ⁸⁷ for details of the acquisition and processing pipeline.

Neonates

Mothers were recruited during the second or third trimester from two obstetrics clinics at Washington University as part of the Early Life Adversity, Biological Embedding, and Risk for Developmental Precursors of Mental Disorders (eLABE) study. This study was approved by the Human Studies Committees at Washington University in St. Louis and written informed consent was obtained from mothers. Neuroimaging was conducted in full-term, healthy neonate offspring shortly after birth (average post-menstrual age of included participants 41.4 weeks, range 38–45 weeks). Of the 385 participants scanned for eLABE, 262 were included in the current analyses (121 female). See ref. ⁸⁸ for additional details of the participants, criteria for exclusion, scanning acquisition protocol and parameters, and processing pipeline.

Analyses

Functional connectivity

For each single-participant dataset, a vertex or voxelwise functional connectivity matrix was calculated from the resting-state fMRI data as the Fisher-transformed pairwise correlation of the time series of all vertices/voxels in the brain. In the ABCD, Washington University 120, eLABE and PRIME-DE datasets, vertex and voxelwise group-averaged functional connectivity matrices were constructed by first calculating the vertex or voxelwise functional connectivity within each participant as the Fisher-transformed pairwise correlation of the time series of all vertices or voxels in the brain, and then averaging these values across participants at each vertex or voxel.

Seed-based functional connectivity

We defined a continuous line of seeds down the left precentral gyrus by selecting every vertex in a continuous straight line on the cortical surface between the most ventral aspect of the medial motor area (approximate MNI coordinates (−4, −31, 54)) and the ventral lip of the precentral gyrus right above the operculum (approximate MNI coordinates (−58, 4, 8)). For each seed, we examined its map of functional connectivity as the Fisher-transformed correlation between that vertex’s time course and that of every other vertex or voxel in the brain.

Network detection in somatomotor cortex

To define the somatomotor regions that were visually identified from the seed-based connectivity analysis in an unbiased fashion for further exploration, we entered each individual adult human participant’s data into a data-driven network detection algorithm designed to identify network subdivisions that are hierarchically below the level of classic large-scale networks (for example, those that produce hand/foot divisions in somatomotor cortex;^23,37). We have previously described how this approach identifies sub-network structures that converge with task-activated regions⁸⁹ and with known neuroanatomical systems⁹⁰.

In each adult participant, this analysis clearly identified network structures corresponding to motor representation of the foot, hand and mouth; and it additionally identified network structures corresponding exactly to the previously unknown connectivity pattern identified from the seed-based connectivity exploration as the inter-effector regions. For simplicity, we manually grouped all inter-effector subnetworks together as a single putative network structure (labelled as inter-effector) for further analysis.

Finally, to identify classic large-scale networks in each participant, we repeated the Infomap algorithm on matrices thresholded at a series of denser thresholds (ranging from 0.2% to 5%), and additionally identified individual-specific networks corresponding to the default, medial and lateral visual, cingulo-opercular, fronto-parietal, dorsal attention, language, salience, parietal memory, and contextual association networks following procedures described in ref. ²⁴. See Supplementary Fig. 6 for these individual-specific networks.

Differences in functional connectivity between inter-effector and foot, hand or mouth regions

Within each adult human participant, we calculated an inter-effector connectivity map as the Fisher-transformed correlation between the average time course of all cortical inter-effector vertices and the time course of every other vertex or voxel in the brain. We then repeated this procedure to calculate a connectivity map for the foot, hand and mouth areas.

To identify brain regions more strongly connected to inter-effector regions than to other motor regions, we computed the smallest positive difference in each voxel or vertex between inter-effector connectivity and any foot, hand or mouth connectivity. That is, we calculated (inter-effector – max[foot, hand, mouth]). This represents a conservative approach that only identifies regions of the brain for which the inter-effector regions are more strongly connected than any of the other motor areas.

Differences in functional connectivity among inter-effector regions

Within each adult human participant, as well as in the HCP group-averaged data, we computed a connectivity map for each of the three distinct inter-effector regions (superior, middle and inferior) as the Fisher-transformed correlation between the average time course of all cortical vertices in the two bilateral regions in each position and the time course of every other vertex or voxel in the brain.

To identify brain regions more strongly connected to one of the inter-effector regions than the other two, we computed the smallest positive difference in each voxel or vertex between that region’s connectivity and either of the other two regions’ connectivity. That is, we calculated (superior inter-effector – max[middle, inferior inter-effector]), (middle inter-effector – max[superior, inferior inter-effector]), and (inferior inter-effector – max[superior, middle inter-effector]). This represents a conservative approach that only identifies regions of the brain for which one inter-effector region is more strongly connected than either of the other two regions.

Functional connectivity with CON

Within each adult human participant, we calculated the functional connectivity between each of the foot, hand, mouth, and inter-effector regions and the CON. This was computed as the Fisher-transformed correlation between (1) the average time course across all vertices in the motor region and (2) the average time course across all vertices in the CON. We conducted paired t-tests across subjects comparing the inter-effector connectivity with CON against each of the foot, hand and mouth connectivity strengths, FDR-correcting for the three tests conducted.

We then calculated the functional connectivity between the inter-effector regions and every other large-scale cortical network in the brain (visual, auditory, salience, premotor, fronto-parietal, default mode, dorsal attention, language, contextual association, and parietal memory). The strength of connectivity between the inter-effector network and the CON was compared against the strength of its connectivity to each of these other networks using paired t-tests, FDR-correcting for the ten tests conducted.

Motor and CON network visualization

Visualization of network relationships was conducted using spring-embedded plots²³, as implemented in Gephi (https://gephi.org/). In each individual adult human participant, nodes were defined as congruent clusters of foot, hand, mouth, inter-effector, and CON networks larger than 20 mm². Pairwise connectivity between nodes was calculated as the Fisher-transformed correlation of their mean time courses. For visualization purposes, graphs were constructed by thresholding the pairwise node-to-node connectivity matrices at 40% density (the general appearance of the graphs did not change across a range of densities).

Functional connectivity with adjacent postcentral gyrus

In each adult human participant, we defined the pre- and postcentral gyri based on the individual-specific Brodmann areal parcellation produced by Freesurfer, which was deformed into fs_LR_32k space to match the functional data. Precentral gyrus was considered to be the vertices labelled as BA 4a and 4p, and postcentral gyrus was the vertices labelled as BA 3b and 2. BA 3a (fundus of central sulcus) was not considered for this analysis. Because the medial aspect of somatomotor cortex (corresponding to representation of the leg and foot) was always classified by Freesurfer as BA 4a, we defined the medial postcentral gyrus as the cortical vertices with y-coordinates farther posterior than the median y-coordinate of the foot region (from the network mapping above).

Within the participant’s precentral gyrus, we labelled vertices as representing foot, hand, mouth or inter-effector according to their labels from the network mapping procedure. We then partitioned the postcentral gyrus into foot, hand, mouth and inter-effector areas depending on which precentral region each vertex was physically closest to. Finally, within each partition (foot, hand, mouth and inter-effector) we calculated the average connectivity between the pre and postcentral gyrus as the Fisher-transformed correlation between the average time courses of all vertices in each area. We then conducted paired t-tests across subjects comparing the inter-effector connectivity with adjacent S1 against each of the foot/hand/mouth connectivity strengths with S1, FDR-correcting for the three tests conducted.

Functional connectivity with middle insula

In each adult human participant, we defined the middle insula based on the individual-specific Freesurfer gyral parcellation using the Destrieux atlas⁹¹, which was deformed into fs_LR_32k space to match the functional data. Middle insula was considered to be the vertices labelled as the superior segment of the circular sulcus of the insula or as the short insular gyrus. We then calculated the functional connectivity between each of the bilateral foot, hand, mouth, and inter-effector regions and the bilateral middle insula. We conducted paired t-tests across subjects comparing the inter-effector connectivity with middle insula against each of the foot, hand and mouth connectivity strengths, FDR-correcting for the number of tests conducted.

Functional connectivity with cerebellum

In each adult human participant, we calculated the functional connectivity between each of the foot, hand, mouth and inter-effector regions with each voxel of the cerebellum. Cerebellar connectivity strengths calculated this way were then mapped onto a cerebellar flat map using the SUIT toolbox⁹². Connectivity strengths were averaged within each of 27 atlas regions⁹³. For each region, we conducted three paired t-tests comparing inter-effector connectivity strength against foot, hand and mouth connectivity strength, FDR-correcting for the total number of tests conducted. Regions were reported if the inter-effector connectivity strength was significantly higher than the connectivity strength of all other motor regions.

Functional connectivity with putamen

In each adult human participant, we divided each unilateral putamen in each hemisphere into quarters by splitting it based on the median of its y (anterior-posterior) and z (dorsal–ventral) coordinates. We then calculated the functional connectivity between each of the foot, hand, mouth and inter-effector regions and each putamen quarter.

For each putamen division, we conducted paired t-tests across subjects comparing the inter-effector connectivity with that putamen division against each of the foot, hand and mouth connectivity strengths, FDR-correcting for the number of tests conducted. We reported divisions in which the inter-effector connectivity was significantly different from all three effector-specific connectivities.

Functional connectivity with thalamus

To investigate subregions of thalamus, we employed the DISTAL atlas v1.1⁹⁴, which contains a number of histological thalamic subregions identified by⁹⁵. This atlas was down-sampled into the 2-mm isotropic space of the functional data. Functional connectivity maps seeded from the foot, hand, mouth, and inter-effector regions in each adult human participant were computed, and mean connectivity values were calculated within each atlas region. The atlas specifies multiple subregions for many nuclei; these subregions were combined and treated as single nuclei for the purposes of connectivity calculation.

For each adult human participant, we averaged the connectivity seeded from the inter-effector regions and from each of the foot, hand and mouth regions across all voxels within each thalamic nucleus. For each thalamic nucleus, we conducted paired t-tests across subjects comparing the inter-effector with the mean of the foot, hand and mouth connectivity strengths, FDR-correcting for the number of thalamic nuclei tested.

Lag structure of RSFC

We used a previously published method for estimating relative time delays (lags) in fMRI data^34,96. In brief, for each session in each adult human participant, we computed a lagged cross-covariance function (CCF) between each pair of vertex or voxel time courses within the motor system and CON in the cortex. Lags were more precisely determined by estimating the cross-covariance extremum of the session-level CCF using three-point parabolic interpolation. The resulting set of lags was assembled into an antisymmetric matrix capturing all possible pairwise time delays (TD matrix) for each session, which was averaged across sessions to yield participant-level TD matrices. Finally, each participant’s TD matrix was averaged across rows to summarize the average time-shift from one vertex to all other vertices. Average time lag was then averaged across all vertices with each of the precentral gyrus foot, hand, mouth and inter-effector regions, and the CON.

We then conducted paired t-tests across subjects comparing (1) the mean lag in inter-effector regions against the mean lags in each of the foot, hand and mouth regions, and (2) the mean lag in CON regions against the mean lags in each of the foot, hand and mouth regions.

Macaque RSFC

We placed connectivity seeds continuously along area 4p in the left hemisphere of each macaque, as well as continuously running from the dorsal cingulate motor area to the rostral cingulate motor area in the dACC (area 24). See Supplementary Fig. 8a for medial cortex seed locations and Supplementary Fig. 8b for all functional connectivity maps from all medial seeds.

Structural MRI

Cortical thickness

Within each adult human participant, the map of cortical thickness generated by the Freesurfer segmentation was deformed into fs_LR_32k space to match the functional data. Precentral gyrus foot, hand, mouth and inter-effector regions were defined as above, and mean cortical thickness was calculated within each region. We then conducted paired t-tests across subjects comparing the inter-effector thickness against each of the foot, hand and mouth thicknesses, correcting for the three tests conducted.

Fractional anisotropy

White matter fibres tracked from separate areas of M1 using diffusion imaging quickly converge into the internal capsule and become difficult to dissociate. As such, we tested for fractional anisotropy differences in the white matter immediately below the precentral gyrus.

To calculate fractional anisotropy beneath the cortex, we first constructed fs_LR_32k-space surfaces 2 mm below each grey-white surface in adult human participants 1–3. To accomplish this, for each vertex on the surface, we computed the 3D vector between corresponding points on the fs_LR_32k pial and the grey-white surfaces, and we extended that vector an additional 2 mm beyond the grey-white surface in order to create a lower surface. We then mapped the fractional anisotropy values using the using the ribbon-constrained method, mapping between the grey-white and the 2 mm-under surfaces. The result is fractional anisotropy values mapped to a lower surface within white matter that is in register to the existing fs_LR_32k surfaces on which the functional data is mapped and the motor regions defined.

Precentral gyrus foot, hand, mouth and inter-effector regions were defined as above, and we calculated mean fractional anisotropy beneath each cortical region.

We then conducted paired t-tests across subjects comparing the mean fractional anisotropy beneath the inter-effector regions against mean fractional anisotropy beneath each of the foot, hand and mouth regions.

Myelin density

Within each adult human participant, we created vertexwise maps of intracortical myelin content following methods described in refs. ^81,97. Precentral gyrus was defined as above. Across participants, we found that baseline myelin density values (both in precentral gyrus and in the whole-brain myelin density map) varied wildly across participants in different datasets, likely based on differences in the T1- and T2-weighted sequences employed. Thus, for optimal visualization of results, in each participant we normalized the myelin density values by dividing the calculated vertexwise myelin densities in precentral gyrus by the mean myelin density across the whole precentral gyrus. Finally, precentral gyrus foot, hand, mouth and inter-effector regions were defined as above, and mean normalized myelin density was calculated within each region. We then conducted paired t-tests across subjects comparing the inter-effector myelin density against each of the foot, hand and mouth myelin densities, correcting for the three tests conducted.

Task fMRI

Movement task battery analysis

Basic analysis of the movement task battery data was conducted using within-participant block designs. To compute the overall degree of activation in response to each motion, data from each run was entered into a first-level analysis within FSL’s FEAT⁹⁸ in which each motion block was modelled as an event of duration 15.4 s, and the combined block waveform for each motion condition was convolved with a haemodynamic response function to form a separate regressor in a generalized linear model (GLM) analysis testing for the effect of the multiple condition regressors on the time course of activity within every vertex or voxel in the brain. Beta value maps for each condition were extracted for each run and entered into a second-level analysis, in which run-level condition betas were tested against a null hypothesis of zero activation in a one-sample t-test across runs (within participant). The resulting t-values from each motion condition tested in this second-level analysis were converted to Z-scores. Z-score activation maps were smoothed with a geodesic 2D (for surface data) or Euclidean 3D (for volumetric data) Gaussian kernel of σ = 2.55 mm.

Movement task battery winner take all

For each vertex within the broad central sulcus area, we identified the movement that produced the greatest activation strength (Z-score from second-level analysis, above) in that vertex, and we assigned that motion to that vertex.

Movement task battery curve fitting

For each vertex within precentral gyrus, we first computed its position along the dorsal–ventral axis of left hemisphere M1. This was done by identifying the closest point within the continuous line of points running down precentral gyrus (defined in ‘Seed-based functional connectivity’), and assigning that closest point’s ordered position within the line to the vertex.

For every movement, we then plotted that dorsal–ventral M1 position against Z-score activation in each vertex. We then fit two curves to each of these relationships. The first curve was a single-Gaussian model of the form:

Activation = a₁ × exp(−((position-b₁)/c₁)²).

The second curve was a double-Gaussian model of the form:

Activation = a₁ × exp(−((position-b₁)/c₁)²) + a₂ × exp(−((position-b₂)/c₂)²).

The a₁ and a₂ parameters in each model were constrained to be positive (to enforce positive-going peaks). Curve fitting was constrained to be conducted within the general vicinity of the activated area in order to avoid fitting negative activations observed in distant portions of M1. For lower extremity movements, this meant excluding the bottom third of M1; for upper extremity movements, the bottom third of M1 plus the medial wall; for face movements, the top third of M1.

Finally, we tested whether the one- or two-peak models better fit the data. This was done by conducting an F-test between the models, computed as:

where SSE represents the sum of squared errors from the model and df represents the degrees of freedom in the model.

The P value was computed from this F by employing the F-statistic continuous distribution function (fcdf.m) in Matlab and using (df_1peak – df_2peaks) and df_2peaks as the numerator and denominator degrees of freedom, respectively.

Movement task battery curve visualization

For each movement, the complete dorsal–ventral M1 position versus Z-score activation profile (from above) was visualized more clearly by fitting a LOWESS curve. These LOWESS curves recapitulated the two-peak activation fits while also revealing additional task responsive cortex.

Movement selectivity

Based on results from the above winner-take-all analysis, we identified the movement that was most preferred at the centre of each the three effector-specific (toe movement, hand movement and tongue movement) and inter-effector regions (abdominal movement, eyelid movement and swallowing). The centre-most movements were selected to avoid issues with spreading, overlapping activation near the borders of effector-specific and inter-effector regions. For every vertex within the precentral gyrus, we compared the strength of activation between the most preferred of the six movements at that vertex against the activation of the second-most preferred movements. The differences between these activation strengths was taken to be the movement selectivity of that vertex.

Movement coactivation

For each region among the six resting-state-derived foot, hand, mouth, and inter-effector regions in the precentral gyrus, we calculated the average activation within that region for each movement, producing a profile of motor activation strengths for that region. We also calculated the average activation within all CON vertices for each movement. To determine the degree to which various regions were coactive across movements, we then correlated each foot, hand, mouth and inter-effector cluster’s profile of activation strengths with that of all other clusters, and with that of the CON. Note: visualization of activation maps revealed some striping, suggesting that the Open and close mouth and the Bend L knee conditions were partially distorted by head motion; therefore, these conditions were excluded from analysis, although their inclusion did not change results.

Laryngeal motor mapping task analysis

As with the movement task battery, analysis of the laryngeal mapping task data was conducted using within-participant block designs. To compute the overall degree of activation in response to each motion, data from each run was entered into a first-level analysis within FSL’s FEAT⁹⁸ in which each motion block was modelled as an event of duration 15.0 s, and the combined block waveform for each motion condition was convolved with a hemodynamic response function to form a separate regressor in a GLM analysis testing for the effect of the multiple condition regressors on the time course of activity within every vertex or voxel in the brain. Beta value maps for each condition were extracted for each run and entered into a second-level analysis, in which run-level condition betas were tested against a null hypothesis of zero activation in a one-sample t-test across runs (within participant). The resulting t-values from each motion condition tested in this second-level analysis were converted to Z-scores.

Note that the laryngeal mapping data was not included in the movement task battery analysis because it was collected on a different scanner with a different sequence, and so would not be directly comparable.

Action control and coordination task analyses

Analysis of the action control task was conducted using within-participant event-related designs. For each separate run, a GLM model was constructed in FEAT⁹⁸ in which separate regressors described the initiation of (1) planning and (2) execution of each type of movement (4 movements × 4 limbs). Each regressor was constructed as a 0-length event convolved with a canonical haemodynamic response, and beta values for each regressor were estimated for every voxel in the brain. These beta value maps for each condition were thus computed for each run and entered into a second-level analysis, in which a t-test across runs contrasted the run-level planning betas against the run-level execution betas.

Human direct electrocortical stimulation site mapping

Each stimulation location reported in ref. ⁵¹ was separately mapped into the MNI-space Conte69 atlas pial cortical surface⁸⁰ by identifying the vertex with the minimal Euclidean distance to the stimulation site’s MNI coordinates. Movements resulting from each site were classified as ‘lower extremity’, ‘upper extremity’ or ‘face’ and coloured accordingly (although no lower extremity movements were reported in the displayed left hemisphere).

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.