Overview: The Path to Whole-Brain Mean Field

In the Mean Field Approximation Tutorial, we learned the core idea: replace microscopic interactions with macroscopic averages and demand self-consistency. We saw this in binary neurons (La Camera) and variational inference.

What You'll Learn

The Architecture: One Mean-Field, Then Spatial Coupling

A common misconception is that this model has "two nested mean-field approximations." This is not quite right. Here's the actual structure:

Layer 1: Single Neurons → Transfer Function (SEMI-ANALYTICAL)
  - AdEx neuron dynamics → F(νₑ, νᵢ)
  - Analytical: membrane stats (μᵥ, σᵥ, τᵥ) from input rates
  - Fitted: effective threshold Vₑff via polynomial (numerical)
  - This is NOT a mean-field approximation — it's deriving the transfer function

Layer 2: Neuronal Population → Mean-Field (THIS IS THE MEAN-FIELD)
  - 10,000 AdEx neurons → E/I population rates (νₑ, νᵢ)
  - Uses the Markovian Master Equation formalism
  - Tracks means AND variances (cₑₑ, cᵢᵢ, cₑᵢ)
  - El Boustani & Destexhe (2009)

Layer 3: Many Regions → Spatial Coupling (NOT ANOTHER MEAN-FIELD)
  - 68 cortical regions connected by DTI tractography
  - Simple weighted sum: νₑᵢₙ(k) = νₑ(k) + G Σⱼ Cₖⱼ νₑ(j, t-Dₖⱼ)
  - This is just adding terms to the ODE — no additional approximation!

The Bridge: From Binary to Biophysical

La Camera's binary neurons are a beautiful abstraction—each neuron is just 0 or 1, spiking with probability σ(I). But real neurons have:

• Continuous membrane potential that integrates input over time
• Voltage-dependent conductances (synapses depend on V)
• Adaptation currents that change firing properties
• Refractory periods after each spike

The Model Hierarchy

                    ABSTRACTION
         High ←─────────────────────────→ Low
         
Binary Neuron   →   LIF   →   AdEx   →   Hodgkin-Huxley
   xᵢ ∈ {0,1}       dV/dt     + exp spike    full ion channels
   P = σ(I)         + reset   + adaptation   + detailed kinetics
   
            ↓           ↓          ↓              ↓
Mean Field:  Easy       Harder     Semi-Analytical   Very Hard
             (La Camera) (Brunel)   (Zerlaut)         (rare)

Why AdEx?

The Adaptive Exponential Integrate-and-Fire (AdEx) model hits a sweet spot:

Leaky Integrate-and-Fire (LIF)

Before AdEx, let's understand its simpler ancestor: the Leaky Integrate-and-Fire model. This bridges binary neurons to continuous voltage dynamics.

The LIF Equations

Membrane potential dynamics:
  τₘ dV/dt = -(V - Eₗ) + Rₘ × I(t)

Spike mechanism:
  if V ≥ Vₜₕᵣₑₛₕ:
    emit spike
    V → Vᵣₑₛₑₜ for duration τᵣₑf

Parameters:
  τₘ = Cₘ/gₗ     membrane time constant (~15-20 ms)
  Eₗ            leak reversal potential (~-65 mV)
  Vₜₕᵣₑₛₕ       spike threshold (~-50 mV)
  Vᵣₑₛₑₜ        reset voltage (~-65 mV)
  Rₘ = 1/gₗ     membrane resistance

Variable Reference

Key Insight: The f-I Curve

The relationship between input current I and firing rate f is the neuron's transfer function. For LIF with constant input:

This is a closed-form solution! But for AdEx with fluctuating input... we won't be so lucky.

The AdEx Model (Brette & Gerstner, 2005)

The Adaptive Exponential Integrate-and-Fire extends LIF with two key features: an exponential spike mechanism and an adaptation current.

The AdEx Equations

