Realistic coupling enables flexible macroscopic traveling waves in the mouse cortex

eLife Assessment

This important work presents a novel computational framework for modeling macroscopic traveling waves in the mouse cortex by integrating open-source connectomic and transcriptomic data into a spiking network model. This approach allows the computational model to assign excitatory/inhibitory connections based on neurotransmitter profiles and extends simulations to the 3D domain. The authors present results that demonstrate how spatiotemporal dynamics such as slow oscillations (0.5-4 Hz) emerge and self-organize at the whole-brain scale. This study provides convincing initial insights into the structural basis of traveling waves at the whole-brain scale, and allows future links between connectome-driven spiking neural networks and whole-hemisphere imaging in the mouse.

https://doi.org/10.7554/eLife.108208.3.sa0

Significance of the findings:

Important: Findings that have theoretical or practical implications beyond a single subfield

Landmark
Fundamental
Important
Valuable
Useful

Strength of evidence:

Convincing: Appropriate and validated methodology in line with current state-of-the-art

Exceptional
Compelling
Convincing
Solid
Incomplete
Inadequate

During the peer-review process the editor and reviewers write an eLife Assessment that summarises the significance of the findings reported in the article (on a scale ranging from landmark to useful) and the strength of the evidence (on a scale ranging from exceptional to inadequate). Learn more about eLife Assessments

Abstract
Introduction
Results
Discussion
Materials and methods
Appendix 1
Data availability
References
Article and author information
Metrics

Abstract

Traveling waves are ubiquitous in neuronal systems across different spatial scales. While microscopic and mesoscopic waves are relatively well studied, the emergence of macroscopic traveling waves remains less understood. Here, by modeling the mouse cortex using spatial transcriptomic and connectivity data, we show that realistic cortical connectivity can generate a significantly higher level of macroscopic traveling waves than artificial local and uniform connectivity across multiple oscillation frequency bands, with the strongest advantage appearing in the theta, alpha, and beta frequency bands. By probing the model in different dynamic regimes, we find that macroscopic wave activity depends on both network connectivity and excitatory coupling strength, with a non-monotonic dependence on coupling. Together, our work shows how flexible macroscopic traveling waves can emerge in the mouse cortex and offers a computational framework to further study traveling waves in the mouse brain at the single-cell level.

Introduction

Traveling waves have emerged as a central topic among different large-scale spatiotemporal dynamics observed in neuronal systems (Muller et al., 2018). These waves are observed at different frequencies during various brain activities, ranging from slow-wave activity (Massimini et al., 2004; Liang et al., 2021), sleep spindles (Muller et al., 2016), to faster oscillations in alpha (Zhang et al., 2018; van Kerkoerle et al., 2014), beta (Bhattacharya et al., 2022; Rubino et al., 2006), and gamma (van Kerkoerle et al., 2014; Aggarwal et al., 2022) frequency bands, as well as in broadband (5–40 Hz) activity (Davis et al., 2020). These waves also appear in various brain regions across species, including the hippocampus (Lubenov and Siapas, 2009; Patel et al., 2012; Zhang and Jacobs, 2015), visual cortex (Sato et al., 2012), prefrontal cortex (Bhattacharya et al., 2022), motor cortex (Rubino et al., 2006), or span across multiple brain regions (Liang et al., 2023). More importantly, these waves have long been conjectured to regulate brain functions through mechanisms such as synaptic plasticity (Ermentrout and Kleinfeld, 2001; Muller et al., 2018). And researchers have recently found that traveling waves can modulate important brain functions such as memory, visual processing, and sleep (Mohan et al., 2024; Aggarwal et al., 2022; Miyamoto et al., 2017).

The prevalence and important implications of traveling waves have prompted researchers to investigate their causes and organization. In these studies, computational modeling has been useful due to its flexibility. Researchers have developed both microscopic models (Davis et al., 2020; Naze et al., 2015; Ermentrout and Kleinfeld, 2001) and mesoscopic models (Roberts et al., 2019; Liang et al., 2023) to explore the mechanical properties and computational implications of traveling waves. However, most of these studies focus on microscopic or mesoscopic waves that appear in a single brain region. Large-scale, macroscopic waves across multiple brain regions have been less studied. Moreover, these models are usually constructed with local or uniform connectivity and lack the realistic geometry and connectivity of the brain. It has been shown that in a coupled neuronal network, the coexistence of global wave activity with locally asynchronous states only arises when the number of oscillators is high enough, where local and global activities can coexist (Davis et al., 2021). At the same time, it has been argued that modeling neuronal networks at single-cell resolution will yield more precise discoveries than the mean-field model (Dura-Bernal et al., 2024; Deco et al., 2008). Therefore, it is ideal to study traveling waves in a large domain such as the cortex at single-cell resolution.

However, there are two significant challenges: First, a microscopic wiring diagram of the mouse brain at the whole-brain level is not available. But we do have mesoscopic connectivity data provided by the Allen Brain Atlas (Oh et al., 2014; Knox et al., 2018; Harris et al., 2019), which will be used in this work. Second, simulating and measuring wave activity in a large 3D computational domain with realistic brain geometry is computationally expensive. Simulations of large-scale neuronal networks with realistic neuronal models are usually conducted on supercomputers (Yamaura et al., 2020; Markram et al., 2015).

In this work, we overcome these two challenges in order to study the emergence and propagation of macroscopic traveling waves in the mouse cortex. First, we propose a network-generating algorithm that integrates a spatial transcriptomic dataset that provides the detailed spatial and molecular profile of neurons (Zhang et al., 2023), with voxelized connectivity data (Knox et al., 2018) of the mouse brain. We use this algorithm to create a neuron-to-neuron connectome for one hemisphere of the mouse cortex, which we call Allen connectivity. We then constructed a 3D computational model with both excitatory and inhibitory neurons on top of the Allen connectivity. To overcome the second difficulty, we built our own parallel computational framework that can simulate the constructed model with a desktop GPU.

We simulate this model under stochastic layer-4 input, and we quantify macroscopic wave activity with a three-dimensional phase-gradient analysis. Across a broad parameter sweep, Allen connectivity produces stronger macroscopic waves than local or uniform connectivity, especially in the theta, alpha, and beta bands. We further find that wave activity depends non-monotonically on excitatory coupling strength and is closely linked to network synchrony. Together, these results show how realistic cortical connectivity can support flexible macroscopic traveling waves and provide a computational framework for studying cortex-wide wave dynamics at single-cell resolution.

Results

Construction of the cortical connectivity

A recent spatial transcriptomic dataset (Zhuang-ABCA1) (Zhang et al., 2023) provides detailed positions of neurons in the mouse brain along with their molecular properties, including the neurotransmitter information (Figure 1a, Use of the spatial transcriptomic data). To build a computational model using this data, we need to connect neurons in Zhuang-ABCA1 based on the realistic connectivity of the mouse brain that has been experimentally measured. Previously, Oh et al., 2014 and Harris et al., 2019 have provided a mesoscopic connectivity matrix, known as the ‘projectome’, which is derived from thousands of viral tracing experiments to measure the strength of connections between different brain regions mesoscopically. This projectome was later refined to a finer level using a statistical model introduced by Knox et al., 2018, where the strength of connections is provided on a voxel (100 µm³) resolution. Because both neuronal (Zhuang-ABCA1) and the voxelized connectivity data are registered with the Allen Brain Atlas Common Coordinate Framework (ABA-CCF) (Wang et al., 2020), we now introduce an algorithm that integrates the two datasets to construct a neuron-to-neuron connectivity.

Figure 1 with 1 supplement see all

Download asset Open asset

Construction of the cortical connectivity.

(a) Visualization of all neurons in the hemispherical dataset we use to build the cortical connectivity. The hemisphere contains about 1 million neurons in total, including about 300,000 cortical neurons. (b) Algorithm of building the connectome: (i) First, we identify neuronal excitatory/inhibitory information based on the spatial transcriptomic data. (ii) We then identify the voxelized position for all the neurons. (iii) Based on the voxelized projection data, we calculate the number of connections between every pair of voxels (see Use of the voxelized connectivity data). (iv) Connections are created randomly between the neurons in each pair of voxels. (v) Repeat for all voxels in regions of interest. (vi) Complete the connectivity. (c) The connectivity matrix shows the number of connections between different cortical regions (in log 10 scale). A full version of the connectivity matrix with annotated regions is shown in Figure 1—figure supplement 1. (d) Visualization of outgoing connections from sampled neurons in the primary visual cortex (VISp, left), primary auditory cortex (AUDp, middle), and the ventral anterior cingulate areas (ACAv, right).

Consider constructing a connectivity for N neurons that are contained in M voxels, and let $X_{i}$ represent the physical position of the ith neuron. We have access to the voxelized projection matrix W, where $W_{a b}$ denotes the projection strength from voxel a to voxel b. Our connectivity algorithm requires only a single density parameter, $ρ$ , which determines the average number of connections per neuron, thereby controlling the network’s sparsity (Use of the voxelized connectivity data). Using this information, we propose the following algorithm to construct the connectivity for the target neurons (Figure 1b and Methods):

Algorithm 1: Construction of cortical connectivity
1: Choose the density parameter $ρ$ 2: for each voxel $i$ do 3: Identify the neurons in voxel $i$ 4: for each voxel $j$ do 5: Identify the neurons in voxel $j$ 6: Decide the number of connections from voxel $i$ to voxel $j, N_{i j}$ , based on the projection strength $W_{i j}$ and the density parameter $ρ$ (see Use of the voxelized connectivity data). 7: Randomly draw $N_{i j}$ connections from the neurons in voxel $i$ to voxel $j$ . 8: end for 9: end for

