Wilson-Cowan Tutorial

Our Journey

Binary Neurons

Micro → Macro

→

E/I Populations

Two coupled eqns

→

Continuous Time

Wilson-Cowan ODEs

→

Stochastic

Micro → Langevin

★

Criticality

Power laws

Introduction: Why This Tutorial?

You've learned the Mean Field Approximation from La Camera's paper. The key insight was simple but powerful:

La Camera's Core Move:

Replace actual spikes x_j with expected rates ⟨x_j⟩ = r_j

This gave us a self-consistent equation for the firing rate:

r = σ(J × r + I_ext)

Now we want to understand the stochastic Wilson-Cowan model used by de Candia et al. to study criticality. But their equations look different! This tutorial bridges the gap.

The Core Pattern: Throughout this tutorial, we follow a microscopic → macroscopic structure:
• Microscopic: Individual neurons with discrete states (aᵢ ∈ {0,1})
• Macroscopic: Population averages and firing rates (r, S)

This is the essence of mean field theory: start with individual neurons, derive equations for collective behavior.

The Questions We'll Answer

How does the microscopic → macroscopic derivation work?
Why do we need separate E and I populations?
What changes when we go from discrete to continuous time?
Where does the noise term come from (finite-size effects)?
What is the critical point and why does it matter?
How does this relate to co-transmission (your ultimate goal)?

Binary Neurons — La Camera's Foundation

La Camera (2022), Section 2-4

The Model: What We Start With

We have a network of N binary neurons. Each neuron i is either spiking (x_i = 1) or silent (x_i = 0).

Variables Introduced

x_i	State of neuron i: 0 (silent) or 1 (spiking)
J_ij	Synaptic weight from neuron j to neuron i
I_ext	External input current (from outside the network)
I_i	Total input to neuron i: I_i = Σ_j J_ijx_j + I_ext

The Spiking Rule (Transfer Function)

Each neuron spikes with a probability that depends on its input via a sigmoid transfer function:

P(x_i = 1 | I_i) = σ(I_i) = 1 / (1 + e^{−2β(I_i − θ)})

New Variables

β	Inverse temperature or gain. Controls how “sharp” the decision is. High β → nearly deterministic. Low β → nearly random.
θ	Threshold. The input value where σ(I) = 0.5. Often set to 0 for simplicity, but in La Camera's formulation it appears inside the sigmoid.

The Threshold θ: When I = θ, the sigmoid outputs exactly 0.5 (maximum uncertainty). For I > θ, spiking becomes more likely; for I < θ, it becomes less likely. Setting θ = 0 is common but not required.

The Mean Field Approximation

The key move: instead of tracking individual spikes, we track the expected firing rate:

r_i = ⟨x_i⟩ = P(x_i = 1) = σ(I_i)

New Variable

r_i	Firing rate of neuron i. This is what mean field theory tracks. It equals the probability that the neuron is spiking at any given moment.

For a homogeneous network (all neurons identical, uniform weights J/N):

Homogeneous Mean Field:

r = σ(J × r + I_ext − θ)

Self-Consistency: The rate r appears on both sides of the equation! We need to find the value where the output equals the input. This is a fixed point.

Try It: Find the Fixed Point

J (coupling):I_ext:β (gain):θ (threshold):

Stage 1 Summary

Binary neurons with probabilistic spiking: P(spike) = σ(input − θ)
The sigmoid has two key parameters: β (gain) and θ (threshold)
Mean field replaces spikes with rates: x_j → r_j
Self-consistent equation: r = σ(J × r + I_ext − θ)
This is a single homogeneous population

E/I Populations — The Wilson-Cowan Step

Chow & Karimipanah (2020)

Why Separate E and I?

La Camera's homogeneous network is a simplification. Real brains have two fundamentally different types of neurons:

Excitatory (~80%)

Release glutamate
Increase postsynaptic activity
Positive effect on target

Inhibitory (~20%)

Release GABA
Decrease postsynaptic activity
Negative effect on target

Dale's Law: A neuron is either excitatory OR inhibitory — it can't switch. This is why we need separate populations!

The Weight Matrix

New Variables: The Weight Matrix

w_EE	E → E: Excitatory neurons exciting each other (recurrent excitation)
w_EI	I → E: Inhibitory neurons inhibiting excitatory neurons
w_IE	E → I: Excitatory neurons driving inhibitory neurons
w_II	I → I: Inhibitory neurons inhibiting each other

Weight Matrix: W =

+w_EE−w_EI

+w_IE−w_II

Two Coupled Mean Field Equations

New Variables

r_E	Firing rate of the excitatory population (fraction of E neurons active)
r_I	Firing rate of the inhibitory population (fraction of I neurons active)

