In this work, we present a co-simulation framework that integrates a microscopically detailed spiking neural network model of the human CA1, simulated using the NEST (based on NEST 3.8.0-post0.dev0; see Supplementary materials) framework (Gewaltig & Diesmann, 2007; Graber et al., n.d.), with macroscopic whole-brain dynamics simulated using The Virtual Brain (TVB; tvb-library 2.3; See supplementary materials) platform (Sanz-Leon et al., 2015). The two main constraints of the framework are scalability and computational efficiency, both of which are essential to support increasing complexity in neural population modeling and whole brain simulations. To tackle the first aspect, the communication protocol transmits rate-coded information for each region in both simulations (see Fig. SM1). To ensure the computational efficiency, the system has been implemented in a containerized environment, enabling straightforward portability and execution across High Performance Computing (HPC) platforms, and has been validated on such systems.

The CA1 region of the human hippocampus was simulated in NEST as in Gandolfi et al. (2023). Briefly, the point-neuron model incorporates data-driven 3D soma positioning (Fig. 1) with the corresponding connectivity matrix generated through a realistic morpho-anatomical connection strategy (Gandolfi et al., 2022). Each neuron is simulated using the Hill-Tononi (HT) model (Hill & Tononi, 2005) which, although originally developed for thalamic neurons, can reproduce hippocampal bursting activity and provides a consistent number of degrees of freedom with a balanced compromise between computational efficiency and detailed dynamics. Moreover, HT neurons support intrinsic fast and slow excitatory and inhibitory input dynamics, enabling synaptic behavior consistent with AMPA, NMDA and GABAergic transmission. Synaptic communication between point neurons is governed by the Tsodyks-Markram model (Tsodyks & Markram, 1997). Such a model captures both the availability of synaptic resources and the probability of neurotransmitter release, enabling short term plasticity. Equations of the HT and Tsodyks-Markram models are detailed in the Supplementary Materials. The network dynamics are simulated at a sub-millisecond time and micrometer spatial resolution, necessitating the use of HPC resources for networks with billions of synapses.

To enable coupling between TVB-based simulation and the NEST-based model, custom input/output devices have been designed in NEST (Kusch et al., 2024) using the Message Passing Interface (MPI) protocol to receive and send data at every simulation time step. The NEST input device receives a rate-coded signal and generates Poissonian spike trains with the received firing rate.

An input device is created for each CA1 mesh vertex, with each device driving a corresponding subpopulation of point neurons, as described below. At each time step, a rate vector is received, containing one value per neuronal subpopulation. The order of the vectors is consistent with the order of creation of the input devices. The NEST MPI-based output device follows a similar design. It records the total number of spike events for each subpopulation and transmits a vector of spike counts back to the TVB process, enabling feedback-driven co-simulation. This will later be used as input for the TVB model. The input/output devices are compatible with multi-node parallel execution, allowing the framework to scale efficiently with increasing network size and complexity, a critical requirement given the growth in computational cost associated with large-scale spiking networks. In the presented study, a down-sampled version of the right CA1 network model was used. The original model, comprising about 5.28 million neurons and 40 billion synapses, was scaled down to 110 000 neurons and 4.3 million connections (Fig. 1B; see Methods). The downscaling procedure has been performed by randomly sampling neurons and, accordingly, preserving only consistent connections. The integrity of the reduced network has been assessed by analyzing the neuronal density distribution (Fig. SM3 A-C), the indegree and outdegree distribution (Fig. SM3 D-E) and the distribution of the connection lengths (Fig. SM3 F) (See Supplementary materials).