Membrane potential:
  Cₘ dV/dt = gₗ(Eₗ - V) + gₗΔₜ exp((V - Vₜ)/Δₜ) + Iₛᵧₙ(V,t) - Iw

Adaptation current:
  τw dIw/dt = a(V - Eₗ) - Iw + Σₛₚᵢₖₑₛ b δ(t - tₛ)

Spike mechanism:
  when V → Vₜ + 5Δₜ (runaway):
    V → Eₗ (reset)
    Iw → Iw + b (adaptation increment)
    refractory for τᵣₑf

What's New Compared to LIF?

Key Parameters

Adaptation: The Key to Sleep/Wake Dynamics

The adaptation current Iw is not just a detail—it's central to brain state dynamics. Different adaptation levels produce dramatically different network behavior.

What is Spike-Frequency Adaptation?

Each spike increments Iw by b:
  Iw → Iw + b  (after spike)

Iw then decays with time constant τw:
  τw dIw/dt = -Iw + a(V - Eₗ)

Result: More spikes → larger Iw → harder to spike → rate decreases

Biological Basis

Adaptation and Brain States

In the whole-brain model (TVB), the parameter b is varied to simulate different states:

Neuromodulatory Control

In reality, acetylcholine (ACh), norepinephrine, and other neuromodulators control adaptation:

• High ACh (wake/REM): Suppresses adaptation → fast, irregular activity
• Low ACh (NREM): Strong adaptation → slow oscillations

The model captures this by varying b as a proxy for neuromodulatory tone.

The Transfer Function: Heart of Mean Field

The transfer function F(νₑ, νᵢ) maps presynaptic firing rates to output firing rate. It's the core ingredient for any mean-field model.

Definition

Given excitatory input rate νₑ, inhibitory input rate νᵢ, and adaptation level W, what is the neuron's output firing rate?

The Problem: No Closed Form for the Threshold

For binary neurons (La Camera), we had a simple sigmoid: P = σ(I). For LIF with constant input, we had the exact f-I curve.

For AdEx with fluctuating conductance-based input? The membrane potential statistics (μᵥ, σᵥ, τᵥ) can be computed analytically, but the effective threshold Vₑff cannot. Why?

• Exponential spike mechanism — nonlinear voltage dependence
• Adaptation current w — builds up with spikes, affects excitability
• Reset mechanism — discontinuous jump after spike
• Subthreshold adaptation (a term) — voltage-dependent current

The Transfer Function Formula

F(νₑ, νᵢ) = 1/(2τᵥ) × erfc((Vₜₕᵣᵉᶠᶠ - μᵥ) / (√2 σᵥ))

where:
  erfc(x) = complementary error function
  μᵥ     = mean membrane potential
  σᵥ     = standard deviation of V fluctuations
  τᵥ     = autocorrelation time of fluctuations
  Vₜₕᵣᵉᶠᶠ = effective (phenomenological) threshold

This formula comes from the diffusion approximation: treating synaptic input as Gaussian noise and computing first-passage time to threshold.

Membrane Potential Fluctuations: μᵥ, σᵥ, τᵥ

The transfer function depends on the statistics of membrane potential fluctuations. These can be computed analytically from the input rates!

The Three Key Statistics (All About Membrane Potential!)

The Calculation Flow

INPUT              →  INTERMEDIATE (Analytical)  →  OUTPUT
(νₑ, νᵢ) Hz           (μᵥ, σᵥ, τᵥ, μG)              νₒᵤₜ Hz
Input firing rates    Membrane potential stats      Output firing rate

Derivation from Conductance Fluctuations (All Analytical!)