E/I Mean Field Equations:

r_E = σ(w_EE × r_E − w_EI × r_I + I_E)

r_I = σ(w_IE × r_E − w_II × r_I + I_I)

Same structure as La Camera! Just two coupled equations instead of one. The firing rate of each population depends only on the current firing rates of both populations. We iterate until self-consistent.

E/I Balance

The brain maintains a delicate balance between excitation and inhibition:

Too much E → Seizures, runaway activity

Too much I → Silence, no computation

Balanced → Rich dynamics, optimal processing

Try It: E/I Fixed Point

w_EE: w_EI: w_IE: w_II: I_E: I_I:

Stage 2 Summary

Separate E and I populations (Dale's Law)
Four connection weights: w_EE, w_EI, w_IE, w_II
Two coupled self-consistent equations
Firing rates depend on current state only (no memory/hysteresis yet)
E/I balance is crucial for brain function

Continuous Time — The Wilson-Cowan ODEs

Chow & Karimipanah (2020), Wilson & Cowan (1972)

From Discrete to Continuous Time

In Stage 2, the firing rate at time n+1 depends only on the rates at time n. But real neurons have temporal dynamics — their activity doesn't instantly snap to new values. Wilson and Cowan introduced differential equations to capture this:

Key Change from Stage 2: In Stage 2, the new firing rate immediately equals σ(input). In Stage 3, the firing rate approaches σ(input) gradually, with a time constant τ. This gives the system memory/hysteresis — the firing rate depends on recent history, not just the current input.

New Variables

τ_E, τ_I	Time constants (in ms). How fast each population responds. Typically τ_E ≈ 10ms, τ_I ≈ 20ms (I neurons are often slower).
dr/dt	Rate of change of firing rate. This is the derivative — how quickly r changes.

Wilson-Cowan Rate Equations:

τ_E × dr_E/dt = −r_E + σ(w_EEr_E − w_EIr_I + I_E)

τ_I × dr_I/dt = −r_I + σ(w_IEr_E − w_IIr_I + I_I)

Understanding the Terms

−r_E (decay term)This represents exponential decay. Without any input, the firing rate decays toward zero. This is NOT just “negative firing rate” — it's the mathematical form of relaxation. The equation dr/dt = −r/τ has solution r(t) = r₀ × e^−t/τ.

+σ(input) (drive term)Input pushes the firing rate toward σ(input). Combined with the decay term, the rate approaches σ(input) exponentially.

τ (time constant)Sets the speed of dynamics. Larger τ → slower response. With a differential equation, you can choose your integration time step (dt) explicitly — smaller dt gives more accurate simulation.

At Equilibrium (Fixed Point): When dr/dt = 0, we get:

0 = −r + σ(input) → r = σ(input)

This is exactly the same fixed point as Stage 2! The difference is the dynamics — how the system gets there. The τ parameter determines how many “previous time steps” influence the current state (the larger τ, the more history matters).

Try It: Watch the Dynamics

w_EE=1.5w_EI=2w_IE=1.5w_II=0.5

τ_E (ms): τ_I (ms): I_E: I_I:

Stage 3 Summary

Continuous time via differential equations: τ × dr/dt = −r + σ(I)
The −r term is decay/relaxation, not “negative rate”
Time constants τ set the dynamics speed and introduce memory
At equilibrium, we get the same fixed points as Stage 2
The difference: how the system approaches equilibrium matters for transient behavior

Stochastic Wilson-Cowan — de Candia's Model

de Candia et al. (2021)

The Two-Level Structure (Just Like Stage 1!)

Remember in Stage 1, we started with individual binary neurons (microscopic), then derived the mean field equations (macroscopic)? de Candia does the same thing, but with a crucial difference: they keep track of finite-size fluctuations.

The Two Levels:
1. Microscopic: N individual neurons, each with state aᵢ ∈ {0, 1}
2. Macroscopic: Population activity S = (1/N) × Σ aᵢ (the mean field)

The deterministic models (Stages 1-3) assume N → ∞ so fluctuations vanish. de Candia keeps N finite, so we get stochastic dynamics!

Part A: The Microscopic Model (Individual Neurons)

Just like La Camera's binary neurons, de Candia starts with N individual neurons. Each neuron i has a binary state:

Variables Introduced (Microscopic Level)

a_i	State of neuron i: 0 (quiescent) or 1 (active/firing)
s_i	Total synaptic input to neuron i: s_i = Σ_j w_ija_j + h
w_ij	Synaptic weight from neuron j to neuron i
h	Small external input (keeps system from getting stuck at all-silent)

de Candia Eq. 1 (Synaptic Input):

s_i = Σ_j w_ij × a_j + h

The Transition Rates (Stochastic Dynamics)

Unlike La Camera's instantaneous updates, de Candia uses a continuous-time Markov process. Neurons don't flip instantly — they transition at specific rates:

Active → Quiescent (1 → 0)Rate = α (constant). An active neuron spontaneously becomes quiescent. This represents spike termination and refractory effects.In time dt: P(deactivate) = α × dt

Quiescent → Active (0 → 1)Rate = f(s_i) (input-dependent). A quiescent neuron becomes active based on its synaptic input. This is where the “activation function” enters.In time dt: P(activate) = f(s_i) × dt

Activation Function (de Candia Eq. 2)

f(s)	Activation rate: f(s) = β·tanh(s) for s > 0, else 0 Unlike the sigmoid in Stage 1, this is zero for negative input — no spontaneous firing without positive drive!
α	Deactivation rate (1/τ). Typically α = 0.1 ms⁻¹ (τ = 10 ms)
β	Activation gain. Typically β = 1.0 ms⁻¹

Key Difference from La Camera: In La Camera's model (Stage 1), neurons have probability σ(I) of being in state 1 at any instant. In de Candia's model, neurons have persistent states that change via stochastic transitions. An active neuron stays active until it deactivates (at rate α). This creates temporal correlations!

The Microscopic Simulation (What the Code Does)

For each neuron i at each timestep dt:

if (neuron[i] is ACTIVE) {
    // Deactivate with probability α × dt
    if (random() < α × dt) {
        neuron[i] = QUIESCENT;
    }
} else {
    // QUIESCENT neuron: compute input
    s_i = Σⱼ wᵢⱼ × aⱼ + h;
    
    // Activate with probability f(s_i) × dt
    if (random() < f(s_i) × dt) {
        neuron[i] = ACTIVE;
    }
}

This is exactly what our SpikingWilsonCowan class does in the CriticalityProject!

Part B: The Mean Field Approximation

Now comes the key step — just like in Stage 1! We go from individual neurons to population averages. But unlike Stage 1, we keep the finite-size fluctuations.

Population Activity (Definition):

S = (1/N) × Σ_i a_i

S is the fraction of active neurons — exactly like r in La Camera's model! But how does S evolve? Let's count the transitions:

Rate of decrease (active → quiescent)There are N×S active neurons, each deactivating at rate α.
Expected decrease: α × S (per unit time)

Rate of increase (quiescent → active)There are N×(1−S) quiescent neurons, each activating at rate f(s).
Expected increase: (1−S) × f(s) (per unit time)

Deterministic Mean Field (Expected Value):

dS/dt = −αS + (1−S)f(s)

The Isomorphism with Stage 1: Compare this to La Camera!
La Camera: r = σ(J×r + I_ext) [self-consistent fixed point]
de Candia: dS/dt = −αS + (1−S)f(s) [time evolution]

At equilibrium (dS/dt = 0): αS = (1−S)f(s), which gives a self-consistent equation! The difference is that La Camera is instantaneous while de Candia has explicit time dynamics.

Part C: Where Does the Noise Come From?

Here's where de Candia differs from La Camera's Stage 1: instead of taking N → ∞ and discarding fluctuations, we keep N finite and derive the noise term.

Finite-Size Fluctuations

Consider N neurons, each firing or quiescent. In time dt:

Some quiescent neurons fire (Poisson process with rate f(s))
Some firing neurons stop (Poisson process with rate α)
Both are random processes — individual neurons flip stochastically
The actual number of transitions fluctuates around the expected value

The Central Limit Theorem Derivation

Let a_i be the state of neuron i (0 or 1). The population average is:

S = (1/N) × Σ a_i

Each a_i is a Bernoulli random variable. By the CLT, for large N:

Var(S) = Var((1/N) × Σ a_i) = (1/N²) × Σ Var(a_i) = (1/N²) × N × σ² = σ²/N

where σ² is the variance of a single neuron's state.

Therefore: Std(S) = σ/√N, so fluctuations scale as 1/√N

For N = 10,000 neurons: fluctuations ~ 1/100 = 0.01 of the mean
For N = 100 neurons: fluctuations ~ 1/10 = 0.1 of the mean (much larger!)

The Variance of Transitions: The noise amplitude isn't just 1/√N — it depends on how many transitions are happening! More transitions = more variance.

Variance = αS + f(s)(1−S)
This is the sum of: (rate of deactivations) + (rate of activations)

The Complete Isomorphism: Microscopic ↔ Macroscopic

Here's the precise correspondence between the two levels of description.They use the same physics — just described at different scales:

One-to-One Correspondence

Aspect	Microscopic (SpikingWilsonCowan)	Macroscopic (StochasticMeanFieldWC)
State variable	a_i ∈ {0, 1} for each neuron i	S ∈ [0, 1] = fraction of active neurons
Interpretation	Binary: active (1) or quiescent (0)	Continuous: population activity rate
Deactivation	P(1→0) = α × dt per neuron	Drift term: −αS
Activation	P(0→1) = f(s) × dt per neuron	Drift term: +(1−S)f(s)
f(s)	β×tanh(s) for s>0, else 0	Same: β×tanh(s) for s>0, else 0
Noise source	Stochastic transitions (Bernoulli)	Langevin: √[αS + f(s)(1−S)] × η(t)
Finite-size	Implicit (N discrete neurons)	Explicit: noise ∝ 1/√N

Why Show Both?
• Microscopic: Where criticality actually emerges — real individual neuron dynamics
• Macroscopic: Faster to compute, analytically tractable, captures essential behavior
• For large N, they produce statistically equivalent results (de Candia Fig. 12)

Complete Variable Reference

S	Fraction of active neurons = (1/N)×Σa_i. Same as firing rate r.
α	Deactivation rate (1→0). Typically α = 0.1 ms⁻¹ (so τ = 10 ms).
β	Activation gain in f(s). Typically β = 1.0 ms⁻¹
w₀	Net coupling strength (see criticality section for full definition).
h	Small external input (prevents absorbing state S=0).
N	Number of neurons. Noise scales as 1/√N.
η(t)	White noise: Gaussian, zero mean, unit variance.

Part D: The Full Langevin Equation (de Candia Eq. 4)

Combining the mean field drift with the finite-size noise gives us the complete Langevin equation:

What is a Langevin Equation?

A Langevin equation is a stochastic differential equation (SDE) of the form:

dX/dt = drift(X) + noise amplitude(X) × η(t)

It combines:

Drift: The deterministic part — where the system would go without randomness
Noise: Random fluctuations that perturb the system. η(t) is “white noise” (Gaussian, uncorrelated in time)

Langevin equations are fundamental in physics (Brownian motion) and are the continuous-time limit of discrete random walks. They're essential for studying systems near criticality!

de Candia Eq. 4 (Full Langevin Form):

dS/dt = −αS + (1−S)f(s) + (1/√N)·√[αS + f(s)(1−S)]·η(t)

−αS (drift: decay)Firing neurons become quiescent at rate α. This pulls S toward 0.

(1−S)f(s) (drift: activation)Quiescent neurons (fraction 1−S) fire at rate f(s). This pushes S up.

(1/√N)·√[...]·η(t)The noise term. η(t) is white noise. The √[αS + f(s)(1−S)] is the noiseamplitude — it depends on both transition rates because fluctuations come from both firing and stopping.

Activation Function

f(s)	Activation rate: f(s) = β·tanh(s) for s > 0, else 0 where s = w₀·S + h (total synaptic input). Unlike the smooth sigmoid, this is zero for negative input, which creates an absorbing state at S=0.

Try It: Stochastic Simulation

α = 0.1 ms⁻¹β = 1.0 ms⁻¹h = 10⁻⁶Critical w₀ᶜ = α/β = 0.1

w₀:0.10N: Duration (ms):

Stage 4 Summary: The Two-Level Structure

Microscopic Model: N neurons with states a_i ∈ {0,1}, transitioning at rates α and f(s)
Macroscopic Model: Population activity S, evolving via Langevin equation
The derivation: Mean field (S = ⟨a⟩) + Central Limit Theorem (noise ∝ 1/√N)
Key isomorphism: P(deactivate) = α×dt → drift term −αS
Key isomorphism: P(activate) = f(s)×dt → drift term (1−S)f(s)
Why both? Microscopic is the “ground truth”; macroscopic is the approximation
Same physics, different levels of description — just like La Camera's Stage 1!

★

Criticality — The Main Event

de Candia et al. (2021), Section “Critical point”

Understanding the Control Parameter w₀

Clarifying w₀: In the full E/I model, we have separate weights w_E (excitatory) and w_I (inhibitory). The effective coupling is their difference:

w₀ = w_E − w_I

In de Candia's simplified model, we work directly with this single parameter w₀. Positive w₀ means excitation dominates; negative w₀ means inhibition dominates.

The Critical Point

At the critical point, the system is perfectly balanced between activity dying out (subcritical) and runaway activity (supercritical).

Critical Condition:

w₀^c = α/β = 0.1 ms⁻¹ / 1 ms⁻¹ = 0.1

Why α/β? At steady state with no noise, dS/dt = 0 means:
αS = (1−S)f(s) = (1−S)β·tanh(w₀S)
For small S near the transition: αS ≈ (1−S)β·(w₀S) ≈ β·w₀·S
The transition between stable S=0 and unstable S=0 occurs when α = β·w₀, i.e., w₀ = α/β.

So: if your activation/deactivation rates satisfy w₀ = α/β, AND your net E/I coupling (w_E − w_I) equals this w₀, then you're at criticality!

The Three Regimes

Subcritical (w₀ < 0.1)

Inhibition dominates

Activity decays to zero
Small, short-lived fluctuations
No sustained dynamics

Critical (w₀ = 0.1)

Perfect E/I balance

Power-law avalanches
Diverging correlation time
Maximal sensitivity

Supercritical (w₀ > 0.1)

Excitation dominates

Runaway activity
Saturation at high S
Seizure-like dynamics

Signatures of Criticality

1. Power-Law Avalanches

Size distribution: P(S) ~ S^−3/2
Duration distribution: P(T) ~ T⁻²

These are the exponents of the mean-field branching process!

2. Critical Slowing Down

Correlation time τ → ∞ at the critical point. The system takes longer to relax.

3. Diverging Fluctuations

The Fano factor (variance/mean) diverges. Fluctuations become much larger than the mean.

4. Finite-Size Scaling

Avalanche cutoff scales with N. Power laws extend further as network grows.

Compare All Three Regimes

Why Does Criticality Matter?
At criticality, the brain has:

Maximum dynamic range (sensitivity to inputs)
Maximum information storage (memory)
Maximum information transmission
Optimal computational capabilities

This is the “critical brain hypothesis” — the brain may operate near criticality!

Summary: The Complete Picture

Stage	Model	Equation	Key Insight
1	La Camera Binary	r = σ(Jr + I − θ)	Self-consistency, threshold
2	E/I Populations	r_E, r_I = σ(Wr + I)	E/I balance, no memory
3	Wilson-Cowan ODE	τ dr/dt = −r + σ(I)	Memory via τ, same fixed points
4	Stochastic WC	dS = drift + √[...]·η	Finite-size fluctuations ~1/√N
★	Criticality	w₀ᶜ = α/β = 0.1	Power laws, optimal computation

Wilson-Cowan Tutorial • PCTvsFEP Project

Based on La Camera (2022), Chow & Karimipanah (2020), de Candia et al. (2021)

Loading tutorial…

Our Journey

Binary Neurons

Micro → Macro

→

E/I Populations

Two coupled eqns

→

Continuous Time

Wilson-Cowan ODEs

→

Stochastic

Micro → Langevin

★

Criticality

Power laws

Introduction: Why This Tutorial?

You've learned the Mean Field Approximation from La Camera's paper. The key insight was simple but powerful:

La Camera's Core Move:

Replace actual spikes x_j with expected rates ⟨x_j⟩ = r_j

This gave us a self-consistent equation for the firing rate:

r = σ(J × r + I_ext)

Now we want to understand the stochastic Wilson-Cowan model used by de Candia et al. to study criticality. But their equations look different! This tutorial bridges the gap.

The Questions We'll Answer

How does the microscopic → macroscopic derivation work?
Why do we need separate E and I populations?
What changes when we go from discrete to continuous time?
Where does the noise term come from (finite-size effects)?
What is the critical point and why does it matter?
How does this relate to co-transmission (your ultimate goal)?

Binary Neurons — La Camera's Foundation

La Camera (2022), Section 2-4

The Model: What We Start With

We have a network of N binary neurons. Each neuron i is either spiking (x_i = 1) or silent (x_i = 0).

Variables Introduced

x_i	State of neuron i: 0 (silent) or 1 (spiking)
J_ij	Synaptic weight from neuron j to neuron i
I_ext	External input current (from outside the network)
I_i	Total input to neuron i: I_i = Σ_j J_ijx_j + I_ext

The Spiking Rule (Transfer Function)

Each neuron spikes with a probability that depends on its input via a sigmoid transfer function:

P(x_i = 1 | I_i) = σ(I_i) = 1 / (1 + e^{−2β(I_i − θ)})

New Variables

β	Inverse temperature or gain. Controls how “sharp” the decision is. High β → nearly deterministic. Low β → nearly random.
θ	Threshold. The input value where σ(I) = 0.5. Often set to 0 for simplicity, but in La Camera's formulation it appears inside the sigmoid.

The Mean Field Approximation

The key move: instead of tracking individual spikes, we track the expected firing rate:

r_i = ⟨x_i⟩ = P(x_i = 1) = σ(I_i)

New Variable

r_i	Firing rate of neuron i. This is what mean field theory tracks. It equals the probability that the neuron is spiking at any given moment.

For a homogeneous network (all neurons identical, uniform weights J/N):

Homogeneous Mean Field:

r = σ(J × r + I_ext − θ)

Self-Consistency: The rate r appears on both sides of the equation! We need to find the value where the output equals the input. This is a fixed point.

Try It: Find the Fixed Point

J (coupling):I_ext:β (gain):θ (threshold):

Stage 1 Summary

Binary neurons with probabilistic spiking: P(spike) = σ(input − θ)
The sigmoid has two key parameters: β (gain) and θ (threshold)
Mean field replaces spikes with rates: x_j → r_j
Self-consistent equation: r = σ(J × r + I_ext − θ)
This is a single homogeneous population

E/I Populations — The Wilson-Cowan Step

Chow & Karimipanah (2020)

Why Separate E and I?

La Camera's homogeneous network is a simplification. Real brains have two fundamentally different types of neurons:

Excitatory (~80%)

Release glutamate
Increase postsynaptic activity
Positive effect on target

Inhibitory (~20%)

Release GABA
Decrease postsynaptic activity
Negative effect on target

Dale's Law: A neuron is either excitatory OR inhibitory — it can't switch. This is why we need separate populations!

The Weight Matrix

New Variables: The Weight Matrix

w_EE	E → E: Excitatory neurons exciting each other (recurrent excitation)
w_EI	I → E: Inhibitory neurons inhibiting excitatory neurons
w_IE	E → I: Excitatory neurons driving inhibitory neurons
w_II	I → I: Inhibitory neurons inhibiting each other

Weight Matrix: W =

+w_EE−w_EI

+w_IE−w_II

Two Coupled Mean Field Equations

New Variables

r_E	Firing rate of the excitatory population (fraction of E neurons active)
r_I	Firing rate of the inhibitory population (fraction of I neurons active)

E/I Mean Field Equations:

r_E = σ(w_EE × r_E − w_EI × r_I + I_E)

r_I = σ(w_IE × r_E − w_II × r_I + I_I)

E/I Balance

The brain maintains a delicate balance between excitation and inhibition:

Too much E → Seizures, runaway activity

Too much I → Silence, no computation

Balanced → Rich dynamics, optimal processing

Try It: E/I Fixed Point

w_EE: w_EI: w_IE: w_II: I_E: I_I:

Stage 2 Summary

Separate E and I populations (Dale's Law)
Four connection weights: w_EE, w_EI, w_IE, w_II
Two coupled self-consistent equations
Firing rates depend on current state only (no memory/hysteresis yet)
E/I balance is crucial for brain function

Continuous Time — The Wilson-Cowan ODEs

Chow & Karimipanah (2020), Wilson & Cowan (1972)

From Discrete to Continuous Time

New Variables

τ_E, τ_I	Time constants (in ms). How fast each population responds. Typically τ_E ≈ 10ms, τ_I ≈ 20ms (I neurons are often slower).
dr/dt	Rate of change of firing rate. This is the derivative — how quickly r changes.

Wilson-Cowan Rate Equations:

τ_E × dr_E/dt = −r_E + σ(w_EEr_E − w_EIr_I + I_E)

τ_I × dr_I/dt = −r_I + σ(w_IEr_E − w_IIr_I + I_I)

Understanding the Terms

+σ(input) (drive term)Input pushes the firing rate toward σ(input). Combined with the decay term, the rate approaches σ(input) exponentially.

At Equilibrium (Fixed Point): When dr/dt = 0, we get:

0 = −r + σ(input) → r = σ(input)

Try It: Watch the Dynamics

w_EE=1.5w_EI=2w_IE=1.5w_II=0.5

τ_E (ms): τ_I (ms): I_E: I_I:

Stage 3 Summary

Continuous time via differential equations: τ × dr/dt = −r + σ(I)
The −r term is decay/relaxation, not “negative rate”
Time constants τ set the dynamics speed and introduce memory
At equilibrium, we get the same fixed points as Stage 2
The difference: how the system approaches equilibrium matters for transient behavior

Stochastic Wilson-Cowan — de Candia's Model

de Candia et al. (2021)

The Two-Level Structure (Just Like Stage 1!)

Part A: The Microscopic Model (Individual Neurons)

Just like La Camera's binary neurons, de Candia starts with N individual neurons. Each neuron i has a binary state:

Variables Introduced (Microscopic Level)

a_i	State of neuron i: 0 (quiescent) or 1 (active/firing)
s_i	Total synaptic input to neuron i: s_i = Σ_j w_ija_j + h
w_ij	Synaptic weight from neuron j to neuron i
h	Small external input (keeps system from getting stuck at all-silent)

de Candia Eq. 1 (Synaptic Input):

s_i = Σ_j w_ij × a_j + h

The Transition Rates (Stochastic Dynamics)

Unlike La Camera's instantaneous updates, de Candia uses a continuous-time Markov process. Neurons don't flip instantly — they transition at specific rates:

Active → Quiescent (1 → 0)Rate = α (constant). An active neuron spontaneously becomes quiescent. This represents spike termination and refractory effects.In time dt: P(deactivate) = α × dt

Activation Function (de Candia Eq. 2)

f(s)	Activation rate: f(s) = β·tanh(s) for s > 0, else 0 Unlike the sigmoid in Stage 1, this is zero for negative input — no spontaneous firing without positive drive!
α	Deactivation rate (1/τ). Typically α = 0.1 ms⁻¹ (τ = 10 ms)
β	Activation gain. Typically β = 1.0 ms⁻¹

The Microscopic Simulation (What the Code Does)

For each neuron i at each timestep dt:

if (neuron[i] is ACTIVE) {
    // Deactivate with probability α × dt
    if (random() < α × dt) {
        neuron[i] = QUIESCENT;
    }
} else {
    // QUIESCENT neuron: compute input
    s_i = Σⱼ wᵢⱼ × aⱼ + h;
    
    // Activate with probability f(s_i) × dt
    if (random() < f(s_i) × dt) {
        neuron[i] = ACTIVE;
    }
}

This is exactly what our SpikingWilsonCowan class does in the CriticalityProject!

Part B: The Mean Field Approximation

Now comes the key step — just like in Stage 1! We go from individual neurons to population averages. But unlike Stage 1, we keep the finite-size fluctuations.

Population Activity (Definition):

S = (1/N) × Σ_i a_i

S is the fraction of active neurons — exactly like r in La Camera's model! But how does S evolve? Let's count the transitions:

Rate of decrease (active → quiescent)There are N×S active neurons, each deactivating at rate α.
Expected decrease: α × S (per unit time)

Rate of increase (quiescent → active)There are N×(1−S) quiescent neurons, each activating at rate f(s).
Expected increase: (1−S) × f(s) (per unit time)

Deterministic Mean Field (Expected Value):

dS/dt = −αS + (1−S)f(s)

Part C: Where Does the Noise Come From?

Here's where de Candia differs from La Camera's Stage 1: instead of taking N → ∞ and discarding fluctuations, we keep N finite and derive the noise term.

Finite-Size Fluctuations

Consider N neurons, each firing or quiescent. In time dt:

Some quiescent neurons fire (Poisson process with rate f(s))
Some firing neurons stop (Poisson process with rate α)
Both are random processes — individual neurons flip stochastically
The actual number of transitions fluctuates around the expected value

The Central Limit Theorem Derivation

Let a_i be the state of neuron i (0 or 1). The population average is:

S = (1/N) × Σ a_i

Each a_i is a Bernoulli random variable. By the CLT, for large N:

Var(S) = Var((1/N) × Σ a_i) = (1/N²) × Σ Var(a_i) = (1/N²) × N × σ² = σ²/N

where σ² is the variance of a single neuron's state.

Therefore: Std(S) = σ/√N, so fluctuations scale as 1/√N

For N = 10,000 neurons: fluctuations ~ 1/100 = 0.01 of the mean
For N = 100 neurons: fluctuations ~ 1/10 = 0.1 of the mean (much larger!)

The Complete Isomorphism: Microscopic ↔ Macroscopic

Here's the precise correspondence between the two levels of description.They use the same physics — just described at different scales:

One-to-One Correspondence

Aspect	Microscopic (SpikingWilsonCowan)	Macroscopic (StochasticMeanFieldWC)
State variable	a_i ∈ {0, 1} for each neuron i	S ∈ [0, 1] = fraction of active neurons
Interpretation	Binary: active (1) or quiescent (0)	Continuous: population activity rate
Deactivation	P(1→0) = α × dt per neuron	Drift term: −αS
Activation	P(0→1) = f(s) × dt per neuron	Drift term: +(1−S)f(s)
f(s)	β×tanh(s) for s>0, else 0	Same: β×tanh(s) for s>0, else 0
Noise source	Stochastic transitions (Bernoulli)	Langevin: √[αS + f(s)(1−S)] × η(t)
Finite-size	Implicit (N discrete neurons)	Explicit: noise ∝ 1/√N

Complete Variable Reference

S	Fraction of active neurons = (1/N)×Σa_i. Same as firing rate r.
α	Deactivation rate (1→0). Typically α = 0.1 ms⁻¹ (so τ = 10 ms).
β	Activation gain in f(s). Typically β = 1.0 ms⁻¹
w₀	Net coupling strength (see criticality section for full definition).
h	Small external input (prevents absorbing state S=0).
N	Number of neurons. Noise scales as 1/√N.
η(t)	White noise: Gaussian, zero mean, unit variance.

Part D: The Full Langevin Equation (de Candia Eq. 4)

Combining the mean field drift with the finite-size noise gives us the complete Langevin equation:

What is a Langevin Equation?

A Langevin equation is a stochastic differential equation (SDE) of the form:

dX/dt = drift(X) + noise amplitude(X) × η(t)

It combines:

Drift: The deterministic part — where the system would go without randomness
Noise: Random fluctuations that perturb the system. η(t) is “white noise” (Gaussian, uncorrelated in time)

Langevin equations are fundamental in physics (Brownian motion) and are the continuous-time limit of discrete random walks. They're essential for studying systems near criticality!

de Candia Eq. 4 (Full Langevin Form):

dS/dt = −αS + (1−S)f(s) + (1/√N)·√[αS + f(s)(1−S)]·η(t)

−αS (drift: decay)Firing neurons become quiescent at rate α. This pulls S toward 0.

(1−S)f(s) (drift: activation)Quiescent neurons (fraction 1−S) fire at rate f(s). This pushes S up.

Activation Function

f(s)	Activation rate: f(s) = β·tanh(s) for s > 0, else 0 where s = w₀·S + h (total synaptic input). Unlike the smooth sigmoid, this is zero for negative input, which creates an absorbing state at S=0.

Try It: Stochastic Simulation

α = 0.1 ms⁻¹β = 1.0 ms⁻¹h = 10⁻⁶Critical w₀ᶜ = α/β = 0.1

w₀:0.10N: Duration (ms):

Stage 4 Summary: The Two-Level Structure

Microscopic Model: N neurons with states a_i ∈ {0,1}, transitioning at rates α and f(s)
Macroscopic Model: Population activity S, evolving via Langevin equation
The derivation: Mean field (S = ⟨a⟩) + Central Limit Theorem (noise ∝ 1/√N)
Key isomorphism: P(deactivate) = α×dt → drift term −αS
Key isomorphism: P(activate) = f(s)×dt → drift term (1−S)f(s)
Why both? Microscopic is the “ground truth”; macroscopic is the approximation
Same physics, different levels of description — just like La Camera's Stage 1!

★

Criticality — The Main Event

de Candia et al. (2021), Section “Critical point”

Understanding the Control Parameter w₀

Clarifying w₀: In the full E/I model, we have separate weights w_E (excitatory) and w_I (inhibitory). The effective coupling is their difference:

w₀ = w_E − w_I

In de Candia's simplified model, we work directly with this single parameter w₀. Positive w₀ means excitation dominates; negative w₀ means inhibition dominates.

The Critical Point

At the critical point, the system is perfectly balanced between activity dying out (subcritical) and runaway activity (supercritical).

Critical Condition:

w₀^c = α/β = 0.1 ms⁻¹ / 1 ms⁻¹ = 0.1

The Three Regimes

Subcritical (w₀ < 0.1)

Inhibition dominates

Activity decays to zero
Small, short-lived fluctuations
No sustained dynamics

Critical (w₀ = 0.1)

Perfect E/I balance

Power-law avalanches
Diverging correlation time
Maximal sensitivity

Supercritical (w₀ > 0.1)

Excitation dominates

Runaway activity
Saturation at high S
Seizure-like dynamics

Signatures of Criticality

1. Power-Law Avalanches

Size distribution: P(S) ~ S^−3/2
Duration distribution: P(T) ~ T⁻²

These are the exponents of the mean-field branching process!

2. Critical Slowing Down

Correlation time τ → ∞ at the critical point. The system takes longer to relax.

3. Diverging Fluctuations

The Fano factor (variance/mean) diverges. Fluctuations become much larger than the mean.

4. Finite-Size Scaling

Avalanche cutoff scales with N. Power laws extend further as network grows.

Compare All Three Regimes

Why Does Criticality Matter?
At criticality, the brain has:

Maximum dynamic range (sensitivity to inputs)
Maximum information storage (memory)
Maximum information transmission
Optimal computational capabilities

This is the “critical brain hypothesis” — the brain may operate near criticality!

Summary: The Complete Picture

Stage	Model	Equation	Key Insight
1	La Camera Binary	r = σ(Jr + I − θ)	Self-consistency, threshold
2	E/I Populations	r_E, r_I = σ(Wr + I)	E/I balance, no memory
3	Wilson-Cowan ODE	τ dr/dt = −r + σ(I)	Memory via τ, same fixed points
4	Stochastic WC	dS = drift + √[...]·η	Finite-size fluctuations ~1/√N
★	Criticality	w₀ᶜ = α/β = 0.1	Power laws, optimal computation

Wilson-Cowan Tutorial • PCTvsFEP Project

Based on La Camera (2022), Chow & Karimipanah (2020), de Candia et al. (2021)