Arbor is an open-source library for building simulations of biophysically detailed neuron models (Abi Akar et al., 2019). It provides an alternative to software like NEURON (Carnevale and Hines, 2006), but with a strong emphasis on modern hardware and scalability to large-scale systems (Hines, 1984). Its overall set of capabilities allows Arbor to model neural networks at a level of resolution beyond point models to explore phenomena like dendritic computation. Thus, support for bulk-synchronous parallelism via MPI, shared memory parallelism by utilizing a thread-pool and job system is central to Arbor, and certain cell types—primarily cable cells—can further leverage SIMD and GPU hardware. Arbor is written in C++, though most users interface with it through an intuitive, high-level Python interface built on top of the lower level implementation. The underlying numerical model of Arbor is the cable equation:

where the membrane potential U is computed over the morphological structure of the neural tree; the spatial coordinate x and the derivative are to be understood within this structure (Von Helmholtz, 1850; Thompson and Kelvin, 1855; Hodgkin et al., 1952; Hodgkin and Huxley, 1952a,b; Loligo, 1952; Hodgkin and Huxley, 1952c; Scott, 1975). The parameters c and σ define the membrane capacitance and longitudinal resistance. The trans-membrane current density i models the entirety of ionic and non-ionic currents. In both NEURON and Arbor, these are calculated from user-specified sets of differential equations, potentially varying along the morphology. The equations for i and U are solved in alternation (Lie-Trotter splitting) using a first-order implicit method.

We begin by selecting a dynamical model that allows for relatively easy yet realistic simulation of a broad spectrum of neural activity at the single-neuron level, governed by a small set of biophysical parameters. The neural model from Depannemaecker et al. (2022) was formulated for Arbor in the Neuron MODeling Language (NMODL). The following equations form the slow part of the system, describing the evolution of ion concentrations due to voltage-gated channels, active pumps, and buffering by an external bath, see Figure 1 for a schematic of the dynamics. It describes the ionic exchanges between the intracellular and extracellular spaces (ICS, ECS) of a neuron immersed within an external bath, acting as a potassium buffer of concentration K_bath. Ions flow between the ICS and ECS through a sodium-potassium pump and the sodium, potassium and chloride voltage-gated channels, driving changes in the internal (K_i, Na_i, Cl_i) and external (K_o, Na_o, Cl_o) ionic concentrations. By gradually increasing the external bath concentration of potassium ions K_bath, the model sequentially presents these patterns: resting state (RS), spike train (ST), tonic spiking (TS), bursting, seizure-like events (SLE), sustained ictal activity (SIA) and depolarization block (DB), see Figure 2. The fast dynamics of the membrane potential V is modeled in Arbor via the cable equations, see above, which require computing the ion current densities i_X used in the simulator update as:

with the ion species X = {K, Na, Cl} and a non-ion current density:

These currents enter the cable equation Equation 1 as the trans-membrane current i via:

Following the Hodgkin-Huxley model, conductivities g_X are written as:

The—internal i and external o—ion concentrations are modeled as:

The variables {ΔK_i, K_g} evolve as:

where γ converts the current density i_X to molar flux, summarizing the effect of the ion pump in Figure 1A and the external buffer. Finally, fast dynamics were reduced and adjusted to mammalian neurons:

based on the observations that the reaction of the sodium gating variable to changes in V is nigh instantaneous and h(t) + n(t) = const.

The resulting ion channel was added to a basic, spherical, single-compartment neuron. After implementing this biophysical model, we reproduced the firing patterns using the parameters of the reported model (Figures 2A–F), see also Depannemaecker et al. (2022) for more details and motivation for model parameter choices. Note that despite the values given in the original publication, neither the Arbor nor the published reference implementation produces the depolarization block pattern at K_bath = 20mM but only at around K_bath = 22.5mM. From here, a simple model network was developed, comprising N total cells, with a mixture of tonic f · N and SLE (1 − f) · N cells, where both sub-populations are assigned individual values for K_bath, sketched in Figure 1. Cells are connected using exponential synapses with an internal weight of w = 0.5 chosen to produce an activity similar to Rabuffo et al. (2025) which uses delta synapses.