Step 1: Conductance statistics (from Campbell's theorem)
  μGₑ = νₑ × Kₑ × τₑ × Qₑ    (mean excitatory conductance, in nS)
  μGᵢ = νᵢ × Kᵢ × τᵢ × Qᵢ    (mean inhibitory conductance, in nS)
  μG  = μGₑ + μGᵢ + gₗ       (total conductance)
  
Step 2: Mean membrane potential (weighted average of reversal potentials)
  μᵥ = (μGₑ×Eₑ + μGᵢ×Eᵢ + gₗ×Eₗ) / μG    (in mV)

Step 3: Variance (from power spectrum integration)
  σᵥ² = Σₛ Kₛνₛ(Uₛτₛ)² / (2(τₘᵉᶠᶠ + τₛ))
  where Uₛ = Qₛ(Eₛ - μᵥ)/μG is the driving force per synapse

Step 4: Autocorrelation time
  τᵥ = [Σₛ Kₛνₛ(Uₛτₛ)²] / [Σₛ Kₛνₛ(Uₛτₛ)²/(τₘᵉᶠᶠ + τₛ)]

The Fluctuation-Driven Regime

Real cortical neurons operate in the fluctuation-driven regime:

• μᵥ is below threshold (no guaranteed spiking)
• σᵥ is large enough that fluctuations occasionally reach threshold
• Spikes are triggered by random coincidences of excitatory inputs

The Effective Threshold: A Phenomenological Trick

The transfer function formula uses an effective threshold Vₜₕᵣᵉᶠᶠ rather than the biophysical threshold Vₜₕ. Why?

The Problem

The simple formula F = erfc((Vₜₕ - μᵥ)/(√2 σᵥ))/(2τᵥ) assumes:

• Linear subthreshold dynamics
• Sharp threshold
• No adaptation
• No reset effects

AdEx violates all of these! The exponential spike, reset, and adaptation all affect the apparent firing rate.

The Solution: Fit an Effective Threshold

Instead of trying to account for each effect analytically, we let the data speak:

1. Simulate AdEx at many (νₑ, νᵢ, W) configurations
2. Compute (μᵥ, σᵥ, τᵥ) analytically
3. Find Vₜₕᵣᵉᶠᶠ that makes F = erfc(...) match the simulated rate
4. Fit Vₜₕᵣᵉᶠᶠ as a polynomial in (μᵥ, σᵥ, τᵥ, μG)

The Polynomial Threshold (Zerlaut 2018)

Vₜₕᵣᵉᶠᶠ(μᵥ, σᵥ, τᵥᴺ, μG) = P₀ 
  + Σₓ Pₓ × (x - x₀)/δx₀
  + PμG × log(μG/gₗ)
  + Σₓ,ᵧ Pₓᵧ × (x-x₀)/δx₀ × (y-y₀)/δy₀

where x, y ∈ {μᵥ, σᵥ, τᵥᴺ}

Fitted Coefficients (from Table 2 of Zerlaut 2018)

Consequences and Limitations of the Effective Threshold

What It Captures Well:

• All the "messy" nonlinear dynamics compressed into one state-dependent value
• Effects of spike-frequency adaptation (at steady state)
• Effects of the exponential spike mechanism
• Conductance-dependent threshold changes

What It Misses (Key Limitations):

Population Dynamics: E/I Networks

Now we can build a population-level mean-field model. Instead of 10,000 individual neurons, we track just a few variables: the excitatory rate νₑ, the inhibitory rate νᵢ, and the adaptation Wₑ.

The Mean-Field Population Equations

First-order (rate) dynamics:
  T dνₑ/dt = Fₑ(νₑ, νᵢ, Wₑ) - νₑ
  T dνᵢ/dt = Fᵢ(νₑ, νᵢ, Wᵢ) - νᵢ

Population adaptation:
  dWₑ/dt = (-Wₑ + bₑνₑ + aₑ(μᵥ - Eₗ)) / τw

where:
  T   ≈ 5 ms (characteristic time scale)
  Fμ  = transfer function for population μ
  Wₑ  = mean adaptation current in E population

Network Architecture

The model assumes a random sparse network:

• N = 10,000 neurons (8,000 E, 2,000 I)
• Connection probability ε = 5%
• Ke = ε × 8000 = 400 excitatory inputs per neuron
• Ki = ε × 2000 = 100 inhibitory inputs per neuron

The Markovian Formalism (El Boustani & Destexhe, 2009)

The first-order equations (just rates) miss something important: finite-size fluctuations. With only 10,000 neurons, there's significant noise in population activity.

Second-Order Description

The Markovian formalism tracks not just means but also variances and covariances:

Five quantities (not six — see explanation below):
  νₑ(t), νᵢ(t)           - mean population rates
  cₑₑ(t), cᵢᵢ(t)         - variances of E and I population activity
  cₑᵢ(t)                 - covariance between E and I

Evolution:
  T ∂νμ/∂t = (Fμ - νμ) + ½ cλη ∂²Fμ/∂νλ∂νη
  
  T ∂cλη/∂t = Aλη + (Fλ - νλ)(Fη - νη)
              + cλμ ∂Fη/∂νμ + cμη ∂Fλ/∂νμ - 2cλη

where Aλη = Fλ(1/T - Fλ)/Nλ  (finite-size noise)

Understanding cₑₑ, cᵢᵢ, cₑᵢ (And Why No cᵢₑ?)

⎛ cₑₑ  cₑᵢ ⎞
⎝ cₑᵢ  cᵢᵢ ⎠

Why Second Order?

Whole-Brain Dynamics: TVB and the Connectome

The final layer: embedding mean-field nodes in a spatial network based on real human brain connectivity.

The Whole-Brain Equations

For each brain region k:

  νₑᵢₙ(k,t) = νₑ(k,t) + νᵃᶠᶠ(k,t) + G Σⱼ Cₖⱼ νₑ(j, t-Dₖⱼ)
  
  T dνₑ(k)/dt = Fₑ(νₑᵢₙ, νᵢ, Wₑ) - νₑ(k)
  T dνᵢ(k)/dt = Fᵢ(νₑᵢₙ, νᵢ, Wᵢ) - νᵢ(k)

where:
  Cₖⱼ = connection strength from region j to k (from DTI)
  Dₖⱼ = propagation delay (from tract length)
  G   = global coupling strength
  νᵃᶠᶠ = afferent input (external drive + noise)

Components of the Whole-Brain Model

The Virtual Brain (TVB)

The Virtual Brain is an open-source platform that implements this whole-brain modeling framework. It provides:

• Pre-computed connectomes from human subjects
• Various neural mass/mean-field models (including AdEx)
• Tools for comparing simulated and empirical data
• GPU acceleration for large-scale simulations

Simulating Brain States

By varying the adaptation parameter b across regions, researchers can simulate:

Sleep Stage   →   Adaptation b   →   Network Dynamics
───────────────────────────────────────────────────────
Wake          →   ~40 pA         →   Desynchronized, high entropy
REM           →   ~60 pA         →   Irregular, dream-like
N2            →   ~80 pA         →   Spindles, K-complexes
N3            →   ~95 pA         →   Slow oscillations, low entropy

Connection to Your Data (Chennu)

The Chennu dataset provides iEEG recordings during sleep. The whole-brain model can generate synthetic signals that match the complexity-entropy characteristics of different sleep stages. By fitting the adaptation parameter, researchers can:

• Reproduce the WPE × SCM plane structure
• Predict transitions between states
• Test mechanistic hypotheses about sleep

Quick Reference

Key Equations

Key Parameters (Table from Zerlaut 2018)

Key Papers

• Brette & Gerstner (2005) - Original AdEx model
• El Boustani & Destexhe (2009) - Markovian formalism
• Zerlaut et al. (2018) - Semi-analytical mean field for AdEx
• Di Volo et al. (2019) - Mean field with adaptation
• Goldman et al. (2023) - TVB + AdEx for sleep

The Three-Layer Architecture (Corrected Understanding)

The Journey Summarized