Mathematical Methods for Neuroscience

Contents

• Introduction and overview
• 0. Review: differential equations, Fourier analysis and probability
• 1. Neuron models
• 2. Dynamical systems
• 3. Bifurcation theory
• 4. The role of noise
• 5. The visual cortex and neural field models
• 6. Symmetry in V1: orientation tuning, cortical patterns and hallucinations
• 7. Mean-field theory I: random networks of rate neurons
• 8. Mean-field theory II: interacting particles, spiking networks and neural masses
• 9. Research directions and open problems
• Exercises
• Projects and assessment
• References

Introduction and overview

Audience: a student with a good M1 in mathematics, physics or computer science who wants to see how dynamical systems, bifurcation theory and stochastic processes become quantitative models of neurons and of the cortex.

Empirical first: inject a constant current into the Hodgkin–Huxley model of the squid giant axon. Up to about 6 µA/cm² the membrane stays at rest; just above, it fires repetitively at about 55 Hz, and never at 5 Hz. Why can the firing rate not start from zero, why do many cortical neurons do the opposite, and why do the geometric hallucinations reported under mescaline so often look like tunnels, funnels, spirals and lattices?

Each of these questions has a precise mathematical answer, and this course builds the tools that give it. Single neurons are low-dimensional dynamical systems whose behaviour changes at bifurcations; noise turns them into diffusion processes with first-passage times; large networks are described by neural fields and by mean-field limits, in which symmetry dictates which spatial patterns appear. The course is taught by Romain Veltz and Quentin Cormier (INRIA, per master-mva.com). It shares its stochastic tools with Stochastic Calculus in Machine Learning and its random-matrix tools with Large Random Matrices and Applications to Machine Learning.

What you will be able to do

• Write, simulate and analyse the space-clamped Hodgkin–Huxley model and its reductions (FitzHugh–Nagumo, integrate-and-fire, quadratic integrate-and-fire).
• Draw the phase portrait of a planar system, find its equilibria, invariant manifolds and limit cycles, and decide their stability.
• State and use Hartman–Grobman, the topological classification of hyperbolic equilibria, structural stability and the centre manifold theorem.
• Recognise saddle-node, transcritical, pitchfork and Hopf bifurcations, reduce a system to its normal form and compute the first Lyapunov coefficient.
• Explain the two classical excitability classes of neurons through the bifurcation that creates firing.
• Solve linear stochastic differential equations, write the Fokker–Planck equation and compute mean first-passage times and the firing rate of a noisy integrate-and-fire neuron.
• Describe the anatomy and functional architecture of the primary visual cortex and model it with neural field equations.
• Use symmetry (O(2) and the Euclidean group) to predict orientation tuning, cortical patterns and the geometry of visual hallucinations.
• Derive the dynamic mean-field equations of random rate networks and the transition to chaos.
• Prove propagation of chaos for interacting diffusions by coupling, and derive the mean-field limits of spiking networks and of neural masses.

Official objectives and syllabus (2026–2027)

The course description on master-mva.com states the objective as presenting some mathematical tools that intervene systematically in many modelling problems in neuroscience. The description is a list of themes, not a session-by-session programme; the Université Paris-Saclay page of the course gives the volume (24 h of lectures, 16 h of tutorials) and the assessment type, but no programme. The table below takes the themes item by item, in the official order.

Every theme of the official list is taught on this page. Beyond the syllabus, the page adds the review chapter (ch. 0), the excitability classes of neurons (ch. 3.5) and a self-study chapter on research directions (ch. 9).

Where this course fits in the programme

Official prerequisites: mathematics at the level of L3 and M1 courses. In practice, the course uses ordinary differential equations, linear algebra, Fourier series and measure-theoretic probability.

• Preparation — the refreshers Linear Algebra, Analysis and Optimization (eigenvalues, matrix exponential, implicit function theorem), Probability and Statistics (Gaussian vectors, conditional expectation, limit theorems) and Scientific Python, NumPy and PyTorch (the simulations of this page). Chapter 0 reviews what is used.
• Related courses in the same semester — Stochastic Calculus in Machine Learning develops Itô calculus and the links between diffusions and partial differential equations in depth; Optimal Transport for Machine Learning studies the Wasserstein distances behind propagation of chaos; Theoretical Foundations of Deep Learning studies the non-convex optimisation of neural networks, to which the mean-field tools of ch. 8 apply (ch. 9.5); Medical Image Analysis Based on Generative, Geometric and Biophysical Models is another Santé-track course of the same semester.
• Courses that build on it or share its tools (semester 2) — Functional Brain Imaging and Brain–Computer Interfaces (the signals that neural masses model), Large Random Matrices and Applications to Machine Learning (the random-matrix tools behind ch. 7), Interactions: Complexity, Stochasticity, Extrema and Rare Events (interacting particle systems, chaos and rare events) and Stopping Times and Online Algorithms (stopping times, of which the first-passage times of ch. 4 are an example).

Course outline

• Review: differential equations, Fourier analysis and probability — existence and uniqueness, Grönwall, linear systems, convolution and Fourier modes, Gaussian tools.
• Neuron models — the membrane circuit, the Hodgkin–Huxley equations, their reductions, synapses and the cable equation.
• Dynamical systems — phase portraits, invariant manifolds, equivalence, Hartman–Grobman, structural stability and the centre manifold.
• Bifurcation theory — saddle-node, transcritical, pitchfork and Hopf bifurcations, normal forms, excitability classes and equivariant bifurcations.
• The role of noise — Brownian motion, Itô calculus, Fokker–Planck equations and the firing rate of noisy neurons.
• The visual cortex and neural field models — the anatomy and functional architecture of V1, the Amari equation, bumps, Turing patterns and fronts.
• Symmetry in V1 — orientation tuning in the ring model, cortical patterns on lattices and the geometry of visual hallucinations.
• Mean-field theory I — random rate networks, the circular law, dynamic mean-field theory and the transition to chaos, the Ben Arous–Guionnet theorem.
• Mean-field theory II — McKean–Vlasov processes, propagation of chaos, mean-field limits of spiking networks and neural masses.
• Research directions and open problems — mean-field models with spikes, numerical bifurcation analysis, structured random networks (self-study).

Session plan

40 h in class = 8 sessions of 5 h (3 h lecture + 2 h tutorial). The volume is the Université Paris-Saclay programme table (24 h of lectures, 16 h of tutorials). Master-mva.com organises the course on Thursdays at ENS Paris-Saclay, with a lecture from 13:30 to 16:30 followed by a tutorial from 16:30 to 18:30, which gives exactly eight such afternoons. The split of the content into sessions is a proposal.

The order differs from the official list in one respect: neuron models (official item 5) come first, because they supply the examples of dynamical systems and bifurcations, and the visual cortex (official items 1 and 4) comes after the bifurcation theory it uses. Level 1 = foundations, 2 = core results, 3 = advanced and research-level material; it never decreases.

Self-study, outside the 40 h (about 4 h in total): ch. 9 Research directions and open problems (≈ 4 h, to be read with the paper-study projects).

How to study this module

• Read ch. 0 before the first session, even quickly: Grönwall’s lemma and the Fourier modes of convolutions return in almost every chapter.
• Run every code block. Each one checks a formula of the text against a simulation (the Hodgkin–Huxley f–I curve, the Siegert formula, the Turing wavelength, the width of the ring-model bump, the variance predicted by dynamic mean-field theory).
• After each session, do the exercises listed in the session plan before reading on: ch. 3 relies on ch. 2, ch. 6 on ch. 3 and ch. 5, ch. 8 on ch. 4.
• Keep a notebook of the normal forms of ch. 3: the same five equations organise ch. 1, ch. 5 and ch. 6.
• The Python environment is NumPy and SciPy only; every simulation runs in a few seconds on a laptop.
• The resit is a paper reading: the paper-study projects below are good training for it.

0. Review: differential equations, Fourier analysis and probability

Neuron and network models are differential equations: sometimes in finite dimension, sometimes on a space of functions, sometimes driven by noise. This chapter collects the M1-level results that the later chapters use without stopping: well-posedness, Grönwall’s lemma, linear stability, the Fourier modes of convolution operators and a few Gaussian identities.

0.1 Existence, uniqueness and Grönwall’s lemma

0.2 Linear systems and linearised stability

0.3 Fourier modes of convolution operators

The neural fields of ch. 5 and the ring model of ch. 6 are built on convolution operators, and their linear stability always reduces to the following computation.

0.4 Probability toolkit

Chapters 4, 7 and 8 use Gaussian vectors, the law of large numbers and the central limit theorem, and the Poisson process. Two identities deserve a statement because they are the engine of the mean-field computations.

1. Neuron models

A neuron is a cell whose membrane potential can make a brief, stereotyped excursion, the action potential or spike, and pass it on to thousands of other cells. This chapter writes the models that reproduce this behaviour: the four equations of Hodgkin and Huxley, their two- and one-dimensional reductions, the synapses that connect neurons, and the cable equation that adds space.

1.1 The membrane as an electrical circuit

The lipid bilayer separating the inside of the cell from the outside behaves as a capacitor, with a specific capacitance of about 1 µF/cm². Ion pumps keep potassium concentrated inside the cell and sodium outside.

Ion channels are pores selective to one ion species. Each population of open channels acts as a conductance in series with a battery, whose voltage is the reversal potential of that ion.

At body temperature RT/F \approx 26.7 mV. With 140 mM of potassium inside and 5 mM outside, E_K \approx -89 mV; sodium, concentrated outside, has a positive reversal potential (Exercise 1).

1.2 The space-clamped Hodgkin–Huxley model

In 1952 Hodgkin and Huxley inserted a wire along the squid giant axon, which makes the potential uniform along the axon (the space clamp), and imposed V by feedback (the voltage clamp). From the measured sodium and potassium currents they built the following model.

The gating variables read as fractions of open gates. The variable m activates the sodium current within a fraction of a millisecond; h inactivates it, and n activates the potassium current, both over several milliseconds.

The spike is a fast positive feedback followed by a slow negative one. A depolarisation opens sodium gates, the inward sodium current depolarises further, and V shoots towards E_{Na}.

Then h closes the sodium channels and n opens the potassium channels, which bring V back below rest. While h and n recover, a second spike is hard to trigger: this is the refractory period.

import numpy as np
from scipy.integrate import solve_ivp

C, gNa, gK, gL = 1.0, 120.0, 36.0, 0.3          # uF/cm^2, mS/cm^2
ENa, EK, EL = 50.0, -77.0, -54.387              # mV

def rates(V):
    an = 0.01 * (V + 55) / (1 - np.exp(-(V + 55) / 10))
    bn = 0.125 * np.exp(-(V + 65) / 80)
    am = 0.1 * (V + 40) / (1 - np.exp(-(V + 40) / 10))
    bm = 4.0 * np.exp(-(V + 65) / 18)
    ah = 0.07 * np.exp(-(V + 65) / 20)
    bh = 1.0 / (1 + np.exp(-(V + 35) / 10))
    return an, bn, am, bm, ah, bh

def hh(t, y, I):
    V, n, m, h = y
    an, bn, am, bm, ah, bh = rates(V)
    dV = (I - gNa * m**3 * h * (V - ENa) - gK * n**4 * (V - EK) - gL * (V - EL)) / C
    return [dV, an * (1 - n) - bn * n, am * (1 - m) - bm * m, ah * (1 - h) - bh * h]

def rate_hz(I, T=300.0, y0=(-65.0, 0.317, 0.053, 0.596)):
    sol = solve_ivp(hh, (0, T), y0, args=(I,), max_step=0.1, rtol=1e-6)
    V, t = sol.y[0], sol.t
    up = np.where((V[:-1] < 0) & (V[1:] >= 0))[0]       # upward crossings of 0 mV
    late = t[up][t[up] > 100.0]                         # discard the transient
    return 0.0 if len(late) < 2 else 1000.0 / np.mean(np.diff(late))

for I in [0, 5, 6, 6.5, 7, 10, 20, 40]:
    print(f"I = {I:5.1f} uA/cm^2  ->  {rate_hz(I):6.1f} Hz")

1.3 Simplified models

Four variables are too many for a phase plane. Two observations reduce the model to two dimensions (Rinzel’s reduction).