The macroscopic brain model in TVB was built using anatomical data derived from T1-weighted and diffusion-weighted magnetic resonance imaging (DW-MRI) (Wang et al., 2024) provided by the human connectome project (Van Essen et al., 2013; Glasser et al., 2013, subject 100206). From these imaging modalities, a cortical surface mesh was extracted (Fig. 1C) and parcellated using the Virtual Epileptic Patient atlas (Wang et al., 2021), generating a map between each mesh vertex and its corresponding brain region. Then, using DW-MRI tractography, a sparse connectivity matrix has been generated between the vertices of the surface mesh, where each weight is proportional to the cross-sectional area of the streamlines connecting the vertices (Triebkorn et al., 2025). Each node of the surface mesh, representing a neural mass of the TVB network, was simulated using the Spatial Epileptor Model (SEM) (Proix et al., 2018) which captures the onset, time course, and termination of ictal-like discharges, as well as their recurrence (Jirsa et al., 2014). Neural masses were coupled via two types of connections: (i) long-range global connectivity derived from tractography (Fig. 2A), and local connections (Fig. 2B), based on a Laplacian kernel spanning over the geodesic distance between nodes along the mesh surface. The resulting global connectivity exhibits a block-structured organization, revealing distinct inter-hemispheric, intra-hemispheric, and subcortical connections, as well as projections from subcortical to cortical regions. Due to the high dimensionality of the system, the visualizations of Fig. 2 A, B display only the sparsity pattern of the connectivity matrix, omitting connection strength to avoid visual clutter. Moreover, the projection of this inherently three-dimensional connectivity structure onto two-dimensional plane introduces artifacts, which stem from mesh node indexing and the high dimensionality of the dataset.

The 6D Spatial Epileptor Model is described by a system of 6 differential equations (see supplementary material) representing distinct dynamics states: two fast discharge variables (\(x_1\), \(y_1\)), two spike and wave events (SWE) states (\(x_2\), \(y_2\)), one permittivity state (z), and an integral coupling function (g). Specifically, \(x_1\) models fast excitatory activity during seizures and \(y_1\) acts as the associated recovery variable. The variables \(x_2\) and \(y_2\) represent slower oscillatory SWE activity and its corresponding recovery dynamics, respectively. The variable z modulates seizure susceptibility and governs transitions between ictal and non-ictal states. The integral coupling function g links the two subpopulations, driving the SWE states by integrating past fast discharge activity via temporal convolution. The external input currents to the rapid discharge population \(I_\) and the SWE population \(I_\), are constant bias terms that dictate the stable points of the corresponding subsystems. The parameter m controls the dynamics of SEM during the ictal period, while the excitability parameter \(x_0\) controls the equilibrium point of the SEM, thus the ability, for each node, to switch between a non-ictal and an ictal state (Fig. 2C). Moreover, the SEM includes global coupling applied to the variable \(x_1\), local coupling influences both the dynamics of \(x_1\) and \(x_2\) and the integral coupling function g. The global coupling is based on a weighted global connectivity signal propagation, where a history buffer containing all neural masses’ previous states is sampled according to a delay matrix based on tract length between regions. Local coupling models the influence of the activity of a node on its neighbors by: (i) rapid discharge coupling (\(loc_\)), (ii) an integral coupling term (\(loc_\)), and (iii) the effects of coupling through SWE events (\(loc_\)), where \(loc_\), \(loc_\), \(loc_\) are the corresponding local coupling strengths. In particular, the g state directly modulates the SWE activity, thus creating a bidirectional link between the fast discharge and the SWE states in neighboring regions. The whole brain mesh consists of 337,981 cortical nodes and 97,509 subcortical nodes, totaling 435,490 total locations. Since both global connectivity, tract length mapping and geodesic distances would result in extremely large and sparse matrices, custom TVB classes have been used. The custom SimulatorSparse class includes a specialized SparseHistory class attribute that keeps track of states for delayed communications. Local coupling is handled by the HeavisideSparse function, optimized via a Numba-accelerated sparse matrix routine. Global connectivity is represented by the ConnectivitySparse class, which stores weights, tract lengths, and delays using sparse matrices.