Algorithm 1: Construction of cortical connectivity

1: Choose the density parameter

ρ

2: for each voxel

i

do
3: Identify the neurons in voxel

i

4: for each voxel

j

do
5: Identify the neurons in voxel

j

6: Decide the number of connections from voxel

i

to voxel

j, N_{i j}

, based on the projection strength

W_{i j}

and the density parameter

ρ

(see Use of the voxelized connectivity data).
7: Randomly draw

N_{i j}

connections from the neurons in voxel

i

to voxel

j

.
8: end for
9: end for

Applying this algorithm to all cortical neurons identified in Zhuang-ABCA1, we obtain a cortical connectivity that consists of 297,812 neurons, and we choose each neuron to have 300 connections on average, resulting in a network sparsity of around 0.1% (Construction and storage of the connectivity). The cortical connectivity we constructed is highly heterogeneous, where the number of connections between different regions can vary across four orders of magnitude (Figure 1c). At the same time, both the in-degree and out-degree distribution of the connectome follow a realistic log-normal distribution (Buzsáki and Mizuseki, 2014; Figure 1—figure supplement 1). Since the complete connectome is difficult to visualize, we show some connections from different cortical regions, including the primary visual, auditory cortex, and the ventral anterior cingulate area (Figure 1d).

Macroscopic traveling waves emerge from random stimulation through realistic connectivity

Equipped with the Allen connectivity, we built a single-cell-resolution cortical model to test whether realistic anatomical wiring is sufficient to convert a generic stochastic input into a coordinated, macroscopic traveling wave. Each neuron is a single-compartment Hodgkin–Huxley unit (Pospischil et al., 2008), with its excitatory (glutamatergic, AMPA) or inhibitory (GABAergic) identity inherited from the Zhuang-ABCA-1 transcriptomic labels, and synapses are conductance-based with first-order kinetics (Neuronal and synaptic model). The full hemispheric network of 297,812 neurons is simulated on a single GPU (Simulation).

Spiking networks of this scale can also generate self-sustained activity in similar regimes (Vogels and Abbott, 2005; Kumar et al., 2008), which differs fundamentally from externally imposed noise (Destexhe and Contreras, 2006), and traveling waves have been reported under such conditions (Davis et al., 2021). Therefore, to mimic the stochastic drive that cortex receives from subcortical regions like thalamus, we deliver a 10-Hz Poisson spike train to layer-4 excitatory neurons (Figure 2a), since layer 4 is the canonical thalamocortical input layer (Rockland, 2019); each Poisson event applies a fixed voltage bump $V_{stim}$ , and no other external input is applied.

Figure 2 with 2 supplements see all

Download asset Open asset

Macroscopic traveling waves emerge in response to random layer-4 stimulation through Allen connectivity.

(a) Stimulation protocol: a 10-Hz random Poisson spike train delivers a fixed-amplitude voltage bump $V_{stim}$ to all layer-4 excitatory neurons; no other external input is applied. (b) Two rows of five spatial snapshots of the local-mean (KNN = 300) intracellular voltage, sampled every 5 ms across the window highlighted in (c). A coherent wavefront propagates along the anterior-to-posterior axis. (c) Global mean intracellular voltage of all neurons (dark blue) overlaid on per-region mean voltages (light traces). The red dashed lines mark the snapshot window visualized in (b). (d) Region-sorted raster plot of ~40,000 neurons sampled evenly from major cortical regions, anterior-to-posterior. See Videos for the movie of the simulation.

Under this protocol, we immediately observe macroscopic traveling waves emerge across the cortex (Figure 2 and Figure 2—video 1). The global mean voltage and the region-sorted raster (Figure 2c, d) reveal oscillatory activity that is well synchronized across regions, while the local-mean intracellular voltage maps over a representative 50 ms window (Figure 2b) reveal a coherent wavefront sweeping along the anterior–posterior axis, consistent with previously reported cortex-wide waves (Massimini et al., 2004; Liang et al., 2021; Aggarwal et al., 2024). The corresponding single-neuron-resolution view of the same simulation, with no spatial averaging, is shown in Figure 2—figure supplement 1.

Realistic connectivity produces stronger macroscopic waves than local or uniform connectivity

In past modeling studies of traveling waves, it is common to adopt uniform-random, local, or a combination of the two kinds of connectivity. Therefore, we now simulate three models differing in connectivity in order to compare and investigate whether such macroscopic waves are presented in Figure 2 only emerge from the Allen connectivity. For comparison, we constructed two additional types of connectivity that are common in past studies of traveling waves:

Uniform connectivity: Neurons are connected across the entire cortex randomly (Figure 3b).
Local connectivity: Neurons are connected to a fixed number of their closest neighbors (Figure 3c).

Figure 3 with 1 supplement see all

Download asset Open asset

Top: schematics of Allen (a), local (b), and uniform connectivity (c); bottom: visualization of outgoing connections from the same 50 randomly selected neurons in the primary visual cortex (VISp) for the three connectivity.

Here we visualize the outgoing connections from 50 randomly selected neurons in the primary visual cortex to show the difference between the three connectivity. Under Allen connectivity, the outgoing connections project to different regions heterogeneously (Figure 3a); with local connectivity, the outgoing connections only reach around the adjacent areas where the sampled neurons are located (Figure 3b); and with uniform connectivity, the outgoing connections tend to cover the whole cortex evenly (Figure 3c). Although the topology of three connectivity is different, they share the same out-degree distribution based on the way they are constructed (see Construction and storage of the connectivity).

To properly quantify the traveling waves observed in simulations, we applied a phase-based method from image processing (Fleet and Jepson, 1990), which has been adopted to analyze neuronal dynamics (Rubino et al., 2006), and extended it to 3D data. We first interpolate the simulated voltage onto a regular Cartesian grid that follows the Allen Brain Institute common coordinate framework (CCFv3); the voltage signal is then band-pass filtered into a narrow frequency band, and the Hilbert transform is applied to the narrowband-filtered grid signal to extract the band-specific generalized phase. We then compute the gradient field of the phase signal on the grid. The phase gradient directionality (PGD), an indicator of planar wave activity (Rubino et al., 2006), is then used as our quantitative measure of macroscopic wave activity (see Quantitative measurement of neuronal activity for the full pipeline).

In Figure 4 and Figure 4—video 1, we compare the PGD produced by three representative simulations with different connectivity but the same parameters otherwise. The phase snapshots in Figure 4 span a 20-ms window centered on the PGD peak. Although both Allen and local connectivity produce visible wave activity, Allen connectivity produces waves that are much more macroscopic, while local connectivity produces locally organized planar waves, the form most frequently reported in previous experimental studies (Lubenov and Siapas, 2009; Bhattacharya et al., 2022) and computational models with local connectivity (Keane and Gong, 2015). Uniform connectivity does not generate any apparent structure. This qualitative observation is reflected in the PGD time series (Figure 4b; analogous comparisons in the delta, alpha, beta, and gamma bands are shown in Figure 4—figure supplements 2–5, respectively): Allen connectivity reaches the highest PGD in this representative simulation, with local and uniform remaining markedly lower.

Figure 4 with 7 supplements see all

Download asset Open asset

Allen connectivity produces a higher level of macroscopic wave activity than local and uniform connectivity.

(a) Theta-band (4–8 Hz) phase snapshots from three representative simulations using Allen (top), local (middle), and uniform (bottom) connectivity under matched Poisson stimulation, sampled at −10, −5, 0, +5, and +10 ms relative to the peak of phase gradient directionality (PGD). Colors represent the theta-band generalized phase of each neuron, normalized between 0 and 1. (b) Theta-band PGD of the three simulations versus time, aligned to the PGD peak. (c) Per-band comparison of the maximum PGD of the three representative simulations shown in (a), one bar per band per connectivity. Please see Videos for animations and Quantitative measurement of neuronal activity (illustrated in Figure 4—figure supplement 1) for the calculation of PGD.

Network coupling and synchrony jointly shape macroscopic waves across the parameter plane

Traveling waves are known to propagate differently under different dynamical states, even when the underlying connectivity is fixed (Muller et al., 2018). To assess whether this connectivity-dependent advantage holds beyond a single set of parameters, we systematically swept the coupling strength $g_{AMPA}$ and the Poisson stimulus magnitude $V_{stim}$ on a 10 × 12 grid for each connectivity (Figure 5a, Simulation). Averaged over the full grid, Allen connectivity continues to produce the highest per-band PGD relative to local and uniform connectivity (Figure 5b), showing that the advantage seen in Figure 4c is not specific to one set of parameters. We do note that for the delta band, local and Allen connectivity are able to generate a similar level of strong macroscopic wave activities. The same parameter plane also reveals three qualitatively distinct dynamic regimes for Allen connectivity (Figure 5a): coherent macroscopic waves at weak coupling, a disordered (asynchronous irregular) state at medium coupling in which wave structure collapses, and a globally synchronous regime at strong coupling that again supports waves.

Figure 5 with 2 supplements see all

Download asset Open asset

Coupling strength, connectivity, and network synchrony jointly shape macroscopic-wave activity.