• The sodium activation is much faster than the rest, so one may set m = m_\infty(V).
• Along trajectories, h and n are strongly anticorrelated. In the simulation of Example 1.4 at 10 µA/cm², n + h stays between 0.78 and 0.85, and a linear fit, on samples equally spaced in time, gives h \approx 0.87 - 1.07\, n.

The result is a planar system in (V, n) with a cubic-shaped V-nullcline. FitzHugh and, independently, Nagumo and coauthors kept only this geometry.

The most drastic simplification keeps only the subthreshold dynamics and replaces the spike by a rule. This is the integrate-and-fire neuron, which goes back to Lapicque (1907).

The quadratic integrate-and-fire neuron keeps the spike initiation as a smooth instability rather than a hard threshold. It is the normal form of the saddle-node bifurcation of ch. 3, and the building block of the neural mass of ch. 8.

Izhikevich’s model (2003) adds a slow adaptation variable to a quadratic neuron. With four parameters it reproduces most firing patterns recorded in cortex (regular spiking, bursting, fast spiking), at the cost of one reset rule.

1.4 Synapse models

At a chemical synapse, a presynaptic spike releases neurotransmitter, which opens receptor channels on the postsynaptic membrane. The postsynaptic neuron sees a transient conductance.

The exponential synapse is a linear filter. If the presynaptic train is Poisson with rate \nu, Campbell’s theorem gives the mean and variance of s in the stationary regime: \mathbb{E}[s] = \nu \tau_s, \qquad \mathrm{Var}[s] = \frac{\nu \tau_s}{2}.

These two numbers are the drift and the noise of the diffusion approximation of ch. 4.5.

1.5 Spatial models: the cable equation

Dendrites and axons are thin cylinders, and the potential varies along them. Charge conservation in a cylinder of radius a, with axial resistivity r_L and membrane capacitance c_m per unit area, gives the cable equation c_m \frac{\partial V}{\partial t} = \frac{a}{2 r_L} \frac{\partial^2 V}{\partial x^2} - i_{\text{ion}} + i_{\text{ext}}.

For a passive membrane, i_{\text{ion}} = (V - E_L)/r_m, where r_m is the specific membrane resistance. With u = V - E_L it reads \tau_m \frac{\partial u}{\partial t} = \lambda^2 \frac{\partial^2 u}{\partial x^2} - u + r_m\, i_{\text{ext}}, \qquad \tau_m = r_m c_m, \quad \lambda = \sqrt{\frac{a\, r_m}{2 r_L}} .

Hodgkin and Huxley solved the active cable equation numerically and obtained a conduction velocity close to the measured one. In practice, a dendritic tree is cut into short compartments coupled by axial conductances, each one a copy of Definition 1.1: these compartmental models are what simulators such as NEURON integrate.

2. Dynamical systems

The models of ch. 1 are too nonlinear to solve in closed form, but their qualitative behaviour can be read from the geometry of their solutions. This chapter gives the vocabulary of that geometry: orbits and phase portraits, invariant manifolds, the equivalence of two systems, structural stability and the centre manifold, the tool that reduces a high-dimensional system to the few directions where something happens.

2.1 Flows, orbits and phase portraits

In the plane, the topology of the plane severely limits what orbits can do.

import numpy as np
from scipy.optimize import brentq

a, b, eps = 0.7, 0.8, 0.08                        # FitzHugh's classical values

def equilibrium(I):
    # v - v^3/3 - (v + a)/b + I = 0 has a single root since b < 1
    v = brentq(lambda v: v - v**3 / 3 - (v + a) / b + I, -3, 3)
    return v, (v + a) / b

def jacobian(v):
    return np.array([[1 - v**2, -1.0], [eps, -eps * b]])

for I in [0.0, 0.3, 0.33, 0.34, 0.5, 1.0, 1.4, 1.5]:
    v, w = equilibrium(I)
    J = jacobian(v)
    tr, det = np.trace(J), np.linalg.det(J)
    kind = ("focus" if tr**2 < 4 * det else "node") if det > 0 else "saddle"
    print(f"I={I:4.2f}  v*={v:+.4f}  tr={tr:+.4f}  det={det:.4f}  "
          f"{'unstable' if tr > 0 else 'stable'} {kind}")

# Hopf points: tr J = 0  <=>  v*^2 = 1 - eps*b
vH = np.sqrt(1 - eps * b)
for v in (-vH, vH):
    print(f"Hopf at v*={v:+.4f}, I_H={(v + a) / b - v + v**3 / 3:.4f}, "
          f"omega={np.sqrt(np.linalg.det(jacobian(v))):.4f}")

2.2 Invariant manifolds

An equilibrium x^* is hyperbolic if no eigenvalue of A = Df(x^*) lies on the imaginary axis. The generalised eigenspaces for eigenvalues of negative and positive real part are the stable and unstable subspaces E^s and E^u.

Invariant manifolds are the thresholds of neuroscience. In a model where rest coexists with a saddle (the Morris–Lecar model in its first-class regime, ch. 3.5), the stable manifold of the saddle separates the initial conditions that return to rest from those that produce a spike: it is a true, sharp threshold.

FitzHugh–Nagumo has no saddle at small I, and its threshold is only a narrow region of rapid change. FitzHugh called it a quasi-threshold: responses grow steeply but continuously with the stimulus.

2.3 Equivalence and the topological classification of equilibria

When are two dynamical systems “the same”? Requiring identical solutions is too strict; the right notion keeps the orbits and forgets the parametrisation.

Smooth conjugacy is too fine for classification: it preserves the eigenvalues of the linearisation, which vary continuously. Topological equivalence forgets them and keeps only signs.

In the plane there are therefore only three topological types of hyperbolic equilibria: sink, source and saddle. A stable node and a stable focus are topologically equivalent; the spiralling, visible as damped oscillations of the membrane potential, is a smooth invariant, not a topological one.

2.4 Structural stability

A model is useful only if its qualitative predictions survive small errors in its equations. Structural stability makes this precise.

The centre \dot x = -y, \dot y = x is not structurally stable. The perturbation \dot x = \varepsilon x - y, \dot y = x + \varepsilon y turns it into a focus, stable or unstable according to the sign of \varepsilon.

In a generic one-parameter family of planar models, structurally unstable systems occur at isolated parameter values. These are the bifurcation points of ch. 3, where the conclusions of a model change.

2.5 The centre manifold in finite dimension

At a non-hyperbolic equilibrium the linearisation does not decide stability, and something happens on the eigenvalues of zero real part. The centre manifold isolates these directions.

To compute h, write that the graph is invariant. With E^h = E^s \oplus E^u and the projections P_c, P_h: Dh(u)\,\big[A_c u + P_c F(u + h(u))\big] = A_h\, h(u) + P_h F\big(u + h(u)\big).

This equation is solved order by order in the Taylor expansion of h.

3. Bifurcation theory

When a parameter of a model (an injected current, a coupling strength) crosses a critical value, the phase portrait can change: equilibria appear or disappear, lose stability, give birth to oscillations. This chapter classifies the generic ways this happens, reduces each to a normal form, and uses them to explain the excitability of neurons and, through symmetry, the patterns of ch. 6.

3.1 What is a bifurcation?

A local bifurcation of an equilibrium requires a non-hyperbolic equilibrium. Indeed, if D_x f(x_0, \mu_0) is invertible, the implicit function theorem continues the equilibrium as a smooth curve x^*(\mu). Its eigenvalues move continuously, so the number of stable ones does not change, and Theorem 2.8 gives local equivalence.

For equilibria, the two codimension-one cases are therefore a zero eigenvalue (saddle-node) and a pair of purely imaginary eigenvalues (Hopf). Symmetry (ch. 3.6) or a trivial solution that persists for all parameters changes the generic picture.

3.2 One-dimensional bifurcations

3.3 The Hopf bifurcation

For a planar system already written with linear part \begin{pmatrix} 0 & -\omega \\ \omega & 0 \end{pmatrix} and nonlinear terms f, g, the classical formula of Guckenheimer and Holmes (1983, §3.4) gives \mathrm{Re}\,c_1 = a with 16 a = f_{xxx} + f_{xyy} + g_{xxy} + g_{yyy} + \frac{1}{\omega}\big[f_{xy}(f_{xx} + f_{yy}) - g_{xy}(g_{xx} + g_{yy}) - f_{xx} g_{xx} + f_{yy} g_{yy}\big].

Kuznetsov (2004, §3.5) gives a coordinate-free formula with the right and left eigenvectors of the Jacobian. Written with the bilinear and trilinear forms of the vector field, it holds in any dimension (same book, ch. 5). It is the one used below.

import numpy as np
a, b, eps = 0.7, 0.8, 0.08
vH = -np.sqrt(1 - eps * b)                        # lower Hopf point
A = np.array([[1 - vH**2, -1.0], [eps, -eps * b]])
B = lambda x, y: np.array([-2 * vH * x[0] * y[0], 0.0])   # quadratic part at the equilibrium
C = lambda x, y, z: np.array([-2 * x[0] * y[0] * z[0], 0.0])

lam, Q = np.linalg.eig(A)
k = np.argmax(lam.imag); om = lam[k].imag; q = Q[:, k]
lamT, P = np.linalg.eig(A.T)
p = P[:, np.argmin(np.abs(lamT + 1j * om))]       # A^T p = -i om p
p = p / np.conj(np.vdot(p, q))                    # normalise <p, q> = 1 (<x,y> = conj(x).y)
dot = lambda u, v: np.vdot(u, v)
l1 = np.real(dot(p, C(q, q, q.conj()))
             - 2 * dot(p, B(q, np.linalg.solve(A, B(q, q.conj()))))
             + dot(p, B(q.conj(), np.linalg.solve(2j * om * np.eye(2) - A, B(q, q))))) / (2 * om)
print(f"omega = {om:.4f}, first Lyapunov coefficient l1 = {l1:.4f}")

3.4 Centre manifold reduction and normal forms

A bifurcation of an equilibrium in \mathbb{R}^n involves only the eigenvalues on the imaginary axis. The parameter-dependent centre manifold makes this exact.

Append \dot\mu = 0 to the system. The parameters become centre directions, and Theorem 2.11 gives a manifold of dimension \dim E^c + p that contains every small bounded solution for every small \mu.

For the Hopf bifurcation of the four-dimensional Hodgkin–Huxley model, the reduced system is two-dimensional, and Theorem 3.5 applies to it. For a fold, it is one-dimensional, and Theorem 3.2 applies.

On the centre manifold, polynomial changes of variables simplify the equations further.

For a Hopf point, use the coordinates (z, \bar z), with eigenvalues i\omega and -i\omega. A monomial z^p \bar z^q in the \dot z equation is resonant when i\omega = (p - q)\, i\omega, \quad \text{that is} \quad p = q + 1 .

The resonant terms are z|z|^2, z|z|^4, \dots: this is why the normal form of Theorem 3.5 keeps only c_1 |z|^2 z at third order.

At a fold, the single eigenvalue 0 makes every monomial resonant. The quadratic term survives, and the normal form is \dot x = \beta + x^2, after scaling by f_{xx} \ne 0.

3.5 Excitability classes of neurons

In 1948 Hodgkin sorted the axons of the crab Carcinus by their response to a constant current. Rinzel and Ermentrout (1989) explained his classes by the bifurcation through which the resting state gives way to repetitive firing.

The quadratic integrate-and-fire neuron (Proposition 1.6) is exactly the SNIC normal form, with f = \sqrt{I}/(\pi\tau). Hodgkin–Huxley and FitzHugh–Nagumo with standard parameters are class 2 (Examples 1.4 and 3.6). The Morris–Lecar model (1981) is a planar conductance-based model, with an instantaneous calcium current and a slower potassium current. It is class 1 or class 2 depending on its parameters, and many cortical pyramidal neurons behave as class 1.

A second distinction follows from the linearisation at rest. Close to a Hopf bifurcation the resting state is a stable focus: it produces damped subthreshold oscillations and responds best to inputs at its resonant frequency. Close to a saddle-node the resting state is a stable node, which simply integrates its inputs.

Izhikevich (2007) calls these neurons resonators and integrators. The difference is smooth, not topological (Theorem 2.8), yet it matters for timing.

