Neural Mass Models

Background reading: what neural mass models are and how the Jansen-Rit model works

The Problem: Simulating Brain Activity

A human brain has roughly 86 billion neurons. Each neuron connects to thousands of others through synapses, producing electrical signals that we can measure with electrodes. When we record SEEG (stereotactic EEG), we are measuring the combined electrical activity of millions of neurons near each electrode contact.

We would like to simulate this activity — to build a mathematical model that produces signals that look like real brain recordings. Why? Because if we can reproduce what we observe, we can understand the mechanisms behind it. For epilepsy, this means understanding what makes certain brain regions generate seizures.

But simulating 86 billion individual neurons is computationally impossible for most purposes. We need a shortcut.

The Shortcut: Populations Instead of Individual Neurons

A neural mass model takes a radical simplification: instead of modeling each neuron individually, it models the average behavior of a large population of similar neurons. Think of it like weather forecasting — you do not track every air molecule; you track temperature and pressure as average properties of large volumes of air.

In the brain, neurons of the same type that are physically close together tend to behave similarly. So we group them into populations and describe each population with just a few numbers:

• Mean membrane potential — the average electrical state of all neurons in the group
• Mean firing rate — how many action potentials (spikes) the group produces per second, on average

This reduces millions of neurons to a handful of variables. A single brain region might be described by 3 populations and 6 variables — instead of millions of individual neuron states.

Two Key Ingredients

Every neural mass model is built from two fundamental operations:

1. The Post-Synaptic Potential (PSP) Operator

When a neuron fires, it sends a signal to connected neurons. The receiving neuron does not respond instantly — the incoming spike is transformed into a smooth, gradual voltage change called a post-synaptic potential. This is modeled as a simple second-order linear system (like a spring that gets pushed and then wobbles back to rest):

Second-order linear system (one per population):

  d²y/dt² = A · a · input - 2a · dy/dt - a² · y

  A = maximum amplitude of the response (how big the voltage bump is)
  a = time constant (how fast it rises and decays)
  y = the resulting post-synaptic potential

An incoming spike goes in; a smooth voltage curve comes out. Excitatory PSPs push the voltage up (making neurons more likely to fire), while inhibitory PSPs push it down (making them less likely to fire).

2. The Sigmoid Function

The other direction: given a population's average membrane potential, how fast does it fire? The relationship is not linear — at very low voltages, neurons barely fire at all. At very high voltages, they fire at their maximum rate. In between, there is a smooth S-shaped transition. This is the sigmoid function:

S(v) = 2 · e0 / (1 + exp(r · (v0 - v)))

  e0 = half the maximum firing rate
  v0 = voltage at which firing is half-maximal
  r  = steepness of the transition

Low voltage in → near-zero firing rate out. High voltage in → maximum firing rate out. The sigmoid converts between the two "languages" of the model: voltage and firing rate.

The Jansen-Rit Model

The Jansen-Rit model (1995) is the most widely used neural mass model. It describes a single cortical column — a small vertical cylinder of brain tissue about 2 mm wide — using three interconnected populations:

How They Connect

The three populations form a loop:

graph TD
  INPUT(["External input p(t)"])
  EXC["Excitatory Interneurons\n(amplify signal)"]
  PYR["**Pyramidal Cells**\n(main output: y1 − y2)"]
  INH["Inhibitory Interneurons\n(dampen signal)"]

  INPUT -->|"p(t)"| EXC
  PYR -->|"C1"| EXC
  EXC -->|"C2 · excitatory"| PYR
  PYR -->|"C3"| INH
  INH -->|"C4 · inhibitory"| PYR

  style EXC fill:#0a3d1a,stroke:#4ade80,color:#4ade80
  style PYR fill:#162450,stroke:#5b8dee,color:#5b8dee
  style INH fill:#3d0a0a,stroke:#f87171,color:#f87171
  style INPUT fill:#2a1850,stroke:#a78bfa,color:#a78bfa

The C values (C1, C2, C3, C4) control how strongly the populations talk to each other. They are all derived from a single connectivity constant C:

C1 = C          # pyramidal → excitatory interneurons
C2 = 0.8 · C    # excitatory interneurons → pyramidal
C3 = 0.25 · C   # pyramidal → inhibitory interneurons
C4 = 0.25 · C   # inhibitory interneurons → pyramidal

The External Input

The excitatory interneurons receive an external input p(t) representing activity arriving from other brain regions and sensory pathways. This input is modeled as random noise with a mean value and some variability:

p(t) = mean_input + noise

# Typically: mean = 220 Hz, noise std = 22 Hz
# This represents the average background bombardment
# from the rest of the brain.

The noise is what keeps the model "alive" — without it, the system would settle to a fixed point and stop oscillating.

The Output

What we measure with an electrode is the difference between excitatory and inhibitory post-synaptic potentials arriving at the pyramidal cells:

output(t) = y1(t) - y2(t)

# y1 = excitatory PSP on pyramidal cells
# y2 = inhibitory PSP on pyramidal cells
# The difference is what the electrode picks up

This is a direct analog of real EEG/SEEG: the signal we record at electrodes is generated primarily by the summed post-synaptic potentials of pyramidal cells.

The Full Equations

Each population has a PSP operator (second-order ODE), which we split into two first-order ODEs. Three populations × 2 variables each = 6 state variables:

# State vector: [y0, y1, y2, y3, y4, y5]
# y0..y2 are PSP outputs, y3..y5 are their derivatives

# Position derivatives (just linking to velocities)
dy0/dt = y3
dy1/dt = y4
dy2/dt = y5

# Velocity derivatives (the actual dynamics)
dy3/dt = A·a·S(y1 - y2) - 2a·y3 - a²·y0
dy4/dt = A·a·(p(t) + C2·S(C1·y0)) - 2a·y4 - a²·y1
dy5/dt = B·b·C4·S(C3·y0) - 2b·y5 - b²·y2

# Output (simulated SEEG signal)
output = y1 - y2

The Parameters

What It Produces

At the default parameter values, the Jansen-Rit model produces a signal with a dominant frequency around 10 Hz — the alpha rhythm, the most prominent rhythm in human EEG. This matches what we observe when a person is relaxed with eyes closed.

By changing the mean input p, the model transitions through different regimes:

Parameter Regime Map

This plot sweeps the mean input p from 0 to 500 Hz and shows the RMS amplitude of the model output at each value. The three behavioral regimes are visible as distinct regions of the curve:

How It Is Solved Numerically

The 6 ODEs cannot be solved with a closed-form formula — they must be integrated numerically. Starting from an initial state (all zeros), the solver takes tiny time steps (typically 0.1 ms) and computes the state at each step using the derivatives.