(a) Schematic of the three dynamic regimes that emerge across the $(g_{AMPA}, V_{stim})$ parameter plane for Allen connectivity. At weak-to-intermediate coupling, the network supports *coherent* macroscopic traveling waves; at medium coupling, increased synaptic drive pushes the network into a *disordered* (asynchronous irregular) state where wave structure collapses; at strong coupling, the network reorganizes into a *globally synchronous* regime that again supports wave propagation. Insets show representative phase snapshots from each regime. (b) Per-band max PGD (mean ± SEM across the Poisson $g_{AMPA} \times V_{stim}$ grid), plus the across-bands aggregate (‘All’), for Allen, local, and uniform connectivity, with significance stars; full per-band $(g_{AMPA}, V_{stim})$ heatmaps in Figure 4—figure supplement 6. (c) Max PGD (alpha band) versus stimulus magnitude, error bars across $g_{AMPA}$ . (d) Max PGD (alpha band) versus excitatory conductance $g_{AMPA}$ , error bars across $V_{stim}$ . (e) Kuramoto order parameter versus $g_{AMPA}$ in the alpha band. (f) Synchrony-PGD scatter, one point per $(g_{AMPA}, V_{stim})$ combination, with Pearson $r$ per connectivity (Allen $r = 0.88$ , local $r = 0.82$ , uniform $r = 0.71$ in alpha; all bands positive and significant at $p < 10^{- 3}$ ). Per-band versions of (c)–(e) are shown in Figure 5—figure supplement 1, and per-band synchrony-PGD scatters in Figure 5—figure supplement 2.

To dissect the non-monotonic coupling response that distinguishes the three regimes, we examine how PGD varies along each axis of the $(g_{AMPA}, V_{stim})$ plane in the alpha band – the band that carries a large advantage for Allen connectivity (Figure 5b). Macroscopic wave activity is relatively insensitive to the stimulus magnitude (Figure 5c) but depends strongly and non-monotonically on the coupling strength (Figure 5d): for Allen connectivity, PGD peaks at weak coupling, collapses through a medium-coupling window ( $g_{AMPA} \approx 0.020 - 0.035 n S$ ), and recovers at strong coupling. The uniform connectivity also follows this non-monotonic relationship, but has much lower max PGD at peaks compared to Allen connectivity. Local connectivity follows a different trajectory, peaking at weak coupling, then decaying almost monotonically with $g_{AMPA}$ . These trends hold band by band (Figure 5—figure supplement 1): the peak–dip–recovery profile is visible in every band for Allen and uniform, while local decays almost monotonically across all bands.

To diagnose the mechanism behind this profile, we measured the Kuramoto order parameter $R (t) = | ⟨ e^{i ϕ (x, t)} ⟩_{x} |$ from the generalized phase field $ϕ (x, t)$ of each frequency band and recorded its maximum over each simulation window (Quantitative measurement of neuronal activity). The resulting synchrony curve (Figure 5e) resembles the trend of the maximum PGD well (Figure 5d). For Allen connectivity in the alpha band, mean synchrony peaks at weak coupling, collapses in the medium-coupling window, and recovers at strong coupling. Across the full 120-point $(g_{AMPA}, V_{stim})$ grid, synchrony and PGD are positively correlated in every band and every connectivity (Pearson $r = 0.30 - 0.92$ , all $p < 10^{- 3}$ ; Figure 5f, with per-band scatters in Figure 5—figure supplement 2). Therefore, the PGD trough at medium coupling may be a synchrony trough: increased synaptic drive pushes the network into an asynchronous irregular state, while strong coupling re-establishes rhythmic synchronization that supports wave propagation. This synchrony-to-wave bottleneck seems more significant to the networks with long-range connectivity (Allen, uniform), partly because long-range connections can augment synchrony in the coupled neuronal network (Bazhenov et al., 2008). We also note that although uniform connectivity is able to achieve almost an identical level of synchrony to that of Allen connectivity, the PGD remains much lower due to the loss of spatial organization within.

The macroscopic-wave findings extend to the local field potential estimate

Across the paper, our analyses use the intracellular voltage as the primary measurement for the PGD. To verify that the macroscopic wave structure we identify is also visible in a quantity comparable to experimentally recorded local field potentials (LFPs), we followed a post hoc pipeline (Telenczuk et al., 2020) that reconstructs the net synaptic current received by every neuron directly from the voltage traces and the simulation’s connectivity (LFP estimation from intracellular voltage). The same 3D phase-gradient pipeline is then applied to the reconstructed LFP. The propagating wavefront seen in the local-average voltage is recovered, with cleaner contrast, in the LFP estimate (Figure 6a; Figure 6—video 1). Voltage and LFP max PGD are tightly correlated across the parameter grid (Pearson $r = 0.78 - 0.89$ per band, $p < 10^{- 25}$ for every band; Figure 6b), and the LFP version carries roughly 2–3 × more PGD than the underlying voltage signal (paired $p < 10^{- 13}$ in every band; Figure 6c). Most importantly, Allen connectivity remains the highest-PGD connectivity in every band when measured on the LFP signal (Figure 6d), and the non-monotonic coupling profile from Figure 6 is preserved (Figure 6e). These results suggest that the macroscopic-wave findings are not specific to intracellular voltage and can also be recovered in an LFP-like signal.

Figure 6 with 1 supplement see all

Download asset Open asset

Macroscopic waves are maintained by the local field potential (LFP) estimate.

(a) Paired snapshots of the local-average voltage (top, clipped to [−65,−50] mV) and the synaptic-LFP estimate (bottom, $z$ -scored) across five consecutive 5 ms frames in a representative Allen Poisson simulation. The same propagating wavefront sweeps across the cortex in both modalities. (b) Voltage versus LFP max PGD scatter for Allen connectivity, one point per $(g_{AMPA}, V_{stim}, band)$ ; Pearson $r = 0.78 - 0.89$ per band, all $p < 10^{- 25}$ . (c) Mean max PGD per band, voltage versus LFP, for Allen connectivity (paired $p < 10^{- 13}$ for every band; LFP carries ≈2–3 × higher PGD than the intracellular voltage). (d) LFP-based mean max PGD per band, three connectivities; Allen wins in every band against both local and uniform. (e) LFP max PGD versus $g_{AMPA}$ for Allen connectivity, one line per band; the same three-regime profile (peak, dip, recovery) seen in voltage is preserved in LFP. See LFP estimation from intracellular voltage for the LFP estimation pipeline.

Discussion

Construction of the Allen connectivity

In this work, we present an algorithm that integrates the spatial transcriptomic dataset, which provides the spatial distribution and molecular properties of neurons across the mouse brain, and the mesoscopic connectivity dataset that offers the connectivity information between brain regions on a voxelized level. The algorithm is efficient and can be easily adapted for similar datasets of the mouse brain or for even other species. In essence, our algorithm belongs to the family of stochastic block models (Holland et al., 1983), where the block structure is given by the voxelization of the Allen Brain Atlas and the inter-block connection probabilities are set by the voxelized projection strengths. The projectome (Knox et al., 2018) we use to build the Allen connectivity is derived from a series of viral-tracing experiments that reflect the connection strength between different regions (Oh et al., 2014). However, the relationship between the actual connectome and the mesoscopic projectome remains unclear, and our current model inherits several important simplifications from this dataset.

First, the voxelized projection data is not cell-type specific. In our model, we distinguish only two neuronal populations: glutamatergic (excitatory) and GABAergic (inhibitory), based on the Zhuang-ABCA-1 transcriptomic dataset. However, real cortical connectivity differs dramatically between more refined cell types: pyramidal neurons and interneurons have distinct connection targets, and connectivity is strongly layer-specific. At the same time, we use a single global density parameter $ρ$ to set the average number of connections per neuron across the entire cortex. A region-specific $ρ$ could better capture the known variation in local synaptic density across cortical areas, for example the higher synaptic density of primary sensory regions relative to higher-order association areas. Future versions of the model could allow a region- and cell-type-specific $ρ$ derived from region-level synapse-density atlases together with cell-type-resolved connectivity (Harris et al., 2019), which would simultaneously address the limitations noted above.

Also, our model uses uniform synaptic weights within each synapse type: every AMPA synapse has conductance $g_{AMPA}$ and every GABA synapse has conductance $g_{GABA}$ . In reality, synaptic strength varies with presynaptic and postsynaptic cell type, anatomical distance, and the distribution of synaptic weights is typically heavy-tailed (e.g. log-normal; Buzsáki and Mizuseki, 2014). Extending our algorithm to sample synaptic weights from realistic distributions would be a natural next step and is likely necessary for quantitative comparisons with electrophysiological recordings.

Cortical model and simulation

Compared to previous large-scale models (Capone et al., 2023; Roberts et al., 2019), our model advances in incorporating the detailed spatial distribution of neurons in the 3D computational domain. Without this information, spatiotemporal dynamics cannot be accurately analyzed and visualized. Recent research (Pang et al., 2023; Gulbinaite et al., 2024) has highlighted the significant role that brain geometry plays in neuronal functions. In this work, we did not further use the spatial information of the cortex, such as the geometry of the cortical manifold, which would be an interesting direction for future work.

At the same time, although we are using one of the largest spatial transcriptomic datasets up to date, 300,000 neurons are still only a small portion of the total neuronal population in the mouse cortex. The number of neurons in a network will have important effects on aspects such as synaptic scaling. Among the total population, we only modeled two kinds of neurons, which is also a simplification. Therefore, enriching our model by including a higher number of different types of neurons and studying how they affect the conclusion in this paper is highly meaningful.

For simulations, we choose to randomly stimulate the total population of layer-4 neurons as a way to mimic subcortical input and generate traveling waves, which can be unrealistic. Subcortical structures, such as the thalamus, are vital to cortical dynamics like slow-wave activity (Steriade et al., 1993) and are known to regulate traveling waves (Ye et al., 2023). Therefore, a direct and important future improvement would be adding subcortical structures to the model. In this work, we choose two parameters that affect the dynamic regime of the network, which are the level of the stimuli controlled by the Poisson voltage bump and the coupling strength between neurons. Other parameters, including excitation-inhibition balance, are also likely to shape the dynamical regime (Ahmadian and Miller, 2021) and should be explored in future work.