In addition to the elementary spherical morphology, we also investigated two multi-compartmental neuron models. The first model included a single dendritic segment of 25μm, subdivided into 5μm compartments. The second model extended the dendrite into a randomly generated tree composed of 5μm compartments. To introduce variability across cells, the random number generator was seeded with each cell's unique identifier (see Figure 3 for examples). In future work, these synthetic morphologies will be replaced with reconstructions derived from neural imaging data available in public databases. Supporting such models in Arbor will require only a simple command to load per-cell morphology data from disk.

Following Deco et al. (2013), we used the reduced Wong-Wang model (Wang, 2002) to simulate resting-state activity and to investigate the dynamics of local brain regions embedded within a large-scale brain network. The mean firing rate H(x_I) and mean synaptic gating variable S_I of region I are described by:

where x_I is the synaptic input to the I-th region, ω = 1 denotes the local excitatory recurrence, c_IK is the strength of the structural connection from the local area I to K, and G is a global coupling strength. The parameters are set to the same values as those used in the TVB implementation. J_N = 0.2609nA is the synaptic coupling of NMDA receptors and I₀ = 0.33nA is the baseline external input. The kinetic parameters are τ_s = 100ms and γ = 0.641. The parameters of the input-output function H are a = 0.27nC⁻¹, b = 0.108kHz, and d = 154ms. Depending on the tuning of G, the system exhibits a multi-stable regime, with steady states of high and low spiking activity. Here, we set G = 0.096.

Both Arbor and TVB offer support for attaching a second simulator to perform co-simulation, potentially at different scales. Co-simulation from TVB's viewpoint is the simpler technology of the two frameworks, since TVB is designed to execute as a single process. TVB allows for exchange of any variable relevant to the region models and any number of variables. The co-simulation partner is encapsulated in one or more TVB regions, called proxy nodes, see Figure 4. These proxies present a conforming interface to TVB, exchanging the salient variables as a table, one row per time-step, one column per variable. As TVB advances in lockstep on a global time-step, this is almost identical to normal operation. However, co-simulation introduces the concept of an ‘epoch’ to TVB, i.e., the length of time that conforms to the smallest delay τ_min in the set of inter-region connections delays τ_IJ, with I and J referring to two connected regions. These delays are part of the connectome data used to construct a TVB simulation. In the case that a connectome contains zero-valued delays these must be replaced with a pre-defined finite value. Further, it is required that the time-step evenly divides τ_min. Co-simulation thus can integrate all nodes' state, including the proxy, for one epoch τ_min without exchanging data. This is correct as an event emanating from any region I at time t influences any other region J at time t + τ_IJ ≥ t + τ_min. Only after an epoch, data need to be exchanged between the proxy and the rest of the TVB regions. A TVB–NEST demonstration has been published to showcase the interaction between a local network of spiking neurons and the whole-brain network dynamics (Kusch et al., 2024).

Arbor has a different design in terms of connectivity. Interaction between physically separate cells is mediated by action potentials, e.g., when the membrane potential triggered by dedicated sources crosses a configurable threshold. Cells are connected by wiring these sources to corresponding sinks like synapses via an abstract connection object comprising a delay and weight, modeling transmission and attenuation via an axon. In contrast to TVB, Arbor is fundamentally a distributed system and internally employs the same approach to decoupling via the minimum network delay as explained above. To initiate co-simulation, Arbor utilizes an additional interface to manage the spike exchange from external connections that are originated from outside but terminate at cells simulated in Arbor. On the technical side, the latter part leverages MPI_Allgatherv through an inter-communicator and effects that the concatenation of all spikes sent from all MPI ranks running TVB arrive on all ranks running Arbor and vice versa (Figure 4B). This allows co-simulation in conjunction with arbitrary numbers of ranks on both sides and even in compounds with more than two simulators.

Finally, bi-directional translation between TVB's variable concept and Arbor's representation of action potentials is required. As the former depends on the region models used, we chose to bundle this with the remaining TVB functionality as part of the Arbor proxy node. For the TVB models used in this study, the main variable is the per-region mean activity rate ν_I which is conceptually compatible with the concept of spike generation. For each region I connected to the proxy node P, i.e., with connectome weight c_IP > 0, a set of synthetic events must be generated such that the mean activity conforms to ν_I. This is an ambiguous process, even if we prescribe a population (list of cell identifiers) and a per-cell distribution, e.g., a Poisson point process, from which to draw events, which likely must be resolved by ensembles of simulations. In general, this is both model dependent and mathematically intractable, so we leave the general case as a customization point in the framework. For our running example, however, we make the following choice: Event timings for the current step k will be drawn from a uniform distribution and dispatched to all cells in the Arbor network. Note that while these events are created at given time a per-connection delay is applied and thus delivery occurs at a later time.