For validation, a reduced mesh was used, including only regions from the right hemisphere: fusiform gyrus, para-hippocampal gyrus, entorhinal cortex, subiculum, CA1, CA2, CA3, CA4 and dentate-gyrus, as extracted from MRI data. The excitability parameter (\(x_0\)) was set to model CA1 as the epileptogenic zone (\(x_0\)= -1.6), while the surrounding areas were defined as propagation zones (\(x_0\) = -1.9). The system was initialized at a stable point (\(x_1\) = -1.5, \(y_1\) = -11, z = 3, \(x_2\) = -0.9, \(y_2\) = 0.3, g = -0.1), except for a localized subset of adjacent CA1 nodes, designated as the seizure onset zone, where the epileptic oscillations triggered pathological propagation (\(x_1\) = 0, \(y_1\) = -5, z = 3, \(x_2\) = 0, \(y_2\) = 0, g = 0). The local coupling parameters were calibrated to ensure coherent propagation of both dynamics \(x_1\) and \(x_2\) throughout the network with the following values (k = 0.636, \(\gamma _\) = 0.34, \(\gamma _\) = 0.064, \(\gamma _\) = 0.032, \(\gamma _\) = 1, \(\theta _\) = -1, \(\theta _\) = -1, \(\theta _\) = -0.5, \(I_\) = 3.1, \(I_\) = 0.45, \(\tau _0\) = 2857, \(\tau _1\) = 1, \(\tau _2\) = 10, tt = 0.17)

High Resolution version of TVB, using subject’s mesh points as nodes. (A) Global Connectivity Matrix obtained from Tractography, (B) Local Connectivity based on Geodesic distance between the nodes. (C) SEM activity (\(x_2-x_1\)) timeseries of a single TVB node and corresponding \(x_1\), \(y_1\) subpopulation dynamic driven by permittivity state z reproducing seizure like oscillations. (D) Assignment between TVB mesh nodes and NEST model neurons using the HippUnfold protocol (DeKraker et al., 2022, 2023) to bring the two sets of points in a common coordinate system. (E) Histogram showing the number of vertices in the mesh assigned to a certain number of neurons

The two simulators were coupled through MPI based communication protocol that enabled bidirectional exchange of activity across the two spatio-temporal scales. To have topologically coherent co-simulation across models, neurons were assigned to vertices of the TVB mesh by mapping each vertex to the spatially nearest NEST neurons. This mapping was performed using geodesic Laplace coordinates along the Anterior-Posterior (AP) and Proximal-Distal (PD) axes of the hippocampal surfaces, as computed by the HippUnfold toolbox (DeKraker et al., 2022, 2023). The HippUnfold toolbox provides a surface-based topology-preserving representation of the hippocampus based on thickness, curvature, gyrification, texture and other morphological and laminar features (DeKraker et al., 2022). This surface-based representation is better suited for anatomically meaningful alignment between subjects than conventional volumetric registrations (DeKraker et al., 2023).

After projecting both the 3D coordinates of the point neurons and the mesh vertices into a common 2D intrinsic coordinate system, a KDTree clustering algorithm was used to identify the closest neurons to each mesh vertex (Fig. 2D). The assignment was limited to mesh nodes labeled as CA1, CA2 or Subiculum. To refine anatomical labeling, mesh nodes originally labeled CA2 or subiculum but assigned to any point neuron were relabeled CA1. In contrast, mesh nodes initially labeled CA1 that remained unassigned were relabeled CA2 or subiculum, depending on the identity of the nearest non-CA1 mesh node.

