Heterochronous maturation of superficial and deep prefrontal cortex myelin balances developmental plasticity and stability

Project Lead

Valerie J. Sydnor

Project Mentor

Beatriz Luna

Collaborators

Daniel Petrie, Shane D. McKeon, Alyssa Famalette, Will Foran, Hoby Heatherington, Finnegan J. Calabro

Project Start Date

September 2023

Current Project Status

Manuscript in preparation

Datasets

Luna-7T

Github Repository

https://github.com/LabNeuroCogDevel/corticalmyelin_maturation

Project Directories

The project directory on rhea is: /Volumes/Hera/Projects/corticalmyelin_development

The directory structure on rhea is as follows below:

BIDS: BIDS-compliant MP2RAGE data
BIDS/derivatives: Processed image derivatives including unicort-corrected UNI and T1 maps, volumetric R1 maps, and surface measures
Dicoms: Organized dicoms for BIDSifying
EEG: EEG aperiodic activity spreadsheets
Figures: Manuscript figure panels
Maps: Brain maps galore! HCP-MPP atlas, S-A axis, Neurosynth z-score maps, AHBA dense gene expression, cytoarchitectural variation, MWF imaging, and templateflow templates
code: local clone of corticalmyelin_maturation repo
gam_outputs: Gam-derived statistics (developmental_effects, sex_effects, sensitivity_analyses, eeg_associations, cognitive_associations)
sample_info: sample demographics, behavioral data, imaging QC spreadsheets, and final participant lists
software/: Directory with project software installs and containers

The project directory on the PSC is: /ocean/projects/soc230004p/shared/datasets/7TBrainMech

Conference Presentations

Sydnor, V.J., Petrie, D., McKeon, S.D., Famalette, A., Foran, W., Calabro, F.J., Luna, B. Gradients of myelin growth across layers of human prefrontal cortex: Balancing adolescent circuit plasticity and stability. Poster presented at the Frontal Cortex Gordon Research Conference, Holderness, NH, Aug 2024 best poster award

CODE DOCUMENTATION

Below is a detailed walk through of all of the code run for this project, including code for image processing, statistical analysis, results derivation, and visualization! The entire analytic workflow is described, and links to the corresponding code on github are provided. Please feel free to reach out to Valerie Sydnor (sydnorvj@upmc.edu) with any questions.

Image Processing Pipeline: Structural and R1 Data Processing

Organization and Curation of BIDS

Structural data dicoms from the Luna-7T Brain Mechanisms project were first organized into {subject}/{session} organization (required for heudiconv) with the script /BIDS/BIDS_heudiconv/organize_dicoms_forBIDS.sh. These dicoms were then converted to niftis according to the Brain Imaging Data Structure via heudiconv (version 0.13.1). Specifically, data were BIDSifyed using the /BIDS/BIDS_heudiconv/7TBrainMech_MP2RAGE_heuristic.py heuristic and the /BIDS/BIDS_heudiconv/BIDSify_7T.sh script.

In order to check the acquisition parameters for the structural MP2RAGE data acquired for this study, Curation of BIDS on Disk (CuBIDS) was run with /BIDS/CuBIDS/cubids_acquisition_group.sh.

Processing of MP2RAGE-Derived UNI Images

MP2RAGE-derived T1w Uniform (UNI) images were denoised off of the scanner and then corrected for residual B1+ transmit field biases using SPM’s UNICORT algorithm. UNICORT was applied in /image_processing/B1+_transmitfield_correction/unicort/unicort_inhomogeneitycorrection_UNI.m using the configuration parameters specified in /image_processing/B1+_transmitfield_correction/unicort/unicort_configparams_UNI.m. UNICORT outputs were organized and the B1-corrected image was put into the BIDS directory for use with BIDS apps via /image_processing/B1+_transmitfield_correction/unicort/organize_unicort_files.sh.

Following UNICORT, corrected UNI images were run through FreeSurfer to generate person-specific cortical surface reconstructions to map R1 data to. Images were processed with FreeSurfer’s longitudinal analysis stream (FS version 7.4.1) using a containerized FreeSurfer BIDS App. FreeSurfer was run on the Pittsburgh Super Computer by submitting subject-specific jobs via slurm/sbatch. Jobs were submitted with the script /image_processing/longitudinal_freesurfer/freesurfer_submitjobs_PSC.sh, which calls /image_processing/longitudinal_freesurfer/longitudinal_freesurfer_call.sh. Cross-sectional, template, and longitudinal steps of recon-all were run for all subjects (and sessions). Job completion was checked with /image_processing/longitudinal_freesurfer/check_freesurferjobs_complete.sh.