3.6 Equivariant bifurcations

Neural fields on a ring or on the plane (ch. 5–6) are invariant under rotations, translations and reflections. Symmetry forces multiple eigenvalues, which are non-generic otherwise, and organises the solutions that bifurcate.

The case used in ch. 6 is \mathrm{O}(2) acting on \mathbb{C} \simeq \mathbb{R}^2 by rotations z \mapsto e^{i\phi} z and the reflection z \mapsto \bar z.

4. The role of noise

Recorded in vivo, cortical neurons fire irregularly: the coefficient of variation of their interspike intervals is often close to 1, the value of a Poisson process (Softky and Koch, 1993). Channels open at random, synapses release transmitter at random, and each neuron is bombarded by thousands of inputs. This chapter builds the stochastic calculus needed to model this variability and computes the firing rate of a noisy neuron in closed form.

4.1 Brownian motion

A continuous path with non-zero quadratic variation has infinite total variation on every interval. Therefore \int H\, dB cannot be defined path by path as a Stieltjes integral; Itô’s construction uses the L^2 structure instead.

The heuristic rule (dB_t)^2 = dt summarises Proposition 4.2.

4.2 The Itô integral and stochastic differential equations

Simulations use the Euler–Maruyama scheme, which replaces dB by a Gaussian increment of variance \Delta t: X_{k+1} = X_k + b(X_k)\, \Delta t + \sigma(X_k)\, \sqrt{\Delta t}\; \xi_k, \qquad \xi_k \sim \mathcal{N}(0, 1) \text{ i.i.d.}

Its strong error is of order \sqrt{\Delta t} in general (Kloeden and Platen, 1992).

4.3 Itô’s formula and the Ornstein–Uhlenbeck process

The extra term \frac12 \sigma^2 \partial_{xx} F is the whole difference with ordinary calculus. For example, d(B_t^2) = 2B_t\, dB_t + dt, so B_t^2 - t is a martingale.

4.4 The Fokker–Planck equation and first-passage times

The operator L\varphi = b\varphi' + \frac12\sigma^2\varphi'' is the generator of X. In the stationary regime with zero flux, the Fokker–Planck equation integrates once and gives p_\infty(x) \propto \frac{1}{\sigma^2(x)}\, \exp\Big(\int^x \frac{2 b(y)}{\sigma^2(y)}\, dy\Big).

For the Ornstein–Uhlenbeck process this is the Gaussian of Example 4.6.

A neuron fires when its potential first reaches the threshold. The mean time to get there solves an ordinary differential equation.

4.5 Noisy integrate-and-fire neurons

Suppose a leaky integrate-and-fire neuron receives C independent Poisson inputs of rate \nu, each making the potential jump by J. Then dV = -V\,dt/\tau + J\, dN_t, where N is Poisson with rate C\nu.

If J is small and C\nu large, the input is approximated by its mean and its Gaussian fluctuations (the diffusion approximation). The result is \tau\, dV = (\mu - V)\, dt + \sigma\sqrt{\tau}\, dB below threshold, with \mu = C J \nu \tau, \qquad \sigma^2 = C J^2 \nu \tau .

import numpy as np
from scipy.integrate import quad
from scipy.special import erfcx

tau, theta, Vr, tref = 20.0, 20.0, 10.0, 2.0      # ms, mV (values of Brunel, 2000)
mu, sigma = 18.0, 5.0                              # mean input and noise amplitude (mV)

# Siegert formula; e^{u^2} (1 + erf u) = erfcx(-u) avoids overflow
integral = quad(lambda u: erfcx(-u), (Vr - mu) / sigma, (theta - mu) / sigma)[0]
nu_theory = 1000.0 / (tref + tau * np.sqrt(np.pi) * integral)

# Euler-Maruyama for tau dV = (mu - V) dt + sigma sqrt(tau) dB, N neurons in parallel
rng = np.random.default_rng(0)
N, dt, T, T0 = 1000, 0.02, 600.0, 100.0           # T0: discarded transient (ms)
V = np.full(N, Vr); ref = np.zeros(N); spikes = 0
for k in range(int(T / dt)):
    active = ref <= 0
    V[active] += (mu - V[active]) * dt / tau \
                 + sigma * np.sqrt(dt / tau) * rng.standard_normal(active.sum())
    ref -= dt
    fired = V >= theta
    if k * dt > T0:
        spikes += fired.sum()
    V[fired] = Vr; ref[fired] = tref
nu_sim = 1000 * spikes / (N * (T - T0))
print(f"Siegert: {nu_theory:.2f} Hz   Euler-Maruyama: {nu_sim:.2f} Hz")

The density of the potential of a noisy integrate-and-fire neuron solves the Fokker–Planck equation of Proposition 4.7 on (-\infty, \theta), with an absorbing boundary at \theta. The outgoing flux at \theta is the firing rate, and it is reinjected at V_r.

When the input rate \nu is itself the output rate of the network, this becomes a nonlinear equation: the starting point of the mean-field theory of spiking networks (Brunel and Hakim, 1999; ch. 8).

Noise can also create order. In an excitable system driven by noise, the regularity of the noise-induced spikes is maximal at an intermediate noise level: this is coherence resonance (Pikovsky and Kurths, 1997).

5. The visual cortex and neural field models

The primary visual cortex contains more than a hundred million neurons in each human hemisphere, far too many to model one by one. At the scale of a millimetre, however, its organisation is regular enough to be described by continuous fields of activity. This chapter describes that organisation, then introduces the neural field equations, their localised states (bumps), their spatial instabilities (Turing patterns) and their travelling fronts.

5.1 Anatomy of the primary visual cortex

Light is transduced in the retina, whose ganglion cells send their axons through the optic nerve to the lateral geniculate nucleus (LGN) of the thalamus. The LGN projects to the primary visual cortex V1, also called striate cortex or Brodmann area 17, at the back of the occipital lobe.

V1 is a folded sheet about 2 mm thick, organised in six layers:

• layer 4 receives most of the thalamic input;
• layers 2/3 project to other cortical areas and carry long-range horizontal connections within V1;
• layers 5 and 6 project to subcortical targets, layer 6 back to the LGN.

About four neurons in five are excitatory (pyramidal and spiny stellate cells), the others inhibitory interneurons.

The map from the visual field to V1 is retinotopic: neighbouring points of the visual field excite neighbouring points of the cortex. Each hemisphere represents the opposite half of the visual field, and the centre of gaze, the fovea, occupies a disproportionate area.

In human V1, Horton and Hoyt (1991) estimated the linear magnification at eccentricity E (degrees) as about 17.3/(E + 0.75) mm per degree. The following model captures this geometry.

5.2 Functional architecture of V1

Hubel and Wiesel (1959, 1962) found that most V1 neurons respond to bars or edges of a preferred orientation in a small region of the visual field, their receptive field. Neurons met along a vertical electrode penetration share the same preferred orientation: they form an orientation column. Along a tangential penetration, the preferred orientation rotates smoothly.

Optical imaging (Bonhoeffer and Grinvald, 1991) revealed the whole orientation preference map. It is smooth except at isolated points, the pinwheels, around which every orientation is represented exactly once.

A hypercolumn is a patch of about 1 mm² that contains all orientations, and both eyes, for one location in the visual field. Kaschube and coauthors (2010) measured a pinwheel density close to \pi per hypercolumn area in three distantly related species.

The horizontal connections of layers 2/3 extend over several millimetres. They link preferentially columns with similar orientation preference, and in the tree shrew they run preferentially along the axis, in the visual field, of the preferred orientation (Bosking and coauthors, 1997).

5.3 From populations to neural fields

Wilson and Cowan (1972) described a local population by the fractions E and I of active excitatory and inhibitory cells. In a simplified form, their equations read \tau_E \dot E = -E + S_E(w_{EE} E - w_{EI} I + h_E), \qquad \tau_I \dot I = -I + S_I(w_{IE} E - w_{II} I + h_I).

In 1973 they extended the model to a continuum of populations along the cortex. Amari (1977) analysed the one-population version that is now standard.

5.4 Bumps: localised persistent activity

Take a one-dimensional field with a Heaviside rate function S(u) = H(u - \kappa), where \kappa is a threshold, a constant input h and an even kernel w. A bump is a stationary solution that exceeds \kappa exactly on an interval.

For a Mexican hat, W rises to a maximum then falls towards its limit W(\infty). A narrow and a wide bump coexist when W(\infty) < \kappa - h < \max W; the wide one is stable (Exercise 8). Every translate of a bump is again a bump, so a family of bumps forms a continuous attractor: the bump stores a position.

Such persistent, localised activity models spatial working memory and head-direction cells.

5.5 Turing instability and travelling fronts

Two cases arise when the gain S'(u_0) increases.

• If \hat w is maximal at k = 0, as for a purely excitatory kernel, the first unstable mode is homogeneous: a fold as in Example 3.4.
• If \hat w is maximal at some k_c > 0, as for a Mexican hat, a periodic pattern of wavelength 2\pi/k_c appears when the gain crosses 1/\hat w(k_c): a Turing instability.

The one-dimensional problem is invariant under translations and reflections. On the critical modes e^{\pm i k_c x} these act as \mathrm{O}(2) on \mathbb{C}, so Proposition 3.13 applies: the patterns form a circle, all translates of one another.

With a_e = a_i = 1, \sigma_e = 1 and \sigma_i = 2: k_c = 1.3596, \qquad 2\pi/k_c = 4.621, \qquad \hat w(k_c) = 0.4725 , so the homogeneous state u_0 = 0 of \tau u_t = -u + w * \tanh(\beta u) loses stability at \beta_c = 1/\hat w(k_c) = 2.117.

import numpy as np

ae, se, ai, si = 1.0, 1.0, 1.0, 2.0                # Mexican hat: short excitation, wide inhibition
w_hat = lambda k: ae * np.exp(-se**2 * k**2 / 4) - ai * np.exp(-si**2 * k**2 / 4)
kc = np.sqrt(4 * np.log(ai * si**2 / (ae * se**2)) / (si**2 - se**2))
beta_c = 1 / w_hat(kc)
print(f"k_c = {kc:.4f}, wavelength = {2 * np.pi / kc:.3f}, critical gain = {beta_c:.3f}")

# tau u_t = -u + w * tanh(beta u) on a ring of length L, pseudo-spectral convolution
L, n, beta, dt = 20 * np.pi, 512, 2.5, 0.05
k = 2 * np.pi * np.fft.fftfreq(n, d=L / n)
rng = np.random.default_rng(1)
u = 1e-3 * rng.standard_normal(n)
for _ in range(int(300 / dt)):
    u += dt * (-u + np.real(np.fft.ifft(w_hat(k) * np.fft.fft(np.tanh(beta * u)))))