Montage Concord uses Runge-Kutta 4 (RK4) for this. RK4 evaluates the derivatives at four points within each time step and takes a weighted average, giving much better accuracy than simply stepping forward by the derivative (Euler's method). The noise enters through the external input p(t), which is pre-generated as a random time series and passed to the derivative function at each step.

After integration at the fine time step (0.1 ms = 10,000 Hz), the output is downsampled to the requested sampling rate (e.g. 256 Hz or 1024 Hz) to match typical EEG recording rates.

Where It Fits in the Concord Pipeline

In Montage Concord, the Jansen-Rit model implements the Model ABC defined in concord-core. The simulation pipeline is:

The ModelOutput contains a simulated signal in the same format as a real Recording, so the same metrics and visualization tools work on both. This is how model fitting will work: simulate, compute metrics on the simulation, compare to metrics on real data, adjust parameters, repeat.

Beyond Jansen-Rit

Jansen-Rit is the foundation. The planned model packages extend it in different directions:

All share the same infrastructure: the Model ABC, models-utils (integrators, sigmoids), and the ParameterVector → ModelOutput data flow.

The Gap Between Model and Reality

The Jansen-Rit model is a powerful explanatory tool, but it is not a literal replica of real EEG. Understanding where the model diverges from measured data is essential — both for interpreting model output honestly and for designing fitting procedures that know what they can and cannot recover.

1/f Noise: The Missing Background

Real EEG power spectra follow a characteristic 1/f shape: power decreases roughly as the inverse of frequency, forming a straight line on a log-log plot. This "pink noise" floor underlies all oscillatory peaks. It is thought to arise from the superposition of many independent processes at different timescales — synaptic noise, dendritic filtering, volume conduction through tissue.

The Jansen-Rit model does not produce 1/f noise. Its output has oscillatory peaks (driven by the excitatory-inhibitory loop) sitting on top of the white noise injected via p(t). White noise has flat power across all frequencies — the opposite of the 1/f slope seen in real data.

This matters in practice because many EEG metrics — band power, spectral edge frequency, 1/f exponent itself — are shaped as much by the aperiodic background as by the oscillatory peaks. A model that matches the alpha peak but ignores the 1/f floor will produce misleading metric values. Common strategies for dealing with this include:

• Fitting metrics that are robust to the background — relative band power, peak frequency, and spectral ratios are less sensitive to the 1/f floor than absolute power.
• Adding 1/f noise post-hoc — after simulation, synthetic 1/f noise can be mixed into the model output before metric comparison. The 1/f exponent and amplitude become additional parameters to fit.
• Separating periodic and aperiodic components — spectral parameterization (e.g., FOOOF) isolates oscillatory peaks from the 1/f background, allowing fitting to target the periodic component only.

Non-Identifiability: Multiple Solutions to the Same Data

The Jansen-Rit model has 10 free parameters. The measured data it must match is typically summarized by a handful of features: dominant frequency, relative band powers, waveform statistics. There are far more degrees of freedom in the model than constraints in the data. This means the problem is non-identifiable — multiple different parameter combinations can produce nearly identical output features.

Non-identifiability is not a bug — it is an inherent property of fitting a complex model to summary statistics. Strategies for dealing with it include:

A well-designed fitting engine should account for this reality. Rather than returning a single "answer," it is valuable to expose the loss landscape and use multi-start optimization to map out the space of viable solutions.

What the Model Cannot Represent

Some aspects of real SEEG data are structurally absent from the Jansen-Rit model. These are not parameter-tuning problems — they are limitations of the model's mathematical structure:

Implications for Model Fitting

These gaps have practical implications for any fitting procedure built on top of the Jansen-Rit model:

1. Metric selection matters more than optimizer choice. Choosing which features to match (PSD peaks vs. broadband power vs. time-domain statistics) implicitly decides which aspects of the data the model is asked to explain — and which are ignored.
2. Residual structure is informative. When the best-fit model still does not match certain features, that mismatch tells you something — the signal contains components (1/f noise, HFOs, artifacts) outside the model's scope.
3. Confidence requires exploration. A single optimized parameter set is a point estimate. To know how much to trust it, you need to explore the neighborhood — are there other parameter sets that fit nearly as well? How sensitive is the fit to small perturbations?
4. Pre-processing is part of the model. Re-referencing scheme, filtering, artifact rejection, and epoch selection all shape the data the model sees. Two different pre-processing pipelines applied to the same recording can yield different "best-fit" parameters.

Further Reading

• Jansen, B.H. & Rit, V.G. (1995). Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns. Biological Cybernetics, 73(4), 357-366.
• David, O. & Friston, K.J. (2003). A neural mass model for MEG/EEG: coupling and neuronal dynamics. NeuroImage, 20(3), 1743-1755.
• Wendling, F., et al. (2002). Epileptic fast activity can be explained by a model of impaired GABAergic dendritic inhibition. European Journal of Neuroscience, 15(9), 1499-1508.
• Donoghue, T., et al. (2020). Parameterizing neural power spectra into periodic and aperiodic components. Nature Neuroscience, 23(12), 1655-1665. [Spectral parameterization / FOOOF]
• He, B.J. (2014). Scale-free brain activity: past, present, and future. Trends in Cognitive Sciences, 18(9), 480-487. [1/f noise in neural signals]
• Raue, A., et al. (2009). Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood. Bioinformatics, 25(15), 1923-1929. [Non-identifiability in ODE models]