The framework employs a Singularity container (Kurtzer et al., 2017) for portability across HPC systems and is organized into three main Python-based processes. The launcher process manages the sequential initialization of the two simulations. It first starts the TVB simulation, which generates the model, and then opens an MPI communication port via a custom model class. This port expects as many inputs/outputs as there are surface vertices labeled CA1, the region targeted for co-simulation. When the MPI port is active, the launcher starts the NEST simulation, which connects to the port using custom developed streaming I/O devices: spike_detector_mpi and spike_generator_mpi. For each CA1 mesh node, a pair of these devices is instantiated and connected to the corresponding neurons, based on the AP-PD mapping. Spike_generator_mpi receives input rates from the MPI port and generates the corresponding spike trains using an inhomogeneous Poisson process. These spikes are relayed to NEST point neuron models through parrot neurons and static synapses. In contrast, spike_detector_mpi is designed to transmit the number of spikes detected within each simulation time step back to the MPI port. However, in this preliminary version of the framework, this feedback is disabled. TVB influences CA1 spiking activity, but not vice versa.

To convert the absolute state variables of the SEM into meaningful spiking input, a Rectified Linear Unit (ReLU)-like transformation is applied to a combination of \(x_2\) slow wave events and \(x_1\) fast discharges (Eq. (1)-(3)). This produces a firing rate which, scaled by the simulation timestep, defines the Poisson parameter used to probabilistically generate spikes at each time step for each subpopulation.

$$\begin p_= (x_1 + o_)\cdot g_\cdot dt \end$$

(1)

$$\begin p_= (x_2 + o_)\cdot g_\cdot dt \end$$

(2)

$$\begin spike\_generation = Poisson(p_ + p_) \end$$

(3)

Where \(o_\), \(o_\) are the conversion offsets and \(g_\), \(g_\) are the conversion gains. The conversion parameters were configured such that at least one NEST spike is generated for each SWE Complex, modeled through the \(x_2\) state, within the connected subpopulation. Additionally, the probability of generating in response to the rapid discharge component (\(x_1\)) denoted as \(p_\), was set to be half the probability associated with the SWE complex (\(p_\)). This parameterization results in Poisson-distributed firing rates ranging from 0 to 0.31 spikes/s for \(x_2\)-driven activity and from 0 to 0.16 spikes/s for \(x_1\)-driven activity. These values were derived using the following \(o_\)= 0.5, \(o_\)= 2, \(g_\)= 1.25, \(g_\)= 0.375 and a simulation step dt = 0.1 ms for both NEST and TVB simulations. This conversion scheme is grounded in the observation of spiking synchronization driven by slow SWE dynamics, as described in Truccolo et al. (2014), where neuronal firing showed temporal clustering at the onset of the spike phase during recorded seizures. A secondary weaver spiking activity was also observed during the wave phase, which in the present model is captured via the \(x_1\) component, modulating the contribution of rapid discharges in the SEM-to-spike-rate transformation. At each simulation time step, the received values \(p_\) and \(p_\) drive a Poisson process that determines whether each subpopulation generates a spike (i.e., \(spike\_generation\)). To validate the co-simulation framework, an initial test was performed with SNN connectivity disabled, in order to confirm the synchronization between SEM-derived \(x_1\) and \(x_2\) activity and subpopulation firing rates. Subsequent simulations with NEST connectivity enabled the investigation of the dynamics of network-level signal propagation. Extensive tests were conducted to assess the correct interaction between TVB and NEST simulators (see Supplementary materials for a full description).

(A-C) Temporal evolution of the TVB SEM, showing the propagation of seizure-like events at three different time points. The onset site of the seizure comprises a subset of CA1 vertices and is pointed by the red arrow. (D-F) Temporal evolution of the NEST spiking neural network showing coherent activity (red dots) between TVB vertices and CA1 point neurons. Circled regions in C and F highlight the same regions in the two models. (G) Propagation of TVB node activity (\(x_2-x_1\)) as a time series across 16,636 vertices. (H) Timeseries of a single TVB node. (I) Comparison between TVB node \(x_2\) driving signal (upper panel, blue), the activity converted spike probability sent to the NEST model (upper panel, yellow) and generated spike activity in the NEST model (lower panel). (J) The phase of the spike probability sent from TVB to NEST computed using the Hilbert transform and (K) average phase locked spike histogram across CA1 vertices and corresponding spiking subpopulations