m = np.argmax(np.abs(np.fft.fft(u))[1:n // 2]) + 1          # dominant Fourier mode
print(f"pattern: {m} periods on the ring, wavenumber {2 * np.pi * m / L:.2f}, "
      f"amplitude {u.max():.3f}")

6. Symmetry in V1: orientation tuning, cortical patterns and hallucinations

The three applications announced by the syllabus share one mechanism: a homogeneous state of a symmetric neural field loses stability, and symmetry dictates the shape of what replaces it. On a ring of orientations, the result is a tuning curve; on the cortical sheet, a periodic pattern; seen through the retinotopic map, a geometric hallucination.

6.1 The ring model of orientation selectivity

The input that a V1 neuron receives from the LGN is only weakly tuned to orientation, yet its output is sharply tuned. Moreover the width of its tuning curve barely changes with stimulus contrast. Ben-Yishai, Bar-Or and Sompolinsky (1995) proposed that recurrent connections within a hypercolumn sharpen the tuning, and analysed the following model.

The hypercolumn is a ring of populations indexed by their preferred orientation \vartheta \in [0, \pi), with rates m(\vartheta, t): \tau\, \partial_t m(\vartheta) = -m(\vartheta) + \Big[\frac{1}{\pi}\int_0^\pi \big(J_0 + J_2 \cos 2(\vartheta - \vartheta')\big)\, m(\vartheta')\, d\vartheta' + h(\vartheta)\Big]_+ .

Here [x]_+ = \max(x, 0), and the input is h(\vartheta) = c\,\big(1 - \varepsilon + \varepsilon\cos 2(\vartheta - \vartheta_0)\big), with contrast c, anisotropy \varepsilon \in [0, \tfrac12] and stimulus orientation \vartheta_0. Uniform inhibition corresponds to J_0 < 0; J_2 > 0 favours similar orientations.

The model is equivariant under the rotations \vartheta \mapsto \vartheta + \phi and the reflection \vartheta \mapsto -\vartheta when \varepsilon = 0: the group \mathrm{O}(2) of ch. 3.6.

For J_2 > 2 the linear regime is unstable, and the rates vanish on part of the ring. This marginal phase is where the model explains sharp, contrast-invariant tuning.

import numpy as np
from scipy.optimize import brentq

n = 180
th = np.arange(n) * np.pi / n                       # preferred orientations in [0, pi)
def ring(J0, J2, eps, c=1.0, T=200.0, dt=0.05):
    W = (J0 + J2 * np.cos(2 * (th[:, None] - th[None, :]))) / n   # (1/pi) * integral -> mean
    h = c * (1 - eps + eps * np.cos(2 * (th - np.pi / 2)))         # weakly tuned input at 90 deg
    m = np.zeros(n)
    for _ in range(int(T / dt)):
        m += dt * (-m + np.maximum(W @ m + h, 0.0))
    return m

for J2, eps in [(1.0, 0.1), (3.0, 0.1), (3.0, 0.02)]:
    m = ring(J0=-2.0, J2=J2, eps=eps)
    half_width = 0.5 * np.degrees(np.pi * np.mean(m > 1e-9))
    print(f"J2={J2}, eps={eps}: max/min rate = {m.max():.3f}/{m.min():.3f}, "
          f"half-width of the active zone = {half_width:.1f} deg")

# marginal phase: half-width theta_c solves theta - sin(4 theta)/4 = pi/J2 (eps -> 0)
tc = brentq(lambda t: t - np.sin(4 * t) / 4 - np.pi / 3.0, 0.1, np.pi / 2)
print(f"theory for J2 = 3: theta_c = {np.degrees(tc):.1f} deg")

6.2 Cortical pattern formation and Euclidean symmetry

Consider now a neural field on the plane, \tau\partial_t u = -u + w * S(u) + h, with an isotropic kernel w(|x - y|). It is equivariant under the Euclidean group \mathrm{E}(2) of translations, rotations and reflections.

By Proposition 5.6, the homogeneous state destabilises on all wave vectors of length k_c at once: a whole circle of critical modes. To obtain a finite-dimensional problem, one looks for patterns that are periodic on a lattice.

• On a square lattice the critical modes are e^{i k_1 \cdot x} and e^{i k_2 \cdot x} with k_1 \perp k_2, |k_j| = k_c, and their conjugates: four real dimensions.
• On a hexagonal lattice they are e^{i k_j \cdot x} for three vectors at 120° with k_1 + k_2 + k_3 = 0: six real dimensions.

Write u \approx \sum_j z_j\, e^{i k_j \cdot x} + \text{c.c.}; the symmetry group of the lattice acts on the amplitudes z_j.

6.3 Visual hallucinations

Klüver classified the geometric hallucinations reported under mescaline into four form constants: gratings and lattices (honeycombs, chessboards); cobwebs; tunnels and funnels; spirals. Similar images are reported with flickering light, on falling asleep and in migraine auras.

Ermentrout and Cowan (1979) proposed that a drug raises the excitability of V1 until its homogeneous state undergoes the Turing instability of ch. 5. The cortex then shows stripes, squares or hexagons, and the subject perceives their preimage under the retinotopic map.

Tunnels, funnels and spirals are thus three faces of one cortical object, stripes in three directions. The model explains why these shapes recur across subjects and causes: they are images of the generic planforms of a Euclidean-invariant field.

Bressloff, Cowan, Golubitsky, Thomas and Wiener (2001) added orientation. Activity lives on the feature space \mathbb{R}^2 \times S^1 of Definition 5.2. The lateral connections link similar orientations along their common axis, so the Euclidean group acts by the shift-twist action (y, \vartheta) \mapsto (R_\phi\, y + d,\ \vartheta + \phi), \qquad (y, \vartheta) \mapsto (\kappa y,\ -\vartheta), where R_\phi is a rotation and \kappa the reflection (y_1, y_2) \mapsto (y_1, -y_2).

The bifurcating planforms are then patterns of oriented contours. Their preimages reproduce the cobwebs and the contoured lattices that a theory without orientation cannot produce.

7. Mean-field theory I: random networks of rate neurons

Cortical connectivity is sparse and heterogeneous, and the simplest model of a large network draws its synaptic weights at random. This chapter studies such random networks of rate neurons. Their quiescent state is governed by the circular law of random matrices; beyond it, dynamic mean-field theory, borrowed from the physics of spin glasses, predicts a transition to chaos; Ben Arous and Guionnet made that theory rigorous.

7.1 Random rate networks

Sompolinsky, Crisanti and Sommers (1988) introduced this model. It is the neural counterpart of the Langevin dynamics of soft spins with random couplings studied by Sompolinsky and Zippelius (1982) for spin glasses; there, the couplings are symmetric.

7.2 The circular law and the quiescent state

The origin x = 0 is always an equilibrium, with Jacobian -I + J since \phi'(0) = 1. Its stability is a question about the spectrum of a large non-symmetric random matrix.

7.3 Dynamic mean-field theory

For g > 1 the network does not settle: each unit fluctuates irregularly forever. Dynamic mean-field theory (DMFT) describes these fluctuations exactly in the limit N \to \infty.

The heuristic. The input of unit i, \eta_i(t) = \sum_j J_{ij}\phi(x_j(t)), is a sum of N weakly dependent terms with random signs. One treats it as a centred Gaussian process, independent across units and of J_{ij}.

Its covariance follows from \mathbb{E}[J_{ij}^2] = g^2/N: C_\eta(t, s) = \sum_j \mathbb{E}[J_{ij}^2]\; \mathbb{E}\big[\phi(x_j(t))\,\phi(x_j(s))\big] = g^2\, C_\phi(t, s).

Each unit then obeys the self-consistent single-unit problem \dot x = -x + \eta, \qquad \eta \sim \mathcal{GP}\big(0,\ g^2\, C_\phi\big), \qquad C_\phi(t, s) = \mathbb{E}\big[\phi(x(t))\,\phi(x(s))\big].

The asymmetry of J is what makes the input independent of the unit’s own past. With symmetric couplings, a reaction term appears: this is the harder spin-glass case.

In a stationary regime, let \Delta(\tau) = \mathbb{E}[x(t)\, x(t + \tau)] and \Delta_0 = \Delta(0). Applying (1 + \partial_t)(1 + \partial_s) to the covariance of x gives \Delta''(\tau) = \Delta(\tau) - g^2\, C_\phi\big(\Delta(\tau); \Delta_0\big), where C_\phi(\Delta; \Delta_0) is the average of \phi(X_1)\phi(X_2) over a centred Gaussian pair with variances \Delta_0 and covariance \Delta.

The static solution of part 2 exists but is unstable for g > 1. The network instead follows the solution of part 3, whose autocorrelation decays to zero: the activity is chaotic, and Sompolinsky, Crisanti and Sommers showed that its largest Lyapunov exponent is positive.

import numpy as np
from scipy.optimize import brentq

# --- simulation of dx/dt = -x + J tanh(x), J_ij ~ N(0, g^2/N) ---
rng = np.random.default_rng(0)
N, dt, T, T0 = 400, 0.1, 200.0, 100.0
G = rng.standard_normal((N, N)) / np.sqrt(N)
print(f"spectral radius of G: {np.abs(np.linalg.eigvals(G)).max():.3f} (circular law: 1)")

# --- mean-field predictions, Gauss-Hermite quadrature for E over Z ~ N(0,1) ---
z, wz = np.polynomial.hermite_e.hermegauss(80); wz = wz / wz.sum()
E = lambda f: np.sum(wz * f(z))
def static_q(g):        # fixed point: q = g^2 E tanh^2(sqrt(q) Z)
    return brentq(lambda q: g**2 * E(lambda z: np.tanh(np.sqrt(q) * z)**2) - q, 1e-6, 50)
def dmft_delta0(g):     # chaos: Delta0^2 / 2 = g^2 Var[ln cosh(sqrt(Delta0) Z)]
    Phi = lambda D: np.log(np.cosh(np.sqrt(D) * z))
    return brentq(lambda D: g**2 * (np.sum(wz * Phi(D)**2) - np.sum(wz * Phi(D))**2) - D**2 / 2,
                  1e-3, 50)

for g in [0.5, 0.9, 1.5, 2.0]:
    J, x, q = g * G, rng.standard_normal(N), []
    for k in range(int(T / dt)):
        x = x + dt * (-x + J @ np.tanh(x))
        if k * dt > T0:
            q.append(np.mean(x**2))
    line = f"g = {g}: simulated mean of x_i^2 = {np.mean(q):.3f}"
    if g > 1:
        line += f", DMFT Delta0 = {dmft_delta0(g):.3f}, static q = {static_q(g):.3f}"
    print(line)

7.4 The rigorous theory: Ben Arous and Guionnet

The derivation of ch. 7.3 assumed that the input of a unit is Gaussian and independent of the couplings. Ben Arous and Guionnet proved the corresponding statement for the Langevin dynamics of soft spins introduced by Sompolinsky and Zippelius, with independent (asymmetric) couplings.

Part 3 is the dynamic mean-field equation of ch. 7.3, with a linear interaction instead of \phi. For rate networks, Moynot and Samuelides (2002) proved the analogous result in discrete time, and Faugeras, Touboul and Cessac (2009) derived the mean-field equations of multi-population networks with random weights and noisy inputs.

The symmetric case, where the reaction term appears, is harder; Ben Arous and Guionnet (1997) treated it. Symmetric spin-glass dynamics show aging, studied in physics by Cugliandolo and Kurchan (1993) and rigorously for spherical spin glasses by Ben Arous, Dembo and Guionnet (2001).

8. Mean-field theory II: interacting particles, spiking networks and neural masses

Ch. 7 dealt with random couplings. Here every neuron interacts with every other through an average, as in a population of statistically identical cells. The mathematics is that of interacting particles: as the number of neurons grows, each one behaves like an independent copy of a single process driven by its own law. This chapter proves this propagation of chaos, then applies it to networks of spiking neurons and to the neural masses used to model EEG signals.

8.1 McKean–Vlasov processes

Consider N neurons with states X_i \in \mathbb{R} (one dimension for simplicity) and dX_i = \Big[f(X_i) + \frac1N \sum_{j=1}^N K(X_i, X_j)\Big]\, dt + \sigma\, dB_i, \qquad X_i(0) \sim \mu_0 \text{ i.i.d.}

The interaction depends on the other neurons only through their empirical measure: the particles are exchangeable.

8.2 Propagation of chaos

The initial conditions are i.i.d., hence chaotic. Propagation of chaos means that this asymptotic independence survives the interaction at every later time.

The proof is the coupling method developed by Tanaka (1978), for the Boltzmann equation of Maxwellian molecules, and by Sznitman (1991). It needs the interaction to be Lipschitz; the spiking models of §8.3 do not satisfy this, which makes them harder.

import numpy as np

# dX_i = (-X_i + I + (J/N) sum_j tanh(X_j)) dt + sigma dB_i, independent N(0,1) initial data
J, I, sigma, dt, T, R = 0.8, 0.2, 0.5, 0.02, 4.0, 100      # R independent replicas
rng = np.random.default_rng(0)
for N in [50, 200, 800, 3200]:
    X = rng.standard_normal((R, N))
    for _ in range(int(T / dt)):
        m = np.tanh(X).mean(axis=1, keepdims=True)          # empirical mean field
        X += (-X + I + J * m) * dt + sigma * np.sqrt(dt) * rng.standard_normal((R, N))
    mN = np.tanh(X).mean(axis=1)                            # one value per replica
    print(f"N = {N:5d}: mean field {mN.mean():.4f}, std over replicas {mN.std():.4f}, "
          f"std * sqrt(N) = {mN.std() * np.sqrt(N):.3f}")

8.3 Mean-field limits of spiking networks

Three families of spiking models have been studied.

• Smooth conductance-based neurons. Baladron, Fasoli, Faugeras and Touboul (2012) proved mean-field limits and propagation of chaos for noisy networks of FitzHugh–Nagumo and Hodgkin–Huxley neurons, whose coefficients are only locally Lipschitz.
• Noisy integrate-and-fire neurons with mean-field coupling. The density solves the Fokker–Planck equation of ch. 4.5, with the firing rate \nu(t), the flux at threshold, added to the drift. For strong excitation, solutions can blow up in finite time (Cáceres, Carrillo and Perthame, 2011), a signature of synchronisation. For weak coupling, Delarue, Inglis, Rubenthaler and Tanré (2015) proved global solvability and the particle approximation.
• Neurons with stochastic spiking (escape noise). This is the model below.

Each neuron has a potential that follows \dot X_i = b(X_i) between spikes and spikes at the random rate f(X_i(t)). At a spike, its potential is reset to 0 and the potential of every other neuron jumps by J/N.

This model was studied by De Masi, Galves, Löcherbach and Presutti (2015), by Fournier and Löcherbach (2016), and by Cormier, Tanré and Veltz (2020, 2021). Its mean-field limit is the McKean–Vlasov equation with jumps dX_t = b(X_t)\, dt + J\, \mathbb{E}\big[f(X_t)\big]\, dt - X_{t^-} \int_0^\infty \mathbb{1}_{\{u \le f(X_{t^-})\}}\, \mathcal{N}(dt, du), where \mathcal{N} is a Poisson random measure of intensity dt\, du; the last term implements spikes at rate f(X) and the reset to 0.

8.4 Neural masses

A neural mass describes the average activity of a cortical column with a handful of ordinary differential equations. Such models are fitted to EEG and MEG recordings.

The model of Jansen and Rit (1995) couples three populations: pyramidal cells, excitatory interneurons and inhibitory interneurons. Each synapse is a second-order linear filter, \ddot y = A a\, S(\cdot) - 2a\dot y - a^2 y, whose impulse response is A a\, t\, e^{-at}. A sigmoid converts the mean potential into a rate. The resulting six equations produce rhythms in the alpha band.

Wilson–Cowan and Jansen–Rit equations are phenomenological: they are not derived from a network of spiking neurons. For one class of networks an exact derivation exists.

The Lorentzian assumption on the initial data can be relaxed: Ott and Antonsen (2008) showed, for oscillator networks, that such a manifold of states attracts the dynamics. Unlike Wilson–Cowan models, these equations contain the mean potential v, and they produce damped oscillations caused by transient spike synchrony.

import numpy as np
from scipy.integrate import solve_ivp

Delta, J = 1.0, 15.0                                  # Lorentzian width, coupling (tau = 1)

def fixed_points(eta):
    # r > 0 solves -pi^2 r^4 + J r^3 + eta r^2 + Delta^2/(4 pi^2) = 0, then v = -Delta/(2 pi r)
    roots = np.roots([-np.pi**2, J, eta, 0.0, Delta**2 / (4 * np.pi**2)])
    out = []
    for r in sorted(roots[(abs(roots.imag) < 1e-9) & (roots.real > 0)].real):
        v = -Delta / (2 * np.pi * r)
        A = np.array([[2 * v, 2 * r], [J - 2 * np.pi**2 * r, 2 * v]])
        out.append((r, v, np.linalg.eigvals(A)))
    return out

for eta in [-6.0, -5.0, -4.0]:
    for r, v, ev in fixed_points(eta):
        print(f"eta = {eta}: r* = {r:.4f}, v* = {v:+.4f}, eigenvalues {np.round(ev, 3)}")

# response of the firing-rate equations to a step of input current
def mpr(t, y, eta, I):
    r, v = y
    return [Delta / np.pi + 2 * r * v, v**2 + eta + J * r + I(t) - (np.pi * r)**2]
sol = solve_ivp(mpr, (0, 40), [0.1, -2.0], args=(-5.0, lambda t: 3.0 if 10 < t < 30 else 0.0),
                max_step=0.01)
r = sol.y[0]
print(f"rate before the step {r[sol.t < 10][-1]:.3f}, peak {r[(sol.t > 10) & (sol.t < 30)].max():.3f}, "
      f"plateau {r[sol.t < 30][-1]:.3f}, after the step {r[-1]:.3f}")

9. Research directions and open problems

This chapter is for self-study. It points to active questions that extend each part of the course, with a few entry papers. It is a map, not a survey: the selection follows the themes of the syllabus.

9.1 Long-time behaviour of mean-field spiking models

Propagation of chaos (Theorem 8.4) holds on a fixed time interval [0, T], with constants that grow with T. What happens as t \to \infty is a different question, because the limit equation is nonlinear and may have several invariant laws or periodic solutions.

For the model with stochastic spiking of ch. 8.3, Cormier, Tanré and Veltz proved that, for small coupling J, the invariant probability is unique and attracts the solutions started from suitable initial laws (2020). They also gave conditions under which a Hopf bifurcation of the McKean–Vlasov equation produces a periodic solution (2021).

Open questions include:

• the global picture for large J;
• the relation between these mean-field cycles and the metastable oscillations of finite networks, whose law is invariant but whose trajectories oscillate for very long times;
• finite-size corrections, for which mesoscopic population equations have been proposed (Schwalger, Deger and Gerstner, 2017).

Hawkes processes offer another route from spikes to fields. Chevallier, Duarte, Löcherbach and Ost (2019) derived neural field equations as mean-field limits of spatially extended nonlinear Hawkes processes, which connects ch. 5 and ch. 8.

9.2 Blow-up and synchronisation in integrate-and-fire populations

The nonlinear Fokker–Planck equation of excitatory integrate-and-fire populations can blow up in finite time (Cáceres, Carrillo and Perthame, 2011). The blow-up corresponds to a macroscopic fraction of neurons firing at once.

Delarue, Inglis, Rubenthaler and Tanré (2015) proved that solutions exist globally when the coupling is weak. Defining and characterising the continuation of solutions after a synchronisation event, and its particle approximation, is an active topic in probability and partial differential equations.

9.3 Numerical bifurcation analysis in large dimension

Discretised neural fields and Fokker–Planck equations have 10^4 to 10^6 unknowns. Continuing their equilibria and periodic orbits requires matrix-free Newton–Krylov solvers and eigenvalue methods that exploit the structure of the problem.

Veltz and Faugeras (2010) analysed the stationary solutions of neural field equations, locally and globally. The Julia package BifurcationKit.jl, developed by Veltz, implements numerical continuation for such large systems.

Computing bifurcation diagrams of McKean–Vlasov equations, where the state is a probability density, remains a demanding test case.

9.4 Structured random networks

The i.i.d. Gaussian couplings of ch. 7 are a null model. Real circuits have cell types with signed outputs (Dale’s law), excitation balanced by inhibition, and low-dimensional structure.

• Van Vreeswijk and Sompolinsky (1996) showed that balanced excitation and inhibition produce chaotic, irregular firing.
• Kadmon and Sompolinsky (2015) extended dynamic mean-field theory to networks with several cell types.
• Mastrogiuseppe and Ostojic (2018) related low-rank connectivity to low-dimensional dynamics and computations.
• Helias and Dahmen (2020) give a field-theoretic account of these methods.

A rigorous theory of dynamic mean-field equations for nonlinear networks with symmetric, sparse or structured couplings is largely open. Trained recurrent networks, such as those obtained with the FORCE method of Sussillo and Abbott (2009), raise the same questions.

9.5 Mean-field limits of learning

The tools of ch. 8 also describe the training of wide neural networks. Mei, Montanari and Nguyen (2018) and Chizat and Bach (2018) showed that gradient descent on a two-layer network with many neurons behaves like a McKean–Vlasov dynamics on the distribution of the neurons’ parameters, a Wasserstein gradient flow. Propagation of chaos then quantifies the gap between finite and infinite width.

This is the meeting point of the course with Theoretical Foundations of Deep Learning and Optimal Transport for Machine Learning.

9.6 Cortical geometry and perception

The theory of ch. 6 predicts the form constants of hallucinations, but not which one a given subject sees, nor how patterns move. Rule, Stoffregen and Ermentrout (2011) modelled the geometric patterns induced by flickering light, where the stimulus frequency selects the pattern.

The universality of the pinwheel density (Kaschube and coauthors, 2010) raises the question of which developmental dynamics produce it. Models of the feature space \mathbb{R}^2 \times S^1 with its shift-twist symmetry also underlie geometric models of contour perception.

Exercises

The exercises are sorted by level: ★ applies the course, ★★ combines results or asks for a short proof, ★★★ is open-ended or research-flavoured. The session plan says when to do each one.

1. ★ Nernst potentials and the resting potential. At 37 °C, RT/F = 26.7 mV. A mammalian neuron has 140 mM of potassium inside and 5 mM outside, 15 mM of sodium inside and 145 mM outside. Compute E_K and E_{Na}. Then show from Definition 1.1 that, with constant conductances and no injected current, the membrane settles at a weighted mean of the reversal potentials, and explain why rest is close to E_K.

   Solution: the Nernst formula gives E_K = 26.7 \ln\frac{5}{140} \approx -89.0 \text{ mV}, \qquad E_{Na} = 26.7 \ln\frac{145}{15} \approx +60.6 \text{ mV}.

   Setting \dot V = 0 gives the conductance-weighted mean V_{\text{rest}} = \frac{\sum g_{\text{ion}} E_{\text{ion}}}{\sum g_{\text{ion}}} .

   At rest the potassium (and leak) conductances are much larger than the sodium one, so the weighted mean lies near E_K, a few millivolts above it.

2. ★ The f–I curve of the leaky integrate-and-fire neuron. Take \tau = 20 ms, E_L = V_r = -70 mV, \theta = -55 mV and \tau_{\text{ref}} = 2 ms. Compute the firing rate for RI = 20 mV, the smallest current that produces firing, and the behaviour of the rate for large I.

   Solution: by Proposition 1.5, V_\infty = -50 mV and T = 2 + 20 \ln\frac{20}{5} = 29.7 \text{ ms}, \qquad f \approx 33.6 \text{ Hz}.

   Firing requires RI > \theta - E_L = 15 mV. For large I, the logarithm is about (\theta - V_r)/RI, so T \approx \tau_{\text{ref}} + 300/(RI) ms, and the rate saturates at 1/\tau_{\text{ref}} = 500 Hz.

3. ★ Planar linear systems and the resting state of FitzHugh–Nagumo. Classify the origin for A_1 = \begin{pmatrix} -1 & 2 \\ -2 & -1 \end{pmatrix}, \quad A_2 = \begin{pmatrix} 1 & 3 \\ 1 & -1 \end{pmatrix}, \quad A_3 = \begin{pmatrix} -3 & 1 \\ 0 & -1 \end{pmatrix}. Which of them are topologically equivalent? Then check, with the formulas of Example 2.4, the eigenvalues of FitzHugh–Nagumo at I = 0.

   Solution: by Proposition 0.3:

   • A_1 has \tau = -2, \Delta = 5 and \tau^2 < 4\Delta: a stable focus.
   • A_2 has \Delta = -4 < 0: a saddle.
   • A_3 has \tau = -4, \Delta = 3 and \tau^2 \ge 4\Delta: a stable node.

   By Theorem 2.8, A_1 and A_3 are topologically conjugate; A_2 is not.

   For FitzHugh–Nagumo, v^* = -1.1994, \mathrm{tr}\,J = -0.5026 and \det J = 0.1081, so the eigenvalues are -0.251 \pm 0.212\, i.

4. ★ The bottleneck of a saddle-node. For \dot x = \mu + x^2 with \mu > 0, compute the time T(\mu) to go from x = -1 to x = 1, and its behaviour as \mu \to 0.

   Evaluate it for \mu = 0.01, and say how the firing frequency of a class 1 neuron changes when I - I_c is divided by 4.

   Solution: T(\mu) = \frac{2}{\sqrt\mu}\arctan\frac{1}{\sqrt\mu}, which behaves like \pi/\sqrt\mu. For \mu = 0.01, T = 29.42 against \pi/\sqrt\mu = 31.42.

   Dividing \mu by 4 doubles the passage time: by Proposition 3.9, the frequency of a class 1 neuron is divided by 2 when I - I_c is divided by 4.

5. ★ Itô calculus warm-up. (a) Use Itô’s formula to compute \mathbb{E}[B_t^4]. (b) Show that X_t = \exp(\sigma B_t - \sigma^2 t/2) solves dX = \sigma X\, dB. (c) For the Ornstein–Uhlenbeck process of Example 4.6 with \mu = 0, recover the stationary variance by applying Itô’s formula to X^2.

   Solution:

   • (a) d(B^4) = 4B^3\, dB + 6B^2\, dt, so \mathbb{E}B_t^4 = 6\int_0^t s\, ds = 3t^2.
   • (b) With F(t, x) = e^{\sigma x - \sigma^2 t/2}, the dt terms -\frac{\sigma^2}{2}F + \frac{\sigma^2}{2}F cancel, leaving dX = \sigma X\, dB.
   • (c) d(X^2) = 2X\, dX + \frac{\sigma^2}{\tau}\, dt. In the stationary regime the expectation of the left-hand side vanishes, so 0 = -\frac{2}{\tau}\,\mathbb{E}X^2 + \frac{\sigma^2}{\tau}, \qquad \mathbb{E}X^2 = \frac{\sigma^2}{2}.

6. ★★ A centre manifold that reverses the naive conclusion. Study the stability of the origin for \dot x = x^2 y - x^5, \qquad \dot y = -y + x^2 . Compare with the conclusion obtained by setting y = 0.

   Solution: the linearisation has eigenvalues 0 and -1. Write h(x) = \alpha x^2 + \gamma x^3 + \beta x^4 + \delta x^5 + \dots and use the invariance equation h'(x)(x^2 h - x^5) = -h + x^2.

   The left-hand side is 2\alpha^2 x^5 + O(x^6). Matching the powers x^2 to x^5 gives \alpha = 1, \gamma = \beta = 0 and \delta = -2, so h(x) = x^2 - 2x^5 + O(x^6).

   The reduced equation is \dot x = x^2 h(x) - x^5 = x^4 - x^5 + O(x^6).

   Points with x > 0 on the centre manifold move away, so the origin is unstable (by Theorem 2.11). Setting y = 0 would give \dot x = -x^5 and the wrong conclusion, stability.

7. ★★ First Lyapunov coefficients. (a) Apply the formula of ch. 3.3 to \dot x = \mu x - y - x(x^2 + y^2), \dot y = x + \mu y - y(x^2 + y^2), and describe the cycle. (b) Same question with the cubic terms reversed in sign. (c) Show that \dot x = -y + x^2 + xy, \dot y = x + \mu y undergoes a subcritical Hopf bifurcation at \mu = 0.

   Solution:

   • (a) With f = -x^3 - xy^2 and g = -x^2 y - y^3, the third derivatives are f_{xxx} = -6, \quad f_{xyy} = -2, \quad g_{xxy} = -2, \quad g_{yyy} = -6, and the quadratic bracket vanishes. So 16a = -16 and a = -1: supercritical, with a stable cycle of radius \sqrt\mu for \mu > 0, since \dot\rho = \rho(\mu - \rho^2).

   • (b) Now a = +1: subcritical, with an unstable cycle of radius \sqrt{-\mu} for \mu < 0.

   • (c) At \mu = 0 the linear part is the rotation with \omega = 1. The nonlinear terms are f = x^2 + xy and g = 0, so f_{xx} = 2 and f_{xy} = 1.

     Then 16a = f_{xy}(f_{xx} + f_{yy}) = 2, so a = 1/8 > 0. The eigenvalues \mu/2 \pm i\sqrt{1 - \mu^2/4} cross the axis with speed 1/2: the bifurcation is subcritical.

8. ★★ Amari bumps for a Mexican hat. Let w(x) = (1 - |x|)\, e^{-|x|}. Compute W, find the bumps of Theorem 5.5 when \kappa - h = 0.2, check the profile conditions and decide their stability.

   Solution: W(a) = a e^{-a}, which increases to its maximum 1/e at a = 1 and then decreases. Bumps exist when 0 < \kappa - h < 1/e; for \kappa - h = 0.2, a e^{-a} = 0.2 gives a_1 = 0.2592, \qquad a_2 = 2.5426 .

   Stability: w(a_1) > 0, so the narrow bump is unstable, and w(a_2) < 0, so the wide bump is stable.

   Profile of the wide bump: at its centre u = a_2 e^{-a_2/2} + h, which exceeds \kappa by 0.51; far away u \to h < \kappa.

9. ★★ Turing instability of a difference of Gaussians. Take a_e = 2, \sigma_e = 1, a_i = 1.5 and \sigma_i = 3 in Proposition 5.7, and S(u) = \tanh(\beta u) around u_0 = 0.

   Find k_c, the wavelength and the critical gain, and check that the instability is not homogeneous.

   Solution: a_i\sigma_i^2 = 13.5 > a_e \sigma_e^2 = 2, so k_c^2 = \frac{4 \ln 6.75}{8} = 0.9548, \qquad k_c = 0.977, \qquad 2\pi/k_c = 6.43 .

   Then \hat w(k_c) = 1.400 while \hat w(0) = 0.5. By Proposition 5.6, the state loses stability at \beta_c = 1/1.400 = 0.714, on the mode k_c first since \hat w(k_c) > \hat w(0): a Turing instability.

10. ★★ Mean first-passage times. (a) Show with Proposition 4.8 that a Brownian motion \sigma B started at x < \theta reaches \theta almost surely but with infinite mean time. (b) With a reflecting barrier at 0, compute T(x) on [0, \theta].

    (c) Show that for \mu far below \theta the Siegert formula gives \frac{1}{\nu} \approx \tau\sqrt\pi\; \frac{e^{U^2}}{U}, \qquad U = \frac{\theta - \mu}{\sigma}, with \tau_{\text{ref}} negligible.

    Solution:

    • (a) With b = 0, \psi = 1 and the inner integral \int_{-\infty}^y 2\sigma^{-2}\, dz diverges, so T = \infty; the hitting itself is almost sure by recurrence of one-dimensional Brownian motion.
    • (b) Solve \frac{\sigma^2}{2}T'' = -1 with T'(0) = 0 and T(\theta) = 0: T(x) = (\theta^2 - x^2)/\sigma^2.
    • (c) For large U, the integrand e^{u^2}(1 + \mathrm{erf}\,u) is dominated by the upper end, where 1 + \mathrm{erf}\,u \approx 2, and \int^U e^{u^2} du \approx e^{U^2}/(2U). The rate is exponentially small in U^2: an Arrhenius law.

11. ★★ Breaking the O(2) symmetry. (a) Show that a smooth map f : \mathbb{C} \to \mathbb{C} that commutes only with rotations has the form p(|z|^2)\, z with p complex. (b) For \dot z = \mu z - |z|^2 z with \mu > 0, find the equilibria and the eigenvalues at z^* = \sqrt\mu.

    (c) Add a small real input: \dot z = \mu z - |z|^2 z + \varepsilon with \varepsilon > 0. Show that only two equilibria remain near the circle, and decide their stability. Relate this to the ring model of ch. 6.1.

    Solution:

    • (a) The first step of the proof of Proposition 3.13 holds without the reflection, except that f(r) need not be real.
    • (b) The circle |z| = \sqrt\mu, plus the origin. At z^* the radial eigenvalue is -2\mu and the tangential one is 0.
    • (c) In polar coordinates, \dot\rho = \mu\rho - \rho^3 + \varepsilon\cos\phi and \dot\phi = -\varepsilon\sin\phi/\rho. Equilibria near the circle need \sin\phi = 0. At \phi = 0 the angular perturbation decays (stable); at \phi = \pi it grows (saddle).

    In the ring model, the anisotropy \varepsilon of the input pins the bump at the stimulus orientation, and the opposite orientation is unstable.

12. ★★ The static mean-field equation (Python). (a) Expand \mathbb{E}[\tanh^2(\sqrt q Z)] to second order in q and deduce the behaviour of the static solution of Proposition 7.4 as g \to 1^+. (b) Compute it numerically for g \in \{1.05, 1.1, 1.5, 2\}. (c) Compare with the dynamic \Delta_0 of Example 7.5 and explain why the network does not follow the static solution.

    Solution: (a) \tanh^2 x = x^2 - \frac23 x^4 + O(x^6) and \mathbb{E}Z^4 = 3, so \mathbb{E}\tanh^2(\sqrt q Z) = q - 2q^2 + O(q^3). The equation q = g^2(q - 2q^2) gives q \approx \frac{g^2 - 1}{2 g^2} \approx g - 1 .

    (b) The code below gives q = 0.0533, 0.1127, 0.7934 and 2.1215 for the four values of g. The expansion gives 0.0465, 0.0868, 0.278 and 0.375: it is accurate only near g = 1.

    import numpy as np
    from scipy.optimize import brentq
    
    z, wz = np.polynomial.hermite_e.hermegauss(80)
    wz = wz / wz.sum()                                   # weights of E over Z ~ N(0, 1)
    for g in [1.05, 1.1, 1.5, 2.0]:
        F = lambda q: g**2 * np.sum(wz * np.tanh(np.sqrt(q) * z)**2) - q
        q = brentq(F, 1e-8, 50)
        print(f"g = {g}: static q = {q:.4f}, expansion {(g**2 - 1) / (2 * g**2):.4f}")

    (c) For g = 1.5 the static value is 0.793, while the network fluctuates with variance about 0.72 to 0.75. The static solution describes a network sitting at a fixed point x^*.

    Heuristically, the Jacobian -I + J\,\mathrm{diag}(\tanh'(x^*)) has its eigenvalues in a disc of centre -1 and radius \rho = g\, \big(\mathbb{E}[\tanh'(\sqrt q Z)^2]\big)^{1/2}.

    Numerically \rho = 1.067 for g = 1.5 and \rho = 1.167 for g = 2. Such fixed points are therefore unstable, and the network follows the chaotic solution of Proposition 7.4 (part 3).

13. ★★★ When the mean field is bistable (Python). In the particle system of Example 8.5, set I = 0, J = 2 and \sigma = 1.2.

    (a) Show that the stationary McKean–Vlasov solutions have \mu = \mathcal{N}(Jm, \sigma^2/2) with m = \mathbb{E}\tanh(Jm + \sigma Z/\sqrt2), and that there are three of them.

    (b) Simulate the particle system for N = 2, 5, 10, 20 up to T = 4000, and count the switches of the empirical mean field between the wells \pm m^*.

    A switch is a passage from above m^*/2 to below -m^*/2, or back. (c) Explain why this does not contradict Theorem 8.4.

    Hints and expected result:

    • (a) In a stationary state the mean field is a constant m, and each particle is an Ornstein–Uhlenbeck process (Example 4.6). The slope of the right-hand side at m = 0 is J\,\mathbb{E}\big[\cosh^{-2}(\sigma Z/\sqrt2)\big] \approx 1.33 > 1, so besides m = 0 there are two symmetric solutions \pm m^* with m^* \approx 0.82.
    • (b) With time step 0.02 and a fixed seed, one run gives about 230, 70 and 20 switches for N = 2, 5 and 10, and none for N = 20. With the stricter thresholds \pm m^* the counts drop to about 140, 45 and 15. The mean time between switches grows roughly exponentially with N, as for an escape over a barrier.
    • (c) The bound of Theorem 8.4 grows like e^{3LT}, so it is uninformative on the switching time scales. The limits N \to \infty and t \to \infty do not commute.

14. ★★★ Hexagons before the instability. Use the hexagonal amplitude equations of Proposition 6.4 with real coefficients, a < 0 and a + 2b < 0. (a) Compute the stripes. (b) Look for hexagons z_1 = z_2 = z_3 = R and show that they exist for some \mu < 0, before the homogeneous state loses stability. (c) What does this predict for the onset of geometric hallucinations as excitability increases, and which form constants should appear first? Discuss what a full answer would require.

    Hints and expected result:

    • (a) Stripes have |z_1|^2 = -\mu/a, for \mu > 0.
    • (b) Hexagons solve \mu + \alpha R + (a + 2b) R^2 = 0, with real roots when \mu \ge \mu_- = \frac{\alpha^2}{4(a + 2b)} , which is negative. Hexagons therefore exist before the instability; they bifurcate subcritically, which implies hysteresis.
    • (c) With \alpha \ne 0, the first patterns are hexagonal, whose log-polar images are lattice-like (first form constant), and they appear abruptly.

    A full answer needs the stability of each branch (eigenvalues in the six-dimensional amplitude space, Golubitsky, Stewart and Schaeffer, 1988), the coefficients computed from a specific kernel and rate function, and the contoured version of Bressloff and coauthors (2001).

15. ★★★ Several stationary states of a spiking mean-field model (Python). For the model of Proposition 8.6 with b(x) = 1 - x, J = 0.5 and f(x) = x^2: (a) solve the self-consistency equation; (b) simulate the particle system with N = 2000 and compare the firing rates.

    (c) Then take I = 0.5 and f(x) = x^8, and explore J. Find parameters with several stationary rates, and discuss what the simulation of a finite network shows there.

    Solution for (a) and (b): the code below gives \alpha^* = 0.6241 for the limit equation and a rate of 0.6227 for the particle system.

    import numpy as np
    from scipy.integrate import quad
    from scipy.optimize import brentq
    
    I, J = 1.0, 0.5
    f = lambda x: x**2                                   # spiking intensity
    
    def mean_isi(alpha):
        xinf = I + J * alpha                             # the flow -x + xinf never reaches xinf
        b = lambda x: xinf - x
        H = lambda x: np.exp(-quad(lambda y: f(y) / b(y), 0, x)[0])   # survival up to level x
        return quad(lambda x: H(x) / b(x), 0, xinf, limit=200)[0]
    
    alpha = brentq(lambda a: a * mean_isi(a) - 1, 1e-3, 10)
    print(f"self-consistent rate alpha* = {alpha:.4f}")
    
    # particle system: each spike resets the neuron to 0 and kicks every neuron by J/N
    rng = np.random.default_rng(0)
    N, dt, T, T0 = 2000, 0.005, 60.0, 10.0
    X, count = np.zeros(N), 0
    for k in range(int(T / dt)):
        fired = rng.random(N) < f(X) * dt
        X += (I - X) * dt + J * fired.sum() / N
        X[fired] = 0.0
        if k * dt > T0:
            count += fired.sum()
    print(f"particle system, N = {N}: rate = {count / (N * (T - T0)):.4f}")

    Expected result for (c): for J = 1.5 the equation has three roots, \alpha \approx 0.0043, 0.411 and 3.46. A low state, where neurons wait near x_\alpha and fire at rate about f(x_\alpha), coexists with a high state.

    Finite networks switch between them at rare random times. The middle root is expected to be unstable; proving it requires the linearisation of the McKean–Vlasov equation, as in the work of Cormier, Tanré and Veltz.

Projects and assessment

The official assessment is an exam for the first session, with a resit in the form of a paper reading (master-mva.com, “Mode de validation”). The Université Paris-Saclay page classifies the course as “Evaluation Continue non Intégrale” (partial continuous assessment), and no weighting is displayed. The proposals below are not official.

The written-exam proposal prepares the first session; the paper-study projects train the skills of the resit, reading a paper, reproducing its results and presenting them, and can serve as homework. Every project is evaluated with the rubric at the end of this section, and its success criterion is the target of the “Numerical results” line.

Solo projects

1. [Solo · ≈ 20 h] The bifurcation diagram of the Hodgkin–Huxley model. Paper study of Rinzel and Miller (1980). Question: where does repetitive firing appear, and why do rest and firing coexist? Implement:

   • continuation of the equilibria in I, with their eigenvalues;
   • the first Lyapunov coefficient at the Hopf point, with the formula of Example 3.6 in dimension 4;
   • stable cycles by simulation and unstable cycles by shooting (Newton’s method on the return map), to locate the fold of cycles.

   Success criterion: Hopf point at I \approx 9.78 µA/cm², \ell_1 > 0, fold of cycles near 6.3 µA/cm². Deliverables: code, a 4–6 page report, a 10-minute talk. Chapters 1, 2, 3.

2. [Solo · ≈ 20 h] Noise-driven firing: Siegert formula and Fokker–Planck solver. Question: how does noise shape the f–I curve and the irregularity of firing? Implement:

   • Euler–Maruyama simulations of the noisy integrate-and-fire neuron for several \sigma;
   • a finite-volume solver of the Fokker–Planck equation with absorbing threshold and reinjection at V_r;
   • the coefficient of variation of the interspike intervals.

   Success criterion: the firing rates of the simulations (with a Brownian-bridge test for crossings between two steps, or a time step halved until convergence) and of the Fokker–Planck solver agree within 2 % with Theorem 4.9 over a grid of (\mu, \sigma), and the CV tends to 1 for strongly subthreshold input. Deliverables: code, a 4–6 page report, a 10-minute talk. Chapter 4.

3. [Solo · ≈ 20 h] Orientation tuning in the ring model. Paper study of Ben-Yishai, Bar-Or and Sompolinsky (1995). Question: when do recurrent connections, rather than the input, set the tuning width? Reproduce the phase diagram in the plane (J_0, J_2) (homogeneous regime, marginal phase, unbounded activity), the tuning width against contrast, and the response to a rotating stimulus. Success criterion: the measured half-width matches Proposition 6.2 within 2° as \varepsilon \to 0. Deliverables: code, a 4–6 page report, a 10-minute talk. Chapters 3, 6.

Group projects

1. [Group of 2–3 · ≈ 30 h each] The transition to chaos in random networks. Paper study of Sompolinsky, Crisanti and Sommers (1988). The work splits as follows:

   • member A simulates networks up to N = 2000 and computes the largest Lyapunov exponent (Benettin’s algorithm);
   • member B solves the dynamic mean-field equation of Proposition 7.4 for the whole autocorrelation function \Delta(\tau);
   • member C studies the spectra of the Jacobians, at rest and along trajectories, against the circular law.

   Success criterion: \Delta_0 within 3 % of the simulations for g \in [1.2, 3] at N = 2000, and a Lyapunov exponent that changes sign near g = 1. Deliverables: shared code, an 8–10 page report, a 20-minute talk. Chapters 0, 7.

2. [Group of 2–3 · ≈ 30 h each] Geometric visual hallucinations. Paper studies of Ermentrout and Cowan (1979) and Bressloff and coauthors (2001). The work splits as follows:

   • member A simulates a two-dimensional neural field on a torus with FFT convolutions and checks the selected wavelength;
   • member B derives the amplitude equations of Proposition 6.4 for the chosen kernel and rate function, and predicts stripes or hexagons;
   • member C maps the patterns back to the visual field with the inverse log-polar map and compares them with Klüver’s form constants.

   Success criterion: wavelength within 5 % of 2\pi/k_c, and both stripes and hexagons obtained by changing the offset of the rate function. Deliverables: shared code, an 8–10 page report with figures, a 20-minute talk. Chapters 5, 6.

3. [Group of 2–3 · ≈ 30 h each] From spiking neurons to an exact neural mass. Paper study of Montbrió, Pazó and Roxin (2015). The work splits as follows:

   • member A simulates 10^4 theta neurons with Lorentzian excitabilities, with a time step adapted to the fastest neurons;
   • member B computes the fixed points and the folds of the firing-rate equations of Theorem 8.7 as functions of \bar\eta;
   • member C compares responses to current steps and quantifies finite-size effects.

   Success criterion: stationary rates within 5 % of Theorem 8.7, and the bistable range in \bar\eta located on both sides. Deliverables: shared code, an 8–10 page report, a 20-minute talk. Chapters 1, 3, 8.

Written exam (proposal)

Three hours, no documents except one handwritten A4 sheet, four exercises. The official page says only that the first session is an exam: this format is a proposal.

1. Question (4 points). State the centre manifold theorem (hypotheses, tangency, reduction principle). Apply it to \dot x = xy, \dot y = -y - x^2: why does the linearisation decide nothing, and what is the stability of the origin?

   Expected: the statement of Theorem 2.11; the eigenvalues 0 and -1; h(x) = -x^2 + O(x^4) from the invariance equation; the reduced equation \dot x = -x^3 + O(x^5), hence asymptotic stability, with algebraic decay.

2. Question (5 points). Consider \dot x = \mu x - y - x^3, \dot y = x + \mu y. Find the eigenvalues at the origin, show that a Hopf bifurcation occurs at \mu = 0, compute the sign of the first Lyapunov coefficient, and give the amplitude and period of the cycle for small \mu.

   Expected: eigenvalues \mu \pm i, crossing the axis with speed 1. With f = -x^3 and g = 0, the formula gives 16a = f_{xxx} = -6, so a = -3/8: supercritical.

   A stable cycle of radius \sqrt{8\mu/3} exists for \mu > 0, with period close to 2\pi.

3. Question (5 points). Let f : \mathbb{C} \to \mathbb{C} be smooth and \mathrm{O}(2)-equivariant. Prove that f(z) = p(|z|^2)\, z with p real. Deduce that a bifurcation from z = 0 produces a circle of equilibria, each with a zero eigenvalue. Explain what this means for the ring model with an untuned input.

   Expected: the proof of Proposition 3.13. Rotations give f(re^{i\phi}) = e^{i\phi}f(r), the reflection makes f(r) real, and parity gives smoothness in r^2. The equilibria form the circle p(\rho^2) = 0, and the tangent vector iz^* lies in the kernel. In the ring model, a bump can sit at any orientation and drifts neutrally.

4. Question (6 points). Consider, with i.i.d. initial data, dX_i = \Big(-X_i + \frac{J}{N}\sum_j \tanh X_j\Big)\, dt + \sigma\, dB_i .

   (a) Write the McKean–Vlasov limit and its Fokker–Planck equation (2 points). (b) Show that the stationary laws are \mathcal{N}(Jm, \sigma^2/2), where m solves m = \mathbb{E}\tanh(Jm + \sigma Z/\sqrt2) (2 points). (c) Prove that m = 0 is the only solution when J < 1 (2 points).

   Expected: (a) Definition 8.1 with K(x, y) = J\tanh y. (b) For fixed m, each particle is an Ornstein–Uhlenbeck process whose stationary law is given by Example 4.6, and consistency gives the equation. (c) The right-hand side G(m) has |G'(m)| \le J < 1, so G is a contraction with fixed point 0.

Grading rubric (proposal)

References

Official references (master-mva.com, 2026–2027)

• W. Gerstner and W. M. Kistler. Spiking Neuron Models: Single Neurons, Populations, Plasticity. Cambridge University Press, 2002.
• Yu. A. Kuznetsov. Elements of Applied Bifurcation Theory. Springer, Applied Mathematical Sciences 112, 1995 (3rd edition 2004, cited in the text).
• E. M. Izhikevich. Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. MIT Press, 2006 as listed on master-mva.com (the printed edition is dated 2007).
• G. B. Ermentrout and D. H. Terman. Mathematical Foundations of Neuroscience. Springer, Interdisciplinary Applied Mathematics 35, 2010.

Textbooks

• P. Dayan and L. F. Abbott. Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems. MIT Press, 2001.
• W. Gerstner, W. M. Kistler, R. Naud and L. Paninski. Neuronal Dynamics: From Single Neurons to Networks and Models of Cognition. Cambridge University Press, 2014.
• J. Guckenheimer and P. Holmes. Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields. Springer, Applied Mathematical Sciences 42, 1983.
• L. Perko. Differential Equations and Dynamical Systems. Springer, Texts in Applied Mathematics 7, 3rd edition, 2001.
• J. Carr. Applications of Centre Manifold Theory. Springer, Applied Mathematical Sciences 35, 1981.
• M. Haragus and G. Iooss. Local Bifurcations, Center Manifolds, and Normal Forms in Infinite-Dimensional Dynamical Systems. Springer, 2011.
• M. Golubitsky, I. Stewart and D. G. Schaeffer. Singularities and Groups in Bifurcation Theory, Volume II. Springer, Applied Mathematical Sciences 69, 1988.
• I. Karatzas and S. E. Shreve. Brownian Motion and Stochastic Calculus. Springer, Graduate Texts in Mathematics 113, 2nd edition, 1991.
• B. Øksendal. Stochastic Differential Equations: An Introduction with Applications. Springer, 6th edition, 2003.
• P. E. Kloeden and E. Platen. Numerical Solution of Stochastic Differential Equations. Springer, 1992.
• A.-S. Sznitman. Topics in propagation of chaos. In École d’Été de Probabilités de Saint-Flour XIX — 1989, Lecture Notes in Mathematics 1464, Springer, 1991.
• P. C. Bressloff. Waves in Neural Media: From Single Neurons to Neural Fields. Springer, 2014.
• S. Coombes, P. beim Graben, R. Potthast and J. Wright (eds.). Neural Fields: Theory and Applications. Springer, 2014.
• M. Helias and D. Dahmen. Statistical Field Theory for Neural Networks. Springer, Lecture Notes in Physics 970, 2020.
• H. Klüver. Mescal and Mechanisms of Hallucinations. University of Chicago Press, 1966.

Neuron and synapse models (ch. 1, 3)

• L. Lapicque. Recherches quantitatives sur l’excitation électrique des nerfs traitée comme une polarisation. Journal de Physiologie et de Pathologie Générale 9, 1907.
• A. L. Hodgkin. The local electric changes associated with repetitive action in a non-medullated axon. Journal of Physiology 107, 1948.
• A. L. Hodgkin and A. F. Huxley. A quantitative description of membrane current and its application to conduction and excitation in nerve. Journal of Physiology 117, 500–544, 1952.
• R. FitzHugh. Impulses and physiological states in theoretical models of nerve membrane. Biophysical Journal 1, 445–466, 1961.
• J. Nagumo, S. Arimoto and S. Yoshizawa. An active pulse transmission line simulating nerve axon. Proceedings of the IRE 50, 2061–2070, 1962.
• C. Morris and H. Lecar. Voltage oscillations in the barnacle giant muscle fiber. Biophysical Journal 35, 193–213, 1981.
• J. Rinzel and R. N. Miller. Numerical calculation of stable and unstable periodic solutions to the Hodgkin–Huxley equations. Mathematical Biosciences 49, 27–59, 1980.
• J. Rinzel and G. B. Ermentrout. Analysis of neural excitability and oscillations. In C. Koch and I. Segev (eds.), Methods in Neuronal Modeling, MIT Press, 1989.
• G. B. Ermentrout and N. Kopell. Parabolic bursting in an excitable system coupled with a slow oscillation. SIAM Journal on Applied Mathematics 46, 233–253, 1986.
• E. M. Izhikevich. Simple model of spiking neurons. IEEE Transactions on Neural Networks 14, 1569–1572, 2003.
• A. Destexhe, Z. F. Mainen and T. J. Sejnowski. Synthesis of models for excitable membranes, synaptic transmission and neuromodulation using a common kinetic formalism. Journal of Computational Neuroscience 1, 195–230, 1994.
• M. L. Hines and N. T. Carnevale. The NEURON simulation environment. Neural Computation 9, 1179–1209, 1997.
• M. M. Peixoto. Structural stability on two-dimensional manifolds. Topology 1, 101–120, 1962.

Noise (ch. 4)

• N. Wiener. Differential-space. Journal of Mathematics and Physics 2, 131–174, 1923.
• W. R. Softky and C. Koch. The highly irregular firing of cortical cells is inconsistent with temporal integration of random EPSPs. Journal of Neuroscience 13, 334–350, 1993.
• N. Brunel and V. Hakim. Fast global oscillations in networks of integrate-and-fire neurons with low firing rates. Neural Computation 11, 1621–1671, 1999.
• N. Brunel. Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons. Journal of Computational Neuroscience 8, 183–208, 2000.
• A. S. Pikovsky and J. Kurths. Coherence resonance in a noise-driven excitable system. Physical Review Letters 78, 775–778, 1997.
• D. J. Higham. An algorithmic introduction to numerical simulation of stochastic differential equations. SIAM Review 43, 525–546, 2001.

Visual cortex, neural fields and symmetry (ch. 5, 6)

• D. H. Hubel and T. N. Wiesel. Receptive fields of single neurones in the cat’s striate cortex. Journal of Physiology 148, 574–591, 1959.
• D. H. Hubel and T. N. Wiesel. Receptive fields, binocular interaction and functional architecture in the cat’s visual cortex. Journal of Physiology 160, 106–154, 1962.
• E. L. Schwartz. Spatial mapping in the primate sensory projection: analytic structure and relevance to perception. Biological Cybernetics 25, 181–194, 1977.
• J. C. Horton and W. F. Hoyt. The representation of the visual field in human striate cortex: a revision of the classic Holmes map. Archives of Ophthalmology 109, 816–824, 1991.
• T. Bonhoeffer and A. Grinvald. Iso-orientation domains in cat visual cortex are arranged in pinwheel-like patterns. Nature 353, 429–431, 1991.
• W. H. Bosking, Y. Zhang, B. Schofield and D. Fitzpatrick. Orientation selectivity and the arrangement of horizontal connections in tree shrew striate cortex. Journal of Neuroscience 17, 2112–2127, 1997.
• F. Wolf and T. Geisel. Spontaneous pinwheel annihilation during visual development. Nature 395, 73–78, 1998.
• M. Kaschube, M. Schnabel, S. Löwel, D. M. Coppola, L. E. White and F. Wolf. Universality in the evolution of orientation columns in the visual cortex. Science 330, 1113–1116, 2010.
• H. R. Wilson and J. D. Cowan. Excitatory and inhibitory interactions in localized populations of model neurons. Biophysical Journal 12, 1–24, 1972.
• H. R. Wilson and J. D. Cowan. A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue. Kybernetik 13, 55–80, 1973.
• S. Amari. Dynamics of pattern formation in lateral-inhibition type neural fields. Biological Cybernetics 27, 77–87, 1977.
• R. Ben-Yishai, R. L. Bar-Or and H. Sompolinsky. Theory of orientation tuning in visual cortex. Proceedings of the National Academy of Sciences USA 92, 3844–3848, 1995.
• G. B. Ermentrout and J. D. Cowan. A mathematical theory of visual hallucination patterns. Biological Cybernetics 34, 137–150, 1979.
• P. C. Bressloff, J. D. Cowan, M. Golubitsky, P. J. Thomas and M. C. Wiener. Geometric visual hallucinations, Euclidean symmetry and the functional architecture of striate cortex. Philosophical Transactions of the Royal Society of London B 356, 299–330, 2001.
• M. Rule, M. Stoffregen and B. Ermentrout. A model for the origin and properties of flicker-induced geometric phosphenes. PLoS Computational Biology 7, e1002158, 2011.
• R. Veltz and O. Faugeras. Local/global analysis of the stationary solutions of some neural field equations. SIAM Journal on Applied Dynamical Systems 9, 954–998, 2010.

Random networks and spin glasses (ch. 7)

• H. Sompolinsky and A. Zippelius. Relaxational dynamics of the Edwards–Anderson model and the mean-field theory of spin-glasses. Physical Review B 25, 6860–6875, 1982.
• H. Sompolinsky, A. Crisanti and H. J. Sommers. Chaos in random neural networks. Physical Review Letters 61, 259–262, 1988.
• S. Geman. The spectral radius of large random matrices. Annals of Probability 14, 1318–1328, 1986.
• T. Tao, V. Vu and M. Krishnapur. Random matrices: universality of ESDs and the circular law. Annals of Probability 38, 2023–2065, 2010.
• G. Ben Arous and A. Guionnet. Large deviations for Langevin spin glass dynamics. Probability Theory and Related Fields 102, 455–509, 1995.
• G. Ben Arous and A. Guionnet. Symmetric Langevin spin glass dynamics. Annals of Probability 25, 1367–1422, 1997.
• L. F. Cugliandolo and J. Kurchan. Analytical solution of the off-equilibrium dynamics of a long-range spin-glass model. Physical Review Letters 71, 173–176, 1993.
• G. Ben Arous, A. Dembo and A. Guionnet. Aging of spherical spin glasses. Probability Theory and Related Fields 120, 1–67, 2001.
• O. Moynot and M. Samuelides. Large deviations and mean-field theory for asymmetric random recurrent neural networks. Probability Theory and Related Fields 123, 41–75, 2002.
• O. Faugeras, J. Touboul and B. Cessac. A constructive mean-field analysis of multi-population neural networks with random synaptic weights and stochastic inputs. Frontiers in Computational Neuroscience 3:1, 2009.

Interacting particles, spiking networks and neural masses (ch. 8, 9)

• H. P. McKean. A class of Markov processes associated with nonlinear parabolic equations. Proceedings of the National Academy of Sciences USA 56, 1907–1911, 1966.
• H. Tanaka. Probabilistic treatment of the Boltzmann equation of Maxwellian molecules. Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete 46, 67–105, 1978.
• J. Baladron, D. Fasoli, O. Faugeras and J. Touboul. Mean-field description and propagation of chaos in networks of Hodgkin–Huxley and FitzHugh–Nagumo neurons. Journal of Mathematical Neuroscience 2:10, 2012.
• M. J. Cáceres, J. A. Carrillo and B. Perthame. Analysis of nonlinear noisy integrate & fire neuron models: blow-up and steady states. Journal of Mathematical Neuroscience 1:7, 2011.
• F. Delarue, J. Inglis, S. Rubenthaler and E. Tanré. Global solvability of a networked integrate-and-fire model of McKean–Vlasov type. Annals of Applied Probability 25, 2015.
• A. De Masi, A. Galves, E. Löcherbach and E. Presutti. Hydrodynamic limit for interacting neurons. Journal of Statistical Physics 158, 866–902, 2015.
• N. Fournier and E. Löcherbach. On a toy model of interacting neurons. Annales de l’Institut Henri Poincaré, Probabilités et Statistiques 52, 2016.
• Q. Cormier, E. Tanré and R. Veltz. Long time behavior of a mean-field model of interacting neurons. Stochastic Processes and their Applications 130, 2020.
• Q. Cormier, E. Tanré and R. Veltz. Hopf bifurcation in a mean-field model of spiking neurons. Electronic Journal of Probability 26, 2021.
• B. H. Jansen and V. G. Rit. Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns. Biological Cybernetics 73, 357–366, 1995.
• E. Ott and T. M. Antonsen. Low dimensional behavior of large systems of globally coupled oscillators. Chaos 18, 037113, 2008.
• E. Montbrió, D. Pazó and A. Roxin. Macroscopic description for networks of spiking neurons. Physical Review X 5, 021028, 2015.
• T. Schwalger, M. Deger and W. Gerstner. Towards a theory of cortical columns: from spiking neurons to interacting neural populations of finite size. PLoS Computational Biology 13, e1005507, 2017.
• J. Chevallier, A. Duarte, E. Löcherbach and G. Ost. Mean field limits for nonlinear spatially extended Hawkes processes with exponential memory kernels. Stochastic Processes and their Applications 129, 1–27, 2019.
• R. Veltz. BifurcationKit.jl, a Julia package for numerical bifurcation analysis. Software, 2020.
• C. van Vreeswijk and H. Sompolinsky. Chaos in neuronal networks with balanced excitatory and inhibitory activity. Science 274, 1724–1726, 1996.
• J. Kadmon and H. Sompolinsky. Transition to chaos in random neuronal networks. Physical Review X 5, 041030, 2015.
• F. Mastrogiuseppe and S. Ostojic. Linking connectivity, dynamics, and computations in low-rank recurrent neural networks. Neuron 99, 609–623, 2018.
• D. Sussillo and L. F. Abbott. Generating coherent patterns of activity from chaotic neural networks. Neuron 63, 544–557, 2009.
• S. Mei, A. Montanari and P.-M. Nguyen. A mean field view of the landscape of two-layer neural networks. Proceedings of the National Academy of Sciences USA 115, E7665–E7671, 2018.
• L. Chizat and F. Bach. On the global convergence of gradient descent for over-parameterized models using optimal transport. Advances in Neural Information Processing Systems 31, 2018.