Synaptic and conduction delays are essential to cortical dynamics such as traveling waves (Petkoski and Jirsa, 2022; Deco et al., 2009; Golomb and Ermentrout, 2000). In our model, because the synapses are modeled by a first-order kinetic model (Neuronal and synaptic model), both excitatory and inhibitory synapses have a rising time that mimics the synaptic delay. However, we did not include conduction delay in our study. Conduction delay is thought to have a proportional relationship with the white matter path length of the connection between two neurons (Lemaréchal et al., 2022). In our model, although we have the 3D positions of neurons, we do not have geometric information about the cortical manifold. Two neurons can be very close in the Cartesian coordinates measured by the Euclidean distance, but very far in terms of the length of the actual connection in the brain. Therefore, incorporating accurate conduction delay in the model is an important future direction.

Quantitative measurements of 3D traveling waves

There are many techniques available for identifying and measuring large-scale neuronal spatiotemporal patterns (Townsend and Gong, 2018; Gutzen et al., 2024). While methods for detecting wave dynamics in three-dimensional data do exist (Budzinski et al., 2023), most published algorithms are designed to analyze 2D data, so our simulation data, which is intrinsically 3D, presents new challenges for measurement. For instance, we need to identify neuron boundaries to exclude unrealistic boundary effects when calculating the phase gradient field (Quantitative measurement of neuronal activity). Additionally, working with high-resolution 3D data is computationally expensive. Despite these challenges, we extended and verified the previous phase-based method to quantify the level of macroscopic wave activity in our simulations. However, as noted in the original work (Rubino et al., 2006), PGD only indicates planar wave activity and may fail to identify other wave types such as spiral waves, which are observed in the cortex (Muller et al., 2016). We also did not analyze other cortical dynamics, such as sinks and sources, which have appeared in our simulation. Extending more sophisticated methods from 2D to 3D to identify diverse spatiotemporal patterns in experimental or simulated data will be a meaningful direction for future research.

Comparing with experimental data

Across the paper, we have chosen the intracellular voltage of neurons as the primary data for the analysis. The voltage signals, though, do not directly resemble experimental measures such as LFP, provide a quick way to analyze the overall dynamics of the whole network. There is strong experimental evidence of cortical spatiotemporal waves from measurements of LFP and multiunit activity reflecting extracellular electrical potentials in regions across the cortex (Destexhe et al., 1999). Recent experimental studies of cortical wave patterns have achieved significant improvements in spatiotemporal resolution via methods such as wide-field calcium imaging (Manley et al., 2024) and optical voltage imaging (Liang et al., 2021). Our simulation framework allows many opportunities to compare with these experimental data.

However, one must be careful when comparing these simulations to experimental data, especially LFPs or other methods measuring extracellular potential rather than intracellular activity. Further work is needed to connect high-resolution mathematical models with experimental data based on the imaging techniques used. A more detailed quantitative comparison with experimental cortical-wave studies, such as the cortex-wide voltage-imaging data of Mohajerani et al., 2013 or the bidirectional visual-evoked waves reported by Aggarwal et al., 2022, is left as a future direction.

Neural oscillation and traveling wave

Different neural oscillations are known to have different generating mechanisms. For example, the gamma oscillation is conjectured to be caused by coordinated interaction of excitatory and inhibitory neurons in the cortex (Buzsáki and Wang, 2012). And slow-wave activity is conjectured to be caused by thalamocortical interaction (Steriade et al., 1993). In our work, we did not fully investigate how different kinds of oscillations emerge in the network from realistic mechanisms, but only studied how these oscillations produce different levels of macroscopic traveling waves in the cortex a posteriori. In the future, it is essential to combine the established understanding of mechanisms behind different oscillations with experiments to further verify the conclusions drawn in this work.

Materials and methods

Use of the spatial transcriptomic data

Request a detailed protocol

The original spatial transcriptomic study (Zhang et al., 2023) provides four datasets of cell information identified from the mouse brain. Among these, we selected the largest dataset, Zhuang-ABCA-1, which includes 2.8 million cells from one hemisphere, with approximately 2.6 million (N = 2,616,328) cells registered to the Allen Brain Atlas. Of these 2.6 million cells, about half are non-neuronal cells and were removed before constructing the connectivity. Since we are only modeling the glutamatergic and GABAergic neurons, we only kept these two kinds of neuronal populations, totaling approximately 1 million neurons (N = 1,050,440). As a result, our final dataset consists of N = 1,050,440 neurons (Figure 1a) in one hemisphere. Within this population, the isocortex contains 297,812 neurons. The cortical connectivity is constructed by applying our algorithm to the total cortical population.

Use of the voxelized connectivity data

Request a detailed protocol

The projection data we use in this study comes from Knox et al., 2018, which provides the connection strength between any pairwise voxel in the isocortex. The most important part of using the projection data is determining the number of connections between neurons in any pairwise voxels based on the connection strength. Let $W$ represent the voxelized connectivity matrix, where $W_{a b}$ represents the connection strength from voxel $a$ to voxel $b$ . In this work, we assume the number of connections between two voxels is proportional to the connection strength, meaning the number of connections from voxel $a$ to voxel $b$ is $N_{a b} = λ W_{a b}$ , where $λ$ is a factor that will be calculated based on the density parameter.

Before building the connectivity, we first choose the density parameter, $ρ$ , which is the average number of connections per neuron. The density parameter $ρ$ is a free parameter that controls the trade-off between computational feasibility and biological fidelity: higher values of $ρ$ yield more connections per neuron and thus higher biological realism, at the cost of greater memory and runtime (Morrison et al., 2007). Let vector $x$ represent the number of neurons in each voxel. If $λ = 1$ , then the total number of connections in the whole connectivity is $N_{t o t a l} = x^{T} W x$ , and the average number of connections per neuron is $\bar{N} = \frac{x^{T} W x}{M}$ , where $M$ is the total number of neurons. Therefore, if we let $λ = \frac{ρ}{\bar{N}} = \frac{ρ * M}{x^{T} W x}$ , we will get the ideal average number of connections per neuron. Consequently, $N_{a b} = \frac{ρ * M}{x^{T} W x} W_{a b}$ . In the code, after calculating $λ$ beforehand, we can multiply the whole projection matrix $W$ by $λ$ to derive a connection matrix, whose each entry now represents the number of connections between any two voxels.

Construction and storage of the connectivity

Request a detailed protocol

Three types of connectivity are used in our simulations: Allen, local, and uniform connectivity. Allen connectivity is constructed by applying our algorithm to the cortical neurons in Zhuang-ABCA1 based on the method described above. Uniform and local connectivity are constructed in a similar fashion: for each neuron, we first read out the number of outgoing connections in the corresponding Allen connectivity, then we randomly choose the same number of neurons from the total population without repetition for the uniform connectivity, or the same number of nearest neighbors for local connectivity. In this way, the out-degree distribution of all three connectivity will be the same, making them more comparable. In particular, local connectivity is constructed through nearest-neighbor coupling, where we use the K-nearest neighbors (KNN) algorithm (Cover and Hart, 1967) to choose a fixed number of outgoing connections to nearest neighbors for each neuron based on spatial position. For all three connectivity, self-connection is not allowed.

When running our connectivity algorithm to build the Allen connectivity, a large number of neurons can pose some challenges to our connectivity algorithm. To tackle this challenge, we loop through the voxels for a significant speed-up. As a result, our algorithm will take $O (M N)$ time, where M is the number of neurons that contain all the neurons and N is the number of voxels. N has a maximum since the computational domain is finite, and in the end, this algorithm is about two orders of magnitude faster than iterating through every neuron.

In our model, we choose $ρ = 300$ , which means there are about 90 million synapses in the model since we have 297,812 neurons. To store such huge connectivity, we use a compressed row form representation of the connectivity matrix, which is also used in similar packages (Dai et al., 2020).

Neuronal and synaptic model

Request a detailed protocol

We use a Hodgkin–Huxley type neuron model that can produce different firing behaviors of a cortical neuron (Pospischil et al., 2008). Each neuron has four ionic channels: voltage-gated sodium, voltage-gated potassium, slow non-inactivating potassium, and leak channels. The pyramidal neurons will have all four ionic channels, while the inhibitory neurons do not have the slow non-inactivating potassium channels because they are fast-firing and do not usually exhibit firing adaptation. Detailed equations of each current and the change of gating variables are provided in the appendix. The parameters of the neuronal model are also drawn from the original paper and presented in Appendix 1—table 1.

For synapses, we model both glutamatergic (AMPA) and GABAergic synapses in our cortical network. Both synapses are modeled by a conductance-based, first-order kinetic model (Ermentrout and Terman, 2010). Please see the supplement for more details.

Simulation

Request a detailed protocol

To simulate the computational model, we have to solve two sets of equations at every time step: the Hodgkin–Huxley equations for the electrophysiology model and the conductance-coupled synaptic equations. We solve the Hodgkin–Huxley equations using an established implicit leapfrog method (DeWoskin et al., 2015), allowing for large timesteps due to increased stability. Both electrophysiology and synaptic equations are solved in parallel on GPU in our framework. The overall parallelism is based on the neuron, which means each set of Hodgkin–Huxley equations of a neuron is put up on an individual thread of GPU to solve at every timestep. Similarly, at any time step, we first check for all neurons firing and then assign each thread to process each of the firing neurons’ downstream connections. Please see the supplements for a detailed description of the numerical methods.

