Guide to the code behind "Psychedelics disrupt activity propagations in the default mode network of humans and mice"
This document outlines the steps and methods used in the project. Below is a structured guide for image processing, derivations, and analyses. All image processing was run in a Linux environment using a SLURM cluster for high-compute jobs. In this context, sbatch refers to submitting a job to the SLURM job scheduler. Note that fmriprep and xcpd calls utilize their singularity images, which need to be installed locally. In addition to fmriprep and xcpd, a TON of this code leverages tools built by other people. A summary table of prerequisite code/dependencies are available at the end of this markdown, as is a short description of our Demo script.
I'll occasionaly refer to study 1, study 2, and study 3. Study 1 is our MDMA sample, 2 is psilocybin, and 3 is LSD/mice.
If you want to find the code behind specific figures, I would use cntrl-f (apple-f for some of us) for "2A", "3C", etc. within this README.md. It should point you to the script that actually generates the figure from the processed data.
This is non-comprehensive, but you might find sbatch_OpFl.sh (a parent script that sequentially launches individual processing scripts on scan sessions) useful for further orienting yourself to order-of-operations. If you are interested in cortical propagations broadly, you might also find the replication guide from our previous paper useful.
This is the fmriprep call used. These derivatives become the inputs of TSNR (temporal signal to noise ratio) scripts as well as xcpd (eXtensible connectivity pipeline) scripts. It is written to run on a SLURM cluster (Simple Linux Utility for Resource Management): - full_fmriprep.sh
After fmriprep has ran, this script is the first in deriving a TSNR mask. It simply masks out the brain from a nifti. This provides an out-of-brain comparison for the magnitude of signal fluctuations: TSNR_mask_0_antimask.sh
This uses derivatives to calculate SNR for each individual subject: TSNR_mask_1_ExtractSNR_subjectwise.py
This script then combines average SNR maps across participants for a global average map. TSNR_mask_2_Combine_SNRmaps.py
Next is a quick downsampling of this mask. We'll need it in fsaverage5 space for regularized non negative matrix factorization. Here's that script, DS_surf_TSNR_fs5.sh
Next, this script combines the averaged left and right hemisphere TSNR maps. The bottom 15 percent of voxels (lowest TSNR) are masked out for all subsequent analyses. Note that psilocybin TSNR maps indicated that psilocybin study data all met the 15th percentile calculated from the MDMA study, likely due to use of multiecho. The same mask is applied to that dataset for equivalence: TSNR_mask_3_avSNR-mw_to_label
Now all we need to do is downsample the mask to appropriate resolution for optical flow, fsaverage4. Here's that script. Shoutout to Connectome Workbench (wb_command) for making this vanilla and simple. Speaking of connectome workbench, you can use their GUI to visualize the TSNR mask t create the constituent components for Figure S1.
That's it for SNR. Now, let's implement the xcpd sbatching script (useful for head-motion-prone scan protocols, runs on fmriprep outputs): sbatch_xcpd.sh
The WashU crew was kind enough to send their data over already processed. More information on their dataset is available here.For full transparency, the code used to download the images and their associated files are linked below:
Download psilocybin data, text files: DL_psilo_txts.sh
Download psilocybin data, framewise displacement files - DL_psilo_fd.sh
Download psilocybin data: neuroimages - DL_psilo.sh
Now we'll downsample both datasets to optical flow resolution. It's a bit computationally intensive (prohibitively intensive at the moment), so we'll downsample these to fsaverage4 resolution rather than the fsaverage5 used for NMF. It uses the same basic code as downsampling to fsaverage5, but with fs4 as the target surface. Here is DS_surf_ts_mdma.sh for study 1 data. For study 2 (psilocyin data), we'll first need to run RS_mask_psil, which pulls just the resting state images out of the downloaded psilocybin data. After doing so, we can proceed with the equivalent downsampling scrip DS_surf_ts_psil.sh for study 2 data.
OK, now we have all our processed data. To make it suitable for NMF and optical flow, we are going to motion-mask it. This goes a step further than normal temporal censoring: we are only interested in sequences of neuroimages uninterrupted by head motion. So rather than simply dropping high-motion frames, we'll only retain sequences of at least 8 TRs that are uninterrupted-by-motion. We'll save out the indices of these sequences for optical flow, as we're only interested in activity movement within clean segments (i.e., if frames 9 and 10 have high motion and we retain frames 1-8 and 11-20, we don't want to try and estimate activity movement between 8 and 11, just within the clean segments.) The script that takes care of this is MotMask.m. It also converts files to .mgh internally, which are easier for matlab to work with. As will be the case for many scripts below, an equivalent script exists for the psilocybin study. In general, the only differences are the input filepath and output file path. For this script in particular, some additional lines of code exist. As we get further in the pipeline, the file derivatives will be more equivalent, requiring fewer differences between MDMA study and Psilocybin Study scripts.
Nice. Done with fMRI preprocessing!
The first thing to do is to list out sessions to process. This will enable us to loop over sessions to-be-processed within python: sesh_lister.sh. Note that several different sesh listers exist for different drugs and pre vs. post conditions, but all use the same framework. More information on this dataset is available here.
Next, we'll obtain a group-level mask so that no mouse has pixels included that other mice don't. We'll combine this with downsampling to saveout a mask that is applicable to processing in downsampled space. As for human data, there are two different downsampled resolutions, with the greater resolution being provided to NMF and the lesser resolution being provided to optical flow for computational tractability.
Obtain group-level mask at factor of 2 downsample (NMF) - Group_Mask_and_DS_oneHalf.py
Obtain group-level mask at factor of 6 downsample (OpFl) - Group_Mask_and_DS_oneSixth.py
Just for cleanliness, these scripts are ran independently of downsampling the data for use in NMF and optical flow (even though internal calculations are matched). Here are those scripts. Note they also include gaussian smoothing for SNR purposes.
Downsample by factor of 2 for NMF - BP_Smooth_oneHalf.py
Downsample by factor of 6 for Optical flow - DS_Smooth_oneSixth_Drug.py (adapted for each drug with file path and extensions)
By the end of this section, you might notice that the mouse data is generally a little less cumbersome to process. That pattern will continue.
Before evaluating how DMN function is altered, we have to derive our DMN definition. This is done with regularized non-negative matrix factorization (NMF). See this link for the original paper. In case it's useful context, we've used this approach fairly regularly across half a dozen publications/preprints now: Cui et al., 2020 Pines et al. 2022 Shanmugan et al., 2022 Keller et al., 2023 Hermosillo et al., 2024, etc. etc.
Downsampling to fs5 First, the ciftis are downsampled to fsaverage5 to play nicely with NMF. That downsample script, DS_surf_ts_mdma_fs5.sh, is available here
NMF steps Alright, now that we have data in fsaverage5, we'll run it through the NMF pipeline. Basically the steps are to 1) prepare the data for internal NMF processing, 2) create a few candidate solutons using different combinations of brain scans 3) determine the optimal solution of the candidate solutions 4) convert the solution from a matrix to a proper brainmap for further use. Note that both human and mouse NMF utilized extra scans to for greater data volume and reliability. Also note you'll need the helper scripts that come with this code. If you are internal to our lab, you can find the helper scripts at:
/oak/stanford/groups/leanew1/users/apines/scripts/PersonalCircuits/scripts
If you are external to our lab, check out Hongming Li's repository.
NMF script 1: this script will pull some information about the cortical surface you are using (fsaverage5 here) as well as the SNR mask generated in 1A. Here's the link to Step_1_CreatePrepData.m.
NMF script 2: this script will launch a bunch of parallel matrix decomposition jobs (into a SLURM queue). The internal workflow is to concatenate the functional timeseries data into a single giant matrix, and to decompose the matrix into spatiotemporal components with NMF. Regularization parameters were set a priori, in accordance with previous literature. k=4 components was set in accordance with a nice coarse network solution we previously derived in an independent sample. The reason for using a coarse solution is outlined in the manuscript, but essentially it boils down to inclusive masking to capture propagations entering and exiting the core DMN. Here's the link to Step_2nd_ParcellationInitialize.m
NMF script 3: this script picks out the optimal NMF solution generated by the previous script. Optimal is defined as consistency with other solutions and operationalized via normalized cuts. This step runs rather quickly and doesn't require a bunch of job submissions and factorization of giant concatenated matrices. Here's the link to Step_3rd_SelRobustInit.m.
NMF script 4: this script just takes the derived solution and converts it to connectome workbench style file formats (giftis). Simple as. Here's the link to Step_4_Extract_Group_Atlas.m. This is how you generate the files for figure S2, which are visualized in the connectome workbench GUI
Downsampling to fs4 Downsampling the NMF solution to fsaverage4 for optical flow: this step is likely familiar by this point. To downsample fsaverage5 resolution networks to fsaverage4, use DS_surf_Networks_fs5tofs4.sh.
func.gii to .mat The last step is to convert this brainmap back into a plain ol' matrix. This makes matlab less fussy down the road when loading it back in. You can use Netgiis_2_mat.m is pretty straightforward though.
Thankfully, we can use the same code to derive functional networks in mice with some slight adaptations. The adaptations are pretty in-the-weeds but pretty simple at a conceptual level. Basically NMF uses an abutment matrix for spatial regularization, where all points in space are arranged in an adjacency matrix with all other points in space. Neighboring voxels are assigned a "1" to reflect sharing a physical edge; this adjacency matrix is leveraged by regularization ultimately making them more likely to be assigned to the same network. The main shift from the original code is from a volumetric image to a 2D brain surface captured via cranial windows. So essentially we just have to remove a dimension. If you look at the equivalent step 1 script, which we will return to below, it calls constructW_flat.m instead of the original script (the original script has to perform across a z-dimension, which is obviated here). This script also has some visualization checks to make sure everything is running as intended. Everything else is extremely similar to the human-brain-adapted NMF we'vea already done by this point. For compelteness, here's an equivalent set of instructions to the human section above.
NMF script 1: this script will pull some information about the cortical surface you are using (cranial window field-of-view here). You can find it here.
NMF script 2: this script will launch a bunch of parallel matrix decomposition jobs (into a SLURM queue). The internal workflow is to concatenate the functional timeseries data into a single giant matrix, and to decompose the matrix into spatiotemporal components with NMF. Regularization parameters did need some adaptation to Ca2+/mouse data from fMRI, but all optical flow and statistical analyses were ran after crystallizing these paramaters on the basis of correspondence with prior accounts of the mouse DMN (see Whitesell et al for a good example). k=13 components was set because 1) the DMN takes up a smaller portion of cortical area in mice and 2) the medial wall is not readily visible from the cranial window, meaning that in addition to the DMN occupying less area of mouse cortex we also are not capturing a chunk of the area it does occupy on the medial wall. For those reasons, greater granularity is needed to pull out specific DMN components. The script can be found here.
NMF script 3: Again, this script picks out the optimal NMF solution generated by the previous script and itself runs rather quickly. The script can be found here.
NMF script 3.5: This one didn't get it's own arabic numeral because it's just a plotting script. Nice to triple check that everything is working as intended. Script can be found here. Panels for Figure S3 generated within.
Great. Now we just have to downsample the mouse DMN to optical flow resolution. It's in python to mirror the same downsampling used for mouse time series data. The script to use for this step is DS_Smooth_NMFtoOF_resolution.py.
Ok, these are just some simple checks to make sure our derived DMN maps onto previous accounts of the DMN. Specifically, we'll use the spin test to see if our DMN is more aligned with previously published maps than 10,000 spun nulls. The spatial permutation approach allows us to better-account for autocorrelation in brain maps. There's a subfolder in the scripts directory for 1D, titled DMN_validate.
We'll start with the ROIs from Goldstein-Piekarski et al., Biological Psychiatry 2021. More specifically, these ROIs are available on our lab's github page.
First, we have to convert them from volumetric to fs5 so we can use the spin test/have them in the space space as the DMN delineation we just conducted. The script to use for this procedure is biol_psych_rois_to_fs5.sh.
In parallel, we can "spin" our DMN 10,000 times to get our null maps. That script explicitly uses the spin test repo linked. The script is named spin_DMN.m.
After the spins are complete (they might take a while, consider sbatching them with something like sbatch_matlab (replacing the script ran in matlab with spin_DMN) so you don't have to sit and wait. The equivalent of sbatch is qsub if you're on an SGE (sun grid engine) system. You can use this script to formally evaluate above-chance localization of our current DMN boundaries to those derived in biol. psych. 2021, via t-test. There's a baby script in the main .rmd that will remain to generate figures/stats after running these, but this is the majority of the computation. If you want to spot check the tiny chunk of code in R, you can skip ahead to Stats_n_Viz.md: control-F for '###### spin tests' to find where these analyses were run.
The other previously established DMN map we evaluate is from the epic 35-figure Yeo et al., 2011. If you don't already know, his lab runs a great repo, which is where I pulled the .annot files from. Because our NMF-DMN maps are already spun 10,000x from the previous comparison, we don't have to respin to generate our null distributions. Just loop over the spins and run the same t-tests in yeo7_overlap.m.
That's it for step 1! As noted above, data were disimilar upon receipt, but are gradually standardized into similar filenames, extensions, variable names, etc. as they go through more and more shared code. This means that the rest of the scripts should be more equivalent.
Optical flow is relatively similar to image registration tools we already commonly use in neuroimaging fields. Just like a lot of image registrations, optical flow seeks to find the deformation field that explains the displacement of signal between two images. In registration, we might use a similar algorithm to find the transformation matrix needed to map one brain image onto another. In optical flow, we instead use this optimization to estimate how signal moves over time between two temporally adjacent images in a timeseries. BOLD signal for humans, and calcium signal for mice. Just like image co-registrations, this is computationally intensive.
For mice/3D data (x and y over time), this is fairly tractable with some downsampling. Specifically We use NeuroPattToolbox from the Brain Dynamics Group to delineate activity movement trajectories in our mouse data. More info on this software in the table at the end of this markdown.
This gets a little more complicated for human data, because it's acquired in 4D (x y z over time). The workaround for this is to take advantage of freesurfer which maps BOLD, to the cortical surface, and as an intermediate step, to a sphere. Once the data is in spherical space, we can then leverage a nice spherical optical flow codebase that makes this all computationally tractable on the cortical surface (largely through spherical basis functions). This was originally designed to track the movement of progenitor cells on zebrafish gastrulae, but we have sucessfully adapted it for brain data in our prior publication. Citations for freesurfer and spherical optical flow can also be found in the table at the end of this markdown.
Once we have our resulant vector fields, which describe the movement of BOLD/Ca2+, we then compare it to the gradient of the default mode network (NMF) to determine if the direction of activity movement is in the bottom-up or top-down direction.
OpFl MDMA: To run optical flow on data from study 1 (MDMA), you can use this OpFl_mdma.m. This will take the .mgh output from the motion mask script, ensure it's aligned with the fsaverag4 spherical cortical surface, and run optical flow between each temporally adjacent image WITHIN low-motion segments. The latter is ensure by loading in the _ValidSegments_Trunc.txt file generated during the motion masking step. The real meat of this script comes between lines 103 and 114, where optical flow is ran via the of function. This function is directly pulled from the Kirisits et al. paper and repository. We'll populate a matlab struct (us) with the fields vf_left and vf_right to store vector fields describing the motion of BOLD signal between timepoints. Parameters are the same as those provided by default in the code and used in our Neuron paper. These will undergo further processing to get our metrics of interest. There's two scripts to aid in submitting these scripts to the cluster. The first is sbatch_OpFl.sh, which runs the whole pipeline for optical flow data for an individual subject session. It's useful to comment certain chunks out if you don't want to run the whole pipeline. The second helper script is sbatch_all.sh, which submits sbatch_OpFl.sh for each subject for each session as a slurm job. Optional code is in there to NOT submit a job if a certain file (i.e., and output file indicating that the pipeline ran successfully already) exists.
OpFl Psilocybin: This script is identical, but works with file paths/extensions from the psilocybin processing stream. The equivalent helper scripts are sbatch_OpFl_psil.sh for individual subject sessions, and sbatch_psil_all.sh for launching all subjects/sessions.
OpFl Mice: Within the "mice" folder, you'll find Mouse_OpFl.m. This script is also equivalent, but not exactly the same. First, it has to deal with the different file structure the mouse data is in, and it uses the 3D (x y time) optical flow repository provided by Pulin Gong's group (instead of 4D, x y z time). A lot of the code is commented out because we don't utilize the majority of the very cool feature-extraction code they've compiled for this application. There's also some built-in visualization code if you want to double or triple-check that stuff is processing as expected. There's only one helper script for this one because there's only one session (with 5-6 acquisitions) to submit for each mouse for the main analyses. Just sbatch sbatch_mouseOpFl.sh once for each mouse, entering the mouse ID as the only agument. Finally, all mouse optical flow steps are listed for LSD, but equivalent scripts for Diazepam and Dexmedetomidine exist for all of these. They are included in the mouse helper script so they can all be in one place. They run identically, but there are fewer acquisition timepoints for the drugs in some instances. This shouldn't disrupt processing, as scripts are designed to run at the level of individual acquisitions. Individual acquisitions might fail if the data doesn't exist (as one would hope), but the whole processing stream shouldn't be interrupted.
Human and mouse optical flow both use the same underlying algorithm.
Now that we have optical flow estimations, we can pull out the magnitude of activity displacement. We'll just use some good ol' 2,500 year old math for this one. Specifically we apply Pythagorean theorem by calculating the square root of (x displacement^2 + y displacement^2) as our aggregate magnitude of each vector at each point in space over each point in time. We'll just average this over the entire DMN mask over all timepoints to get a single measurement per scan.
Human/spherical optical flow measurements technically start out as 3D (x y and z components for each "between" timepoint), but because we conducted optical flow on the sphere, movement of activity orthogonal to the surface of the sphere is negligible. We take advantage of this redundancy by using cart2sphvec, a matlab-ordained function,to obtain activity movement vectors in a tangent plane (tangential to the spherical surface). We're left with azimuth and elevation, which are equivalent to x and y. This has also been validated previously in our optical flow work in task-fMRI and neurodevelopment.
For the sake of convenience, note that these extractions (2B and 2C) are integrated into sbatch_OpFl.sh and sbatch_OpFl_psil.sh. If you want to run analyses broadly, it might be easier to use these pipeline-level scripts rather than run extractions individually for each instance.
Extract Magnitudes MDMA: We'll run this for individual subjects and sessions once all optical flow runs for that subject have been completed. This can be done serially for each subject/session individually or at the group-level (running all optical flow delineations and then all vector-field-derived extractions). All it does is get the average magnitude of optical flow vectors within the DMN over the ascquisition. The script for this is Extract_DMNMag.m.
Extract Magnitudes Psil: Again, extremely similar to what we ran for study 1. Just differences in the filepath. The script is called Extract_DMNMag_psil.m.
Extract Magnitudes Mice: As usual, the script takes some adaptation to fit with the data structure of the mouse data. It's extremely similar under the hood though. sqrt(x^2+y^2) averaged within the DMN over time for each mouse. That script is called Extract_DMNMag_mice.m.
Aggregating DMN magnitudes for assessment of session-to-session differences
These script just aggregate the derivatives we've calculated. For MDMA, we'll use Extract_DMNMag_dif.m to loop over subjects and sessions. If you are running this off of Sherlock (Stanford's slurm cluster), I'd just request an interactive node with sdev, ml matlab, open up matlab with matlab -nodisplay, and run this single script. Should run in <10 seconds. You can then pull the output with either scp or OnDemand to your local machine.
For psilocybin, we'll use Extract_DMNMag_dif_psil.m. Again, just a tiny bit more complicated because of the variable number of sessions per participant. Same protocol as above for study 1, but note that many files won't be found! This is because the script is designed to index all possible sessions/scans, but subjects all have a different number of sessions and scans per each session. It's a little simpler to script for all possible sessions/scans and remove blank rows after than to prescribe the exact number of sessions per subject and scans for each session a priori.
For LSD, we'll use this Extract_mag_dif_mice.m. For LSD there's one mouse missing a 6th session, but it's not as heterogeneously populated as the psilocybin data.
Extract Relative angles MDMA: This will follow a lot of the same steps as magnitude extractions, and is also called by the same helper scripts. We will again use cart2sphvec to obtain activity movement vectors in a tengent plane, but will instead measure the angle (in degrees) of these movement vectors relative to the gradient of the DMN. The surfaces are loaded in the same way, as are the optical flow vectors. As for magnitude, angles are calculated at each point in the DMN at each timepoint, but are averaged across space and time to yield a single observation per acquisition. Specifically, the percentage of all vectors within the DMN over time that are flowing in the bottom-up direction is saved as our primary measurement. As prior, we define bottom-up as < 90 degrees from the gradient of the DMN (BOLD flowing into the DMN) and top-down as > 90 degrees from the gradient of the DMN (BOLD flowing out of the DMN). That script is called Extract_RelativeAngles.m. Note that we are using "gradient" in the classic calculus definition, not in the sense used in the beautiful work by Paquola, Bernhardt, Valk, Margulies, Misic, Shafiei, Vogel, Smallwood, too many other brilliant scientists to list, etc. For the rest of this guide, I'll refer to "the gradient of the DMN" as "nabla DMN" for clarity.
Extract Relative angles Psil: As prior, this is quite similar to the script for MDMA. The only difference is the input and output filepaths. The script is called Extract_RelativeAngles_psil.m.
Extract Relative angles mice: This script is simpler (because vectors are already in 2 dimensions), but requires different code than the two above because of the difference in input format. There's lots of visual checks in this one in particular (i.e., steps to visualize intermediate calculations to ensure they are working as expected), but they are commented out for ease of use/computational limits (saving out hundreds to thousands of pngs takes up a lot of space and time). Once the data is loaded in, the operations are the same as in the human fmri scripts. The only difference to note is that the DMN cutoff is set to loadings of .6 rather than .3 to account for the mouse DMNs being two merged NMF components. Here's the script: Extract_RelativeAngles_mice.m
Aggregating Bottom-up angle percentages for assessment of session-to-session differences
Again, this uses a parallel structure to what we used for aggregating magnitudes. The script for combining acquisition level measurements of bottom-up percentages for MDMA is here, with a parallel script for psilocybin here. Note it's more complicated for psilocybin because of the non-fixed # of sessions and conditions per participant. For mice, the equivalent script is here
2D.I Humans
Contrast maps are generated for figures 2C, 2E, 2G, 2I, 2M, 2O, 3C, 3E, 3G, 3I, 3M, 3O, S5, and S7.
To calculate the placebo - drug and no drug - drug contrasts, we'll have to run two other scripts first. The first script averages magnitudes for each subject (MDMA here, psilocybin here. After calculating these values, we can find the average across subjects/sessions with the AvgMag_x_L/R.mat output files generated by using this script for MDMA and this script for psilocybin. These last scripts will print out the images found in figure 2 using Vis_FaceVec_Mag
The same procecdure can be applied to generate contrasts maps for % bottom-up. Extract relative angles provides the output needed for Calc_AvgBup (MDMA) and Calc_AvgBup_psil. After calculating these values, we can use scripts equivalent to aggregate_avgMags to produce the same visualizations. The script for doing so is here for MDMA, and here for psilocybin. The embedded visualization script is Vis_FaceVec, which is the same as the magnitude visualization but with different color ranges.
2D.II Mice
The scripts are a little simpler for mice because the data is more homogenous, 2d, etc. All aggregation and plotting of average magnitudes in LSD vs nodrug contrasts is in a script called Aggregate_AvgMag_LSD.m.
Similarly, aggregation and plotting of average % BUP is within a single script that can be found within the script Aggregate_AvgBup_LSD.m.
Note there are some additional drugs in the mouse study. For Diazepam, aggregating average magnitudes and plotting is within Aggregate_AvgMag_Diaz.m, and the equivalent script for bottom-up % is Aggregate_AvgBup_Diaz.m. For Dexmedetomidine, the average magnitudes script is Aggregate_AvgMag_Dex.m, and the average bottom-up % script is Aggregate_AvgBup_Dex.m.
Great, all of section 3 is a bit more straightforward to calculate. We'll start with FC, then temporal autocrrelation.
3A.I Humans Extract_DMNSeg just takes the time series, applies the same masking used in the optical flow pipeline (in terms of the DMN and temporal masking), and calculates the functional connectivity between the DMN and the rest of the cortex. Here is the psilocybin version (different filepaths). The easiest way to use these is to uncomment them in the helper scripts sbatch_OpFl.sh and sbatch_OpFl_psil.sh
The scripts to aggregate the derived FC values is Extract_DMNSeg_dif.m for MDMA, and Extract_DMNSeg_dif_psil.m for psilocybin. The easiest way to run these (on Sherlock if you are at Stanford) is just using an interactive node and running it through command-line matlab. Runs very quickly.
3A.II Mice We use the same approach in mice, but again have to accomodate the different file formats. Use the script Extract_DMNSeg_mice. I would use the helper script sbatch_mouseOpFl.sh, with Extract_DMNSeg_mice uncommented, to run through all of these. That way you'll only have to run it once per mouse and you can run them all in parallel as slurm jobs. Once completed, the script to aggregate extracted DMN segregation values is Extract_DMNSeg_dif_mice.m, which again is most-easily run through an interactive node and command-line matlab (matlab -nodisplay) on the computing cluster you are using.
3B.I Humans We're following the lead of this paper because they put forward a case that a lot of effects can be simplified to changes in temporal autocorrelation. Thankfully that does not seem to be the case for this propagation stuff, but we should be thinking about this paper in network analyses broadly (and papers that have the potential to simplify neuroscience rather than complicate it are extremely rare). The operationalization of autocorrelation essentially comes down to the correlation between signal and signal shifted "+1" in time. The script to calculate temporal autocorrelation in the DMN for MDMA is Extract_TAutoCor.m, and for psilocybin it is Extract_TAutoCor_psil.m. To aggregate the MDMA autocor values, use Extract_TA_dif.m. To aggregate the psilocybin autocor values, use Extract_TA_dif_psil.m script.
3B.II Mice Extract_TAutoCor_mice is pretty simple given the corr(timeseries,timeseries+1) approach, but has a few dozen lines of code just to handle the input filenames. Uses the same approach as in the human data and the same DMN masking as in all other scripts. Will return the average autocorrelation for each run. I would also run this by sbatching the helper script sbatch_mouseOpFl.sh for each mouse.
To aggregate the data across runs in mice (for autocorrelation), use this script. Fairly straightforward given all other code preceding this script.
We are going to run DMN magnitude analyses first, then integration, then autocorrelation, and end with bottom-up because the first three print out metrics used for the AUC curves (ran in the bottom-up scripting). So for a small portion of the bottom-up scripts, magnitude integration and autocorrelation r scripts are a prerequisite. Note that all aggregation scripts should have been run by this point to organize .csvs for r to easily read.
-
4A.I MDMA After running the aggregation scripts, you should be able to download the aggregated files as .csvs to your local machine. I'd reccomend running the r code locally because rstudio is the bomb. The script used to calculate magnitude effects (decreases in DMN magnitude) for study 1 (MDMA) can be found here. There's also a markdown version of the same script available here.I'm going to leave the within-rmarkdown comments as standalone instructions because this readme is getting long. Note this script will also saveout magnitude for unified DMN analyses, i.e., the AUC curves in figure 4. Within this script, Figures for 1A, 2B, and 2D are generated within
-
4A.II Psilocybin This is structured to be parallel to the study 1 (MDMA) analyses, but requires a little extra scripting just to organize all data in an equivalent fashion. The .rmd can be found here. There's also a markdown (.md) version of the same file here As prior, I'm going to leave the within-markdown comments as instructions. Figures for 1A, 2F, and 2H are generated within
-
4A.III LSD In what has become a pattern by this point, the mouse data will be easier to work with. This time it is because of the scan-session structure rather than the 2D vs. 3D data embedding. The script that calculates equivalent statistics for mouse LSD data is here. Markdown available [here]. Figures for 1A, 2J, 2K, 2L, and 2N are generated within (https://github.com/WilliamsPanLab/PsychedProps/blob/master/scripts/Stats_n_viz_mice_Mag.md).
- 4B.I MDMA Post-aggregation DMN integration evaluation in study 1 is within the bottom-up analysis script (4D.I)
- 4B.II Psilocybin Because the psilocybin scripts take a little more organizing code to get everything in stats-friendly format, the integration and autocorrelation scripts are separated out. You can find the script for evaluating DMN integration here. Some variable names and comments are left as they were when copied over from Stats_n_Viz_psil.Rmd, as the script is equivalent but with different input .csvs.
- 4B.III LSD Post-aggregation DMN integration evaluation in study 3 (mice/LSD) is within the bottom-up analysis script (4D.III).
-
4C.I MDMA Post-aggregation DMN autocorrelation evaluation in study 1 is within the bottom-up analysis script (4D.I)
-
4C.II Psilocybin Because the psilocybin scripts take a little more organizing code to get everything in stats-friendly format, the integration and autocorrelation scripts as separated out. You can find the script for evaluating DMN integration here. Some variable names and comments are left as they were when copied over from Stats_n_Viz_psil.Rmd, as the script is equivalent but with different input .csvs.
-
4C.III LSD Post-aggregation DMN autocorrelation evaluation in study 3 (mice/LSD) is within the bottom-up analysis script (4D.III).
-
4D.I MDMA Remember this script encapsulates DMN integration/segregation and autocorrelation scripting. It can be found in markdown format here, or as an .rmd here. Figures for 1A, 3B, 3D, 4A, 4B, 4C, and 4D are generated within. Figures S4 and S9 are also generated in this script
-
4D.II Psilocybin The equivalent script for psilocybin can be found here. Note this also has the lasting effects lil' chunk of code (just after line 1,000). The .rmd version is here. Figures for 1A, 3F, 3H, 4A, 4B and Figure S8 are generated within
-
4D.III LSD The equivalent knit markdown for mice can be found here. The .rmd version can be found here. Figures for 1A, 3J, 3K, 3L, 3N, 4A, and 4B are generated within
You can find the bootstrap and AUC analyses for MDMA (figure 4) further down the same markdown file used above. The .md file doesn't have line numbers, but you can cntrl-F to "library(pROC)" to find where this section stats. If instead you are looking at the .Rmd, this starts at line 1037.
Same goes for bootstrap and AUC/bootstraps for psilocybin. Check out the same markdown as prior. Corresponding .rmd is here, and AUC stuff also starts at line 1037.
Aaaaand same goes for mice. Here's the link. .rmd link. AUC stuff starts at line 238.
You can find the bootstrap and AUC analyses (figure 4) further down the same markdown file used above. The .md file doesn't have line numbers, but you can cntrl-F to "inter-psychedelic-session to find where this section starts. If instead looking at the .Rmd, this starts at line 1350.
| Software | Citation | Use(s) in this project |
|---|---|---|
| Non-negative matrix factorization adapted for brain data | Li, H. et al. (2017): Large-scale sparse functional networks from resting state fMRI. Neuroimage | NMF for human and mouse DMN delineation |
| Freesurfer(7.4.1) | Dale, A. et al. (1999): Cortical surface-baed analysis. I. Segmentation and surface reconstruction. Neuroimage | Cortical surface modeling |
| Connectome Workbench(1.3.1) | Marcus, D. et al., (2011) Informatics and data mining tools and strategies for the Human Connectome Project. Front. Neuroinformatics | Image transformation and visualization |
| g_ls.m | Cui Z., et al. (2013): PANDA: a pipeline toolbox for analyzing brain diffusion images. Front Hum Neurosci | Listing files in a directory when in matlab |
| fmriprep(20.2.3) | Esteban, O., et al. (2018): fMRIPrep: a robust preprocessing pipeline for functional MRI. Nat Methods | Preprocessing of fmri data (prior to xcp-d) |
| xcpd(0.3.0) | Mehta, K., et al. (2024): XCP-D: A Robust Pipeline for the post-processing of fMRI data. Imaging Neurosci | Post-processing of fmri data |
| fsl(5.0.10) | S.M. Smith, et al. (2004): Advances in functional and structural MR image analysis and implementation as FSL. Neuroimage | Basic image manipulation |
| Spherical Optical Flow(1.0) | Kirisits, C., et al. (2013): Decomposition of optical flow on the sphere. International J. on Geomath | Optical flow on spherical surfaces (human data, inflated coritces) |
| "flat" optical flow | Townsend, R. & Gong, P. (2018): Detection and analysis of spatiotemporal patterns in brain activity. PLoS Comp. Biol. | Optical flow on cortical window data (mice, Ca2+) |
| ggplot(3.5.1) | Wickham, H. (2016): ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York | Most figures made in R |
| nlme(3.1.164) | Pinheiro, J. & Bates, D (2000): Mixed-Effects Models in S and S-PLUS. Springer | Statistical testing in R |
| numpy (1.20.3) | Array programming with NumPy (2020) | Basic numeric operations |
| scipy (1.9.1) | SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python (2020) | Basic signal manipulations |
| h5py (3.7.0) | h5py.org | Interface for HDF5 data format |
| Software |
|---|
| Matlab(2022b) |
| R(4.4.1) |
| Python(3.9) |
| GNU Bash(4.2.46) |
There is currently no formal installation command for this software suite. Please use git clone to copy the code from this repository to your local machine: external packages utilized are captured in the tables above, although most packages are not needed for most analyses. Packages needed for core analyses include Freesurfer and Spherical Optical flow.
Demo.m provides a short demonstration of spherical optical flow and includes generation of a test dataset: two “frames” of a gaussian wave traveling forward across a spherical projection of a cortical hemisphere (left). Part 1 of the demo creates the test dataset. Part 2 runs spherical optical flow to estimate signal directionality from the simulated data. Part 3 extracts angular distances between estimated signal directionality and the “gradient” (nabla) of anterior positioning as a reference direction (equivalent to nablaDMN). Part 4 extracts vector magnitudes from estimated signal movement. Part 5 evaluates derivatives of the test dataset. On our system, Demo.m only takes 36.22 seconds to run, but it's probably more useful to run each line individual rather than the script as a whole. Expected output is 1) the average vector magnitudes of optical flow vectors near the peak of the anterior->posterior wave and 2) the average magnitudes of optical flow vectors on the periphery of the "wave". A logical confirmation of the expected Crest>Periphery magnitude values (i.e., "1") is the 3rd output. Output 4) is the average angular distance between estimated signal flows and the gradient of anterior positioning, followed by a logical confirmation that this is <90 degrees. The 5th output is a demonstration of how the average angular distance is ~90 degrees in null circumstances. The final output is the elapsed time of the script running (i.e., tic/toc in Matlab)