Data Processing for Building Digital Brain Twins with The Virtual Brain

Leon Martin

leon.martin@bih-charite.de

Berlin Institute of Health at Charité University Medicine Berlin

2024-05-16

Main Challenges

Goal of neuroimaging: Elucidate brain function in relation to its multiple layers of organization that operate at different spatial and temporal scales

Warning

1. Challenging integration across spatiotemporal scales (Akil, Martone, and Van Essen 2011)
2. Limited reproducibility of imaging studies (Botvinik-Nezer et al. 2020)
3. Reproducible brain imaging studies require huge sample sizes (Marek et al. 2022)

↳ Brain network models help understanding multi-scale brain dynamics by integrating individual structural and functional data with neural population models (Schirner et al. 2018).

 

(Schirner et al. 2018)

• Figure1
  • Construction: Derived from diffusion-weighted MRI tractography and anatomical MRI region parcellations.
  • Node Injection: Nodes receive subject-specific EEG source activity time series instead of noise.
  • fMRI Prediction: Predicted fMRI time series are matched to each subject’s empirical fMRI data, acquired simultaneously with EEG.
  • Network Structure:
    • Local Networks: Comprise excitatory (E) and inhibitory (I) neural population models driven by EEG source activity (red arrows).
    • Global Coupling: Nodes are connected by structural connectomes (green arrows), representing white matter coupling strengths.
  • Analysis:
    • Synaptic Input Currents (Equations 1 and 2)
    • Firing Rates (Equations 3 and 4)
    • Synaptic Activity (Equations 5 and 6)
  • Purpose: To relate neural population activity and network interactions to observable neuroimaging signals.
  • Visualization: See Video 1 for a brain network model construction and simulation results.
• Figure2
  1. Neuron firing decreases during peaks of α-waves and increases during troughs (Haegens et al., 2011).
  2. Simulations show an inverse relationship between neuron firing and α-phase due to neural excitation and inhibition balance (Atallah and Scanziani, 2009).
  3. On a longer timescale (<0.25 Hz), firing rates oscillate inversely to α-band power (Haegens et al., 2011).
  4. Model simulations link α-band power oscillations to fMRI oscillations, predicting resting-state fMRI time series and spatial patterns (de Munck et al., 2008).
  5. Scale-free fMRI power spectra emerge from long-range network input (He, 2011).
  6. Functional connectivity matrices were predicted over varying time windows (Allen et al., 2014).

Current Modeling Studies

Dynamical Brain Model of Decision Making Processes (Schirner, Deco, and Ritter 2023)

Network Structure Shapes Cortical Travelling Waves (Koller, Schirner, and Ritter 2024)

Longitudinal Alzheimer’s Proteinopathologic Disease Progression Model (Cabrera-Álvarez et al. 2024)

The Virtual Epileptic Patient - probabilistic brain modeling in drug-resistant epilepsy (Wang et al. 2022)

Schirner2023 - Tuning curves for a reduced model (two nodes): - Identical to the 379-nodes BNM. - FC (correlation) between the two nodes increased smoothly and monotonically as a function of their E/I-ratio. - The relationship between E/I-ratio and FC persisted when: - Noise strength (upper panel; Eqs. 5 and 6) was modulated. - Structural coupling strength (lower panel; Eqs. 1 and 2) was modulated. - Both parameters are fixed during the fitting of the full 379-nodes model.

• Fitting results for the full 379-nodes model:
  • One exemplary FC.
  • Empirical FC: Upper triangular portion of the matrix.
  • Simulated FC: Lower triangular portion of the matrix.
  • Joint distributions:
    • Without E/I-tuning (upper panel).
    • With E/I-tuning (lower panel).
• Pearson correlations and RMSEs between empirical and simulated FCs:
  • Three model variants:
    • EI-tuning: Tuning algorithm applied on both E and I.
    • E-tuning: Tuning algorithm applied only on E.
    • Original: Tuning of a scalar global coupling scaling factor to rescale.
  • N = 650 empirical and simulated FCs.

Building Blocks of a Brain Network Model

(Martin et al. 2024)

MR-Physics Overview

Larmor Equation: \[ \omega = \gamma B_0 \]

• \(\omega\): Larmor precession frequency
• \(\gamma\): gyromagnetic ratio (\(\gamma_p ≈ 42.577 \text{ MHz/T}\) For protons of a hydrogen nucleus)
• \(B_0\): External magnetic field strength

(Schild 1992)

Analytical Variability in fMRI Data Interpretation