Each simulation lasts for 3 s, and we only record the last 1 s of the simulation to avoid the transient at the beginning of the simulation. The initial conditions of all neurons are randomly chosen at the beginning of each simulation. To be more specific, the initial voltage of all neurons is drawn from a uniform distribution between –65 and –5 mV, and all gate variables are drawn from a uniform distribution between 0.02 and 0.05. For each simulation, we collect the voltage traces of all neurons with a sampling frequency of 2000 Hz.

In all simulations, the total excitatory population of layer-4 neurons is stimulated with a Poisson spike train at a rate of 10 Hz, delivering a fixed-amplitude voltage bump $V_{stim}$ on each Poisson event. To adjust the level of drive we sweep $V_{stim}$ from 5 to 32.5 mV in steps of 2.5 mV. To adjust the coupling strength we sweep the excitatory conductance $g_{AMPA}$ from 0.005 to 0.050 nS in steps of 0.005 nS, with the inhibitory conductance always fixed at 10% of the excitatory conductance. The full $g_{AMPA} \times V_{stim}$ grid (10 × 12 = 120 simulations per connectivity) is used for the Poisson-stimulation analyses in Figures 4—6.

Quantitative measurement of neuronal activity

Request a detailed protocol

To quantify spatiotemporal dynamics, we extend a phase-based method (Rubino et al., 2006) that has primarily been used for 2D data but naturally generalizes to 3D (Figure 4—figure supplement 1 illustrates the full pipeline on a representative snapshot). Because the Hilbert transform is most interpretable when applied to narrowband signals, the simulated voltage traces are first interpolated onto a regular Cartesian grid that follows the Allen Brain Institute common coordinate framework (CCFv3), then band-pass-filtered into a frequency band, and only then converted to a generalized phase via the Hilbert transform. Concretely:

Spatial interpolation. We construct a regular 3D Cartesian grid following the discretization of the Allen Brain Institute common coordinate framework (CCFv3), of size 132 × 80 × 114 at 100 µm spatial resolution, and linearly interpolate the per-neuron voltage time series onto this grid; grid cells that contain no neurons are masked out and excluded from later steps.
Band-pass filtering. The gridded voltage time series are bandpass-filtered into the five canonical bands (delta 0.5–4 Hz, theta 4–8 Hz, alpha 8–12 Hz, beta 12–30 Hz, gamma 30–100 Hz) using a fourth-order Butterworth filter applied with zero-phase filtfilt.
Generalized phase. For each narrowband-filtered gridded signal, we apply the Hilbert transform and extract the band-specific generalized phase $ϕ (x, t)$ for that band.
Phase gradient directionality. We compute the phase gradient $\nabla ϕ$ on the 3D grid using a finite-difference scheme and use it to compute the PGD $P (t)$ (Rubino et al., 2006), our quantitative measure of macroscopic wave activity:

P (t) = \frac{‖ \bar{\nabla ϕ} ‖}{‖ \nabla ϕ ‖}

where the horizontal bar denotes the spatial average across the grid. Higher PGD means the gradient field is more aligned across spatial locations, indicating a more coherent planar wavefront.

In computing the gradient field, the boundary grid cells that contain no neurons are excluded by masking, implemented with the PyVista package Sullivan and Kaszynski, 2019; both interpolation and gradient computation use the intrinsic routines of the 3D Visualization Toolkit (VTK) (Schroeder et al., 1998).

The Kuramoto order parameter $R (t) = | ⟨ e^{i ϕ (x, t)} ⟩_{x} |$ is computed by averaging $e^{i ϕ}$ across the valid (neuron-containing) grid cells of the same per-band phase field, and reduced to its maximum over each 1-s recording window for the synchrony curves in Figure 6.

The same per-band pipeline is run independently for each of the five bands in the per-band analyses (Figures 4c, 5, Figure 4—figure supplements 2 and 5, Figure 4—figure supplement 6, Figure 5—figure supplements 1 and 2). The LFP analyses in Figure 6 use the same pipeline applied to the per-neuron synaptic-LFP estimates (see LFP estimation from intracellular voltage) after interpolation onto the same grid.

LFP estimation from intracellular voltage

Request a detailed protocol

To connect our simulated intracellular voltage traces to experimentally measured LFPs, we used a post hoc pipeline that reconstructs the synaptic current received by every neuron directly from the voltage traces and the connectivity used in the simulation (Telenczuk et al., 2020). Since LFP primarily reflects the net transmembrane synaptic current in a local population (Destexhe et al., 1999), using the reconstructed synaptic current as an LFP proxy is a natural and widely adopted choice in large-scale simulations.

Concretely, for each neuron $i$ at time $t$ we first compute the sigmoidal neurotransmitter output used in the simulation’s synaptic model,

T_{i} (t) = \frac{1}{1 + \exp (- (V_{i} (t) - V_{t c}) / K_{p})},

with $V_{t c} = 2$ mV and $K_{p} = 5$ mV matching the parameters used in the simulation. We then use the CSR connectivity matrix to propagate this presynaptic output to the postsynaptic neurons, keeping excitatory and inhibitory sources separate. Let $A^{E}$ and $A^{I}$ denote the incoming-connection adjacency matrices restricted to excitatory and inhibitory presynaptic neurons, respectively (obtained by transposing the outgoing CSR matrix and masking columns by the presynaptic E/I label); the excitatory and inhibitory synaptic activation at neuron $j$ is

s_{j}^{E} (t) = \sum_{i} A_{j i}^{E} T_{i} (t), s_{j}^{I} (t) = \sum_{i} A_{j i}^{I} T_{i} (t) .

The net synaptic current per neuron is then

I_{j}^{syn} (t) = g_{AMPA} s_{j}^{E} (t) (V_{j} (t) - E_{AMPA}) + g_{GABA} s_{j}^{I} (t) (V_{j} (t) - E_{GABA}),

with reversal potentials $E_{AMPA} = 0$ mV, $E_{GABA} = - 70$ mV, and the same conductances $g_{AMPA}$ , $g_{GABA}$ used in the corresponding simulation. This per-neuron current $I_{j}^{syn} (t)$ plays the role of an LFP-like proxy: spatially averaging $I^{syn}$ over a local population yields a signal analogous to an extracellularly measured LFP. The pipeline is implemented with sparse matrix multiplications and streams the voltage traces in chunks, so it scales to the full 297,812-neuron population across the 1-s recording window.

We then run the same 3D phase-gradient pipeline described in Quantitative measurement of neuronal activity on the reconstructed LFP. The per-neuron synaptic current is interpolated onto the Allen Brain Atlas 3D grid, bandpass-filtered into the delta, theta, alpha, beta, and gamma bands, converted to a phase field via the Hilbert transform, and used to compute PGD.

Visualization

Request a detailed protocol

Throughout the paper, we visualize neurons by coloring them based on their intracellular voltage. The visualization is straightforward for a single neuron such as in Figure 2—figure supplement 1, but in Figures 2, 4, we also need to calculate the average voltage of 300 neurons nearby. This is done by identifying the closest 300 neighbors for each neuron using the KNN algorithm and calculating the average. All visualizations are made using the PyVista package (Sullivan and Kaszynski, 2019).

AI usage

Request a detailed protocol

Generative AI was used during revision to assist with language editing; all scientific content, analyses, interpretations, and final wording were reviewed and approved by the authors.

Appendix 1

Realistic coupling enables flexible macroscopic traveling waves in the mouse cortex

Neuronal and synaptic model

As we described in the main manuscript, each neuron is modeled as a single-compartment Hodgkin–Huxley neuron, whose intracellular voltage is governed by the following equation:

C \cdot A \cdot \frac{d V}{d t} = - A \cdot (I_{L} + I_{N a} + I_{K} + I_{M}) + I_{s y n}

where $C$ and $A$ represent the specific membrane capacitance and membrane area, $V$ is the membrane potential, $I_{L}$ is the leak current, $I_{N a}$ and $I_{K}$ are the voltage-dependent sodium and potassium currents, and $I_{M}$ is a slow non-inactivating potassium current that is responsible for firing frequency adaptation, and $I_{syn}$ is the synaptic current (defined in the Synaptic Currents section). External stimulation is delivered to the total excitatory population of layer-4 neurons as a Poisson spike train at a rate of 10 Hz; on each Poisson event the membrane potential of the stimulated neuron is instantaneously incremented by a fixed-amplitude voltage bump $V_{stim}$ . The ionic currents are based on descriptions in Pospischil et al., 2008. Please see the following for equations that describe each ionic current:

Ionic currents:

I_{L} = g_{L} (V - E_{L})

I_{N a} = g_{N a} m^{3} h (V - E_{N a})

\frac{d m}{d t} = α_{m} (V) (1 - m) - β_{m} (V) m

\frac{d h}{d t} = α_{h} (V) (1 - h) - β_{h} (V) h

α_{m} (V) = \frac{- 0.32 (V - V_{T} - 13)}{\exp (- \frac{V - V_{T} - 13}{4}) - 1} β_{m} (V) = \frac{0.28 (V - V_{T} - 40)}{\exp (\frac{V - V_{T} - 40}{5}) - 1}

α_{h} (V) = 0.128 \cdot \exp (- \frac{V - V_{T} - 17}{18}) β_{h} (V) = \frac{4}{1 + \exp (- \frac{V - V_{T} - 40}{5})}

I_{K} = g_{K} n^{4} (V - E_{K}) \frac{d n}{d t} = α_{n} (V) (1 - n) - β_{n} (V) n

α_{n} (V) = \frac{- 0.032 (V - V_{T} - 15)}{\exp (- \frac{V - V_{T} - 15}{5}) - 1} β_{n} (V) = 0.5 \cdot \exp (- \frac{V - V_{T} - 10}{40})