The inverse direction, converting spike events to mean rates, while being well-defined, is still subject to customization. We explore two options here. First, simple running averages, i.e., all spikes that originate within the Arbor network during the current epoch, are collected and sorted into bins of width Δt. This list is then normalized to the cell count and time step and sent to TVB as the mean activities as a function of time. Although straightforward, this can lead to unrealistically rough activity traces, especially if cell populations are small. Second, as inspired by high-speed calcium imaging experiments (Grewe et al., 2010), a mechanism to track cell activity through calcium level is implemented as:

with per cell p having a decay parameter τ and a weight β. Computing the activity becomes the average:

yielding a smooth trace. This method of converting discrete spiking events into a continuous interval variable is also used in a few plasticity models recruiting a negative feedback control mechanism such as synaptic scaling (Van Rossum et al., 2000) and homeostatic structural plasticity (Butz and Van Ooyen, 2013; Diaz-Pier et al., 2016; Lu et al., 2024) models. The choice of the calcium kernel parameters is based on trial-and-error to match the output of Arbor and the magnitude of activity generated by TVB. Figure 5 compares the impact of this choice on the macro-scale network. In small networks and over short timescales defined by the epoch length as shown in the example, spiking activity occurs in noncontinuous bursts, which is dubious in conjunction with the smooth dynamics of the chosen TVB model. Figure 5A shows the propagation of this noncontinuous activity into the TVB regions, while using the Ca-like model (B) provides smooth dynamics in both the Arbor and TVB models. We thus will use the latter in all simulations from here on out. In general, both methods require scaling by the number of cells in the proxy region to arrive at a scale-free activity measure. A local scaling factor G_A is used to convert between the activity of the detailed network and the activity of the region modeled in TVB. In general, G_A needs to be adjusted to the choice of connectome and TVB model, similar to the choice of the global coupling strength G in the RWW model. In this study, G_A = 100 is used as it produces seizure-like propagation patterns comparable to those found in similar studies, see e.g., Melozzi et al. (2017) and Courson et al. (2024) and references therein.

As an initial step, we examined a network of morphologically detailed cells (comprising of a single-compartment soma and a random dendritic tree, see bottom-right inset panel of Figure 6) without embedding them into a co-simulation framework. The objective was to assess how interactions among neurons exhibiting distinct firing patterns influence network dynamics within a small population of SLE neurons at the local (Arbor) level. We constructed a network comprising 80 tonic-spiking neurons (K_bath = 9.5mM) and 20 SLE neurons (K_bath = 17.0mM), initially without internal synaptic connectivity. The system was simulated for T = 5s, and the resulting membrane potentials are presented in Figure 6 (this initial phase is indicated by the gray-shaded region). During this initial phase, all neurons independently showed their intrinsic firing behavior, consistent with the activity shown in Figure 6. Following this baseline simulation, integration was paused, the network was reconfigured to full connectivity, and the simulation was resumed. The introduction of connectivity produced an immediate response across the network. The SLE neurons progressively transitioned toward tonic-like spiking, mirroring the dominant (80%) tonic-spiking subpopulation. At the same time, tonic-spiking neurons began to exhibit the bursting activity characteristic of the SLE phenotype (Figure 6, T ≥ 5s)). Note that the use of different K_bath values for neurons within the same group is, at this stage, driven by computational considerations and the need to validate the performance of the Arbor-TVB co-simulator. This choice can be adjusted to study more biologically realistic scenarios which involve the interaction of different neural populations or transition of network dynamics, which is further addressed in the Discussion section.

Next, we embedded a network of detailed neurons as a proxy node in a simulation of neural mass models in TVB. The general setup is similar as before, however, the proxy node now consists of only SLE-type point neurons in a fully connected network. Figure 7A illustrates the evolution of the membrane potential of individual neurons within the Arbor network, with all neurons exhibiting identical dynamical behavior. The pattern of SLE activity is modulated by neuronal coupling.

We run simulations over 20s, and investigate the effect of SLE emergence in the network once all brain areas have reached their steady state. The production of these patterns in the Arbor node generates changes in the firing rates of local TVB brain areas. Even though changes in activity occur in most brain areas, these fluctuations occur within different ranges.