Processing of MP2RAGE-Derived R1 Maps

MP2RAGE-derived quantitative T1 maps were corrected for B1+ transmit field biases using UNICORT by running /image_processing/B1+_transmitfield_correction/unicort/unicort_inhomogeneitycorrection_T1map.m, which uses /image_processing/B1+_transmitfield_correction/unicort/unicort_configparams_T1map.m as the config file. Corrected quantitative T1 maps were renamed, organized, and cleaned up with /image_processing/B1+_transmitfield_correction/unicort/organize_unicort_files_T1map.sh.

UNICORT-corrected quantitative T1 maps were next converted into volumetric R1 maps (s^-1). R1 maps were calculated in /image_processing/R1_surfacemaps/R1_volumetricmaps.sh by coverting the T1 readout from ms to s and then calculating R1 = 1/T1. Volumetric R1 maps were then projected to participant-specific cortical surfaces on the PSC by running /image_processing/R1_surfacemaps/R1_vol2surf_submitjobs_PSC.sh. This script calls /image_processing/R1_surfacemaps/voltosurf_projection_nativefsaverage.sh for every session in order to project R1 data to the surface at 11 projection fractions: 0 (GM/WM boundary), 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 (pial boundary).

Calculationof Anatomical and R1-based Measures

After running FreeSurfer and projecting depth-dependent R1 data to the cortical surface, anatomical metrics (cortical thickness, cortical curvature) and depth-dependent R1 metrics could be quantified in individual cortical regions using atlases of interest. This was accomplished by running /surface_metrics/fstabulate/fstabulate_R1anat_submitjobs_PSC.sh, which ultimately uses the shell and python scripts in /surface_metrics/fstabulate. The workhorse script directly called via fstabulate_R1anat_submitjobs_PSC.sh is /surface_metrics/fstabulate/collect_stats_to_tsv.sh; this script parses user-defined inputs and warps fsaverage parcellations to native space, runs FreeSurfer’s qcache + mris_anatomical_stats + mri_segstats to get regional metrics of interest, and converts metric stats files to tsvs. After fstabulate_R1anat_submitjobs_PSC.sh completed, /surface_metrics/surface_measures/collate_study_anatomicalstats.py was run to combine session-specific tsv metric files into study parquet files with all participants’ anatomical and depth-dependent R1 data in atlases of interest. These parquets are used to generate final input dataframes for statistical modeling.

Image Processing Pipeline: EEG Aperiodic Activity

Creation of an EEG Electrode Cortical Surface Atlas

This project examines relationships between cortical R1 and EEG-derived measures of aperiodic activity. In order to most accurately compare aperiodic activity in individual electrodes with R1 in the cortical area covered by the electrodes, EEG electrode positions were mapped to the fsaverage cortical surface and then warped to participant-specific surfaces. First, spherical ROIs were generated from EEG electrode coortdinates in MNI space. Subsequently, ROIs were mapped to the fsaverage surface, turned into a cifti label file, dilated, and converted to freesurfer-compatible annot files. This was done through /surface_metrics/EEGelectrode_atlas/EEGatlas_fsaverage_annots.sh, which internally calls /surface_metrics/EEGelectrode_atlas/EEGatlas_fsaverage_giftis.py.

Quantification of Aperiodic Activity

EEG aperiodic activity measures used in this project were initially derived in McKeon et al., 2024, DCN. Code for preprocessing and parameterization of aperiodic activity are available here: https://github.com/LabNeuroCogDevel/7T_EEG/tree/main/Aperiodic_MRS_Development.

Preprocessing was done through 01_Aperiodic_Preprocessing.sh and the aperiodic component of the power spectrum was parameterized using FOOOF via 02_runFOOOF.py, 03_ExtractFOOOFmeasures.py, and 04_CreateFOOOFdataframes.R. Raw aperiodic activity measures were then cleaned (exlusion of electrodes with high spectral fit errors, low R2, and outlier values) and averaged across nearby electrodes for this project in /image_processing/EEG_aperiodic/fooof_aperiodicactivity.R. The main analysis used the EEG electrode cortical surface atlas described above. A sensitivity analysis was additionally performed using aperiodic measures averaged in Brodmann regions (see /image_processing/EEG_aperiodic/fooof_aperiodicactivity_brodmann.R).

Sample Construction

Construction of Final Sample