I_{M} = g_{M} \cdot p \cdot (V - E_{K}) \frac{d p}{d t} = \frac{p_{\infty} (V) - p}{τ_{p} (V)}

p_{\infty} (V) = \frac{1}{1 + \exp (- \frac{V + 35}{10})}

τ_{p} (V) = \frac{τ_{max}}{3.3 \cdot \exp (\frac{V + 35}{20}) + \exp (- \frac{V + 35}{20})}

Synaptic currents:

In our model, the synaptic current $I_{syn}$ can be excitatory AMPA current or inhibitory GABA current. For the synaptic current, we follow the conductance-based model described in Ermentrout and Terman, 2010. Each synaptic current is calculated as

I_{x} = g_{x} s_{x} (V - V_{x})

\frac{d s_{x}}{d t} = a_{r} T_{in} (1 - s_{x}) - a_{d} s_{x}

where $x$ indicates the synaptic type (AMPA or GABA), $T_{in}$ is the total synaptic input, and $s_{x}$ is the gating variable. To calculate the synaptic input for the ith neuron:

T_{in}^{i} = \frac{1}{N} \sum_{j} w_{i j} T_{out}^{j} and T_{out} = \frac{1}{1 + \exp (- \frac{V - V_{T s}}{K_{T}})}

where $W$ is the connectivity matrix, $T_{out}$ is modeled as a sigmoidal function, and $V_{T s}$ and $K_{T}$ are parameters that control the synaptic output for each neuron.

Based on this model, when a pre-synaptic neuron fires, it will generate a synaptic input that excites the post-synaptic neuron’s synaptic channel and further generates a post-synaptic current. All parameters used in the simulation are presented in Appendix 1—table 1.

Numerical method

Both electrophysiology and synaptic equations are solved with an implicit second-order method, which is based on the trapezoidal rule. We refer to the original manuscript (DeWoskin et al., 2015) for the details of the numerical method.

Appendix 1—table 1

Parameters used in the simulations.

Name	Notation	Excitatory	Inhibitory	Unit
Membrane capacitance	$C$	1	1	μF/cm²
Area of neuron	$A$	0.29	0.22	mm²
Leak reversal potential	$E_{L}$	–70	–70	mV
Sodium reversal potential	$E_{N a}$	50	50	mV
Potassium reversal potential	$E_{K}$	–90	–90	mV
AMPA receptor reversal potential	$E_{A M P A}$	0	0	mV
GABA receptor reversal potential	$E_{G A B A}$	–70	–70	mV
Firing threshold	$θ$	–61.5	–61.5	mV
Rising time of AMPA receptor	$τ_{r}^{A M P A}$	0.94	0.94	ms
Decaying time of AMPA receptor	$τ_{d}^{A M P A}$	0.18	0.18	ms
Rising time of GABA receptor	$τ_{r}^{G A B A}$	10.0	10.0	ms
Decaying time of GABA receptor	$τ_{d}^{G A B A}$	0.25	0.25	ms
Conductance of the leak channels	$g_{L}$	0.02	0.08	mS/cm²
Conductance of Na⁺ channels	$g_{N a}$	50	46	mS/cm²
Conductance of K⁺ channels	$g_{K}$	4.8	5.1	mS/cm²
Conductance of slow K⁺ channels	$g_{M}$	0.13	0	mS/cm²
Time constant for $I_{M}$	$τ_{p}$	1123	824.5	ms
AMPA receptor conductance	$g_{A M P A}$	15–50	15–50	nS
GABA receptor conductance	$g_{G A B A}$	1.5–5	1.5–5	nS
Poisson stimulus voltage bump to layer-4 neurons	$V_{stim}$	5–32.5	0	mV
Poisson stimulation rate to layer-4 neurons	$λ$	10	0	Hz
Firing threshold for synaptic output	$θ_{syn}$	–20	–20	mV
Simulation timestep	$d t$	0.02	0.02	ms

Data availability

Raw data related to this study, including the neuronal and connectivity data, along with the code and simulation data related to the figures, have been uploaded to the following repository: https://doi.org/10.6084/m9.figshare.29410397.

The following data sets were generated

1. Sun G
(2026) figshare
Realistic coupling enables flexible macroscopic traveling waves in the mouse cortex.

https://doi.org/10.6084/m9.figshare.29410397.v2

References

1. Aggarwal A
2. Brennan C
3. Luo J
4. Chung H
5. Contreras D
6. Kelz MB
7. Proekt A
(2022) Visual evoked feedforward-feedback traveling waves organize neural activity across the cortical hierarchy in mice
Nature Communications 13:4754.

https://doi.org/10.1038/s41467-022-32378-x
- PubMed
- Google Scholar
1. Aggarwal A
2. Luo J
3. Chung H
4. Contreras D
5. Kelz MB
6. Proekt A
(2024) Neural assemblies coordinated by cortical waves are associated with waking and hallucinatory brain states
Cell Reports 43:114017.

https://doi.org/10.1016/j.celrep.2024.114017
- PubMed
- Google Scholar
1. Ahmadian Y
2. Miller KD
(2021) What is the dynamical regime of cerebral cortex?
Neuron 109:3373–3391.

https://doi.org/10.1016/j.neuron.2021.07.031
- PubMed
- Google Scholar
(2008) Effect of synaptic connectivity on long-range synchronization of fast cortical oscillations
Journal of Neurophysiology 100:1562–1575.

https://doi.org/10.1152/jn.90613.2008
- PubMed
- Google Scholar
(2022) Traveling waves in the prefrontal cortex during working memory
PLOS Computational Biology 18:e1009827.

https://doi.org/10.1371/journal.pcbi.1009827
- PubMed
- Google Scholar
(2023) Analytical prediction of specific spatiotemporal patterns in nonlinear oscillator networks with distance-dependent time delays
Physical Review Research 5:01.

https://doi.org/10.1103/PhysRevResearch.5.013159
- Google Scholar
1. Buzsáki G
2. Wang XJ
(2012) Mechanisms of gamma oscillations
Annual Review of Neuroscience 35:203–225.

https://doi.org/10.1146/annurev-neuro-062111-150444
- PubMed
- Google Scholar
1. Buzsáki G
2. Mizuseki K
(2014) The log-dynamic brain: how skewed distributions affect network operations
Nature Reviews. Neuroscience 15:264–278.

https://doi.org/10.1038/nrn3687
- PubMed
- Google Scholar
1. Capone C
2. De Luca C
3. De Bonis G
4. Gutzen R
5. Bernava I
6. Pastorelli E
7. Simula F
8. Lupo C
9. Tonielli L
10. Resta F
11. Allegra Mascaro AL
12. Pavone F
13. Denker M
14. Paolucci PS
(2023) Simulations approaching data: cortical slow waves in inferred models of the whole hemisphere of mouse
Communications Biology 6:266.

https://doi.org/10.1038/s42003-023-04580-0
- PubMed
- Google Scholar
1. Cover T
2. Hart P
(1967) Nearest neighbor pattern classification
IEEE Transactions on Information Theory 13:21–27.

https://doi.org/10.1109/TIT.1967.1053964
- Google Scholar
1. Dai K
2. Gratiy SL
3. Billeh YN
4. Xu R
5. Cai B
6. Cain N
7. Rimehaug AE
8. Stasik AJ
9. Einevoll GT
10. Mihalas S
11. Koch C
12. Arkhipov A
(2020) Brain Modeling ToolKit: An open source software suite for multiscale modeling of brain circuits
PLOS Computational Biology 16:e1008386.

https://doi.org/10.1371/journal.pcbi.1008386
- PubMed
- Google Scholar
(2020) Spontaneous travelling cortical waves gate perception in behaving primates
Nature 587:432–436.

https://doi.org/10.1038/s41586-020-2802-y
- PubMed
- Google Scholar
(2021) Spontaneous traveling waves naturally emerge from horizontal fiber time delays and travel through locally asynchronous-irregular states
Nature Communications 12:6057.

https://doi.org/10.1038/s41467-021-26175-1
- PubMed
- Google Scholar
(2008) The dynamic brain: from spiking neurons to neural masses and cortical fields
PLOS Computational Biology 4:e1000092.

https://doi.org/10.1371/journal.pcbi.1000092
- PubMed
- Google Scholar
1. Deco G
2. Jirsa V
3. McIntosh AR
4. Sporns O
5. Kötter R
(2009) Key role of coupling, delay, and noise in resting brain fluctuations
PNAS 106:10302–10307.

https://doi.org/10.1073/pnas.0901831106
- Google Scholar
(1999) Spatiotemporal analysis of local field potentials and unit discharges in cat cerebral cortex during natural wake and sleep states
The Journal of Neuroscience 19:4595–4608.

https://doi.org/10.1523/JNEUROSCI.19-11-04595.1999
- PubMed
- Google Scholar
1. Destexhe A
2. Contreras D
(2006) Neuronal computations with stochastic network states
Science 314:85–90.

https://doi.org/10.1126/science.1127241
- PubMed
- Google Scholar
1. DeWoskin D
2. Myung J
3. Belle MDC
4. Piggins HD
5. Takumi T
6. Forger DB
(2015) Distinct roles for GABA across multiple timescales in mammalian circadian timekeeping
PNAS 112:E3911–E3919.

https://doi.org/10.1073/pnas.1420753112
- PubMed
- Google Scholar
1. Dura-Bernal S
2. Herrera B
3. Lupascu C
4. Marsh BM
5. Gandolfi D
6. Marasco A
7. Neymotin S
8. Romani A
9. Solinas S
10. Bazhenov M
11. Hay E
12. Migliore M
13. Reinmann M
14. Arkhipov A
(2024) Large-scale mechanistic models of brain circuits with biophysically and morphologically detailed neurons
The Journal of Neuroscience 44:e1236242024.