In Figure 7B, we show the time-series of TVB nodes' firing rate throughout the simulation, here with N = 100 in the Arbor detailed network. A transient period is necessary before all nodes reach their steady-state. The emergence of recurrent SLE patterns in the Arbor node triggers a mixture of periodic and seizure-like patterns in the dynamics of the TVB nodes. We highlight traces corresponding to the l CA1 region modeled using detailed cells (orange) and the four regions with the highest activity modeled by the neural mass model. The inset shows the highlighted traces, excluding l CA1, during a single period.

In Figure 8, we present the propagation of SLE originating in l CA1 (star marker), represented as a fully connected network of N = 10000 SLE point neurons. Colors on the brain template show the time distance between seizure emergence in the diseased area and seizure arrival in the different brain regions. The systematic detection of SLE onset in other regions of the whole-brain network (if and when present within the simulation period) is performed as follows. First, the baseline firing activity of each region is defined (first time window of the simulation). The non-SLE nodes initially exhibit a relatively stable firing rate, which gradually becomes influenced by the SLE originating from the Arbor node. When connected to the Arbor node, spiking-tonic brain areas are repeatedly recruited into SLE activity patterns, followed by periods of relaxation. For each region, we define the baseline firing rate as the minimum firing rate observed during the relaxation phase immediately preceding any event in the onset region. When a robust increase in firing rate is detected in a region—indicative of a seizure-like event—the corresponding onset timestamp is recorded. The baseline firing rate serves as a reference, and seizure activity is identified based on a sustained increase in firing rate over a defined duration. We define the onset of a seizure as the beginning of the high-amplitude bursting regime, specifically the point at which the firing rate begins to rise consistently after the first small activity peak. To ensure that only significant deviations are classified as seizures, we require the firing rate to exceed the baseline by at least 10% and remain above this threshold for a minimum of 100 ms. Our seizure onset detection follows a similar (though not identical) approach to that used in related studies—for example, Melozzi et al. (2017) and Courson et al. (2024)—where seizure-like or bursting activity in mouse brain models is identified by monitoring a model variable and applying a threshold to detect the onset. In Figure 8, we also depict the activity time-series of four initially non-SLE nodes of the brain network being recruited in the seizure, namely the ventral part of left Lateral Septal Nucleus (l LSv), l CA3, l ENTCl and r ENTCl. In the main panel (mouse brain template), non-colored regions correspond to areas where seizures either did not occur or had relatively weak effects.

The Arbor network was initialized with fixed connection weights (w = 0.5). To assess robustness, we also tested weights drawn from a normal distribution (μ = 0.5, σ = 0.5), truncated to positive values. Neural activity was averaged over 20 independent realizations (see Supplementary Figures S2, S3). Experiments with inhibitory/excitatory network variants are also currently in progress. Both baseline and seizure-evoked activity differ across brain regions. Due to the symmetrical inter-hemispheric connections in the Allen Mouse Brain SC, l ENTCl and r ENTCl share the same baseline firing rate. As the SLE pattern emerges in the left hippocampus, l ENTCl exhibits higher spiking rates.

We next evaluated the performance of the running example on a single Apple M1 (2021) laptop. Arbor was built with MPI and SIMD (Arm Neon/SVE) support, with cells organized into groups of ten to fully exploit SIMD capabilities. The overall runtime consists of four primary components: (1) Arbor model update, (2) Conversion from spikes to rates, (3) TVB model update, and (4) Conversion from rates to spikes.

The Arbor update runs in parallel with the conversions between rates and spikes, as well as the TVB update. During spike exchange, both simulations synchronize, meaning the slower part must wait in the MPI collective, which accounts for the primary time spent in the collective call. Figure 9 illustrates the total runtime of a 10s simulation for the entire model described above, along with the relative contributions, for system sizes ranging from one to 10,000 cells. Notably, at 10,000 cells, nearly all computational time is spent within the Arbor network model. In future experiments, we plan to leverage additional hardware, including GPUs, to accelerate the Arbor side of the simulation. At this scale, TVB and the spike/rate conversions are potential bottlenecks that will require optimization, potentially through TVB's JIT compilation and GPU acceleration. Additionally, further parallelization and porting of the conversion steps to a more performant programming environment remain promising avenues for improvement.