• 70 independent teams analized the same fMRI dataset.
• Every team chose a unique workflow to analyze the data.
• Choice of analysis pipeline had substantial effects on scientific results conclusions

(Botvinik-Nezer et al. 2020; Lindquist 2020)

• Previous studies with the mixed gambles task suggested that activity in the ventromedial prefrontal cortex and ventral striatum (and other) is related to the magnitude of the potential gain.
• 9 Predefined Hypotheses
  • Each hypothesis predicted fMRI activations in a specific brain region, in relation to a specific aspect of the task (gain or loss amount) and a specific group.
  • analysis teams were instructed to report their results on the basis of a whole-brain analysis
• Need for validation and sharing of complex image processing workflows!!

Figure2: Investigates statistical variability across 64 different research teams.

Part A:

• Spearman correlations calculated between team’s whole-brain maps.
• Maps clustered based on similarity (Ward clustering on Euclidean distances). Color coding:
• Row colors indicate cluster membership (Purple for Cluster 1, Blue for Cluster 2, Grey for Cluster 3).
• Column colors represent hypothesis decisions (Green for “yes”, Red for “no”).

Part B:

• Shows average statistical maps for each cluster, thresholded at uncorrected z > 2.0.
• Displays probability of reporting a positive outcome (P_yes) for each cluster.
• Notes: Maps for Hypothesis 1 and 3 are identical, only differing in the brain regions analyzed.

Key Takeaway:

• Highlights significant variability in data interpretation among different research teams.
• Emphasizes the need for standardized methods in neuroimaging analysis.

Network concepts and terminology

(Sporns 2014, 2013; Koller, Schirner, and Ritter 2024)

• Connection Matrix (a):
  • Summarizes binary pairwise relations among 20 nodes.
  • Black indicates presence, white indicates absence of a symmetric connection.
• Reordered Matrix (b):
  • Nodes reordered for optimal modularity partition.
  • Displays three network modules (red, green, blue).
• Spring-Embedded 2D Network Diagram (c):
  • Visual representation of the network from (a) and (b).
  • Nodes colored by module assignments.
  • Highlighted nodes based on various network measures:
    • Node 1:
      • High degree (many connections).
      • High betweenness centrality (on many short paths).
      • High participation coefficient (broadly distributed connections across modules).
      • Low clustering (neighbors are not mutually connected).
    • Node 15: High clustering.
    • Node 8: Low betweenness.
    • Node 14: Low degree.
    • Node 7: Low participation coefficient.

Simulating Whole-Brain Dynamics

• Multimodal data integration
• Models simulating local neural activity are coupled based on a whole-brain network
• Empirical data features can be used for parameter optimization and model fitting

(Deco et al. 2018)

Building Blocks of a Brain Network Model

(Martin et al. 2024)

Empirical data features

MRI

Figure 1: Average parcellated structural and functional connectomes (84 Areas)

• Structural connectivity (SC): Anatomical connections
• Functional connectivity (FC): Statistical dependencies

(Fischl et al. 2002; Desikan et al. 2006)

Connectome: Weighted Graph with nodes as brain areas and edges characterizing the

Empirical data features

EEG

• Empirical: Data is projected to source space to fit model
• Simulated: Forward-modelling via leadfield matrix

\[ \mathbf{V} = \mathbf{L} \cdot \mathbf{J} \]

(Zorzos et al. 2021)

Source Imaging: Programmatic interpretation of EEG results (What we are intuitively doing when looking at EEG data)

Challenges: - Oblique currents - Electric resistance (of different tissues)

V: Voltages on the scalp L: Leadfield matrix J: Currents in the source

• Source Localization: Based on the anatomical model of a subject’s cortical surface.
• Sensor Registration: Registering sensor locations to the source space.
• Lead-Field Matrix Computation: Using a sphere-head model fit to meshes representing:
  • Boundaries between cortical surface and skull.
  • Boundaries between skull and skin.
  • Boundaries between skin and air.
• Customization: Head model and gain matrix customized to each subject’s anatomy.

Inverse Problem:

• In EEG source localization, the inverse problem involves estimating the distribution of neural sources J from the measured scalp potentials V.
• Considered ill-posed because no unique solution.
• Many different source configurations can produce similar scalp potentials.
• small errors in the measurements can lead to large errors in the estimated sources.

MRI Processing Pipeline

(Schirner et al. 2015; Glasser et al. 2013)