https://doi.org/10.1523/JNEUROSCI.1236-24.2024
- PubMed
- Google Scholar
1. Ermentrout GB
2. Kleinfeld D
(2001) Traveling electrical waves in cortex: insights from phase dynamics and speculation on a computational role
Neuron 29:33–44.

https://doi.org/10.1016/s0896-6273(01)00178-7
- PubMed
- Google Scholar
Book
1. Ermentrout B
2. Terman DH
(2010) Mathematical Foundations of Neuroscience
Springer.

https://doi.org/10.1007/978-0-387-87708-2
- Google Scholar
1. Fleet DJ
2. Jepson AD
(1990) Computation of component image velocity from local phase information
International Journal of Computer Vision 5:77–104.

https://doi.org/10.1007/BF00056772
- Google Scholar
1. Golomb D
2. Ermentrout GB
(2000) Effects of delay on the type and velocity of travelling pulses in neuronal networks with spatially decaying connectivity
Network 11:221–246.

https://doi.org/10.1088/0954-898X_11_3_304
- Google Scholar
(2024) Spatiotemporal resonance in mouse primary visual cortex
Current Biology 34:4184–4196.

https://doi.org/10.1016/j.cub.2024.07.091
- PubMed
- Google Scholar
1. Gutzen R
2. De Bonis G
3. De Luca C
4. Pastorelli E
5. Capone C
6. Allegra Mascaro AL
7. Resta F
8. Manasanch A
9. Pavone FS
10. Sanchez-Vives MV
11. Mattia M
12. Grün S
13. Paolucci PS
14. Denker M
(2024) A modular and adaptable analysis pipeline to compare slow cerebral rhythms across heterogeneous datasets
Cell Reports Methods 4:100681.

https://doi.org/10.1016/j.crmeth.2023.100681
- PubMed
- Google Scholar
1. Harris JA
2. Mihalas S
3. Hirokawa KE
4. Whitesell JD
5. Choi H
6. Bernard A
7. Bohn P
8. Caldejon S
9. Casal L
10. Cho A
11. Feiner A
12. Feng D
13. Gaudreault N
14. Gerfen CR
15. Graddis N
16. Groblewski PA
17. Henry AM
18. Ho A
19. Howard R
20. Knox JE
21. Kuan L
22. Kuang X
23. Lecoq J
24. Lesnar P
25. Li Y
26. Luviano J
27. McConoughey S
28. Mortrud MT
29. Naeemi M
30. Ng L
31. Oh SW
32. Ouellette B
33. Shen E
34. Sorensen SA
35. Wakeman W
36. Wang Q
37. Wang Y
38. Williford A
39. Phillips JW
40. Jones AR
41. Koch C
42. Zeng H
(2019) Hierarchical organization of cortical and thalamic connectivity
Nature 575:195–202.

https://doi.org/10.1038/s41586-019-1716-z
- PubMed
- Google Scholar
(1983) Stochastic blockmodels: First steps
Social Networks 5:109–137.

https://doi.org/10.1016/0378-8733(83)90021-7
- Google Scholar
1. Keane A
2. Gong P
(2015) Propagating waves can explain irregular neural dynamics
The Journal of Neuroscience 35:1591–1605.

https://doi.org/10.1523/JNEUROSCI.1669-14.2015
- PubMed
- Google Scholar
1. Knox JE
2. Harris KD
3. Graddis N
4. Whitesell JD
5. Zeng HK
6. Harris JA
7. Shea-Brown E
8. Mihalas S
(2018) High-resolution data-driven model of the mouse connectome
Network Neuroscience 3:217–236.

https://doi.org/10.1162/netn_a_00066
- PubMed
- Google Scholar
1. Kumar A
2. Schrader S
3. Aertsen A
4. Rotter S
(2008) The high-conductance state of cortical networks
Neural Computation 20:1–43.

https://doi.org/10.1162/neco.2008.20.1.1
- PubMed
- Google Scholar
1. Lemaréchal J-D
2. Jedynak M
3. Trebaul L
4. Boyer A
5. Tadel F
6. Bhattacharjee M
7. Deman P
8. Tuyisenge V
9. Ayoubian L
10. Hugues E
11. Chanteloup-Forêt B
12. Saubat C
13. Zouglech R
14. Reyes Mejia GC
15. Tourbier S
16. Hagmann P
17. Adam C
18. Barba C
19. Bartolomei F
20. Blauwblomme T
21. Curot J
22. Dubeau F
23. Francione S
24. Garcés M
25. Hirsch E
26. Landré E
27. Liu S
28. Maillard L
29. Metsähonkala E-L
30. Mindruta I
31. Nica A
32. Pail M
33. Petrescu AM
34. Rheims S
35. Rocamora R
36. Schulze-Bonhage A
37. Szurhaj W
38. Taussig D
39. Valentin A
40. Wang H
41. Kahane P
42. George N
43. David O
44. F-TRACT consortium
(2022) A brain atlas of axonal and synaptic delays based on modelling of cortico-cortical evoked potentials
Brain 145:1653–1667.

https://doi.org/10.1093/brain/awab362
- PubMed
- Google Scholar
1. Liang Y
2. Song C
3. Liu M
4. Gong P
5. Zhou C
6. Knöpfel T
(2021) Cortex-wide dynamics of intrinsic electrical activities: propagating waves and their interactions
The Journal of Neuroscience 41:3665–3678.

https://doi.org/10.1523/JNEUROSCI.0623-20.2021
- PubMed
- Google Scholar
1. Liang Y
2. Liang J
3. Song C
4. Liu M
5. Knöpfel T
6. Gong P
7. Zhou C
(2023) Complexity of cortical wave patterns of the wake mouse cortex
Nature Communications 14:1434.

https://doi.org/10.1038/s41467-023-37088-6
- Google Scholar
1. Lubenov EV
2. Siapas AG
(2009) Hippocampal theta oscillations are travelling waves
Nature 459:534–539.

https://doi.org/10.1038/nature08010
- PubMed
- Google Scholar
1. Manley J
2. Lu S
3. Barber K
4. Demas J
5. Kim H
6. Meyer D
7. Traub FM
8. Vaziri A
(2024) Simultaneous, cortex-wide dynamics of up to 1 million neurons reveal unbounded scaling of dimensionality with neuron number
Neuron 112:1694–1709.

https://doi.org/10.1016/j.neuron.2024.02.011
- PubMed
- Google Scholar
1. Markram H
2. Muller E
3. Ramaswamy S
4. Reimann MW
5. Abdellah M
6. Sanchez CA
7. Ailamaki A
8. Alonso-Nanclares L
9. Antille N
10. Arsever S
11. Kahou GAA
12. Berger TK
13. Bilgili A
14. Buncic N
15. Chalimourda A
16. Chindemi G
17. Courcol JD
18. Delalondre F
19. Delattre V
20. Druckmann S
21. Dumusc R
22. Dynes J
23. Eilemann S
24. Gal E
25. Gevaert ME
26. Ghobril JP
27. Gidon A
28. Graham JW
29. Gupta A
30. Haenel V
31. Hay E
32. Heinis T
33. Hernando JB
34. Hines M
35. Kanari L
36. Keller D
37. Kenyon J
38. Khazen G
39. Kim Y
40. King JG
41. Kisvarday Z
42. Kumbhar P
43. Lasserre S
44. Le Bé JV
45. Magalhães BRC
46. Merchán-Pérez A
47. Meystre J
48. Morrice BR
49. Muller J
50. Muñoz-Céspedes A
51. Muralidhar S
52. Muthurasa K
53. Nachbaur D
54. Newton TH
55. Nolte M
56. Ovcharenko A
57. Palacios J
58. Pastor L
59. Perin R
60. Ranjan R
61. Riachi I
62. Rodríguez JR
63. Riquelme JL
64. Rössert C
65. Sfyrakis K
66. Shi Y
67. Shillcock JC
68. Silberberg G
69. Silva R
70. Tauheed F
71. Telefont M
72. Toledo-Rodriguez M
73. Tränkler T
74. Van Geit W
75. Díaz JV
76. Walker R
77. Wang Y
78. Zaninetta SM
79. DeFelipe J
80. Hill SL
81. Segev I
82. Schürmann F
(2015) Reconstruction and simulation of neocortical microcircuitry
Cell 163:456–492.

https://doi.org/10.1016/j.cell.2015.09.029
- PubMed
- Google Scholar
1. Massimini M
2. Huber R
3. Ferrarelli F
4. Hill S
5. Tononi G
(2004) The sleep slow oscillation as a traveling wave
The Journal of Neuroscience 24:6862–6870.

https://doi.org/10.1523/JNEUROSCI.1318-04.2004
- PubMed
- Google Scholar
(2017) The roles of cortical slow waves in synaptic plasticity and memory consolidation
Frontiers in Neural Circuits 11:92.

https://doi.org/10.3389/fncir.2017.00092
- PubMed
- Google Scholar
1. Mohajerani MH
2. Chan AW
3. Mohsenvand M
4. LeDue J
5. Liu R
6. McVea DA
7. Boyd JD
8. Wang YT
9. Reimers M
10. Murphy TH
(2013) Spontaneous cortical activity alternates between motifs defined by regional axonal projections
Nature Neuroscience 16:1426–1435.

https://doi.org/10.1038/nn.3499
- PubMed
- Google Scholar
(2024) The direction of theta and alpha travelling waves modulates human memory processing
Nature Human Behaviour 8:1124–1135.