The initial developmental dataset considered for inclusion in this project was comprised of 264 longitudinal sessions (with correctly acquired, i.e. CuBIDS non-variant, UNI and quantitative T1 images) from 159 healthy participants. Sessions were iteratively excluded for:

failing visual quality control of UNI and R1 images (44 sessions/18 subjects excluded)
being collected within 6 months of another session that passed QC (4 sessions/0 subjects excluded)
having a whole-brain mean R1 > 4 SD above the group mean (1 session/1 subject excluded)

This resulted in a final sample of 215 longitudinal sessions collected from 140 individuals, finalized in /sample_construction/finalsample_7Tmyelin.Rmd.

Generation of Final Input Dataframes

After determining the final study sample (final subject and session lists), /surface_metrics/surface_measures/extract_depthdependent_R1.R was run in order to extract depth-specific R1 measures from cortical atlas regions (HCP-MMP atlas and EEG electrode surface atlas) for the final study sample. The extract_depthdependent_R1.R script parses parquets (generated by collate_study_anatomicalstats.py) using the flexible /surface_metrics/surface_measures/extract_surfacestats.R function and produces RDS files with R1 measures for the final sample that are used for statistical modeling.

In addition to generating input dfs for the final study sample, a final list of cortical regions to include in analyses was generated based on the signal to noise ratio of R1 in each region. This was accomplished in /surface_metrics/surface_measures/extract_R1_SNR.R which computes the average (across-depth) R1 signal to noise ratio in each glasser region and identifies low SNR regions in frontal and temporal poles to exclude from analyses (24 regions excluded). The R1 SNR metric used here was calculated via FreeSurfer’s mri_segstats (in /fstabulate/collect_stats_to_tsv.sh), which computes SNR as mean R1/standard deviation of R1 in a region.

Statistical Analysis: Model Fitting with GAMMs

This project uses generalized additive mixed models (GAMs with random intercepts) to characterize developmental changes in R1 as well as quantify associations between R1 and both EEG and cognitive measures. An array of GAM(M) fitting functions from /gam_models/gam_functions.R were used throughout analyses, all which use gam from the mgcv package. gam_functions.R contains the following functions:

gam.statistics.smooths: A function to fit a GAM(M) with one smooth (s(x)) plus optional id-based random effect terms and save out model statistics, spline-derived characteristics, fitted values for the smooth from a prediction df, zero-centered smooth fits, and the instantaneous first derivative of the smooth function. This was used to fit most developmental models in main and sensitivity analyses
gam.agebysex.interaction: A function to fit a GAM(M) with a factor-smooth interaction (s(x) + s(x, by = y) and save out interaction statistics. This was used to fit age-by-sex interaction models
gam.tensorproduct.interaction: A function to fit a GAM(M) with an explictly modeled tensor product interaction term (s(x) + s(y) + ti(x,y)) and optional random intercepts or random intercepts + slopes. This was used to fit age-by-depth interaction models
gam.linearcovariate.maineffect: A function to fit a GAM(M) with a main smooth plus a linear covariate of interest (s(x) + y) and save out statistics for the covariate term. This was used to test for main effects of sex
gam.covariatesmooth.maineffect: A function to fit a GAM(M) with a main smooth plus a smooth for a covariate of interest (s(x) + s(y)) and save out statistics for the covariate smooth. This was used for quantifying relationships between R1 and EEG/behavior
gam.factorsmooth.interaction: A function to fit a GAM(M) with a main smooth plus a factor-smooth interaction for a second smooth term (s(x) + s(y) + s(y, by = z)) and save out main effects and interaction results. This was used to test for superficial/deep depth interactions

These functions were applied in a series of scripts included in the /gam_models directory to model developmental effects (/gam_models/fit_ageGAMs_glasserregions_bydepth.R and /gam_models/fit_agedepth_interactionGAMs_glasserregions.R), to conduct sensitivity analyses (/gam_models/fit_sensitivityGAMs_glasserregions_bydepth.R), to examine relationships with EEG (/gam_models/fit_eegGAMS_electrodeatlas_bydepth.R), and to characterize cognitive links (gam_models/fit_cognitionGAMs_glasserregions_bydepth).

Statistical Analysis: BrainSmash

This project quantifies across-region correlations between brain maps, including whole-brain correlations between R1 and myelin-related measures (MBP, MWF) and frontal lobe-only correlations between development maps and indexes of cortical hierarchy (S-A axis, cytoarchitectural variation gradient). To determine the statistical significance of these spatial correlations, autocorrelation preserving null models were constructed using BrainSmash. The procedure to compute spatial map correlation p-values with BrainSmash nulls requires 1. a parcellated distance matrix for your cortical atlas of interest, and 2. parcellated brain maps (empirical) to use for creating spatial autocorrelation-preserving surrogate maps (nulls). The distance matrix (1) was created by first deriving a vertex-wise geodesic distance matrix for left and right cortical surfaces with /results/brainsmash_nulls/compute_dense_geodesicdistmat.py and then parcellating it with the HCP-MMP atlas in /results/brainsmash_nulls/compute_parcellated_geodesicdistmat.py; a frontal lobe-only distance matrix was generated in /results/brainsmash_nulls/mask_parcellated_distmat_frontallobe.R. Surrogate maps (2) were created for whole-brain R1 /results/brainsmash_nulls/create_R1_surrogatemaps.py, the frontal lobe S-A axis /results/brainsmash_nulls/create_SAaxis_surrogatemaps.py, and the frontal lobe cytoarchitectural variation gradient /results/brainsmash_nulls/create_bigbrain_surrogatemaps.py.

Results: Quantitative Imaging of Myelin Across Cortical Regions and Depths

Cortical Myeloarchitecture Characterization

R1 was validated in our sample as a quantitative 7T imaging measure capable of capturing variation in myelin levels across cortical regions and cortical depths. Characterization of R1-based myeloarchitecture was conducted in /results/R1_corticalmyelin_anatomy/R1_cortical_myeloarchitecture.Rmd which shows spatial variation in R1 and its alignment with myelin basic protein gene expression, myelin water fraction imaging, and the S-A axis as well as depth-dependent profiles of R1. A rendered version of this R markdown is available here!

Myelin Map Generation

The above code utilizes a myelin basic protein gene expression map, a myelin water fraction imaging map, and the S-A axis. These brain maps were obtained as follows:

Myelin basic protein gene expression was obtained from vertex-level gene expression data generated from the AHBA by Wagstyl et al., 2023, elife. The script /results/R1_corticalmyelin_anatomy/myelin_maps/create_AHBAmagicc_genemap.sh gets MBP gene maps by interfacing with /results/R1_corticalmyelin_anatomy/myelin_maps/create_AHBAmagicc_genemap.py to pull gene expression and gene info csvs and produce a parcellated MBP cifti.
Myelin water fraction imaging was quantified in HCP-MMP regions by projecting the MNI-space myelin imaging template from Liu et al., 2019, Journal of Neuroimaging to the cortical surface and parcellating it to produce a MWF cifti. This was done in /results/R1_corticalmyelin_anatomy/myelin_maps/parcellate_MWFmap.sh, which calls /results/R1_corticalmyelin_anatomy/myelin_maps/parcellate_MWFmap.py.
Sensorimotor-association (S-A) axis ranks for HCP-MMP regions were obtained from the archetypal S-A axis repo

BrainSmash Significance

The significance of correlations between regional R1 and myelin maps was assessed using autocorrelation-preserving R1 surrogates generated (detailed above) via BrainSmash. Visualization of example surrogates and significance testing was executed in /results/brainsmash_nulls/R1_cortical_myeloarchitecture_brainsmash.Rmd which can be viewed here.

Results: Heterochronous and Hierarchical Myelin Maturation in the Human PFC

Charting the Development of Myelin Across Cortical Layers

Gradients of myelin maturation across depths of the human frontal cortex were investigated in /results/R1_depth_development/R1_depthdevelopment.Rmd. This markdown file examines trajectories of myelin development across cortical depths, elucidates regional differences in the rate and timing of myelin maturation across depths and their hierarchical spatial embedding, and demarcates cortical zones that show similar depth-dependent developmental profiles. A rendered version of this R markdown is available here! R1_depthdevelopment.Rmd was paired with /results/R1_depth_development/R1_sexdifferences.Rmd, which is knitted here, to examine whether developmental trajectories differed by sex.

BrainSmash Significance

The significance of correlations between developmental maps and hierarchically-expressed variability in cortical properties and cytoarchitecture was assessed using autocorrelation-preserving S-A axis surrogates and cytoarchitectural variation surrogates (detailed above). Visualization of example surrogates and significance testing was executed in /results/brainsmash_nulls/R1_depthdevelopment_brainsmash.Rmd, which can be viewed here.

Results: Sensitivity Analyses

A series of sensitivity analyses were undertaken to ensure that developmental results were not driven by individual differences in sex, structural data quality, or cortical architecture. Sensitivity GAMMs were all fit in /gam_models/results/fit_sensitivityGAMs_glasserregions_bydepth.R and then examined separately in /results/R1_sensitivity_analyses as explained below.

Biological Sex

Developmental models were covaried for sex and developmental results were examined in /results/R1_sensitivity_analyses/R1_sensitivity_sex.Rmd which is knitted for visualization here.

Structural Data Quality

Developmental models were covaried for euler number and developmental results were examined in /results/R1_sensitivity_analyses/R1_sensitivity_euler.Rmd which is knitted here. The euler number was computed by FreeSurfer’s mris_anatomical_stats, which is run as part of the R1 and Anatomical Measure Extraction pipeline.

Cortical Thickness

Developmental models were covaried for regional cortical thickness and developmental results were examined in /results/R1_sensitivity_analyses/R1_sensitivity_corticalthickness.Rmd which is knitted here. Cortical thickness was computed by FreeSurfer’s mris_anatomical_stats in every HCP-MMP region for every participant.

Cortical Curvature

Developmental models were covaried for regional curvature and developmental results were examined in /results/R1_sensitivity_analyses/R1_sensitivity_curvature.Rmd which is knitted here. Curvature represents FreeSurfer’s mean curvature measure, which varies according to folding properties and gyral versus sulcal location.

Cortex Fraction (PVE)

Developmental models were covaried for regional cortex fraction (quantified at each depth) and developmental results were examined in /results/R1_sensitivity_analyses/R1_sensitivity_cortexfraction.Rmd which is knitted here. The cortex fraction was calculated in each region at every cortical depth by first estimating the cortex tissue fraction in volumetric space (i.e., the percent of each voxel occupied by cortex) using FreeSurfer’s mri_compute_volume_fractions, which was applied via /image_processing/cortexfraction_surfacemaps/cortexfraction_volumetricmaps.sh. Volumetric cortex fraction maps were then projected to participant-specific cortical surfaces at varying depths in /image_processing/cortexfraction_surfacemaps/cortexfraction_vol2surf.sh and the average cortex fraction in each region (at each depth) was computed with /surface_metrics/fstabulate/fstabulate_cortexpve_submitjobs_PSC.sh followed by /surface_metrics/surface_measures/collate_cortexpve_stats.py.

Results: Cortical Myelin Supports Mature Cortical Dynamics

Associations Between Superficial and Deep R1 and EEG Aperiodic Activity

To study the impact of intracortical myelination on the expression of mature cortical dynamics, relationships between superficial and deep cortex R1 and EEG-derived aperiodic parameters were delineated in /results/R1_EEG_relationships/R1_aperiodicactivity_associations.Rmd. This markdown presents developmental changes in the aperiodic slope in four cortical areas, quantifies depth-specific relationships between R1 and the aperiodic slope, and tests for depth interactions in R1-EEG relationships. A rendered version of this R markdown is available here! A sensitivity analysis mapping R1 to EEG signals via Brodmann regions was conducted in /results/R1_EEG_relationships/R1_aperiodicactivity_associations_brodmann.Rmd based on this paper.

Results: Cortical Myelin Impacts Learning Rate and Cognitive Processing Speed

Associations Between Cortical R1 and Cognitive Abilities

After characterizing how cortical myeloarchitecture is refined with age in the adolescent prefrontal cortex, we investigated the impact of cortical myelin levels on cognitive abilities. Specifically, we examined relationships between cortical R1 and measures of learning rate and cognitive processing speed derived from each individual’s performance on a 2-stage sequential learning task (featuring spaceships and aliens!). Learning rates were obtained by fitting a 7-parameter reinforcement learning model to trial-level data on the sequential learning task; a seperate learning rate was obtained for each of the 2 stages of the task. Cognitive processing speed was proxied by the average across-trial response time on each of the 2 stages of the task. Findings regarding the influence of cortical R1 on cognitive measures are presented in /results/R1_cognition_relationships/R1_cognition_associations.Rmd. This markdown explores developmental changes in learning rate and processing speed on both stages of the task and characterizes depth-specific and across-depth relationships between R1 and cognitive measures. It furthermore includes Neurosynth decoding of brain regions where myelin is linked to learning rates and a sensitivity analysis using anti-saccade and visually-guided saccade data that demonstrates that prefrontal myelin is important for efficient processing only when engaging higher-order cognitive abilities. A rendered version of the html is viewable here