• T1w and T2w images are processed through the same steps in parallel up until cross-modal registration.
• A field map is used to remove readout distortion in the T1w and T2w images, which slightly improves their alignment and the resulting myelin maps.
• Bias Field Correction: \(\sqrt{T1w * T2w}\)

Minimal Preprocessing with HCP Pipelines

A surface-based approach

(Glasser et al. 2013, 2016)

Local and Global Connectivity

(Sanz-Leon et al. 2015)

Spatial Scales of Structural Connectivity

• Human Cortex Representation
  • Depicted as a triangular mesh (lower panel).
• Local Connections
  • Represented with a kernel function, \(W_{\xi}(z)\), of the exponential family.
  • Associated with intracortical fibers.
• Anatomical Parcellation
  • Cortex divided into a finite number of patches (upper right).
• Connectome
  • Connections between parcels represented as Uv = 2.
• Time-Delays
  • Introduced via long-range connections {\(u_{jk}\), \(\tau_{jk}\)}.

Scaleable workflows with Snakemake

• Workflow management system for creating reproducible and scalable data analysis pipelines.
• Python-based
• Interoperable with hcp cluster or cloud computing and local execution.
• Flexible and efficient job scheduling
• Supports containerized workflows

1. Example workflow DAG. The greenish area depicts the jobs that are ready for scheduling (because all input files are present) at a given time during the workflow execution. We assume that the red job at the root generates a temporary file, which may be deleted once all blue jobs are finished.
2. Suboptimal scheduling solution: two green jobs are scheduled, such that only one blue job can be scheduled and the temporary file generated by the red job has to remain on disk until all blue jobs are finished in a subsequent scheduling step.
3. Optimal scheduling solution: the three blue jobs are scheduled, such that the temporary file generated by the red job can be deleted afterwards.

Rule-based modular workflows

Note

Rule: Set of instructions that specifies input and output files along with a shell command for transformation.

• Benchmarking for resource management
• Option for caching between workflows
• Reusable standardized wrappers for common tasks

classDiagram
    class Rule {
    }

    class Input {
        +List files
    }

    class Output {
        +List files
    }

    class Params {
        +Dict key_values
    }

    class Threads {
        +Int count
    }

    class Container {
        +String name
    }

    class Shell {
        +String command
    }

    class Log {
        +String file
    }

    Rule --> Input
    Rule --> Output
    Rule --> Params
    Rule --> Threads
    Rule --> Container
    Rule --> Shell
    Rule --> Log

rule create_sc:
    input:
        parc=rules.mmp_vol_to_mrtrix.output.parc,
        tracks=rules.run_tractography.output.tracks,
        sift_weights=rules.run_sift.output.sift_weights,
    output:
        weights=join(MRTrix_out,
                    "sub-{subid}_parc-mmp1_sc_weights.csv"),
        lengths=join(MRTrix_out,
                    "sub-{subid}_parc-mmp1_sc_lengths.csv"),
        assignments=join(MRTrix_out,
                        "assignments.txt"),
    threads: config["nthreads"]
    container:
        config["containers"]["mrtrix"]
    shell:
        """
        tck2connectome {input.tracks} {input.parc} {output.weights} \
            -tck_weights_in {input.sift_weights} \
            -out_assignments {output.assignments} \
            -symmetric -zero_diagonal -nthreads {threads}

        tck2connectome {input.tracks} {input.parc} {output.lengths} \
            -tck_weights_in {input.sift_weights} \
            -symmetric -zero_diagonal -scale_length -stat_edge mean \
            -nthreads {threads}
        echo $(date) > {output.done}
        """

1. Rule Definition:
   • Input: Specifies the input files needed for the task.
   • Output: Defines the expected output files.
   • Threads: Configures the number of threads for parallel processing.
   • Container: Specifies the container to use for running the task.
   • Shell: Contains the actual commands to execute.
2. Key Features Highlighted:
   • Benchmarking for Resource Allocation: Important for optimizing workflow performance.
   • Caching Between Workflows: Useful for reusing intermediate results to save computation time.
   • Standardized Wrappers: Promotes reusability and standardization of common tasks.

This structure makes it easy for your audience to understand how Snakemake rules are constructed and the benefits of using them in workflows.

Writing Workflows

graph TD;
    A[Target Output] --> R2[rule 2];
    R2 --> O1[Output 1]
    O1 --> R1[rule 1]
    R1 --> I1[data/Input 1];
    R1 --> I2[data/Input 2];

Top down applicatio of rules within a Snakefile:

• Application of a rule to generate a set of output fiels is called job
• For each input file of a job, Snakemake recursively determines rules that can be applied to generate it.
• Directed acyclic graph (DAG) of jobs
  • nodes: rules
  • edges: dependencies

Example Workflow

HCP’s Diffusion Pipeline

Brain Imaging Data Structure

Tip

• BIDS stands for Brain Imaging Data Structure.
• It is a standard for organizing and describing neuroimaging and electrophysiological data.
• The goal is to enhance reproducibility, transparency, and data sharing in the neuroscience community.

 

 

Key Features

• Organized Directory Structure: Standardizes how data files and metadata are stored.
• Standard Naming Conventions: Ensures consistency across datasets.
• Comprehensive Metadata: Includes detailed descriptions of the data and acquisition methods.

Resources: BIDS Specification, BIDS Starter Kit

Supported Modalities - MRI (Magnetic Resonance Imaging) - MEG (Magnetoencephalography) - EEG (Electroencephalography) - iEEG (Intracranial Electroencephalography) - PET (Positron Emission Tomography)

BIDS Extension Computational Modeling (BEP034)

• Network Graphs
• Systems of Mathematical Equations
• Computer Code
• File formats: tsv, JSON, and XML

from tvb.simulator.models import Epileptor

mathematical objects

• equations
• parameters
• code

spatial objects

• coordinates
• maps
• sensor locations

temporal objects

• time series
• stimulus patterns

BIDS Folder Structure of a TVB Project

📁 Project/

📁 sub-1/

    📁 net/

        📄 sub-1_nt-distances_desc-SC_net.json

        📈 sub-1_nt-distances_desc-SC_net.tsv

        📄 sub-1_nt-weights_desc-SC_net.json

        📈 sub-1_nt-weights_desc-SC_net.tsv

    📁 eq/

        📄 sub-1_desc-JansenRit_lems.json

        📈 sub-1_desc-JansenRit_lems.xml

    📁 ts/

        📄 sub-1_desc-simVOIy1_ts.json

        📈 sub-1_desc-simVOIy1_ts.tsv

        📄 sub-1_desc-simVOIy2_ts.json

        📈 sub-1_desc-simVOIy2_ts.tsv

Thank you!

References

Akil, Huda, Maryann E. Martone, and David C. Van Essen. 2011. “Challenges and Opportunities in Mining Neuroscience Data.” Science 331 (6018): 708–12. https://doi.org/10.1126/science.1199305.

Botvinik-Nezer, Rotem, Felix Holzmeister, Colin F. Camerer, Anna Dreber, Juergen Huber, Magnus Johannesson, Michael Kirchler, et al. 2020. “Variability in the Analysis of a Single Neuroimaging Dataset by Many Teams.” Nature 582 (7810): 84–88. https://doi.org/10.1038/s41586-020-2314-9.

Cabrera-Álvarez, Jesús, Leon Stefanovski, Leon Martin, Gianluca Susi, Fernando Maestú, and Petra Ritter. 2024. “A Multiscale Closed-Loop Neurotoxicity Model of Alzheimer’s Disease Progression Explains Functional Connectivity Alterations.” Eneuro 11 (4): ENEURO.0345–23.2023. https://doi.org/10.1523/eneuro.0345-23.2023.

Deco, Gustavo, Josephine Cruzat, Joana Cabral, Gitte M. Knudsen, Robin L. Carhart-Harris, Peter C. Whybrow, Nikos K. Logothetis, and Morten L. Kringelbach. 2018. “Whole-Brain Multimodal Neuroimaging Model Using Serotonin Receptor Maps Explains Non-Linear Functional Effects of LSD.” Current Biology 28 (19): 3065–3074.e6. https://doi.org/10.1016/j.cub.2018.07.083.

Desikan, Rahul S., Florent Ségonne, Bruce Fischl, Brian T. Quinn, Bradford C. Dickerson, Deborah Blacker, Randy L. Buckner, et al. 2006. “An Automated Labeling System for Subdividing the Human Cerebral Cortex on MRI Scans into Gyral Based Regions of Interest.” NeuroImage 31 (3): 968–80. https://doi.org/10.1016/j.neuroimage.2006.01.021.

Fischl, Bruce, David H. Salat, Evelina Busa, Marilyn Albert, Megan Dieterich, Christian Haselgrove, Andre van der Kouwe, et al. 2002. “Whole Brain Segmentation.” Neuron 33 (3): 341–55. https://doi.org/10.1016/s0896-6273(02)00569-x.