https://doi.org/10.1038/s41562-024-01838-3
- PubMed
- Google Scholar
(2007) Spike-timing-dependent plasticity in balanced random networks
Neural Computation 19:1437–1467.

https://doi.org/10.1162/neco.2007.19.6.1437
- PubMed
- Google Scholar
1. Muller L
2. Piantoni G
3. Koller D
4. Cash SS
5. Halgren E
6. Sejnowski TJ
(2016) Rotating waves during human sleep spindles organize global patterns of activity that repeat precisely through the night
eLife 5:e17267.

https://doi.org/10.7554/eLife.17267
- PubMed
- Google Scholar
(2018) Cortical travelling waves: mechanisms and computational principles
Nature Reviews. Neuroscience 19:255–268.

https://doi.org/10.1038/nrn.2018.20
- PubMed
- Google Scholar
1. Naze S
2. Bernard C
3. Jirsa V
(2015) Computational modeling of seizure dynamics using coupled neuronal networks: factors shaping epileptiform activity
PLOS Computational Biology 11:e1004209.

https://doi.org/10.1371/journal.pcbi.1004209
- PubMed
- Google Scholar
1. Oh SW
2. Harris JA
3. Ng L
4. Winslow B
5. Cain N
6. Mihalas S
7. Wang Q
8. Lau C
9. Kuan L
10. Henry AM
11. Mortrud MT
12. Ouellette B
13. Nguyen TN
14. Sorensen SA
15. Slaughterbeck CR
16. Wakeman W
17. Li Y
18. Feng D
19. Ho A
20. Nicholas E
21. Hirokawa KE
22. Bohn P
23. Joines KM
24. Peng H
25. Hawrylycz MJ
26. Phillips JW
27. Hohmann JG
28. Wohnoutka P
29. Gerfen CR
30. Koch C
31. Bernard A
32. Dang C
33. Jones AR
34. Zeng H
(2014) A mesoscale connectome of the mouse brain
Nature 508:207–214.

https://doi.org/10.1038/nature13186
- PubMed
- Google Scholar
(2023) Geometric constraints on human brain function
Nature 618:566–574.

https://doi.org/10.1038/s41586-023-06098-1
- PubMed
- Google Scholar
1. Patel J
2. Fujisawa S
3. Berényi A
4. Royer S
5. Buzsáki G
(2012) Traveling theta waves along the entire septotemporal axis of the hippocampus
Neuron 75:410–417.

https://doi.org/10.1016/j.neuron.2012.07.015
- PubMed
- Google Scholar
1. Petkoski S
2. Jirsa VK
(2022) Normalizing the brain connectome for communication through synchronization
Network Neuroscience 6:722–744.

https://doi.org/10.1162/netn_a_00231
- PubMed
- Google Scholar
(2008) Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons
Biological Cybernetics 99:427–441.

https://doi.org/10.1007/s00422-008-0263-8
- PubMed
- Google Scholar
(2019) Metastable brain waves
Nature Communications 10:1056.

https://doi.org/10.1038/s41467-019-08999-0
- PubMed
- Google Scholar
1. Rockland KS
(2019) What do we know about laminar connectivity?
NeuroImage 197:772–784.

https://doi.org/10.1016/j.neuroimage.2017.07.032
- PubMed
- Google Scholar
(2006) Propagating waves mediate information transfer in the motor cortex
Nature Neuroscience 9:1549–1557.

https://doi.org/10.1038/nn1802
- PubMed
- Google Scholar
(2012) Traveling waves in visual cortex
Neuron 75:218–229.

https://doi.org/10.1016/j.neuron.2012.06.029
- PubMed
- Google Scholar
Book
(1998)
The Visualization Toolkit an Object-Oriented Approach to 3D Graphics

Prentice-Hall, Inc.
- Google Scholar
(1993) Thalamocortical oscillations in the sleeping and aroused brain
Science 262:679–685.

https://doi.org/10.1126/science.8235588
- PubMed
- Google Scholar
1. Sullivan C
2. Kaszynski A
(2019) PyVista: 3D plotting and mesh analysis through a streamlined interface for the Visualization Toolkit (VTK)
Journal of Open Source Software 4:1450.

https://doi.org/10.21105/joss.01450
- Google Scholar
(2020) A kernel-based method to calculate local field potentials from networks of spiking neurons
Journal of Neuroscience Methods 344:108871.

https://doi.org/10.1016/j.jneumeth.2020.108871
- PubMed
- Google Scholar
1. Townsend RG
2. Gong P
(2018) Detection and analysis of spatiotemporal patterns in brain activity
PLOS Computational Biology 14:e1006643.

https://doi.org/10.1371/journal.pcbi.1006643
- PubMed
- Google Scholar
(2014) Alpha and gamma oscillations characterize feedback and feedforward processing in monkey visual cortex
PNAS 111:14332–14341.

https://doi.org/10.1073/pnas.1402773111
- PubMed
- Google Scholar
1. Vogels TP
2. Abbott LF
(2005) Signal propagation and logic gating in networks of integrate-and-fire neurons
The Journal of Neuroscience 25:10786–10795.

https://doi.org/10.1523/JNEUROSCI.3508-05.2005
- PubMed
- Google Scholar
1. Wang Q
2. Ding SL
3. Li Y
4. Royall J
5. Feng D
6. Lesnar P
7. Graddis N
8. Naeemi M
9. Facer B
10. Ho A
11. Dolbeare T
12. Blanchard B
13. Dee N
14. Wakeman W
15. Hirokawa KE
16. Szafer A
17. Sunkin SM
18. Oh SW
19. Bernard A
20. Phillips JW
21. Hawrylycz M
22. Koch C
23. Zeng H
24. Harris JA
25. Ng L
(2020) The allen mouse brain common coordinate framework: a 3D reference atlas
Cell 181:936–953.

https://doi.org/10.1016/j.cell.2020.04.007
- PubMed
- Google Scholar
(2020) Simulation of a human-scale cerebellar network model on the k computer
Frontiers in Neuroinformatics 14:16.

https://doi.org/10.3389/fninf.2020.00016
- PubMed
- Google Scholar
Preprint
1. Ye Z
2. Bull MS
3. Li A
4. Birman D
5. Daigle TL
6. Tasic B
7. Zeng H
8. Steinmetz NA
(2023) Brain-Wide Topographic Coordination of Traveling Spiral Waves
bioRxiv.

https://doi.org/10.1101/2023.12.07.570517
- Google Scholar
1. Zhang H
2. Jacobs J
(2015) Traveling theta waves in the human hippocampus
The Journal of Neuroscience 35:12477–12487.

https://doi.org/10.1523/JNEUROSCI.5102-14.2015
- PubMed
- Google Scholar
1. Zhang H
2. Watrous AJ
3. Patel A
4. Jacobs J
(2018) Theta and alpha oscillations are traveling waves in the human neocortex
Neuron 98:1269–1281.

https://doi.org/10.1016/j.neuron.2018.05.019
- PubMed
- Google Scholar
1. Zhang M
2. Pan X
3. Jung W
4. Halpern AR
5. Eichhorn SW
6. Lei Z
7. Cohen L
8. Smith KA
9. Tasic B
10. Yao Z
11. Zeng H
12. Zhuang X
(2023) Molecularly defined and spatially resolved cell atlas of the whole mouse brain
Nature 624:343–354.

https://doi.org/10.1038/s41586-023-06808-9
- PubMed
- Google Scholar

Article and author information

Author details

Guanhua Sun

Department of Mathematics, University of Michigan, Ann Arbor, United States

Contribution
Conceptualization, Resources, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing – original draft, Writing – review and editing

Contributed equally with
James Hazelden

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-3710-2807
James Hazelden

Department of Mathematics, University of Michigan, Ann Arbor, United States

Present address
Department of Applied Mathematics, University of Washington, Seattle, United States

Contribution
Conceptualization, Data curation, Software, Validation, Investigation, Methodology

Contributed equally with
Guanhua Sun

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-5155-3858
Ruby Kim

Department of Mathematics, University of Michigan, Ann Arbor, United States

Contribution
Visualization, Writing – review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-0209-4875
Daniel B Forger
1. Department of Mathematics, University of Michigan, Ann Arbor, United States
2. Gilbert S. Omenn Department of Computational Medicine and Bioinformatics, University of Michigan, Ann Arbor, United States
Contribution
Conceptualization, Supervision, Funding acquisition, Validation, Investigation, Writing – original draft, Project administration, Writing – review and editing

For correspondence
forger@umich.edu

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-7581-4031

Funding

United States Army Research Office (W911NF-22-1-0223)

Daniel B Forger

Human Frontier Science Program

https://doi.org/10.52044/hfsp.rgp00192018.pc.gr.166630

Daniel B Forger

National Science Foundation (DMS 2052499)

Daniel B Forger

The funders had no role in study design, data collection, and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Ningyuan Wang for a helpful discussion about work and coding. This work is dedicated to Steven Brown, who provided many helpful suggestions.

We acknowledge the following funding: ARO MURI (Multidisciplinary University Research Initiatives) W911NF-22-1-0223, HFSP (Human Frontier Science Program) RGP0019/2018, and NSF DMS (National Science Foundation Division of Mathematical Sciences) 2052499. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the Army Research Office or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Government purposes, notwithstanding any copyright notation herein. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.

Version history

Preprint posted: July 5, 2025
Sent for peer review: July 5, 2025
Reviewed Preprint version 1: September 9, 2025
Reviewed Preprint version 2: August 18, 2026
Version of Record published: September 4, 2026

Cite all versions

You can cite all versions using the DOI https://doi.org/10.7554/eLife.108208. This DOI represents all versions, and will always resolve to the latest one.

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.