Glasser, Matthew F., Timothy S. Coalson, Emma C. Robinson, Carl D. Hacker, John Harwell, Essa Yacoub, Kamil Ugurbil, et al. 2016. “A Multi-Modal Parcellation of Human Cerebral Cortex.” Nature 536 (7615): 171–78. https://doi.org/10.1038/nature18933.

Glasser, Matthew F., Stamatios N. Sotiropoulos, J. Anthony Wilson, Timothy S. Coalson, Bruce Fischl, Jesper L. Andersson, Junqian Xu, et al. 2013. “The Minimal Preprocessing Pipelines for the Human Connectome Project.” NeuroImage 80 (October): 105–24. https://doi.org/10.1016/j.neuroimage.2013.04.127.

Koller, Dominik P., Michael Schirner, and Petra Ritter. 2024. “Human Connectome Topology Directs Cortical Traveling Waves and Shapes Frequency Gradients.” Nature Communications 15 (1). https://doi.org/10.1038/s41467-024-47860-x.

Lindquist, Martin. 2020. “Neuroimaging Results Altered by Varying Analysis Pipelines.” Nature 582 (7810): 36–37. https://doi.org/10.1038/d41586-020-01282-z.

Marek, Scott, Brenden Tervo-Clemmens, Finnegan J. Calabro, David F. Montez, Benjamin P. Kay, Alexander S. Hatoum, Meghan Rose Donohue, et al. 2022. “Reproducible Brain-Wide Association Studies Require Thousands of Individuals.” Nature 603 (7902): 654–60. https://doi.org/10.1038/s41586-022-04492-9.

Martin, L., K. Bülau, C. Hüttl, R. A. Schmitt, D. Perdikis, L. Domide, H. Taher, L. Stefanovski, and P. Ritter. 2024. “Semantic Annotation of Computational Neuroscience Knowledge: The Virtual Brain Ontology for Reproducible Whole-Brain Models.”

Sanz-Leon, Paula, Stuart A. Knock, Andreas Spiegler, and Viktor K. Jirsa. 2015. “Mathematical Framework for Large-Scale Brain Network Modeling in the Virtual Brain.” NeuroImage 111 (May): 385–430. https://doi.org/10.1016/j.neuroimage.2015.01.002.

Schild, H. H. 1992. MRI Made Easy (... Well Almost). Berlex Laboratories. https://books.google.de/books?id=mG5rAAAAMAAJ.

Schirner, Michael, Gustavo Deco, and Petra Ritter. 2023. “Learning How Network Structure Shapes Decision-Making for Bio-Inspired Computing.” Nature Communications 14 (1). https://doi.org/10.1038/s41467-023-38626-y.

Schirner, Michael, Anthony Randal McIntosh, Viktor Jirsa, Gustavo Deco, and Petra Ritter. 2018. “Inferring Multi-Scale Neural Mechanisms with Brain Network Modelling.” eLife 7 (January): e28927. https://doi.org/10.7554/elife.28927.

Schirner, Michael, Simon Rothmeier, Viktor K. Jirsa, Anthony Randal McIntosh, and Petra Ritter. 2015. “An Automated Pipeline for Constructing Personalized Virtual Brains from Multimodal Neuroimaging Data.” NeuroImage 117 (August): 343–57. https://doi.org/10.1016/j.neuroimage.2015.03.055.

Sporns, Olaf. 2013. “Making Sense of Brain Network Data.” Nature Methods 10 (6): 491–93. https://doi.org/10.1038/nmeth.2485.

———. 2014. “Contributions and Challenges for Network Models in Cognitive Neuroscience.” Nature Neuroscience 17 (5): 652–60. https://doi.org/10.1038/nn.3690.

Wang, Huifang E, Marmaduke Woodman, Paul Triebkorn, Jean-Didier Lemarechal, Jayant Jha, Borana Dollomaja, Anirudh Nihalani Vattikonda, et al. 2022. “Virtual Epileptic Patient (VEP): Data-Driven Probabilistic Personalized Brain Modeling in Drug-Resistant Epilepsy.” medRxiv, January, 2022–01. https://doi.org/10.1101/2022.01.19.22269404.

Zorzos, Ioannis, Ioannis Kakkos, Errikos M. Ventouras, and George K. Matsopoulos. 2021. “Advances in Electrical Source Imaging: A Review of the Current Approaches, Applications and Challenges.” Signals 2 (3): 378–91. https://doi.org/10.3390/signals2030024.