<figure><img src="image_1.png" /></figure>

HAL
open science

# Nonlinear partial differential equations in neuroscience: from modelling to mathematical theory

José A Carrillo, Pierre Roux

► To cite this version:

José A Carrillo, Pierre Roux. Nonlinear partial differential equations in neuroscience: from modelling to mathematical theory. 2025. hal-04875366

HAL Id: hal-04875366

https://hal.science/hal-04875366v1

Preprint submitted on 8 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Nonlinear partial differential equations in neuroscience: from modelling to mathematical theory

José A. Carrillo* Pierre Roux†

January 8, 2025

## Abstract

Many systems of partial differential equations have been proposed as simplified representations of complex collective behaviours in large networks of neurons. In this survey, we briefly discuss their derivations and then review the mathematical methods developed to handle the unique features of these models, which are often nonlinear and non-local. The first part focuses on parabolic Fokker-Planck equations: the Nonlinear Noisy Leaky Integrate and Fire neuron model, stochastic neural fields in PDE form with applications to grid cells, and rate-based models for decision-making. The second part concerns hyperbolic transport equations, namely the model of the Time Elapsed since the last discharge and the jump-based Leaky Integrate and Fire model. The last part covers some kinetic mesoscopic models, with particular attention to the kinetic Voltage-Conductance model and FitzHugh-Nagumo kinetic Fokker-Planck systems.

**Keywords:** Neuroscience, Partial Differential Equations, Neural Networks, Mathematical Biology.

**Mathematics Subject Classification:** Primary: 35-02, 92-02, 35Q80; Secondary: 35B40, 35B10, 65C30.

<footnote>*Mathematical Institute, University of Oxford, Oxford OX2 6GG, UK (carrillo@maths.ox.ac.uk)</footnote>

<footnote>†Institut Camille Jordan, École Centrale de Lyon, Ecully 69134, France (pierre.roux@ec-lyon.fr)</footnote>
# Contents

I Parabolic equations of Fokker-Planck type 4

1 The Nonlinear Noisy Leaky Integrate & Fire model for neural networks 4

1.1 Derivation of the model à la physicienne 5

1.2 Mathematical setting for the NNLIF type models 6

1.3 Characterisation of the stationary states 8

1.4 Construction of solutions 8

1.5 Finite time blow-up in non-delayed networks (d=0) 11

1.6 Desynchronisation in the weakly connected case 13

1.7 The delayed NNLIF model and self-sustained oscillations 16

1.8 Rigorous derivation of the NNLIF model 20

1.9 Generalised solutions and continuation after blow-up 22

1.10 Asymptotic behaviour for general connectivity strength and delay 24

1.11 Variants of the NNLIF model 29

1.12 Numerical methods for NNLIF type systems 37

2 Fokker-Planck representation of a stochastic neural field and the modelling of grid cells 38

2.1 Modelling the hexagonal firing pattern of grid cells 38

2.2 Rigorous mean-field limit towards a new model 39

2.3 Mathematical setting and definitions 41

2.4 Analytical study of the stationary problem 42

2.5 Global-in-time existence of classical solutions 46

2.6 Asymptotic stability of stationary states 48

3 Rate models for decision making tasks 52

3.1 Analysis of the rate model 54

3.2 Reduction to a diffusion model 56

II Hyperbolic equations of transport type 58

4 The Time Elapsed model for neural networks 59

4.1 Neurons structured by time elapsed since the last discharge 59

4.2 Mathematical setting and the choice of parameters 60

4.3 Global existence and well-posedness 62

4.4 Study of the stationary equation 65

4.5 Asymptotic relaxation to a stationary state for the step function case 66

4.6 Relaxation to a stationary state for a smooth discharge function 70

4.7 Self-sustained oscillations 76

4.8 Rigorous derivation from microscopic models 81

4.9 Variants of the Time Elapsed model 82

5 Jump based Leaky Integrate & Fire model 95

5.1 Macroscopic description of deterministic IF neurons 95

5.2 Mathematical setting 96

5.3 Local and global well-posedness of the jump equation 97

5.4 Finite time blow-up in the strongly-excitatory regime 99

5.5 Long time behaviour 100

5.6 Heuristic link to the NNLIF model 101
<table><tr><td><b>6</b></td><td><b>Some other hyperbolic models</b></td><td><b>102</b></td></tr><tr><td></td><td>6.1 Leaky Integrate and Fire with random firing</td><td>102</td></tr><tr><td></td><td>6.2 Continuity equation for pulse coupled oscillators</td><td>103</td></tr></table>

### III Kinetic Fokker-Planck equations 104

<table><tr><td><b>7</b></td><td><b>The kinetic Voltage-Conductance model</b></td><td><b>104</b></td></tr><tr><td></td><td>7.1 Fluctuations at the level of the conductance</td><td>104</td></tr><tr><td></td><td>7.2 Notation for mathematical analysis</td><td>107</td></tr><tr><td></td><td>7.3 The linear Voltage-Conductance model</td><td>107</td></tr><tr><td></td><td>7.4 <i>A priori</i> estimates on solutions in the nonlinear case</td><td>111</td></tr><tr><td></td><td>7.5 Existence theory in the weak interaction regime</td><td>113</td></tr><tr><td></td><td>7.6 Numerical simulations of the Voltage-Conductance equation</td><td>115</td></tr><tr><td></td><td>7.7 The random firing voltage-conductance model and one variable reductions</td><td>115</td></tr><tr><td></td><td>7.8 Other variants of the voltage-conductance model</td><td>118</td></tr></table>

<table><tr><td><b>8</b></td><td><b>Kinetic FitzHugh-Nagumo models</b></td><td><b>123</b></td></tr><tr><td></td><td>8.1 From individual FitzHugh-Nagumo neurons to a mesoscopic system</td><td>123</td></tr><tr><td></td><td>8.2 Existence of weak solutions</td><td>124</td></tr><tr><td></td><td>8.3 Long time behaviour of the mesoscopic system</td><td>125</td></tr><tr><td></td><td>8.4 Periodic solutions arising through noise and interaction</td><td>127</td></tr><tr><td></td><td>8.5 A spatially extended kinetic FitzHugh-Nagumo model and its macroscopic limit</td><td>127</td></tr></table>

## Introduction

Since the first mathematical models for individual neurons [160, 137, 109, 194], a humongous research effort has been devoted to understanding complex phenomena in the brain and the modelling of large neural assemblies. Contrary to several older fields like theoretical physics, which can for example rely upon quantum field theory and general relativity, there is no abstract paradigm or first physics principles allowing to construct knowledge in neuroscience. Furthermore, animal brains exhibit an astonishing level of complexity and, despite recent advances, precise *in vivo* measurements of entire neural networks are still challenging.

For these reasons, researchers seek to develop simple mathematical models that represent the qualitative properties of neural tissues at the mesoscopic and macroscopic scales. They need to be mathematically and numerically tractable while still encompassing as many key features as possible. A popular approach is to start with a simple model for how a state variable *x* evolves for individual neurons, for example with stochastic differential equations (SDE), to choose interaction rules and to derive a mean-field equation for a network containing an infinite number of neurons [19, 18, 35, 202, 201, 204, 3, 106, 104]. In many cases, the probability density $\rho(x, t)$ of finding in the network, at time $t$, a neuron with state variable *x*, solves a nonlinear partial differential equation (PDE).

These nonlinear equations often present unique structures and non-local terms, making it impossible to straightforwardly apply the standard theories of PDEs to them. In the past two decades, a lot of progress has been made in the mathematical understanding of nonlinear partial differential equations modelling neuronal assemblies. Specific tools have been developed and many methods from mathematical physics have been adapted to their specific difficulties.

Although there are already excellent survey articles [13, 66, 65, 265] and books [122, 68] about mathematical modelling in neuroscience, they are either biology-oriented or focused on older types of equations like neural fields, Wilson-Cowan models [270, 271] and dynamical systems. We propose here a review article focused on the mathematical theory for neuroscience models in the form of nonlinear
PDEs, often rigorously obtained from particle systems, with an emphasis on rigorous proofs of the prop-
erties of the solutions. We have no intention, nor do we pretend, to be exhaustive. Our attempt is rather
to accurately describe the process of constructing methods for a few key types of models and to show
how they are related to each other through asymptotic reductions and approximations, while keeping in
mind the application purpose: getting back from proving theorems to challenging our understanding of
the brain.

The structure of this survey is organized around three large families of PDEs. In the first part, we discuss systems of parabolic type, with Fokker-Planck differential operators. This entails smoothness of solutions and allows to use tools such as Green kernels and entropy methods. We discuss the Nonlinear Noisy Leaky Integrate and Fire (NNLIF) model for neural networks, a Fokker-Planck system representing grid cells, and rate-based equations used to model decision making. In the second part, we focus on first-order hyperbolic type PDEs, whose solutions are not always smooth nor unique. It allows us to introduce other key methods like characteristics, the measure solutions framework, and Doeblin-Harris theorems. We mainly discuss the Time Elapsed neuron model and a jump-based Leaky Integrate and Fire model. In the last part, we focus on kinetic models with degenerate diffusion, for which the regularity theory is challenging and the long-time behaviour still poorly understood. We propose an overview of Voltage-Conductance and FitzHug-Nagumo models, both in the form of kinetic Fokker-Planck equations.

**Notation:** Most conventions will be explained when required, but we want to fix for clarity some choices that are used throughout this document. We denote $\mathbb{R}_+ = [0, +\infty)$ (note that $0 \in \mathbb{R}_+$). We write $C([a, b])$ or $C[a, b]$ for the space of continuous functions from the interval $[a, b]$ to $\mathbb{R}$; when the norm is not ambiguous, $C(X, Y)$ is the space of continuous functions from space $X$ to space $Y$, or just $C(X)$ if $Y = \mathbb{R}$. Similarly, $C^k$ will be used to denote spaces with higher regularity than just continuous. The Lebesgue spaces on a set $\Omega$ are denoted $L^p(\Omega)$ or just $L^p$, $p \in [1, +\infty]$, with the classical norms $\|\cdot\|_{L^p}$. Those norms are also written $\|\cdot\|_p$ when there is no ambiguity. Given a weight function $\omega : \Omega \to \mathbb{R}_+$, we denote $L^p(\Omega, \omega)$ or $L^p(\omega)$ the space of measurable functions $f$ with finite weighted $L^p$ norm $\|f\omega^{1/p}\|_{L^p}$. The Sobolev spaces on $\Omega$ are $W^{k,p}(\Omega)$; for $p = 2$, it is the standard Hilbert space $H^k$. Similar to Lebesgue spaces, we consider weighted versions $W^{k,p}(\Omega, \omega)$, $H^k(\Omega, \omega)$. When the notation $W$ is already taken by a connectivity kernel, we use $W^{k,p}(\Omega)$ for Sobolev spaces. Given an interval $I \subset \mathbb{R}$ and a Banach space $X$ we use the Bochner spaces $L^p(I, X)$, also denoted by $L^p(a, b; X)$ for $I = (a, b)$. The space of probability measures on a domain $\Omega \subset \mathbb{R}^d$ will be denoted by $P(\Omega)$ typically in dimension $d \in \{1, 2, 3\}$. We will write $f \equiv a$ to mean "the function $f$ is constant and equal to $a \in \mathbb{R}$".

# Part I

Parabolic equations of Fokker-Planck type

1 The Nonlinear Noisy Leaky Integrate & Fire model for neural networks

We focus here on the so-called Nonlinear Noisy Leaky Integrate & Fire model for neural networks,
or in short form NNLIF<fnref n="1" />. It has been introduced about twenty-five years ago, in the seminal works of
Brunel and Hakim [19, 13]. Their main motivation was to explain a counter-intuitive phenomenon:
self-sustained oscillations in networks of neurons with strong inhibitory connections and low individual
firing rate.

<footnote n="1">¹In some references, the model has no name, or the first N stands for "Network". The probabilistic community refers to it as a nonlinear McKean-Vlasov equation. The name NNLIF was first proposed in [27].</footnote>
Self-sustained oscillations are ubiquitous in neuroscience and play pivotal roles in many core functions like breathing and memory. They are found in many brain areas, like the visual cortex and the olfactory cortex. Previous studies had been focusing on networks where individual neurons are themselves oscillators, leading to total or partial synchronisation in a few clusters under the effect of excitatory connections. However, the role of inhibition in the emergence of periodic activity was still poorly understood. Even more striking was the fact that in networks of neurons that sporadically fire, oscillations at the population scale can be very fast: gamma frequency range (~30Hz) or even faster.

We start by presenting the historical derivation of the NNLIF model through physics arguments. We then set the framework for mathematical analysis of this type of PDE systems and, in order: discuss the existence and multiplicity of stationary states; explain how to prove existence and regularity of solutions through changes of variables and representation formulae; discuss the occurrence and meaning of finite-time blow-up in the non-delayed case; present the way entropy methods can be used to study long-time asymptotics in the weakly connected case; explain the specific methods and heuristic arguments allowing us to understand the NNLIF model in the delayed case. After all these considerations on the standard theory around this type of equations, we propose a brief overview of the rigorous derivation through a mean-field limit, and we present recent methods developed on the one hand to extend solutions after blow-up via time dilations and on the other hand to study the asymptotic behaviour beyond the weakly-connected case via semigroup methods and a sequence of pseudo-equilibria. Last, we present the numerous variants of the standard model and discuss some numerical methods tailored to this specific problem.

## 1.1 Derivation of the model à la physicienne

Let us propose first a physical derivation of the model as can be found in early references [19, 27]. The rigorous mean-field limit from a particle system was established later in 2015, in [82]. The main idea is to consider a large network where a representative neuron of electric potential $V(t)$ is modelled by the Laplacian Integrate & Fire model [160] given by

$$C_m \frac{dV}{dt} = -g_L(V - V_L) + I(t), \quad (1.1)$$

where $C_m$ is the capacitance of the membrane, $g_L$ the leak conductance and $V_L$ the leak potential. The characteristic time of relaxation $\tau_m = \frac{g_L}{C_m}$ is of the order of 2ms and $V_L \simeq -70mV$. The incoming synaptic current is a stochastic process of the form

$$I(t) = J_E \sum_{i=1}^{C_E} \sum_{j \in \mathbb{N}} \delta(t - t_{E_j}^i - d) + J_I \sum_{i=1}^{C_I} \sum_{j \in \mathbb{N}} \delta(t - t_{I_j}^i - d), \quad (1.2)$$

where $\delta$ is the Dirac mass measure at zero, $J_E$ and $J_I$ are the strength of excitatory and inhibitory synapses, $C_E$ and $C_I$ are the number of excitatory and inhibitory connections received on average, and $t_{E_j}^i, t_{I_j}^i$ are the times of the $j^{th}$ action potential emitted by the $i^{th}$ pre-synaptic excitatory or inhibitory neuron; those times are subject to random fluctuations. The parameter $d$ represents the average time it takes for an action potential to go from a neuron to another and be integrated; as we will see, this parameter can have a crucial effect on the solutions and the model is more singular when $d = 0$.

When a neuron reaches the firing potential $V_F \simeq -50mV$, it emits an action potential (it is also said to fire), and then comes back to the reset potential $V_R \simeq -60mV$. Assuming that each neuron in the network fires independently following a Poisson law of intensity $\nu(t)$, the mean and variance of $I(t)$ are $\mu_C(t - d) = b\nu(t - d)$, where

$$b = C_E J_E - C_I J_I,$$

is the average connectivity of the network, and $\sigma_C^2(t-d) = (C_E J_E^2 + C_I J_I^2)\nu(t-d)$. When $b < 0$ we say that the network is average inhibitory, when $b > 0$ we describe it as average excitatory. When $b = 0$, the network is disconnected and only mathematicians care about it, although it is useful as linearized equation for small connectivity.
The stochastic process (1.2) being difficult to study, a classical approach [18, 19, 226, 175, 248, 202] is to replace it by an Ornstein-Uhlenbeck process with same mean and variance

$$I(t)dt \approx \mu_c(t-d)dt + \sigma_c(t-d)dB_t,$$

where $B_t$ is a standard Brownian motion.

If we rescale the voltage and time variables so to have $C_m = g_L = 1$, we can write the following stochastic differential equation for a neuron of this network:

$$dV = (-V + V_L + \mu_C(t-d)) dt + \sigma_C(t-d) dB_t, \quad V < V_F \quad (1.3)$$

with the jump process $\limsup_{t \to t_0^-} V(t) = V_F \implies \liminf_{t \to t_0^+} V(t) = V_R$, and the discharge intensity $\nu(t) = \nu_{ext} + N(t)$, where $N(t)$ is the flux of neurons in the network that cross the firing potential $V_F$ (and emit action potentials) at time $t$. The term $\nu_{ext} \in \mathbb{R}$ represents an external input from outside the network. We consider it constant here, but it may vary in time in more realistic settings. Although it plays an important role in the modelling, as an interface between a group of neurons, other networks and sensory inputs, we will often rescale the voltage $V$ in order to absorb it, which simplifies the mathematical analysis.

The stochastic differential equation (1.3) is associated to a Fokker-Planck partial differential equation, describing the probability density $(v, t) \mapsto \rho(v, t) = P(V(t) \in dv)$ of finding in the network, at time $t \ge 0$, a neuron with membrane potential $v \in (-\infty, V_F]$, given by

$$\frac{\partial \rho}{\partial t}(v,t) + \frac{\partial}{\partial v} \left[ (V_L - v + b_0 + bN(t-d))\rho(v,t) \right] - a(N(t-d)) \frac{\partial^2 \rho}{\partial v^2}(v,t) = \delta_{v=V_R} N(t), \quad (1.4)$$

with $b_0 = b\nu_{ext}$ and $a(N(t-d)) = \frac{\sigma_C^2(t)}{2} = a_0 + a_1 N(t-d)$, $\rho(V_F, t) = \rho(-\infty, t) = 0$ and

$$N(t) = -a(N(t-d)) \frac{\partial \rho}{\partial v}(V_F, t). \quad (1.5)$$

Since $\rho$ must remain a probability density, integrating (1.4) in $v$ yields the relation (1.5). The singular resetting of neurons is embodied by the Dirac mass at potential $V_R$.

## 1.2 Mathematical setting for the NNLIF type models

Up to appropriate variables re-scaling, the Fokker-Planck equation (1.4)–(1.5) can be rewritten as the Cauchy problem

$$\left\{
\begin{array}{ll}
\displaystyle \frac{\partial \rho}{\partial t} + \frac{\partial}{\partial v} \left[ (-v + bN(t-d))\rho \right] - a(N(t-d)) \frac{\partial^2 \rho}{\partial v^2} = \delta_{v=V_R} N(t), & v \in (-\infty, V_F], \\
\displaystyle N(t) = -a(N(t-d)) \frac{\partial \rho}{\partial v}(V_F, t), & \rho(V_F, t) = 0, \quad \rho(-\infty, t) = 0, \quad t > 0, \\
\displaystyle \rho(v, 0) = \rho^0(v) \ge 0, & \int_{-\infty}^{V_F} \rho^0(v) dv = 1, \\
\displaystyle N(t) = N^0(t), & v \in (-\infty, V_F], \\
\end{array}
\right.
\quad \text{(NNLIF)}$$

We call this model the standard NNLIF model, with ($d > 0$) or without ($d = 0$) synaptic delay. When $d = 0$, the initial condition $N^0$ for the firing rate is not required.

Note that the external input $\nu_{ext}$ has been absorbed in the values $V_F$ and $V_R$ through an affine voltage rescaling in order to reduce the number of parameters and simplify the mathematical description. Every time a result is stated in terms of $V_R, V_F, b$ and the function $a$, the role of external inputs $\nu_{ext}$ is hidden in the values of these parameters. The impact of $\nu_{ext}$ is especially visible when taking negative values for $V_F$, which is equivalent to a highly excitatory external input. Note also that we allow for more
general choices for the function a than the form $a(N) = a_0 + a_1 N$ obtained in the previous section. The rigorous derivation from a particle system yields $a_1 = 0$, in contrast to the physical one.

To better understand the role of the Dirac mass $\delta_{v=V_R}$ in the right-hand side, there are two viewpoints:

* (NNLIF) is satisfied in the classical sense on $(-\infty, V_R) \cup (V_R, V_F]$ and in the sense of distributions on $(-\infty, V_F]$;

* (NNLIF) is satisfied in the classical sense on $(-\infty, V_R) \cup (V_R, V_F]$ and the function $\partial_v \rho(\cdot, t)$ has left and right limits at $V_R$ with a jump of size $\frac{N(t)}{a(N(t-d))}$, i.e.

$$ \lim_{v \to V_R^-} \left( \frac{\partial \rho}{\partial v}(v, t) \right) - \lim_{v \to V_R^+} \left( \frac{\partial \rho}{\partial v}(v, t) \right) = \frac{N(t)}{a(N(t-d))}. $$

In order to avoid technical difficulties, we will consider fast-decaying (at $-\infty$) initial conditions and solutions, since in neuroscience applications voltage distributions are not fat-tailed.

**Definition 1.1.** We say that $(\rho, N)$ is a strong fast-decaying solution, of system (NNLIF) on $[0, T^*)$, $T^* \in (0, +\infty]$ if

* $\rho \in C^0((-\infty, V_F] \times [0, T^*)) \cap C^{2,1}(((-\infty, V_R) \cup (V_R, V_F]) \times (0, T^*))$ and $N \in C^0([-d, T^*));$

* Functions $\rho$ and $N$ are non-negative, satisfy (NNLIF) in the classical sense on $(-\infty, V_R) \cup (V_R, V_F]$ and in the sense of distributions in $(-\infty, V_F]$;

* $\forall t \in [0, T^*), \forall Q \in \mathbb{R}[X], \lim_{v \to -\infty} Q(v)\rho(v, t) = 0$ and $\lim_{v \to -\infty} Q(v) \frac{\partial \rho}{\partial v}(v, t) = 0$.

We can also define weak solutions.

**Definition 1.2.** Let $T^* \in (0, +\infty]$, let $\rho \in L^\infty([0, T^*], L_+^1((-\infty, V_F)))$ and $N \in L_{\text{loc},+}^1([0, T^*])$. The pair $(p, N)$ is said to be a fast-decaying weak solution of (NNLIF) if for every test function $\phi \in C^2((-\infty, V_F])$ with polynomial growth<fnref n="2" /> we have

$$ \frac{d}{dt} \int_{-\infty}^{V_F} \rho(v, t) \phi(v) dv = \int_{-\infty}^{V_F} \left[ (-v + bN(t-d)) \frac{\partial \phi}{\partial v} + a(N(t-d)) \frac{\partial^2 \phi}{\partial v^2} \right] \rho(v, t) dv \\ + N(t) (\phi(V_R) - \phi(V_F)), $$

where all the quantities are finite, the time derivative is taken in the sense of distributions, $N(t-d)$ is given by $N^0$ when $t \le d$ and $\lim_{t \to 0} \|\rho(\cdot, t) - \rho^0(\cdot)\|_{L^1} = 0$.

All strong solutions are weak solutions, which can be checked by integration by parts.

**Remark 1.3.** If we apply the definition of weak solution with the test function $\phi = 1$, we obtain for $t > 0$,

$$ \int_{-\infty}^{V_F} \rho(v, t) dv = \int_{-\infty}^{V_F} \rho^0(v) dv = 1. $$

We will refer to the following regularity assumptions on initial data.

**Definition 1.4 (Initial datum for classical solutions).** We say that $\rho^0$ is a compatible initial condition for (NNLIF) if $\rho^0 \in C^0((-\infty, V_F]) \cap C^1((-\infty, V_R) \cup (V_R, V_F]) \cap \mathcal{P}((-\infty, V_F])$ is fast-decaying at $-\infty$, $\rho^0(V_F) = 0$, $\frac{d\rho^0}{dv}$ admits finite left and right limits at $V_R$ and is fast decreasing at $-\infty$. In the case with delay $d > 0$, we require that $N^0 \in C([-d, 0])$ and $N^0(0) = -a(N(-d)) \frac{d\rho^0}{dv}(V_F)$.

<footnote n="2">²There exists $Q \in \mathbb{R}[X]$ such that $\forall v \in (-\infty, V_F]$, $\phi(v) \le Q(v)$.</footnote>
## 1.3 Characterisation of the stationary states

We call a stationary state, or steady state, of (NNLIF) a solution $(\rho_{\infty}, N_{\infty})$ of the system

$$ 
\left\{ 
\begin{array}{l} 
\displaystyle \frac{\partial}{\partial v} [(-v + bN_{\infty})\rho_{\infty}] - a(N_{\infty}) \frac{\partial^2 \rho_{\infty}}{\partial v^2} = \delta_{v=v_R} N_{\infty}, \quad v \in (-\infty, V_F], \\ \\ 
\displaystyle N_{\infty} = -a(N_{\infty}) \frac{\partial \rho_{\infty}}{\partial v}(V_F), \\ \\ 
\displaystyle \rho_{\infty}(V_F) = \rho_{\infty}(-\infty) = 0, \quad \rho_{\infty} \ge 0, \quad \int_{-\infty}^{V_F} \rho_{\infty}(v)dv = 1. 
\end{array} 
\right. 
$$

They are of course the same with and without delay. Direct integration and rescalings allow to characterise the stationary states in a semi-explicit way.

**Theorem 1.5** (Cáceres, Carrillo, Perthame [27]). Let $(\rho_{\infty}, N_{\infty})$ be a stationary state of (NNLIF). It must be of the form

$$ \rho_{\infty}(v) = \frac{N_{\infty}}{a(N_{\infty})} e^{-\frac{(v-bN_{\infty})^2}{2a(N_{\infty})}} \int_{\max(v,V_R)}^{V_F} e^{\frac{(w-bN_{\infty})^2}{2a(N_{\infty})}} dw. \quad (1.6) $$

where $N_{\infty}$ satisfies the fixed-point equation

$$ \frac{1}{N_{\infty}} = I(N_{\infty}), \quad I(N) = \int_{0}^{+\infty} \frac{e^{-\frac{s^2}{2}}}{s} \left( e^{\frac{sV_F - bN}{\sqrt{a(N)}}} - e^{\frac{sV_R - bN}{\sqrt{a(N)}}} \right) ds. \quad (1.7) $$

Notice that given the stationary firing rate $N_{\infty}$, the stationary density is explicitly given. The existence and number of classical stationary states thus depends solely on the solutions of the fixed-point equation (1.7). A couple pages of real analysis lead to a good picture of what can happen.

**Theorem 1.6** (Cáceres, Carrillo, Perthame, [27]). Assume $a(N) \equiv a > 0$; then, for (NNLIF),

* there is a constant $c(a) > 0$ such that if $b \in (-\infty, c(a)]$, there is a unique stationary state;

* if $0 < b < V_F - V_R$ or $0 < 2ab < (V_F - V_R)^2 V_R$, there is at least one stationary state;

* if $b > V_F - V_R$ and $0 < 2ab < (V_F - V_R)^2 V_R$, there are at least two stationary states;

* if $b > \max(2(V_F - V_R), 2V_F I(0))$, there is no stationary state.

When the diffusion coefficient linearly depends on $N$, the classification remains similar, but without the uniqueness for inhibitory or small connectivity.

**Theorem 1.7** (Cáceres, Carrillo, Perthame [27]). Assume $a(N) = a_0 + a_1 N$, with $a_0, a_1 > 0$; then,

* If $b < V_F - V_R$ or $(2a_0b + 2a_1V_R < (V_F - V_R)^2 V_R$ and $b > 0$), there is at least one steady state;

* If $b > V_F - V_R$ and $2a_0b + 2a_1V_R < (V_F - V_R)^2 V_R$, there are at least two stationary states;

* if $b > \max(2(V_F - V_R), 2V_F I(0))$, there is no stationary state.

## 1.4 Construction of solutions

### 1.4.1 Local-in-time existence via a Stefan-like free boundary problem

Because it was the basis for the construction of solutions in other parameter settings or variants of the NNLIF model, we will present here the main ideas of the first deterministic existence proof, as can be found in [48]. It concerns only the case $d = 0$, $a(N) \equiv a > 0$ constant, but it was generalised later to the case with delay $d > 0$ in [31], to the case with internal noise $a(N) = a_0 + a_1 N$ in [91] and
to other PDEs from neuroscience, see, for example (GC) below. Proofs via stochastic methods were independently achieved in [81, 80] for classical solutions and in [82] for generalised solutions.

Here we will present the main ideas of the deterministic proof for the case $d = 0$. The case with delay proceeds similarly. The idea is to transform the original problem on a fixed domain $(-\infty, V_F]$ into a free-boundary Stefan-like problem, then isolate the firing rate in a closed equation and apply the Banach fixed point theorem.

First, let us rescale the variables, in order to have $a = 1$ and $V_F = 0$, which comes at the price of changing $b$ and putting back an external input $\nu_{\text{ext}}$. It is enough to apply the scaling $\rho_{\text{new}}(v, t) = \sqrt{a}\rho_{\text{old}}(\sqrt{a}v + V_F, t)$, $N_{\text{new}}(t) = -\partial_v\rho_{\text{new}}(0, t)$. If we define $\nu_{\text{ext}} = \frac{V_F}{\sqrt{a}}$ and the new parameters $\bar{b} = \frac{b}{\sqrt{a}}$, $\bar{V}_R = \frac{V_R - V_F}{\sqrt{a}}$, then dropping the bars, the "new" pair $(\rho, N)$ satisfies,

$$ \frac{\partial \rho}{\partial t} + \frac{\partial}{\partial v} \left[ (-v + \nu_{\text{ext}} + bN(t))\rho \right] - \frac{\partial^2 \rho}{\partial v^2} = \delta_{v=V_R} N(t), \quad v \in (-\infty, 0], t > 0. $$

Note that this preserves the sign of $b$ and compensates the change $V_F \to 0$ of the firing potential by an equivalent external force $\nu_{\text{ext}}$. This external force is not the same as the one derived in the first Subsection, as $V_F \neq 0$ in real life; it is here purely mathematical.

Secondly, we apply a classical change of variables that transforms the linear Fokker-Planck equation into a heat equation [55],

$$ y = e^t v, \quad \tau = \frac{1}{2}(e^{2t} - 1), \quad w(y, \tau) = \alpha(\tau)\rho(y\alpha(\tau), -\log \alpha(\tau)), $$

with the notation $\alpha(\tau) = \frac{1}{\sqrt{2\tau + 1}} = e^{-t}$. The firing rate in the new coordinates is $M(\tau) = \alpha(\tau)^2 N(t)$.

Now, we apply the second change of variables

$$ x = y - \int_{0}^{\tau} \left( \nu_{\text{ext}} \alpha(s) + b \frac{M(s)}{\alpha(s)} \right) ds = y - \nu_{\text{ext}} (\sqrt{1+2\tau} - 1) - b \int_{0}^{\tau} \frac{M(s)}{\alpha(s)} ds, $$

and we denote $u(x, \tau) = w(y, \tau)$. The pair $(u, M)$ is a solution to the following Stefan-like free boundary problem with a moving distribution in the right-hand side given by

$$ \left\{ \begin{array}{ll} \displaystyle \frac{\partial u}{\partial \tau}(x, \tau) = \frac{\partial^2 u}{\partial y^2}(x, \tau) + M(\tau)\delta_{s(\tau)+\frac{V_R}{\alpha(\tau)}}(x) & x \in (-\infty, s(\tau)], \tau > 0 \\ \displaystyle M(\tau) = -\frac{\partial u}{\partial x}(s(\tau), \tau) & \tau > 0 \\ \displaystyle s(\tau) = s(0) - \nu_{\text{ext}}(\sqrt{1+2\tau}-1) - b \int_0^\tau \frac{M(s)}{\alpha(s)} ds & \tau > 0 \\ \displaystyle u(-\infty, \tau) = u(s(\tau), \tau) = 0, \quad u(x, 0) = u^0(x) & x \in (-\infty, s(0)], \tau > 0. \end{array} \right. \quad (1.8) $$

Notice that

1. the mass is conserved: for all $\tau$, $\int_{-\infty}^{s(\tau)} u(x, \tau)dx = \int_{-\infty}^{0} u^0(x)dx$;

2. the flux through the point $s_1(\tau) := s(\tau) + \frac{V_R}{\alpha(\tau)}$ is exactly

$$ M(\tau) := -\partial_x u(s(\tau), \tau) = \partial_x u(s_1(\tau)^-, \tau) - \partial_x u(s_1(\tau)^+, \tau). $$

Away from the singular points $s(\tau)$ and $s_1(\tau) = s(\tau) + \frac{V_R}{\alpha(\tau)}$, the problem is a linear heat equation. Hence, given the heat kernel and its space derivative

$$ G(x, \tau, \xi, \eta) = \frac{1}{\sqrt{4\pi(\tau - \eta)}} e^{-\frac{(x-\xi)^2}{4(\tau-\eta)}}, \quad \frac{\partial G}{\partial \xi}(x, \tau, \xi, \eta) = \frac{x-\xi}{4\sqrt{\pi}(\tau-\eta)^{3/2}} e^{-\frac{(x-\xi)^2}{4(\tau-\eta)}}, $$
we have the Green identity

$$ \frac{\partial}{\partial \xi} \left( G \frac{\partial u}{\partial \xi} - u \frac{\partial G}{\partial \xi} \right) - \frac{\partial}{\partial \eta} (Gu) = 0. $$

Fix any $\tau \in \mathbb{R}_+$ and any $x \in (-\infty, s_1(\tau)) \cup (s_1(\tau), s(\tau))$. Integrating the identity and using the regularity and decay of both the solution and the heat kernel $G$ yields

$$ u(x, \tau) = \int_{-\infty}^{0} G(x, \tau, \xi, 0) u^{0}(\xi) d\xi - \int_{0}^{\tau} M(\eta) G(x, \tau, s(\eta), \eta) d\eta \\ + \int_{0}^{\tau} M(\eta) G(x, \tau, s_1(\eta), \eta) d\eta. \quad (1.9) $$

The first term is the solution to the linear heat equation on the whole real line with initial condition $u^0$; the other two terms represent the loss of neurons at the free boundary and their return at a moving reset point $s_1(\tau)$. Now that we have $u$ in terms of only $u^0$ and the firing rate $M$, it is possible to isolate $M$ by taking the partial derivative along $x$ and then the limit $x \to s(\tau)$. It leads to the close formula

$$ M(\tau) = -2 \int_{-\infty}^{0} G(s(\tau), \tau, \xi, 0) \frac{du^0}{dx}(\xi) d\xi + 2 \int_{0}^{\tau} M(\eta) \frac{\partial}{\partial x} G(s(\tau), \tau, s(\eta), \eta) d\eta \\ - 2 \int_{0}^{\tau} M(\eta) \frac{\partial}{\partial x} G(s(\tau), \tau, s(\eta) + \frac{V_R}{\alpha(\eta)}, \eta) d\eta. \quad (1.10) $$

From there, it is possible to apply the Banach fixed point theorem in $\{M \in C([0, \tau_0]) \mid \|M\|_{\infty} \le m\}$ with $\tau_0$ small enough and $m$ large enough. For

$$ \tau_0 < \left\| \frac{du_0}{dx} \right\|_{L^\infty}^{-1}, $$

there exists a unique continuous local-in-time solution $M$ on $[0, \tau_0]$ and the maximal time of existence is $\tau^* = \sup\{\tau \in \mathbb{R}_+ | M(\tau) < +\infty\}$. The existence of the firing rate implies, through the Duhamel formula (1.9), the existence of a solution $(u, M)$ to the Stefan-like problem (1.8). Undoing all changes of variables yields the following result.

**Theorem 1.8** (Carrillo, Gonzalez, Gualdani, Schonbeck [48]). Assume $d = 0$, $a(N) = a > 0$ is constant, and $\rho^0$ is a compatible initial condition. Then there exists a unique strong solution to (NNLIF) whose maximal time of existence $T^*$ satisfies

$$ T^* = \sup\{\tau \in \mathbb{R}_+ | N(t) < +\infty\}. \quad (1.11) $$

### 1.4.2 Global-in-time existence and uniform bounds in the inhibitory case

Consider again the case $d = 0$, $a(N) \equiv a > 0$ constant. Then, in the case $b < 0$ of an inhibitory network, strong solutions are in fact global in time. Two different methods were used to prove it.

First, in [48], global-in-time existence for inhibitory networks is proved using the Stefan problem reformulation. The main idea is that the nonlinear term in the free boundary motion is monotonically increasing, yielding the key estimate

$$ s(\tau) - s_1(\eta) \ge |V_R| - |\nu_{\text{ext}}|(\tau - \eta), \quad \eta < \tau. $$

Then, all terms in the fixed-point formulation (1.10) can be controlled uniformly on a short enough time and a contradiction argument with (1.11) proves global-in-time existence.

Then, in [52], the same result was achieved and improved by using universal super-solutions. The idea is to use the fact that in the inhibitory case, non-increasing super-solutions of the linear problem are super-solutions for the nonlinear problem and can control the firing rate.
**Definition 1.9** (universal super-solutions). Assume $d = 0$, $a(N) = a > 0$ and $b < 0$. Let $V_0 \in (-\infty, V_F)$ and $T > 0$. $(\bar{\rho}, \bar{N})$ is a universal super-solution to (NNLIF) on $[V_0, V_F] \times [0, T]$ if

$$ \frac{\partial \bar{\rho}}{\partial t} - \frac{\partial}{\partial v}[v\bar{\rho}] - a \frac{\partial^2 \bar{\rho}}{\partial v^2} \geq \delta_{v=V_R} \bar{N}(t) $$

in the distributional sense on $[V_0, V_F] \times [0, T]$, in the classical sense away from $V_R$, if $\bar{N}(t) = -a\partial_v\bar{\rho}(V_F, t)$ and if for all $t$, $\bar{\rho}(\cdot, t)$ is non-increasing on $[V_0, V_F]$.

This definition allows the following comparison principle.

**Lemma 1.10.** Assume $d = 0$, $a(N) = a > 0$ and $b < 0$. Let $V_0 \in (-\infty, V_F)$ and $T > 0$. Let $(\rho, N)$ be a classical solution of (NNLIF) defined up to time $T$ and $(\bar{\rho}, \bar{N})$ be a universal super-solution on $[V_0, V_F] \times [0, T]$. Assume that for all $(v, t) \in [V_0, V_F] \times [0, T]$,

$$ \bar{\rho}(v, 0) \geq \rho^0(v), \quad \bar{\rho}(V_0, t) \geq \rho(V_0, t) $$

Then $\bar{\rho} \geq \rho$ on $[V_0, V_F] \times [0, T]$ and if $\bar{\rho}(\cdot, 0) - \rho^0$ is not identically zero, then $\bar{\rho} > \rho$ on $(V_0, V_F) \times (0, T]$.

Using two specific families of super-solutions constructed at hand, and taking advantage of the Stefan reformulation and classical heat equation estimates, [52] proves

**Theorem 1.11** (Carrillo, Perthame, Salort, Smets [52]; generalising [48]). Assume $d = 0$, $a(N) = a > 0$ is constant and $b < 0$. Let $\rho^0$ be a compatible initial condition. Then there exists a unique global-in-time solution $(\rho, N)$ to (NNLIF) and

$$ \rho \in L^\infty((-\infty, V_F] \times [0, +\infty)), \quad N \in L^\infty(0, +\infty). $$

The result can be improved to $\rho^0 \in L_+^1 \cap L^\infty(-\infty, V_F)$ with

$$ \limsup_{v \to V_F} \frac{\rho^0}{V_F - v} < +\infty $$

if there exists a strong solution. The influence of less stringent decays near the boundary has been studied in the context of finance models in [130, 83].

Note the crucial assumption $a(N) = a > 0$ constant; as we shall discuss, in the case with linear internal noise $a(N) = a_0 + a_1 N$, $a_0, a_1 > 0$, finite time blow-up occurs even in the inhibitory case.

## 1.5 Finite time blow-up in non-delayed networks (d=0)

When there is no synaptic delay ($d = 0$), the NNLIF model is prone to finite time blow-up. As we can see in Theorem 1.8, finite-time blow-up is caused by the divergence of the firing rate $N(t)$ at the blow-up time $T^*$:

$$ \limsup_{t \to T^*} N(t) = +\infty. $$

From a biology perspective, it is a rough depiction of a macroscopic part of the neural network synchronising and firing at the same time.

A first proof of this phenomenon was provided in [27], using a contradiction argument on the exponential moment

$$ M_{\mu}(t) = \int_{-\infty}^{V_F} e^{\mu v} \rho(v, t) dv. $$

Indeed, for $\mu$ large enough, it satisfies the differential inequality

$$ \frac{d}{dt} M_{\mu}(t) \geq b\mu N(t) \left[ M_{\mu}(t) - \frac{e^{\mu V_F} - e^{\mu V_R}}{b\mu} \right], $$

which allows to prove that if the right-hand side is initially positive, then $M_\mu$ is increasing and, using another lower-bound, that $M_\mu$ is unbounded, contradicting the natural bound $M_\mu(t) \leq e^{\mu V_F}$.
**Theorem 1.12** (Cáceres, Carrillo, Perthame [27]). Assume $d = 0$, $a(N) \ge a_0 > 0$ and $b > 0$. Let $\mu > \max(\frac{V_F}{a_0}, \frac{1}{b})$. If

$$b\mu \int_{-\infty}^{V_F} e^{\mu v} \rho^0(v) dv > e^{\mu V_F} - e^{\mu V_R}, \qquad (1.12)$$

then there is no global-in-time weak solution of (NNLIF).

The condition (1.12) is satisfied, for example, whenever $\rho^0$ is concentrated close enough to the firing threshold: if $\rho^0$ is fully supported in $[V_F - \varepsilon, V_F]$ with $\varepsilon$ small enough, then $b\mu \ge 1$ and

$$\int_{-\infty}^{V_F} e^{\mu v} \rho^0(v) dv = \int_{V_F - \varepsilon}^{V_F} e^{\mu v} \rho^0(v) dv \\ \ge e^{-\mu \varepsilon} e^{\mu V_F} \int_{V_F - \varepsilon}^{V_F} \rho^0(v) dv = e^{-\mu \varepsilon} e^{\mu V_F} \int_{-\infty}^{V_F} \rho^0(v) dv = e^{-\mu \varepsilon} e^{\mu V_F} \\ > e^{\mu V_F} - e^{\mu V_R}.$$

The smaller $b > 0$, the higher $\mu$, and then the smaller $\varepsilon$ we need. For a low excitatory connectivity $b > 0$, blow-up occurs when the initial density is very concentrated near the firing threshold. In contrast, the stronger the connectivity, the easier it is for the network to synchronise.

In fact, it was proved in [233] that for a strong enough excitatory connectivity, all initial conditions lead to finite-time blow-up.

**Theorem 1.13** (Roux, Salort [233]).

* If $V_F \le 0$, and $b \ge V_F - V_R$, no fast-decaying weak solution of (NNLIF) can be global-in-time.

* If $V_F > 0$, there exists a value $b_c > 0$ such that if $b \ge b_c$, no fast-decaying weak solution of (NNLIF) can be global-in-time.

The main idea is to study the linear moment

$$M_1(t) = \int_{-\infty}^{V_F} (V_F - v) \rho(v, t) dv$$

and to prove that its evolution in time leads to enough concentration of the density to use the contra-diction method on the exponential moment. Note that the case $V_F \le 0$ implies the existence of an excitatory external input $\nu_{\text{ext}}$ which is strong enough to cause neurons to fire even in the absence of noise. It explains why finite-time blow-up, an avatar of synchronisation, occurs more easily in this parameter range.

As seen in Theorem 1.11, constant diffusion $a(N) \equiv a > 0$ and an inhibitory connectivity $b < 0$ prevent blow-up from occurring. However, when there is linear internal noise in the form $a(N) = a_0 + a_1 N$, $a_0, a_1 > 0$, there can be finite time blow-up even in the inhibitory case.

**Theorem 1.14** (Carrillo, Perthame, Salort, Smets [52]). Assume $d = 0$, $a(N) = a_0 + a_1 N$, $a_0, a_1 > 0$. Let $\mu > 0$ be such that $a_1\mu^2 + \mu b > 1$ and $a_0\mu^2 - V_F > 0$. If

$$\int_{-\infty}^{V_F} e^{\mu v} \rho^0(v) dv > e^{\mu V_F} - e^{\mu V_R},$$

then there is no global-in-time weak solution of (NNLIF).

The reason for this discrepancy between the cases $a_1 = 0$ and $a_1 \ne 0$ is that in the latter $N(t)$ is defined by the self-consistency equation

$$N(t) = -(a_0 + a_1 N(t)) \frac{\partial \rho}{\partial v}(V_F, t),$$
which can be rewritten as

$$N(t) = \frac{-a_0 \partial_v \rho(V_F, t)}{1 + a_1 \partial_v \rho(V_F, t)}.$$

Because $N(t)$ must remain finite and non-negative for classical solutions, internal noise adds the well-posedness constraint $a_1 \partial_v \rho(V_F, t) < 1$.

As per the biological interpretation of finite time blow-up as a partial or total synchronisation event, the solutions should not stop at the first blow-up time. In principle, it should be possible to compute which part of the mass blows-up in the density and to reset it at $V_R$ as a proportionate delta Dirac mass. From the mathematical analysis perspective, this problem is challenging and, although some progress have been made, a lot of questions are still open. In what follows, we will discuss some existing probabilistic [82, 236, 252] and deterministic [91, 90, 46] methods. Apart from conjectures [130] and results for related models (see Theorem 1.50 below), the blow-up asymptotics, namely the divergence speed of $N(t)$ close to $T^*$, are poorly understood (see however [83, 136] for the Stefan formulation).

## 1.6 Desynchronisation in the weakly connected case

We have seen that in some cases the NNLIF model synchronises through finite-time blow-up. An opposite phenomenon can occur: convergence of the solution towards a stationary state. A stationary density $\rho_\infty(v)$ does not mean that nothing is happening in the network. At the microscopic level, many neurons fire continuously, but at the global scale, the distribution is not evolving. This is akin to the convergence of a gas towards a state of higher entropy. As we shall see, the relative entropy method that works well for the linear Fokker-Planck equation on the whole space (see the classical review [172]) can be adapted, with enough work, to the weakly nonlinear regime of the NNLIF model, as was the case for many biological models in PDE form [178, 179, 208]. In the case of the NNLIF model, specific difficulties due to the non-local structure require appropriate control of the firing rate $N(t)$.

In all this subsection, we choose $d = 0$ and $a(N) \equiv a > 0$ constant.

### 1.6.1 A priori convergence in relative entropy

**Linear case** $b = 0$. Let us first showcase the method on the linear case $b = 0$ where neurons are disconnected from each other. Theorem 1.6 ensures that there exists a unique stationary state $(\rho_\infty, N_\infty)$. For any convex and $C^2$ function $G: \mathbb{R} \to \mathbb{R}_+$, we have the entropy dissipation

$$\frac{d}{dt} \int_{-\infty}^{V_F} G\left(\frac{\rho}{\rho_\infty}\right) \rho_\infty dv = -a \int_{-\infty}^{V_F} \left[\frac{\partial}{\partial v}\left(\frac{\rho}{\rho_\infty}\right)\right]^2 G''\left(\frac{\rho}{\rho_\infty}\right) \rho_\infty dv - N_\infty \left[G\left(\frac{N}{N_\infty}\right) - G\left(\frac{\rho}{\rho_\infty}\right) - \left(\frac{N}{N_\infty} - \frac{\rho}{\rho_\infty}\right)G'\left(\frac{\rho}{\rho_\infty}\right)\right] \Big|_{v=V_R}.$$

Note that by convexity of $G$, the second part of the right-hand side is always non-positive. We can choose $G(x) = (x - 1)^2$ and apply a Poincaré-like inequality to the first term, which yields

$$\frac{d}{dt} \int_{-\infty}^{V_F} \left( \frac{\rho(v, t) - \rho_\infty(v)}{\rho_\infty(v)} \right)^2 \rho_\infty(v) dv \le -2a\nu \int_{-\infty}^{V_F} \left( \frac{\rho(v, t) - \rho_\infty(v)}{\rho_\infty(v)} \right)^2 \rho_\infty(v) dv,$$

where $\nu$ is the Poincaré constant. By Grönwall's lemma,

$$\int_{-\infty}^{V_F} \left( \frac{\rho(v, t) - \rho_\infty(v)}{\rho_\infty(v)} \right)^2 \rho_\infty(v) dv \le e^{-2a\nu t} \int_{-\infty}^{V_F} \left( \frac{\rho^0(v) - \rho_\infty(v)}{\rho_\infty(v)} \right)^2 \rho_\infty(v) dv,$$

which is the convergence in relative entropy of $\rho(\cdot, t)$ towards $\rho_\infty$ at exponential speed.
**Weakly nonlinear case.** Proceeding like in the linear case, for any convex and $C^2$ function $G$ and a solution with appropriate decay, we compute the entropy dissipation

$$ \frac{d}{dt} \int_{-\infty}^{V_F} G\left(\frac{\rho}{\rho_\infty}\right) \rho_\infty dv = -a \int_{-\infty}^{V_F} \left[\frac{\partial}{\partial v}\left(\frac{\rho}{\rho_\infty}\right)\right]^2 G''\left(\frac{\rho}{\rho_\infty}\right) \rho_\infty dv \\ -N_\infty\left[G\left(\frac{N}{N_\infty}\right) - G\left(\frac{\rho}{\rho_\infty}\right) - \left(\frac{N}{N_\infty} - \frac{\rho}{\rho_\infty}\right)G'\left(\frac{\rho}{\rho_\infty}\right)\right] \Big|_{v=V_R} \\ +b(N(t) - N_\infty) \int_{-\infty}^{V_F} \frac{\partial \rho_\infty}{\partial v} \left[G\left(\frac{\rho}{\rho_\infty}\right) - \frac{\rho}{\rho_\infty}G'\left(\frac{\rho}{\rho_\infty}\right)\right] dv. $$

For $|b|$ small enough, using natural upper bounds, Sobolev injections, a Poincaré-like inequality of constant $\nu$ and Grönwall's lemma, we can get to

$$ \int_{-\infty}^{V_F} G\left(\frac{\rho}{\rho_\infty}\right) \rho_\infty dv \le e^{-\nu at} + \frac{8b^2}{a} \int_0^t |(N(s)-N_\infty)|^2 ds \int_{-\infty}^{V_F} G\left(\frac{\rho^0}{\rho_\infty}\right) \rho_\infty dv. \quad (1.13) $$

In order to conclude, we need a uniform estimate of the form

$$ \int_0^t N(s)^2 ds \le C(1+t), \quad (1.14) $$

where $C$ does not depend on $b$, so we can take $|b|$ small enough with respect to $\nu$, $a$ and $C$.

We can apply again a relative entropy method, but in a neighbourhood of $V_F$, in order to isolate the firing rate. However, it would not yield a $b$-uniform estimate. The trick is to compare the solution $\rho$ not with the stationary state $\rho_\infty$ but rather with $\rho_\infty^1$ obtained with another parameter $b_1$.

Fix $b_1 > 0$ such that there is at least one stationary state $(\rho_\infty^1, N_\infty^1)$. We work near $V_F$ by fixing some $V_M \in (V_R, V_F)$ and using the cutoff function $\gamma(v) = 1_{v>V_M} e^{\frac{1}{(V_F-V_M)^2-(V_F-v)^2}}$. Define

$$ W(v,t) = \left(\frac{p(v,t)}{p_\infty^1(v)}\right)^2 p_\infty^1(v) = \frac{p(v,t)^2}{p_\infty^1(v)} \quad \text{and} \quad I(t) = \int_{-\infty}^{V_F} W(v,t)\gamma(v)dv. $$

The key reason for this quantity $W$ is that applying the boundary conditions for the equations (NNLIF) associated to the parameters $b$ and $b_1$, we can define by continuity

$$ -a \frac{\partial W}{\partial v}(V_F, t) = \frac{N(t)^2}{N_\infty^1}, \quad W(V_F, t) = 0. $$

Then, it all comes down to finding differential inequalities for the modified relative entropy $I(t)$.

* **Inhibitory case:** Assume $b \le 0$. There exist $C_1, C_2, C_3 > 0$ independent of $b$ such that

$$ \frac{d}{dt} I(t) \le -C_1 N(t)^2 + C_2 + (bN(t) - C_3)I(t). $$

* **Excitatory case:** Assume $b > 0$. There exist $C_1, C_2, C_3, C_4 > 0$ independent of $b$ such that

$$ \frac{d}{dt} I(t) \le N(t)^2 (C_1 b^2 - C_2 + C_3 b^2 I(t)) + C_3 - C_4 I(t). \quad (1.15) $$

These differential inequalities allow to get uniform $L^2$ control over $N$. However, because of the multiplicative term $I(t)$ in (1.15), we need $b$ small enough with respect to $I(0)$ in the excitatory case $b > 0$. In both cases, it yields an estimate of the form (1.14) (see Theorem 1.17 for a generalised statement). Applying these $L^2$ estimates to (1.13), and choosing also $b$ small enough yields:
**Theorem 1.15** (Carrillo, Perthame, Salort, Smets [52]). Assume $d = 0$, $a(N) \equiv a > 0$ constant.
Let us fix $b_1 > 0$ such that there exists at least one corresponding stationary state $(\rho_\infty^1, N_\infty^1)$ and
$V_M \in (V_R, V_F)$. Let $(\rho, N)$ a classical solution of (NNLIF) such that,

$$ \int_{V_M}^{V_F} \frac{p^0(v)^2}{p_\infty^1(v)} dv < +\infty, \quad \int_{-\infty}^{V_F} \left( \frac{\rho^0(v) - \rho_\infty(v)}{\rho_\infty(v)} \right)^2 \rho_\infty(v) dv < +\infty. $$

There exist positive constants $b^*, \mu > 0$ such that for all $b \in [-b^*, b^*]$,

$$ \int_{-\infty}^{V_F} \left( \frac{\rho(v, t) - \rho_\infty(v)}{\rho_\infty(v)} \right)^2 \rho_\infty(v) dv \le e^{-\mu t} \int_{-\infty}^{V_F} \left( \frac{\rho^0(v) - \rho_\infty(v)}{\rho_\infty(v)} \right)^2 \rho_\infty(v) dv. $$

When $b \le 0$, $b^*$ and $\mu$ can be chosen independently of $\rho^0$ if we start at a larger time $T = T(\rho^0)$,

$$ \int_{-\infty}^{V_F} \left( \frac{\rho(v, t) - \rho_\infty(v)}{\rho_\infty(v)} \right)^2 \rho_\infty(v) dv \le e^{-\mu(t-T)} \int_{-\infty}^{V_F} \left( \frac{\rho(v, T) - \rho_\infty(v)}{\rho_\infty(v)} \right)^2 \rho_\infty(v) dv. $$

In the excitatory case $b > 0$, $\mu$ and $b^*$ heavily depend on $\rho^0$ and there is nothing we can do about it, since there is blow-up for some initial conditions however small we choose $b > 0$ (Theorem 1.12). The more concentrated $\rho^0$ is around $V_F$, the smaller $b$ has to be. In the inhibitory case, this a priori estimate is true for any small enough $b$ independently of $\rho^0$. Although the $L^2$ estimate on $N$ holds for any inhibitory connectivity parameter $b \in (-\infty, 0]$, we still need $b$ small for $\mu$ to be positive.

### 1.6.2 $L^q$ estimates on the firing rate and global existence for weak interaction

Theorem 1.15 is an a priori convergence result, valid as long as the solution exists; it does not enforce global-in-time existence by itself. Indeed, convergence happens in the relative entropy space $L^2(\rho_\infty^{-1})$ in which the boundary derivative operator is not bounded (see Subsection 1.10.1 below for another approach involving a new stronger space). Global-in-time existence for the weakly nonlinear case was achieved first in [81] using probability methods on the associated stochastic differential equation (see Subsection 1.8.1 below for an introduction of this SDE).

**Theorem 1.16** (Delarue, Inglis, Rubenthaler, Tanré [81]). For any initial condition $V_0 < V_F$, there exists $b^* = b^*(V_0)$ such that for all $b \in [0, b^*$], there is a unique global-in-time strong solution to the stochastic differential equation (1.22).

Later, the method of [52] was extended in [233] to prove a similar result in a deterministic PDE framework. The idea is to use an iteration argument to get from the $L^2$ estimates of [52] to $b$-uniform $L^q$ estimates for any integer $q \ge 2$.

**Theorem 1.17** (Roux, Salort [233], generalising [52]). Assume $d = 0$, $a(N) \equiv a > 0$ is constant. Fix $q \ge 2$ integer and $b_1 > 0$ such that there is at least one stationary state $(\rho_\infty^1, N_\infty^1)$ and $V_M \in (V_R, V_F)$.
For any compatible initial condition $\rho^0$ such that $\int_{V_M}^{V_F} \frac{p^0(v)^q}{p_\infty^{-1}(v)} dv < +\infty$, the solution of (NNLIF) satisfies

* There exists $C_q$ depending only on $b_1$, $q$, $a$, $V_F$ and $T = T(\rho^0)$ such that for all $b \le 0$,

$$ \forall t \ge T, \quad \int_T^t N(s)^q ds \le C_q (1 + t - T). $$

* There exist $C_q$ depending only on $b_1$, $q$, $a$, $V_F$, $\rho^0$ and $b^* = b^*(\rho^0)$ such that for all $0 < b \le b^*$,

$$ \forall t \ge 0, \quad \int_0^t N(s)^q ds \le C_q (1 + t). $$
Then, coming back to the Stefan-like free-boundary problem reformulation exposed above, it is possible to use $L^3$ estimates on the firing rate in the representation formula (1.10) and to repeat the contradiction argument used by [48] for the inhibitory case.

**Theorem 1.18** (Roux, Salort [233]). Assume $d = 0$, $a(N) \equiv a > 0$ is constant. Let $\rho^0$ be a compatible initial condition. There exists $b^* = b^*(\rho^0)$ such that for all $b \in (-\infty, b^*]$, the unique strong solution of (NNLIF) is global-in-time.

The reason $L^2$ estimates are not enough but $L^3$ estimates on $N$ yield global-in-time existence can be seen in the Hölder inequality

$$ \int_{0}^{\tau} M(\eta) \frac{\partial G}{\partial x}(s(\tau), \tau, s(\eta), \eta) d\eta \le \left( \int_{0}^{\tau} M(\eta)^3 d\eta \right)^{\frac{1}{3}} \left( \int_{0}^{\tau} \left| \frac{\partial G}{\partial x}(s(\tau), \tau, s(\eta), \eta) \right|^{\frac{2}{3}} d\eta \right)^{\frac{2}{3}}. $$

Because the free boundary position $s(\tau)$ is Lipschitz, the second integral is convergent. However, if the exponents $\frac{1}{2} + \frac{1}{2} = 1$ are used for the inequality, then the second part is a diverging integral.

## 1.7 The delayed NNLIF model and self-sustained oscillations

As we said above, the goal of Brunel and Hakim in [19] was to explain fast periodic activity in the strongly inhibitory case. In order to observe it in the model, we need to consider the delayed case $d > 0$.

In this subsection, we consider the case $a(N) \equiv a > 0$ constant. We talk about the first wave of results on the delayed NNLIF model. Recent advances are described below in Subsection 1.10

### 1.7.1 Global-in-time existence and relaxation for weak interaction

In [31], the aforementioned techniques for the non-delayed case $d = 0$ have been extended to the case $d > 0$. It is first still possible to use the Stefan free-boundary reformulation to prove local-in-time existence, with similar representation formulae for $\rho$ and $N$ in modified variables. The delayed structure even allows for shortcuts. The maximal time of existence is still characterised by

$$ T^* = \sup\{\tau \in \mathbb{R}_+ | N(t) < +\infty\}. $$

A striking difference with the case $d = 0$ is that there is unconditional global-in-time existence of strong solutions. Because of the delayed impact of the firing rate, the notion of universal super-solution can be extended to include the nonlinear part of the problem. On time intervals of length $d$, nonlinear super-solutions can be constructed using the firing rate at previous times in $[t - d, t)$. These super-solutions prove that $T^* = +\infty$ in all cases.

**Theorem 1.19** (Cáceres, Roux, Salort, Schneider [31]). Assume $d > 0$ and $a(N) \equiv a > 0$. Let $\rho^0$ be a compatible initial condition. Then there exists a unique global-in-time strong solution to (NNLIF).

Then, the work on $b$-uniform $L^2$ estimates and generalised entropy can be extended to the case with delay. Unfortunately, during the process, the iteration argument on time intervals of length $d$ forces the appearance of additional terms depending on $d$, the consequence being that smallness of $d$ needs also to be assumed and not only that of $b$.

**Theorem 1.20** (Cáceres, Roux, Salort, Schneider [31]). Under the assumptions of Theorem 1.15 but with $d > 0$, there exist $b_1 = b_1(d) > 0$, $b_2 = b_2(d, \rho^0) > 0$, $\mu = \mu(d, \rho^0) > 0$ such that for all $b \in [-b_1, b_2]$, there exists $S_0(\rho^0, N^0, b, d)$ satisfying

$$ \int_{-\infty}^{V_F} \left( \frac{\rho(v, t) - \rho_{\infty}(v)}{\rho_{\infty}(v)} \right)^2 \rho_{\infty}(v) dv \le e^{-\mu t} S_0(\rho^0, N^0, b, d). $$

As we shall see in Theorem 1.35, the fact that the smallness assumption depends on $d$ is a technical hypothesis that other methods can remove [25, 26].
### 1.7.2 Numerical insights, heuristics and partial results

**The excitatory delayed NNLIF: convergence versus infinite-time blow-up** We assume for now that we are in the excitatory and delayed case: $b > 0$ and $d > 0$. As we said, all solutions are global-in-time. Two very different regimes were identified (See Subsection 1.10.2 below for another perspective on this dichotomy using the sequence of pseudo-equilibria).

* When $b$ is very large ($b > b^*$): Since the stationary states are the same for $d = 0$ and $d > 0$, Theorem 1.6 tells us that there are no stationary states. In this regime, there are also no periodic solutions.

**Theorem 1.21 (Cáceres, Roux, Salort, Schneider [31].)** If $b > V_F - V_R$ and $V_F \le 0$, then for any $d \ge 0$ there are no classical periodic solutions to (NNLIF).

What then happens? This question was answered with detailed numerical simulations in [30]. For a large excitatory connection, there is infinite-time blow-up of the firing rate $N(t)$:

$$ \lim_{t \to +\infty} N(t) = +\infty \quad (\text{numerical result}). $$

The density then converges towards a plateau state

$$ \rho(\cdot, t) \underset{t \to +\infty}{\longrightarrow} \frac{1}{V_F - V_R} \mathbb{1}_{[V_R, V_F]}(\cdot) \quad (\text{numerical result}). \quad (1.16) $$

A natural interpretation is that neurons fire, return to $V_R$ and are then sent increasingly fast to $V_F$. The plateau is made up of a trail of neurons that are firing and resetting. Asymptotically, all of them join the train, and the density becomes 0 outside of $[V_R, V_F]$ and uniform inside. This feature seems meaningless from the neuroscience point of view. A simple fix is to add a refractory period, which will change the structure of stationary states in the highly excitatory regime [29, 33] (see the system (1.32) below).

* When $b$ is not too large ($b < b^*$): In this case, there is at least one stationary state. Periodic solutions have never been observed numerically in this regime either. Numerical solutions exhibit different behaviours depending on the value of $b$, $a$, $V_R$, $V_F$ and the initial condition:

i) if $b < V_F - V_R$, there is a unique equilibrium; solutions either converge directly towards this stationary state or the firing rate increases sharply, all or part of the density is absorbed by the boundary $V_F$ and after resetting at $V_R$ of a portion of the mass, there is convergence towards a stationary state.

ii) if $V_F - V_R \le b < b^*$ (note that $b^* \ge V_F - V_R$), there are often two equilibria; either the solutions converge directly towards a stationary state or the density converges towards a plateau state as in (1.16).

When the delay $d$ tends to 0, this dichotomy sharpens into convergence versus finite-time blow-up. The authors in [30] argue that the global-in-time blow-up friendly physical solutions of [82], see Definition 1.27, will converge towards the stationary state after blow-up in the non-delayed case $d = 0$. Note that these physical solutions are not defined in the case of large excitatory $b$, consistent with the fact that the plateau state is meaningless from the probabilistic perspective.

**The inhibitory delayed NNLIF and the ongoing quest for periodic solutions** Despite more than twenty years of numerical evidence [19, 33, 139], no analytical existence result was ever achieved for periodic solutions of (NNLIF). In [142], a heuristic approach was proposed to understand the appearance of self-sustained oscillations in the strongly inhibitory case.
<figure><img src="image_3.png" /><figcaption>Figure 1: Simulation of the PDE (NNLIF) with $V_F = 1$, $V_R = 0$, $a = 0.2$, $d = 1$ and $b = -45$; (A.) solution at three different times; (B.) evolution in time of the firing rate $N(t)$. Ikeda, Roux, Salort, Smets [142].</figcaption></figure>

We can look for a solution composed of a periodic wave of unit mass $\varphi(v, t) = \phi(v - c(t))$ defined on the whole real line $\mathbb{R}$ plus a corrective term $R(v, t)$ needed to account for the boundary and reset conditions:

$$p(v, t) = \phi(v - c(t)) + R(v, t),$$

The wave $\varphi : (v, t) \mapsto \phi(v - c(t))$ solves

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial \varphi}{\partial t} + \frac{\partial}{\partial v} [(-v + bN(t-d))\varphi] - a \frac{\partial^2 \varphi}{\partial v^2} = 0, \quad v \in \mathbb{R}, t > 0, \\ \varphi(-\infty, t) = \varphi(+\infty, t) = 0, \quad \int_{-\infty}^{+\infty} \varphi(v,t)dv = 1. \end{array} \right. \qquad (1.17) $$

The remainder term $R : (v, t) \mapsto R(v, t)$ must be a solution of

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial R}{\partial t} + \frac{\partial}{\partial v} [(-v + bN(t-d))R] - a \frac{\partial^2 R}{\partial v^2} = (\delta_{V_R} - \delta_{V_F})N(t), \quad v \in \mathbb{R}, \\ R(-\infty, t) = R(+\infty, t) = 0, \quad R(v, 0) = p^0(v) - \varphi(v, 0), \quad v \in \mathbb{R}. \end{array} \right. \qquad (1.18) $$

<figure><img src="image_4.png" /><figcaption>Figure 2: Simulation of the PDE (NNLIF) with $V_F = 0$, $V_R = -2$, $a = 0.2$, $d = 1$ and $b = -35$; (A.) solution at three different times; (B.) evolution in time of the firing rate $N(t)$. Ikeda, Roux, Salort, Smets [142].</figcaption></figure>

**Theorem 1.22** (Ikeda, Roux, Salort, Smets [142]). Let $\phi : \mathbb{R} \to \mathbb{R}$ be the function defined by $\phi(v) = \frac{1}{\sqrt{2\pi a}} e^{-\frac{v^2}{2a}}$. Let $c$ be a solution of

$$c'(t) + c(t) = bN(t - d).$$

Then the function $\varphi$ defined by $\varphi(v, t) = \phi(v - c(t))$ is a solution of (1.17).
Unfortunately, systems (1.17) and (1.18) are strongly coupled through the firing rate:

$$N(t) = -a\partial_v\varphi(V_F, t) - a\partial_v R(V_F, t).$$

In order to make the problem autonomous and thus theoretically tractable, we assume that $b\partial_v R(V_F, t) \simeq 0$ in an appropriate sense when $b \to -\infty$. Hence, we replace $N(t)$ by the simpler firing rate

$$\mathcal{N}(c(t)) = -a\frac{\partial\varphi}{\partial v}(V_F, t).$$

in equation 1.22. It reduces the problem to the autonomous delayed differential equation (DDE)

$$c'(t) + c(t) = b\mathcal{N}(c(t-d)), \quad \mathcal{N}(c) = \frac{1}{\sqrt{2\pi a}}(V_F - c)e^{-\frac{(V_F-c)^2}{2a}}. \quad (1.19)$$

Exploiting the classical methods for DDEs [129, 249, 126], it is possible to prove existence of a periodic solution $c$, and thus a periodic solution $\varphi(v, t) = \frac{1}{\sqrt{2\pi a}}e^{-\frac{(v-c(t))^2}{2a}}$ of the simplified nonlinear PDE

$$\frac{\partial \varphi}{\partial t} + \frac{\partial}{\partial v} \left[ (-v + b\mathcal{N}(t-d))\varphi \right] - a \frac{\partial^2 \varphi}{\partial v^2} = 0, \quad \mathcal{N}(t) = -\frac{\partial \varphi}{\partial v}(V_F, t) \quad v \in \mathbb{R}, t > 0.$$

**Theorem 1.23** (Ikeda, Roux, Salort, Smets [142]). Assume $V_F^2 \ge a$, then there exist $b_* < 0$ and $d_* > 0$ such that (1.19) has a non-constant periodic solution $c(t)$ for any $b < b_*$ and $d > d_*$.

The way to prove this result is to use the Browder fixed point theorem.

**Definition 1.24** (Ejective fixed-point [17]). Let $E$ a Banach space, $D \subset E$ a closed subset, $\mathcal{F}: D \to D$ a continuous map. A fixed-point $\bar{x}$ of $\mathcal{F}$ is ejective if there is a neighbourhood $U \subset D$ of $\bar{x}$ such that

$$\forall x \in U \setminus \{\bar{x}\}, \exists n \in \mathbb{N}, \mathcal{F}^n(x) \notin U.$$

**Theorem 1.25** (Browder's fixed-point theorem [17]). Let $D$ be a closed, bounded and convex subset of an infinite dimensional Banach space and let $\mathcal{F}: D \to D$ a continuous and compact map. Then $\mathcal{F}$ admits at least one fixed-point which is not ejective.

Spectral study of (1.19) allows to prove that for large enough delay and strong enough inhibition, the stationary state $c_\infty$ is an ejective fixed-point of a carefully chosen functional $\mathcal{F}$. Then, there is a non-ejective fixed-point, which is associated to a periodic solution of (1.19).

In [142], asymptotic information is provided in the limit $b \to -\infty$. Assume for the sake of clarity that $V_F = 0$, $d = 1$ and make the change of variables

$$\beta = \frac{1}{\log(-b)}, \quad u(s) = \beta \log \left( -\frac{1}{\sqrt{a}} c \left( \frac{s}{\beta} \right) \right),$$

then we come to the following equivalent system:

$$u'_{\beta}(t) + 1 = \exp \left( \frac{1}{\beta} (1 - u_{\beta}(t)) + f(u_{\beta}(t - \beta)) \right), \quad 1 \ge u_{\beta}([-\beta, 0]) > 0 \quad (1.20)$$

with $f(x) = \frac{1}{\beta}x - \frac{1}{2}e^{\frac{2}{\beta}x}$. Equation (1.20) has a unique positive stationary state $\bar{u}_\beta = \frac{\beta}{2} \log(\frac{2}{\beta})$. Note that $\lim_{\beta \to 0} \bar{u}_\beta = 0$.

**Theorem 1.26** (Ikeda, Roux, Salort, Smets [142]). Assume that for all $\beta > 0$, $u_\beta([-\beta, 0]) \equiv 1$ and let $T > 2$. Then, there exists a constant $C_M \in (0, 1]$ such that $\lim_{\beta \to 0} u_\beta = P$, in $L^1(0, T)$, where P has the following form:

$$P(t) = 1 - t \text{ on } (0, 1); \quad P(t) = C_M + 1 - t \text{ on } (1, 1 + C_M); \quad P \text{ is } C_M \text{ periodic on } (1, T).$$
The assumption $u_β([-β, 0]) = 1$ is purely technical, as can be checked by numerical simulations. In original variables, the result can be interpreted as follows

* The period $T_b$ of solutions of (1.19) evolves in $O(\log(-b))$ when $b \to -\infty$.

* More precisely, we have in original variables

$$c(t) \approx -\sqrt{a} \exp \left( \log(-b) P \left( \frac{t}{\log(-b)} \right) \right).$$

In decay phases, $c$ behaves like $c(t) \approx -\sqrt{a}e^{C_M \log(-b)-t}$, which constitutes an exponential decay over a spatial area of length $S_b \approx \sqrt{a}(-b)^{C_M}$ during a time $T_b \approx C_M \log(-b)$.

* between decays phases, $c$ has very rapid growth which is discontinuous in the re-scaled limit.

For another asymptotic approach using the limit $d \to +\infty$ with $b$ fixed instead of $b \to -\infty$ with $d$ fixed, see Subsection 1.10.2 below.

# 1.8 Rigorous derivation of the NNLIF model

## 1.8.1 Physical solutions from the limit of a particle system

Although the NNLIF model was introduced in the late 90's, the first rigorous construction from a mean-field limit of interacting particles was published way later, in 2015, in the article [82]. Note that an overview on technical questions related to stochastic processes depending on singular hitting times was provided by the authors in the companion document [80].

For the sake of simplicity, the authors fix $V_R = 0$, $V_F = 1$ and make the assumption of a fully connected network with mean-field scaling. They consider the system of $K$ particles satisfying the stochastic differential equation

$$V_t^{i,K} = V_0^{i,K} - \int_0^t V_s^{i,K} ds + \frac{b}{K} \sum_{j=0}^K M_t^{j,K} + B_t^i - (V_F - V_R) M_t^{i,K}, \quad (1.21)$$

where $V_t^i$ represents the electrical potential of the $i^{th}$ neuron, all $V_0^{i,K}$ are independent and identically distributed with law $X_0$, and $(B_t^i)_{t \ge 0}$ are independent Brownian motions on a probability space $(\Omega, \mathcal{A}, \mathbb{P})$ with respect to some adequate filtration. The processes $(M_t^{i,K})$ count the number of times the $i^{th}$ neuron has spiked during the interval of time $[0, t]$, which can be precisely defined as càdlàg processes:

$$M_t^{i,K} = \sum_{l \ge 1} 1_{[0,t]} (\tau_l^{i,K}), \quad \tau_l^{i,K} = \inf \left\{ t > \tau_{l-1}^{i,K} \mid V_{t-}^{i,K} + \frac{\alpha}{K} \sum_{j=1}^K (M_t^{j,K} - M_{t-}^{j,K}) \ge V_F \right\},$$

and $\tau_0^{i,K} = 0$. This system was previously introduced and studied in a physics framework in [204]. In order to avoid neurons spiking multiple times at once and ensure that spikes happen in a well-ordered physical way, [82] imposes $b < V_F - V_R$ (see also the problem of eternal blow-up in [91], described later in Section 1.9) and defines rigorously the notion of cascade: a first set $\Gamma_0$ of neurons spike, then a second set of neurons $\Gamma_1$ spike given that neurons in $\Gamma_0$ spike, etc. A cautious definition yields uniqueness for the solution of the particle system (1.21).

Since the counting processes $M_t^{i,K}$ are naturally monotone, [82] uses the M1 Skorokhod topology to perform the mean-field limit $K \to +\infty$. The limiting system writes

$$V_t = V_0 - \int_0^t V_s ds + b \mathbb{E}[M_t] + B_t - (V_F - V_R) M_t. \quad (1.22)$$

Because $\mathbb{E}[M_t]$ might be a discontinuous function, defining a notion of solution for (1.22) is challenging. Recall, for any càdlàg function $F$, the notation $\Delta F(t) = F(t) - F(t^-)$.
**Definition 1.27** (Physical solution [82]). We say that the càdlàg adapted process $(V_t, M_t)_{t \ge 0}$ is a physical solution of (1.22) if

• $(M_t)_{t \ge 0}$ has integrable marginal distributions and for all $t \ge 0$, $\mathbb{P}(\Delta M_t \le 1) = 1$;

• (1.22) holds true $\mathbb{P}$-almost surely with

$$M_t = \sum_{l \ge 1} 1_{[0,t]}(\tau_l), \quad \tau_l = \inf \{t > \tau_{l-1} \mid V_{t-} + b \Delta \mathbb{E}[M_t] \ge V_F\}, \quad \tau_0 = 0;$$

• The discontinuity points of $t \mapsto \mathbb{E}[M_t]$ satisfy

$$\Delta \mathbb{E}[M_t] = \inf \{\eta \ge 0 \mid \mathbb{P}(V_{t-} + b\eta \ge V_F) < \eta\}.$$

Through the mean-field limit in M1 Skorokhod topology, [82] proves the existence of a physical solution for any initial condition of law in $L^p, p \ge 1$, and supported away from $V_F$. Note that for a fixed end-time $T > 0$, almost surely

$$\lim_{K \to +\infty} \frac{1}{K} \sum_{i=1}^{K} M_t^{i,K} = \mathbb{E}[M_t]$$

for all $t \in [0, T)$, except at any discontinuity points of $t \mapsto \mathbb{E}[M_t]$, of which there are countably many.

Physical solutions can also be obtained in the vanishing delay limit $d \to 0$, also through convergence in the M1 Skorokhod topology [82]. In [81], with more details in the notes [80], it is also proved that, at least up to the first blow-up time, which corresponds to the first discontinuity point of $\mathbb{E}[M_t]$, the law of a physical solution of (1.22) is a solution in the weak sense of the PDE (NNLIF) with $d = 0$ and $a \equiv a > 0$.

## 1.8.2 Smooth deterministic solutions from an iteration approach in the linear case

The link between the stochastic differential equation (1.22) and the PDE (NNLIF) is investigated in depth in [166] for the simpler linear case $b = 0, d = 0, a(N) \equiv a > 0$ constant. They choose for simplicity $V_R = 0, V_F = 1, a = 1$.

The idea of [166] is to rely on an iteration perspective, considering a collection of processes representing the laws of the successive firing events for a representative neuron. They consider a sequence of independent Ornstein-Uhlenbeck processes $(Y_t^n)_{n \in \mathbb{N}}$, each representing the trajectory from $V_R$ to $V_F$ after a spike event and before the next. They introduce a sequence of hitting times $(\tau_n)_{n \in \mathbb{N}}$ at $V_F = 1$ for the processes $Y^n$ and the jump times for the process $V_t$ are defined as

$$T_n = \sum_{i=0}^{n} \tau_i.$$

The original process can then be recovered with the formula $V_t = Y_{t-T_{k-1}}^k$ on the time interval $[T_{k-1}, T_k)$. The process $n_t = k-1$ keeps track of which process $Y^k$ is being used, resulting in a pair $(V_t, n_t) = (Y_{t-T_{k-1}}^k, k-1)$. Denoting

$$F_k(v, t) = \mathbb{P}(V_t \le v, n_t = k), \qquad (1.23)$$

the cumulative distribution function after $k$ jumps, it is possible to construct an infinite series of PDE systems such that $F_k$ is associated to the solution $\rho_k$ of each of them. The method relies on the key iteration formula

$$F_n(x, t) = \int_0^t F_{n-1}(x, t-s) dG_1(s),$$
where $G_1$ is the cumulative distribution for the first spike time given by $G_1(t) = \mathbb{P}(\tau_1 \le t)$. The collection of systems write for $n \ge 1$

$$ egin{cases} 
\frac{\partial \rho_n}{\partial t} - \frac{\partial}{\partial v}(v\rho_n) - \frac{\partial^2 \rho_n}{\partial v^2} = \delta_{v=v_R} N_{n-1}(t), & v \in (-\infty, 1], t > 0, \\ 
N_n(t) = -\frac{\partial \rho_n}{\partial v}(1, t), & \rho_n(1, t) = \rho_n(-\infty, t) = 0, \quad t > 0, \\ 
\rho_n(v, 0) = 0 \ge 0, & v \in (-\infty, 1]. 
\end{cases} $$

and for the initial density $\rho_0$,

$$ \frac{\partial \rho_0}{\partial t} - \frac{\partial}{\partial v}(v \rho_n) - \frac{\partial^2 \rho_0}{\partial v^2} = 0, \quad N_0(t) = -\frac{\partial \rho_0}{\partial v}(1, t), \quad \rho_0(-\infty, t) = \rho_0(1, t) = 0, $$

and $\rho_0(v, 0) = \rho^0(v) \ge 0$, $\int_{-\infty}^{V_F} \rho^0(v)dv = 1$. The function $\rho^0$ is the probability density of finding a neuron which has never spiked in the network. The density $\rho_n$ is the probability density of neurons which have spiked $n$ times but not $n+1$ times. Each firing rate $N_k$ is a source term in the next equation for $\rho_{k+1}$ whose initial condition is uniformly 0. It is then possible to prove rigorously the decomposition

$$ \rho(v, t) = \sum_{k=0}^{+\infty} \rho_k(v, t). $$

Each of these densities $\rho_k$ is associated to the cumulative distribution function $F_k$ in (1.23). It can be used to construct strong solutions to the original PDE (NNLIF) from the stochastic process (1.22).

**Theorem 1.28** (Liu, Wang, Zhang, Zhou [166]). Assume $d = 0$, $b = 0$, $a(N) \equiv 1$, $V_F = 1$ and $V_R = 0$. If $\rho^0$ is a compatible initial condition compactly supported in $(-\infty, V_F)$ or if $\rho^0 = \delta_{v=v_0}$ is a Dirac mass at a point $V_0 \in (-\infty, V_F)$, then the unique solution to (1.22) has a continuously evolved probability density $\rho(\cdot, t)$ which is for positive times a global-in-time strong solution to (NNLIF).

# 1.9 Generalised solutions and continuation after blow-up

## 1.9.1 Case with internal noise

Most of the results on the NNLIF model were performed in the case of a constant noise $a(N) \equiv a > 0$. The mean-field limit in [82] leads to a constant noise coefficient. However, earlier derivations by physicists often propose a noise strength $a(N) = a_0 + a_1 N$, $a_0, a_1 > 0$. As we have seen in Theorem 1.14, this internal noise case is more ill-posed as even in the inhibitory case $b < 0$, finite time blow-up can occur. Paradoxically, this parameter range is more amenable for the construction of generalised solutions with PDE tools.

In [91], an approach is proposed that allows both the construction of smooth classical solutions and the extension beyond blow-up points if some conditions are satisfied. The idea is to take advantage of the fact that, although in the constant noise case blow-up is equivalent to $\limsup_{t \to T^*} [-\partial_v \rho(V_F, t)] = +\infty$, the internal noise case allows blow-up to happen while $\partial_v \rho(V_F, t)$ remains bounded. It can be seen by solving for $N$ in $N = -(a_0 + a_1 N) \partial_v \rho(V_F, t)$:

$$ N(t) = \frac{-a_0 \partial_v \rho(V_F, t)}{1 + a_1 \partial_v \rho(V_F, t)}. $$

The main idea of [91] is to rescale time in function of the firing rate, in order to turn the divergent speed in the drift term into a locally-bounded one while enjoying uniform parabolicity in the second order term. They propose the change of variable, for some $c > 0$,

$$ d\tau = (c + N(t))dt. $$
Before applying the change of variable, they substitute $N(t)$ for $-(a_0 + a_1 N)\partial_v\rho(V_F, t)$ in the right-hand side, which leads to the $\tau$ time-scale evolution equation

$$ \frac{\partial \tilde{\rho}}{\partial \tau} + \frac{\partial}{\partial v} \left[ (-v\tilde{N} - cb\tilde{N} + b)\tilde{\rho} \right] - (a_c\tilde{N} + a_1) \frac{\partial^2 \tilde{\rho}}{\partial v^2} = -(a_c\tilde{N} + a_1) \frac{\partial \tilde{\rho}}{\partial v} (V_F, \tau) \delta_{v=V_R}, \quad (1.24) $$

with $a_c = a_0 - ca_1$, $\tilde{\rho}(-\infty, t) = \tilde{\rho}(V_F, t) = 0$, and the new firing rate

$$ \tilde{N}(\tau) = \frac{1}{N+c} = \frac{[1+a_1 \partial_v \tilde{\rho}(V_F, \tau)]_+}{-a_0 \partial \tilde{\rho}(V_F, \tau) + c[1+a_1 \partial_v \tilde{\rho}(V_F, \tau)]_+} \quad (1.25) $$

This modified firing rate is uniformly bounded: $0 \le \tilde{N}(\tau) \le \frac{1}{c}$, which implies uniform parabolicity because

$$ \min \left( a_1, \frac{a_0}{c} \right) \le a_c \tilde{N}(\tau) + a_1 \le \max \left( a_1, \frac{a_0}{c} \right). $$

The original time can then be recovered via the formula $t(\tau) = \int_0^\tau \tilde{N}(\eta)d\eta$. Adapting the method of [48], [91] constructs global-in-time solutions to the time-dilated problem.

**Theorem 1.29** (Dou, Zhou [91]). For any compatible initial condition $\rho^0$, any fixed $c$, there exists a unique global-in-time strong solution to (1.24)–(1.25).

It is then possible to recover a $c$-independent unique generalised solution for the original time problem (NNLIF) which is defined up to time

$$ T^* = \int_0^{+\infty} \tilde{N}(\tau) d\tau. $$

The generalised firing rate in original time $t$ is then

$$ N(t) = \frac{1}{\tilde{N}(\tau(t))} - \frac{1}{c} \in [0, +\infty], \quad \tau(t) = \sup \left\{ \tau \in (0, +\infty) \mid \int_0^\tau \tilde{N}(\eta) d\eta = t \right\}. $$

**Theorem 1.30** (Dou, Zhou [91]). Assume $d=0$, $a(N) = a_0 + a_1 N$, $a_0, a_1 > 0$. For any compatible initial condition $\rho^0$, there exists a unique generalised solution to (NNLIF) defined on the maximal interval $[0, T^*)$.

* If $b < 0$ and $V_F \ge 0$, then $T^* = +\infty$.

* If $0 < b < V_F - V_R$, $\rho^0(v)e^{-\frac{b}{a_1}v}$ is bounded and $V_F \ge -\frac{a_0}{a_1}b$, then $T^* = +\infty$.

* If $V_F - V_R < b$, there exist initial conditions such that $T^* < +\infty$.

The case where generalised solutions cannot be global-in-time in the high connectivity case is called an eternal blowup. In the time-dilated problem (1.24), the solution never exits blow-up regime and $\tilde{N}$ stays zero indefinitely. Note that this nonphysical case was excluded by design in the mean-field limit of [82], as they assume $b < V_F - V_R$. In the other cases with weaker or inhibitory interaction, the idea of the proof is to use an entropy method in time-dilated timescale to prove that if the solution enters blow-up, it must exit in finite time. A crucial ingredient is to use, the fact that if $T^*$ is assumed by contradiction to be finite, then $\tilde{N}$ must be small for large times, enabling weakly-nonlinear analysis.

Those generalised solutions come from specific choices: how to rewrite the $N(t)\delta_{v=V_R}$ source term and how to do the dynamics inside blow-up. Other choices could be made, leading to other types of generalised solutions. The approach in [91] enjoys many good mathematical properties, but comes with the caveat of allowing neurons to spike multiple times during a blow-up event.
### 1.9.2 Case without internal noise

When the diffusion is constant $a(N) \equiv a$, the method presented above no longer works because the equation in a dilated timescale is no longer parabolic. Moreover, in the case without internal noise, the Dirichlet boundary condition at $V_F$ is not satisfied during blow-up. To overcome these issues, an other approach that involves time dilation was recently developed in [90]. Before blow-up, the problem (NNLIF) can be recast on the whole real line by putting a $-N(t)\delta_{v=V_F}$ in the right-hand side. In order to generalise this idea after blow-up, a more general measure $S$ has to be considered. More precisely, introduce

$$d\tau = N(t)dt, \quad n(v, \tau) = \rho(v, t), \quad Q(\tau) = \frac{1}{N(t)}.$$

Then, $(n, Q, S)$ satisfies

$$\frac{\partial n}{\partial \tau} + \frac{\partial}{\partial v} [(-Q(\tau) + b)n] - aQ(\tau) \frac{\partial^2 n}{\partial v^2} = \delta_{v=V_R} - S(v, \tau), \quad \tau \ge 0, \ v \in \mathbb{R}. \quad (1.26)$$

In the classical regime when $\int_{V_F}^{+\infty} n(v, \tau)dv = n(V_F, \tau) = 0$, then $S(v, \tau) = \delta_{v=V_F}$ and
$$-a\partial_v n(V_F, \tau)Q(\tau) = 1.$$

In the blow-up regime when $\int_{V_F}^{+\infty} n(v, \tau)dv > 0$, $Q(\tau) = 0$. The other conditions for the generalised boundary condition write

$$S(\cdot, \tau) \ge 0, \quad \int_{\mathbb{R}} S(\cdot, \tau)(dv) = 1, \quad S(\cdot, \tau) \equiv 0 \text{ on } (-\infty, V_F). \quad (1.27)$$

Whenever $Q(\tau) > 0$, a solution to (1.26)–(1.27) and related conditions corresponds to a solution to (NNLIF). In order to obtain a generalised solution in this context, [90] uses the NNLIF model with random firing (1.33) presented below in Subsection 1.11.2, with the choice $\phi_\varepsilon = \frac{1}{\varepsilon}1_{v \ge V_F}$. Working in the dilated timescale $\tau$ allows to obtain $\varepsilon$-uniform estimates and pass to the limit $\varepsilon \to 0$.

Given an interval $I \subset \mathbb{R}$, consider $M(I)$ the space of bounded measures on $I$, $M_+(I) \subset M(I)$ the space of positive bounded measures on $I$ and $P_2(I) \subset M_+(I)$ the space of probability measures with second moment.

**Theorem 1.31** (Dou, Perthame, Salort, Zhou [90]). Assume the initial condition $n^0 \in L^2(\mathbb{R}) \cap P_2(\mathbb{R})$. Then, there is a global-in-time weak solution to (1.26)–(1.27) such that $n \in L^\infty(0, +\infty; L^1(\mathbb{R})) \cap L^\infty_{\text{loc}}(0, +\infty; L^2(\mathbb{R}))$ and for all $\tau_0 > 0$, in $(0, \tau_0)$ we have $Q \in M_+(0, \tau_0)$, $\int_{\mathbb{R}} v^2 n(v, \tau)dv \le C(\tau_0)$. Moreover, for all $\psi \in L^2(\mathbb{R}) \cap L^\infty(\mathbb{R}) \cap C(\mathbb{R})$, $\tau \mapsto \int_{\mathbb{R}} \psi(v)n(v, \tau)dv \in C(0, +\infty)$ and

$$Q(\tau) \int_{V_F}^{+\infty} n(v, \tau)dv = 0.$$

Similar to the generalised solution in the internal noise case described above, the generalised solution $\rho$ defined on $[0, T^*)$ in original timescale can be global-in-time ($T^* = +\infty$) or enter an eternal blow-up ($T^* < +\infty$). Again we have

$$T^* = \int_{0}^{+\infty} Q(\tau)d\tau.$$

When $b < V_F - V_R$, [90] proves that $T^* = +\infty$; when $b > V_F - V_R$, examples of eternal blow-up are found, which correspond to the plateau states identified in [30] that we have described above.

## 1.10 Asymptotic behaviour for general connectivity strength and delay

In most previous study, analytical results on the asymptotic behaviour of (NNLIF) required smallness assumption on both the connectivity parameter $b$ and the delay parameter $d$. In two companion papers
[25, 26], new methods were proposed in order to study the long-time behaviour for large b and/or large d. The authors fix $a(N) \equiv 1$ for convenience, but claim their strategy could be extended to any constant $a(N) \equiv a > 0$ with only technical changes. For simplicity of notations in their approach, they also consider initial conditions of the form

$$ \rho^0(v, t), \quad N^0(t) = \frac{\partial \rho^0}{\partial t}(V_F, t), \quad v \in (-\infty, V_F], t \in [-d, 0]. $$

from which initial conditions in the sense of Definition 1.4 can easily be retrieved provided enough regularity and decay at infinity.

### 1.10.1 Improved convergence estimates with a spectral gap approach

In order to better understand the nonlinear system (NNLIF) in the case of strong nonlinearities, a possibility is to first look at a linearised version of the firing rate impact on the network. Fix an external firing rate $N \ge 0$, a connectivity parameter $b \in \mathbb{R}$ and define the linear operator

$$ L_N \rho = \frac{\partial^2 \rho}{\partial v^2} + \frac{\partial}{\partial v}[(v - bN)\rho] + \delta_{v=V_R} N_{\rho}, \quad N_{\rho} = -\frac{\partial \rho}{\partial v}(V_F, t), $$

naturally associated to the linear PDE $\partial_t \rho = L_N \rho$, $\rho(V_F, t) = 0$, $\rho(\cdot, 0) = \rho^0(\cdot)$. Now write $\rho = \rho_\infty + u$ as a perturbation $u$ from a stationary state $(\rho_\infty, N_\infty)$ and denote $L_\infty = L_{N_\infty}$. Then $u$ solves the equation

$$ \frac{\partial u}{\partial t} = L_{\infty}u - bN_{u}(t-d)\frac{\partial \rho_{\infty}}{\partial v} - bN_{u}(t-d)\frac{\partial u}{\partial v}, $$

which yields, ditching the nonlinear term for small $u$, the linearised problem

$$ \frac{\partial u}{\partial t} = Tu, \quad Tu = L_{\infty}u - bN_{u}(t-d)\frac{\partial \rho_{\infty}}{\partial v}, \quad u(-\infty, t) = u(V_{F}, t) = 0. \quad (1.28) $$

Because of the structure of the stationary states, the need to control the derivative at the boundary and the framework for classical well-posedness, a natural space for the study of strong solutions is

$$ X = \{ \rho \in C(-\infty, V_F] \cap C^1(-\infty, V_R] \cap C^1[V_R, V_F] \mid \rho(V_F) = 0, \|\rho\|_X < +\infty \}, $$

where

$$ \|\rho\|_X = \|\rho\|_\infty + \|\partial_v \rho\|_\infty + \|\rho\|_{L^2(\phi)} + \|\partial_v \rho\|_{L^2(\phi)}, \quad \phi(v) = e^{-\frac{(v-bN)^2}{2}}. $$

It can be readily checked that this space embeds into the natural space for relative entropy estimates: $\|\rho\|_{L^2(\rho_\infty^{-1})} \le C \|\rho\|_X$ and $X \subset L^2(\rho_\infty^{-1})$. Note that the spectrum of the linear problem in $L^2(\rho_\infty^{-1})$ had been studied before [25], although less thoroughly, in [48].

**Definition 1.32.** We say that a stationary state $(\rho_\infty, N_\infty)$ is linearly stable if there exists $C \ge 1$ and $\lambda > 0$ such that all solutions $u$ to (1.28) from an initial condition $u_0 \in C([-d, 0], X)$ such that $\int_{-\infty}^{V_F} u_0(v, \cdot) dv \equiv 0$ on $[-d, 0]$, satisfy

$$ \|u(\cdot, t)\|_X \le C e^{-\lambda t} \sup_{\tau \in [-d, 0]} \|u_0(\cdot, \tau)\|_X. $$

It is linearly unstable if for any $C \ge 1$, there exists $u_0 \in C([-d, 0], X)$ and $T \ge 0$ such that

$$ \|u(\cdot, T)\|_X \ge C \sup_{\tau \in [-d, 0]} \|u_0(\cdot, \tau)\|_X. $$
Given this definition, the next step is to study the evolution along the linear flow of the firing rate of the small perturbation $u$, which can be done by introducing the functions

$$q(v, t) = e^{tL_{\infty}} \frac{\partial \rho_{\infty}}{\partial v}, \qquad N_q(t) = -\frac{\partial q}{\partial v}(V_F, t),$$

and the Laplace transform

$$\hat{N}_q(\xi) = \int_0^{+\infty} e^{-\xi t} N_q(t) dt, \qquad (1.29)$$

which [25] proves to be well-defined for a large enough real part $\mathcal{R}(\xi) > -\lambda$, as there exist $C \ge 1, \lambda > 0$ such that $\hat{N}_q(t) \le Ce^{-\lambda t}$.

The cornerstone of the use of this Laplace transform is the following result.

**Theorem 1.33** (Cáceres, Cañizo, Ramos-Lora [25]). Let $(\rho_{\infty}, N_{\infty})$ be a stationary state and associated $\hat{N}_q$ as in (1.29); it is linearly stable if and only if all zeros of the analytic function

$$\Phi_d(\xi) = 1 + b \hat{N}_q(\xi) e^{-\xi d}$$

are located on the real-negative half plane $\{\xi \in \mathbb{C} | \mathcal{R}(\xi) < 0\}$.

From this characterisation, a surprising consequence is that the linear stability of stationary states of (NNLIF) can be related to the real function $I(N)$ used for the study of stationary states.

**Theorem 1.34** (Cáceres, Cañizo, Ramos-Lora [25]). Let $(\rho_{\infty}, N_{\infty})$ be a stationary state and associated $\hat{N}_q$ as in (1.29);

* $\rho_{\infty}$ is linearly unstable for all delay $d \ge 0$ under the sufficient condition

$$\frac{d}{dN} \left( \frac{1}{I(N)} \right)_{|N=N_{\infty}} > 1;$$

* $\rho_{\infty}$ is linearly stable for all delay $d \ge 0$ under the sufficient condition

$$|b| \int_{0}^{+\infty} |N_q(t)| dt < 1;$$

* $\rho_{\infty}$ is linearly unstable if the delay $d$ is large enough under the sufficient condition

$$\frac{d}{dN} \left( \frac{1}{I(N)} \right)_{|N=N_{\infty}} < -1;$$

where the function $N \mapsto I(N)$ is the one defined in (1.7).

This classification is not exhaustive theoretically, but numerical results show that it is almost comprehensive in practice, because $N_q$ has a fixed sign, leading to

$$|b| \int_{0}^{+\infty} |N_q(t)| dt = \left| \frac{d}{dN} \left( \frac{1}{I(N)} \right)_{|N=N_{\infty}} \right|$$

and the limit case $= \pm 1$ is marginal. From these linear stability results, [25] draws rigorous consequences for the nonlinear system (NNLIF).

First, the weakly-nonlinear stability performed with relative entropy in previous studies [27, 52, 31, 32, 33] can be performed in the smaller space $X$ with improved assumptions on $d$.
**Theorem 1.35** (Cáceres, Cañizo, Ramos-Lora [25]). Assume $a(N) \equiv 1$ and $|b|$ is small enough to have a unique stationary state $(\rho_\infty, N_\infty)$. There exists $C \ge 1$, $\lambda > 0$ depending only on $V_R, V_F$ such that for all initial condition $\rho^0 \in C([-d, 0], X)$, there exists $b^* > 0$ depending only on $V_R, V_F$ and

$$ \sup_{t \in [-d, 0]} \|\rho^0(\cdot, t) - \rho_\infty(\cdot)\|_X, $$

such that for all $b \in (-b^*, b^*)$, the solution $\rho$ to (NNLIF) satisfies

$$ \forall t \ge 0, \quad \|\rho(\cdot, t) - \rho_\infty(\cdot)\|_X \le C e^{-\frac{\lambda}{2+d}t} \|\rho^0(\cdot, 0) - \rho_\infty(\cdot)\|_X. $$

The small connectivity limit $b^*$ depends on the initial condition $\rho^0$, but is independent on the delay value $d$, significantly improving the earlier findings of [31] (see Theorem 1.20 above).

Second, it can be proved that linear stability entails nonlinear stability beyond the weakly-nonlinear case when there is no delay $d = 0$.

**Definition 1.36.** We say that the flow $(e^{tT})_{t \ge 0}$ of the linearised problem (1.28) has a spectral gap $\lambda$ with constant $C \ge 1$ if for all $u_0 \in X$ such that $\int_{-\infty}^{V_F} u_0(v)dv$, the solution $e^{tT}u_0$ to (1.28) satisfies

$$ \|e^{tT}u_0\|_X \le C e^{-\lambda t} \|u_0\|_X. $$

**Theorem 1.37** (Cáceres, Cañizo, Ramos-Lora [25]). Assume $a(N) \equiv 1$ and $d = 0$. Let $(\rho_\infty, N_\infty)$ be a stationary state. Assume the flow of the linearised problem (1.28) has a spectral gap of size $\lambda > 0$ with constant $C \ge 1$. There exists $\varepsilon > 0$ depending on $b, V_F, V_R, \rho_\infty, \lambda$ such that for any initial condition $\rho^0 \in X$, if $\|\rho^0 - \rho_\infty\|_X < \varepsilon$, then the solution $\rho$ to (NNLIF) satisfies

$$ \forall t \ge 0, \quad \|\rho(\cdot, t) - \rho_\infty(\cdot)\|_X \le 2C e^{-\frac{\lambda}{2}t} \|\rho^0 - \rho_\infty\|_X. $$

The authors of [25] claim that the result should hold for non-zero delay $d > 0$ and that an extension of the result is mainly a technical challenge.

### 1.10.2 The discrete sequence of pseudo-equilibria

In a companion paper [26] to [25], the troubling link between the function $I(N)$ defined in (1.7) and the stability of stationary states is studied further with the introduction of a sequence of pseudo-equilibria.

Consider the simple case $a \equiv 1$. Their main idea is to remark that when the delay $d$ is large enough, the system behaves as a linear equation with a given external drift on each time interval of size $d$. Given a fixed firing rate $N \ge 0$ acting on the network as an eternal source, the solution of (NNLIF) would converge exponentially fast to the probability density

$$ \rho_{\text{pseudo}}(v) = \tilde{N} e^{-\frac{(v-bN)^2}{2}} \int_{\max(v, V_R)}^{V_F} e^{\frac{(w-bN)^2}{2}} dw, \quad \tilde{N} = \frac{1}{I(N)}, \quad (1.30) $$

which is a stationary state of the full nonlinear system (NNLIF) if and only if $\tilde{N} = N$ or equivalently if $NI(N) = 1$. If the firing rate was $N_{0,\infty}$ on $[-d, 0]$, then on $[0, d]$ with $d$ large enough, the firing rate $N(t)$ of a solution to (NNLIF) has the time to approach a new equilibrium value $N_{1,\infty} \simeq -\partial_v \rho(V_F, d)$ given by (1.30). One can thus define a pseudo-equilibria sequence in order to approach the behaviour of a solution on a sequence of time intervals $[kd, (k+1)d]$.

**Definition 1.38.** We call the firing rate sequence associated to the initial firing rate $N_{0,\infty}$ the sequence $(N_{k,\infty})_{k \ge 0}$ recursively defined by

$$ N_{k+1,\infty} = \frac{1}{I(N_{k,\infty})}. $$

The associated pseudo-equilibria sequence $(\rho_{k,\infty})_{k \ge 0}$ is defined by

$$ \rho_{k,\infty} = N_{k,\infty} e^{-\frac{(v-bN_{k-1,\infty})^2}{2}} \int_{\max(v, V_R)}^{V_F} e^{\frac{(w-bN_{k-1,\infty})^2}{2}} dw. $$
A careful analysis of the function $\frac{1}{I(N)}$ allows to do an almost comprehensive classification of the asymptotic behaviour of the firing rate sequence.

**Theorem 1.39** (Cáceres, Cañizo, Ramos-Lora [26]). Assume $a(N) \equiv 1$.

* In the excitatory case $b > 0$:

- if there is a unique stationary firing rate $N_\infty$ and $N \mapsto \frac{1}{I(N)}$ crosses the diagonal $N \mapsto N$, $(N_{k,\infty})_{k \ge 0}$ converges monotonically to $N_\infty$;

- if there is a unique stationary firing rate $N_\infty$ and $N \mapsto \frac{1}{I(N)}$ touches tangentially the diagonal, $(N_{k,\infty})_{k \ge 0}$ converges to $N_\infty$ when $N_{0,\infty} < N_\infty$ and diverges when $N_\infty < N_{0,\infty}$;

- if there are two stationary firing rates $N_\infty^1 < N_\infty^2$, $(N_{k,\infty})_{k \ge 0}$ monotonically converges to $N_\infty^1$ when $N_{0,\infty} < N_\infty^2$ and diverges when $N_\infty^2 < N_{0,\infty}$;

- if there is no stationary state, $(N_{k,\infty})_{k \ge 0}$ diverges.

* In the inhibitory case $b < 0$: there exists a critical connectivity $b^* < 0$ such that

- if $b^* < b < 0$, $(N_{k,\infty})_{k \ge 0}$ converges monotonically to the unique stationary firing rate $N_\infty$;

- if $b < b^*$, there exists $N^-$, $N^+$ such that $N_\infty \in (N^-, N^+)$ and $(N_{k,\infty})_{k \ge 0}$ converges to the 2-cycle $\{N^-, N^+\}$.

Across a wider parameter range, numerical simulations allow to check that the behaviour of the nonlinear system (NNLIF) asymptotically follows the sequence of pseudo-equilibria when the delay is large enough, $d = 10$ being enough in many cases.

In the excitatory case $b > 0$, when there is a unique stationary state $(\rho_\infty, N_\infty)$, in most cases the sequence of firing rates converges to $N_\infty$. Starting from a pseudo-equilibrium as an initial condition, the nonlinear system for positive delays (hence, no finite-time blowup) will numerically converge to the stationary state. When $b$ is large enough to have two stationary states $(\rho_\infty^1, N_\infty^1)$ and $(\rho_\infty^2, N_\infty^2)$ with $N_\infty^1 < N_\infty^2$, the asymptotic behaviour of (NNLIF) also agrees with the sequence of pseudo-equilibria:

* if the initial firing rate is above $N_\infty^2$, there is divergence in infinite time $\lim_{t \to +\infty} N(t) = +\infty$ and the density converges to a plateau state, i.e. a uniform distribution between $V_R$ and $V_F$;

* if the initial firing rate is below $N_\infty^2$, there is convergence to the stationary state with lower firing rate $(\rho_\infty^1, N_\infty^1)$.

In the case where there is no stationary state, both the sequence of pseudo-equilibria and the solutions of (NNLIF) converge towards the plateau state. Note that in the excitatory case, the delay does not matter unless it is $d = 0$ (then the technical problem of finite-time blow-up arises) and only the initial condition and the multiplicity of stationary states matter. The behaviour of the sequence is a good predictor of the nonlinear asymptotics.

In the inhibitory case $b < 0$, the initial condition is irrelevant, but the value of the delay matters. Fix, for example, $V_F = 1, V_R = 0, a \equiv 1$. For small enough $b$ and any delay $d \ge 0$, there is convergence of the firing rate sequence to the unique stationary firing rate $N_\infty$ and convergence of the nonlinear solution of (NNLIF) towards the unique stationary state. There is a critical value $b^* \simeq -9.4$ such that when $b < b^*$, the sequence of firing rates converges to a 2-cycle $\{N^-, N^+\}$ and the sequence of pseudo-equilibria to a 2-cycle $\{\rho^-, \rho^+\}$. Then, the solution of (NNLIF) follows the same behaviour only when the delay $d$ is large enough. For example, with the parameters specified, $(\rho, N)$ converges to the unique stationary state $(\rho_\infty, N_\infty)$ when $d = 2$ but to a periodic solution when $d = 10$. The only impact of the initial condition is which state in the two-cycle is reached first.

In the case of a small enough connectivity, the link between the sequence of pseudo-equilibria and the stationary state is rigorously shown in [26]. The idea is to introduce, given a solution $(\rho, N)$ to (NNLIF),

$$ \rho_{k+1}(v, t) = \rho(v, t + kd), \quad N_k(t) = N(t + kd), $$
which allows to rewrite (NNLIF) in the time interval $(kd, (k+1)d)$ as

$$ \frac{\partial \rho_k}{\partial t} + \frac{\partial}{\partial v} [(-v + bN_{k-1}(t-d))\rho_k] - \frac{\partial^2 \rho_k}{\partial v^2} = \delta_{v=V_R} N_k(t). \quad (1.31) $$

The stationary states to (1.31) correspond to the sequence of pseudo-equilibria. Assuming a spectral gap and a regularisation property on the linearised system, [26] proves under the functional analysis artillery of [25] that there exists a connectivity $b^* > 0$, a delay $d^* > 0$ and constants $Q > 0, \lambda > 0$, all depending on $V_R, V_F$ and $\sup_{t \in [0,d]} \|\rho(\cdot, t)\|_X$ such that for all $-b^* < b < b^*$ and $d > d^*$,

$$ \|\rho(\cdot, t)\|_X \le Qe^{-\mu t} \|\rho^0 - \rho_\infty\|_X. $$

Note that this result was already proved in previous literature [52, 31, 25]. However, the relationship with the sequence of pseudo-equilibria sheds new light onto the problem. Moreover, in [31], the smallness condition for $b$ was given by the constraint $N_\infty b^2 e^{\lambda d} \le C$, where $C, \lambda$ are functional analysis constants, but with the pseudo-equilibria method it is replaced by the $d$-independent condition $b^4 \le C\lambda$, where $C, \lambda$ are related to the spectral gap and other functional analysis constraints.

## 1.11 Variants of the NNLIF model

The system (NNLIF) is a simplistic representation of a large network of neurons. Some behaviours of the solutions, such as convergence to a plateau state for a large connectivity $b > 0$ and positive delay $d > 0$, are not realistic features in neuroscience. Many studies have proposed and investigated variants of the original model displaying additional features or modified equations. It allows in particular to find the minimal structures necessary to observe more realistic neuroscience patterns.

### 1.11.1 Adding a refractory state

In [29], the addition to (NNLIF) of a refractory state was rigorously studied. This models an important feature of neurons: after emitting an action potential, they cannot easily spike shortly after and have to enter a state of recovery. Mathematically, we can represent it by replacing the reset term $\delta_{v=V_R} N(t)$ by a source from an external pool of neurons $R(t)$ which are in a refractory state. Note that previous physics-oriented studies like [19, 18] had already proposed this idea and had done preliminary numerical studies. The NNLIF model with refractory period reads

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial \rho}{\partial t} + \frac{\partial}{\partial v} [(-v + bN(t))\rho] - a(N(t)) \frac{\partial^2 \rho}{\partial v^2} = \delta_{v=V_R} M(t), \\ N(t) = -a(N(t)) \frac{\partial p}{\partial v}(V_F, t), \\ \displaystyle \frac{dR}{dt}(t) = N(t) - M(t), \quad R(0) = R^0 \ge 0, \\ \rho(V_F, t) = \rho(-\infty, t) = 0, \quad \int_{-\infty}^{V_F} \rho^0(v)dv + R^0 = 1, \quad \rho(v, 0) = \rho^0(v) \ge 0. \end{array} \right. \quad (1.32) $$

The choice in [18] was for example $M(t) = N(t - \tau_0)$, meaning that neurons firing at time $t$ are ready to fire again after a fixed time $\tau_0$. In [29], the authors opt for the smoother

$$ M(t) = \frac{R(t)}{\tau}, \quad a(N) \equiv a > 0. $$

Weak and strong solutions can be defined similar to Section 1.2 for (NNLIF). The system enjoys the mass conservation property

$$ \int_{-\infty}^{V_F} \rho(v, t)dv + R(t) = \int_{-\infty}^{V_F} \rho^0(v)dv + R^0 = 1. $$

If $b > 0$, Theorem 1.12 still holds with a similar proof, which indicates that a refractory period, unlike a time delay, is not enough to prevent finite-time blow-up. There is finite time blow-up when, assuming $a(N) > a_m > 0$,
* for a fixed $\rho^0$, when $b > 0$ is large enough;

* for a fixed $b > 0$, when

$$b\mu \int_{-\infty}^{V_F} \rho^0(t)e^{\mu v}dv > e^{\mu V_F}, \quad \mu > 2\max\left(\frac{V_F}{a_m}, \frac{1}{b}\right)$$

Stationary states of (1.32) are similar to those of (NNLIF) but with $R_\infty = \tau N_\infty$. However, the classification of their existence is strikingly different.

**Theorem 1.40** (Cáceres, Perthame [29]). Assume $M(t) = \frac{R(t)}{\tau}$, $a(N) \equiv a > 0$; then

* if $b \le 0$ there is a unique stationary state of (1.32);

* if $b > 0$ there is an odd number of stationary states of (1.32).

In contrast with the standard NNLIF model, there is always at least one stationary state, even for very high connectivities. The refractory state removes the possibility of the plateau state, which was a form of singular stationary state for (NNLIF). For system (1.32), there are in practice one or three stationary states (two in special cases; see footnote). The stationary density is still described by (1.6), but the stationary firing rate is now solution of $N(I(N) + \tau) = 1$ with the same function $I$ as in (1.7).

The natural relative entropy for (1.32) given a stationary state ($\rho_\infty, N_\infty, R_\infty$) is

$$\mathcal{H}(\rho|\rho_\infty)(t) = \int_{-\infty}^{V_F} \left( \frac{\rho(v,t) - \rho_\infty(v)}{R_\infty} \right)^2 \rho_\infty(v)dv + \left( \frac{R(t) - R_\infty}{R_\infty} \right)^2 R_\infty.$$

It is used in [29] to prove exponential convergence to the unique stationary state in the linear case.

**Theorem 1.41** (Cáceres, Perthame [29]). Assume $b = 0$, $a(N) \equiv a > 0$. There exists $\mu > 0$ such that if $\mathcal{H}(\rho^0|\rho_\infty) < +\infty$, the solution to (1.32) satisfies

$$\mathcal{H}(\rho|\rho_\infty)(t) \le e^{-\mu t}\mathcal{H}(\rho^0|\rho_\infty).$$

The result was extended to the weakly-nonlinear case in [33] for solutions starting close enough to the stationary state in relative entropy and constant diffusion $a(N) \equiv a > 0$.

**Theorem 1.42** (Cáceres, Schneider [33]). Assume $|b|$ is small enough and $a(N) \equiv a > 0$. There exists $\mu > 0$ such that if $\mathcal{H}(\rho^0|\rho_\infty) \le \frac{1}{2|b|}$, the solution to (1.32) satisfies

$$\mathcal{H}(\rho|\rho_\infty)(t) \le e^{-\mu t}\mathcal{H}(\rho^0|\rho_\infty).$$

### 1.11.2 NNLIF model with random firing

As we have seen above, a positive delay coefficient $d > 0$, representing travel of action potential along axons and synaptic integration, prevents finite-time blow-up and replaces it either by a sharp peak of the firing rate followed by convergence to a stationary state or by an unphysical plateau state. A smoother way to prevent finite-time blow-up is to reduce the deterministic firing of neurons when they reach the firing potential $V_F$ by randomizing firing events. The model is then extended to the whole real line $\mathbb{R} = (-\infty, +\infty)$ and the discharge of neurons is represented by the decay term $\phi_\varepsilon(v)\rho(v, t)$, where the function $\phi_\varepsilon(v)$ is uniformly equal to zero on $(-\infty, V_F)$ and takes positive values on $(V_F, +\infty)$. To our knowledge, the first rigorous study was performed in [29], where the authors propose the model

$$\left\{
\begin{array}{l}
\displaystyle \frac{\partial \rho}{\partial t} + \frac{\partial}{\partial v} [(-v + bN(t))\rho] - a(N(t)) \frac{\partial^2 \rho}{\partial v^2} + \phi_\varepsilon(v) \rho = N(t) \delta_{v=v_R}, \quad v \in \mathbb{R}, \\[2.5ex]
N(t) = \int_{-\infty}^{+\infty} \phi_\varepsilon(v) \rho(v, t) dv, \\[2.5ex]
\displaystyle \rho(+\infty, t) = \rho(-\infty, t) = 0, \quad \int_{-\infty}^{+\infty} \rho^0(v) dv = 1, \quad \rho(v, 0) = \rho^0(v) \ge 0,
\end{array}
\right.
\qquad (1.33)$$

<footnote n="3">³Up to multiplicity. There are a finite number of values $b$ where $1/I(N)$ touches tangentially the diagonal, with only two distinct stationary states, often during transitions from 1 to 3 stationary states.</footnote>
along with a second version with refractory period:

$$ egin{cases} 
\frac{\partial \rho}{\partial t} + \frac{\partial}{\partial v} [(-v + bN(t))\rho] - a(N(t)) \frac{\partial^2 \rho}{\partial v^2} + \phi_{\varepsilon}(v)\rho = \frac{R(t)}{\tau} \delta_{v=v_R}, \quad v \in \mathbb{R}, \\ 
N(t) = \int_{-\infty}^{+\infty} \phi_{\varepsilon}(v)\rho(v,t)dv, \quad \frac{dR}{dt}(t) = N(t) - \frac{R(t)}{\tau}, \quad R(0) = R^0 \ge 0, \\ 
\rho(+\infty, t) = \rho(-\infty, t) = 0, \quad \int_{-\infty}^{+\infty} \rho^0(v)dv + R^0 = 1, \quad \rho(v, 0) = \rho^0(v) \ge 0. 
\end{cases} 
\quad (1.34) $$

Since the obstruction for global-in-time existence of solutions to (NNLIF) is finite-time divergence of the firing rate $N(t)$, the article [29] provides arguments for the absence of blow-up in (1.33) by proving locally uniform bounds for the firing rate.

**Theorem 1.43** (Cáceres, Perthame [29]). Assume $a(N) \le a_0 + a_2 N^2$ for some $a_0, a_2 \ge 0$, $\phi_{\varepsilon}(v) = \frac{1}{\varepsilon}(v - V_F)_+$ and

$$ \int_{-\infty}^{+\infty} (1 + |v|^3) \rho^0(v) dv < +\infty. $$

Then, there exists $C > 0$ such that any weak solution of (1.33) or of (1.34) satisfies

$$ N(t) \le \frac{1}{\varepsilon} \max \left( C, \int_{-\infty}^{+\infty} (1 + |v|^3) \rho^0(v) dv \right)^{\frac{1}{3}} e^{\frac{C}{3\varepsilon^2} t}. $$

The result is proved with the choice of a linear by part firing function

$$ \phi_{\varepsilon}(v) = \frac{1}{\varepsilon}(v - V_F)_+. $$

For the choice $\phi_{\varepsilon}(v) = \frac{1}{\varepsilon}1_{(V_F, +\infty)}(v)$, [29] notices the simpler direct control

$$ N(t) \le \frac{1}{\varepsilon}. $$

Numerical simulations in [29] show that solutions of (1.34) can display convergence towards a time-periodic solution, which was never the case for system (NNLIF) without delay ($d = 0$). We could think that the random firing mechanism effectively plays the role of a delay as some randomly chosen neurons will have a delayed impact on the network upon crossing the threshold $V_F$, but a notable difference is that those periodic oscillations are in the excitatory case, which the delayed NNLIF model with deterministic firing cannot reproduce. This type of periodic activity was rigorously proved to arise through a Hopf bifurcation for a hyperbolic variant of the model [70, 71], where there is no diffusion $a(N) \equiv 0$ (See (6.2) In Subsection 6.1 below).

When the parameter $\varepsilon$ tends to 0, the NNLIF model with random firing formally converges towards the standard model (NNLIF). This limit was rigorously proved in the linear case $b = 0, d = 0, a(N) \equiv a > 0$ in [165], using the iteration scheme method exposed above in Section 1.8.2. The authors also show with numerical simulations that the convergence should hold in the nonlinear case $b \ne 0$. Rigorous convergence in the nonlinear case was then proved rigorously in [90], allowing to construct generalised solutions to (NNLIF) (See Subsection 1.9.2 above).

### 1.11.3 A system of excitatory-inhibitory populations

In [18], a model was proposed for two interacting populations of neurons, one excitatory and the other inhibitory. In many areas of the brain, distinct interlaced populations play excitatory and inhibitory roles, like inhibitory interneurons among pyramidal cells [77]. The excitatory-inhibitory NNLIF system
involves two probability densities $\rho_I(v, t)$ and $\rho_E(v, t)$ solving

$$ egin{cases} 
\frac{\partial \rho_I}{\partial t} + \frac{\partial}{\partial v} [h^I(v, N_E, N_I) \rho_I] - a_I(N_E, N_I) \frac{\partial^2 \rho_I}{\partial v^2} = N_I(t) \delta_{v_R}, \\ 
\frac{\partial \rho_E}{\partial t} + \frac{\partial}{\partial v} [h^E(v, N_E, N_I) \rho_E] - a_E(N_E, N_I) \frac{\partial^2 \rho_E}{\partial v^2} = N_E(t) \delta_{v_R}, \\ 
N_\alpha(t) = -a_\alpha(N_E(t), N_I(t)) \frac{\partial \rho_\alpha}{\partial v}(V_F, t), \quad \alpha = E, I \\ 
\rho_\alpha(V_F, t) = \rho_\alpha(-\infty, t) = 0, \quad \rho_\alpha(v, 0) = \rho_\alpha^0(v), \quad \int_{-\infty}^{V_F} \rho_\alpha^0(v) dv = 1, \quad \alpha = E, I 
\end{cases} 
\quad (1.35) $$

with

$$ 
\begin{aligned} 
h^\alpha(v, N_E, N_I) &= -v + b_E^\alpha N_E - b_I^\alpha N_I, \\ 
a_\alpha(N_E, N_I) &= a_0^\alpha + a_E^\alpha N_E + a_I^\alpha N_I, \quad \alpha = E, I 
\end{aligned} 
$$

where $b_E^E, b_E^I, b_I^E, b_I^I, a_E^E, a_I^E, a_0^E$ and $a_0^I$ are all non-negative parameters. The first rigorous study of such a model was proposed in [32].

The action of the inhibitory population, even when the inhibitory coefficients are large is not necessarily enough to prevent entirely finite time blow-up. Under an *a priori* assumption of the form $\int_0^t N_I(s)ds \le Ct$, [32] proves that as soon as $b_E^E > 0$, finite-time blowup can occur for initial conditions where $\rho_E^0$ is concentrated around $V_F$. In unpublished results one can find in [232], this technical assumption was removed and replaced by assumptions on coefficients only: $a_E$ and $a_I$ positive constant functions, $b_E^E > 0$ and $b_E^I > \frac{1}{V_F-V_R}b_E^Ib_E^E$. However, these results do not offer a complete understanding on finite-time blow-up in excitatory-inhibitory NNLIF systems as many simple at hand counter-example prove that the conditions are not sharp. The interplay between the parameters and the two initial conditions makes any attempt at a general classification challenging.

Regarding stationary states, for $\alpha = E, I$, the stationary densities are of the form

$$ \rho_{\alpha, \infty}(v) = \frac{N_{\alpha, \infty}}{a_{\alpha}(N_{E, \infty}, N_{I, \infty})} e^{-\frac{(v-b_E^\alpha N_E + b_I^\alpha N_I)^2}{2a_\alpha(N_{E, \infty}, N_{I, \infty})}} \int_{\max(v, V_R)}^{V_F} e^{-\frac{(w-b_E^\alpha N_E + b_I^\alpha N_I)^2}{2a_\alpha(N_{E, \infty}, N_{I, \infty})}} dw. $$

where $N_{\alpha, \infty}$ satisfies

$$ \frac{1}{N_{\alpha, \infty}} = \int_{0}^{+\infty} \frac{e^{-\frac{s^2}{2}}}{s} \left( e^{\frac{s \sqrt{V_F - b_E^\alpha N_E + b_I^\alpha N_I}}{\sqrt{a_\alpha(N_{E, \infty}, N_{I, \infty})}}} - e^{\frac{s \sqrt{V_R - b_E^\alpha N_E + b_I^\alpha N_I}}{\sqrt{a_\alpha(N_{E, \infty}, N_{I, \infty})}}} \right) ds. $$

It is possible to provide a partial classification.

**Theorem 1.44** (Cáceres, Schneider [32]). Assume $b_I^E, b_E^I > 0$, and $a_E, a_I$ are constant positive functions. Then for $\alpha = E, I$,

* if $(V_F - V_R)^2 < (V_F - V_R)(b_E^E - b_I^I) + b_E^E b_I^I - b_I^E b_E^I$, either there is no stationary state of (1.35) or an even number of stationary states; if $b_E^E$ is large enough with respect to other parameters, there is no stationary state;

* if $(V_F - V_R)(b_E^E - b_I^I) + b_E^E b_I^I - b_I^E b_E^I < (V_F - V_R)^2$, there is an odd number of stationary states of (1.35); if $b_E^E$ is small enough, there is a unique stationary state.

As this result and numerical examples in [32] show, there can be one, two, three or no stationary states; the two population system exhibits a more complicated picture than the standard model (NNLIF).

Given a stationary state $(\rho_{E,\infty}, \rho_{I,\infty}, N_{E,\infty}, N_{I,\infty})$, a natural relative entropy for asymptotic convergence is

$$ \mathcal{H}[\rho_E, \rho_I](t) = \int_{-\infty}^{V_F} \left( \frac{\rho_E(v, t) - \rho_{E,\infty}(v)}{\rho_{E,\infty}(v)} \right)^2 \rho_{E,\infty}(v) dv + \int_{-\infty}^{V_F} \left( \frac{\rho_I(v, t) - \rho_{I,\infty}(v)}{\rho_{I,\infty}(v)} \right)^2 \rho_{I,\infty}(v) dv, $$
**Theorem 1.45** (Cáceres, Schneider [32]). Assume $b_E^E, b_E^I, b_I^E, b_I^I$ are all small enough, $a_E, a_I$ are constant positive functions. Then, if

$$ \mathcal{H}[\rho_E, \rho_I](0) \le \frac{1}{2 \max(b_E^E + b_I^E, b_E^I + b_I^I)}, $$

there exists $\mu > 0$ such that the solution of (1.35) satisfies,

$$ \forall t \ge 0, \quad \mathcal{H}[\rho_E, \rho_I](t) \le e^{-\mu t} \mathcal{H}[\rho_E, \rho_I](0). $$

### 1.11.4 Combination of excitatory-inhibitory populations, refractory state and delay

In [33], a more realistic NNLIF-type model is studied theoretically and numerically. It combines separate excitatory and inhibitory populations, a refractory state and transmission delays. The model writes, for $\alpha = E, I$,

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial \rho_{\alpha}}{\partial t} + \frac{\partial}{\partial v} [h^{\alpha} \rho_{\alpha}] - a_{\alpha} \frac{\partial^2 \rho_I}{\partial v^2} = \frac{R_{\alpha}(t)}{\tau_{\alpha}} \delta_{v_R}, \\[1.5em] \displaystyle N_{\alpha}(t) = -a_{\alpha} \frac{\partial \rho_{\alpha}}{\partial v}(v_F, t), \\[1.5em] \displaystyle \frac{d R_{\alpha}}{dt}(t) = N_{\alpha}(t) - \frac{R_{\alpha}(t)}{\tau_{\alpha}}, \qquad R_{\alpha}(0) = R_{\alpha}^0 \\[1.5em] \displaystyle \rho_{\alpha}(v_F, t) = \rho_{\alpha}(-\infty, t) = 0, \quad \rho_{\alpha}(v, 0) = \rho_{\alpha}^0(v), \qquad \int_{-\infty}^{v_F} \rho_{\alpha}^0(v) dv + R_{\alpha}^0 = 1, \end{array} \right. \quad (1.36) $$

with

$$ \begin{aligned} h^{\alpha}(v,t) &= -v + b_E^{\alpha} N_E(t - d_E^{\alpha}) - b_I^{\alpha} N_I(t - d_I^{\alpha}), \\ a_{\alpha}(t) &= a_0^{\alpha} + a_E^{\alpha} N_E(t - d_E^{\alpha}) + a_I^{\alpha} N_I(t - d_I^{\alpha}), \quad \alpha = E, I \end{aligned} $$

where $b_E^E, b_I^E, b_I^E, a_E^E, a_I^E, a_0^E, a_0^I, d_E^E, d_I^E, d_I^E$ and $d_I^I$ are all non-negative parameters, and $\tau_E, \tau_I > 0$ are positive parameters. It is of course also possible to replace the smooth refractory input $\frac{R_{\alpha}(t)}{\tau_{\alpha}}$ by $N_{\alpha}(t - \tau_{\alpha})$. The mass of the density plus the refractory population is conserved for $\alpha = E, I$.

The study of stationary states proceeds similarly to the other NNLIF-type models. Adding a refractory state to the excitatory-inhibitory population prevents the case where no stationary state can exist.

**Theorem 1.46** (Cáceres, Schneider [33]). Assume $b_E^I, b_E^E > 0$ and $a_E, a_I$ are positive constant functions. Then, there is always an odd number$^4$ of stationary states to (1.36). If $b_E^E$ is small enough, or $\tau_E$ is large enough, there is a unique stationary state.

After proving Theorem 1.42, [33] proves an equivalent result for the general system (1.36) in the non-delayed case $d_E^E = d_I^E = d_I^I = 0$, using an appropriate relative entropy functional and adequate smallness assumptions on coefficients and initial conditions. The finite time blow-up result proved for system (1.35) is also extended to (1.36) in the case $d_E^E = 0$.

Note that in the case where all delay coefficients are positive, it is possible to rewrite (1.36) in a form that allows to extend the existence proof of [31], thus proving global-in-time existence of smooth solutions. More generally, the construction of local-in-time solutions of (1.36), although not done yet, would be a (very) technical extension of [48] whose main benefit would be the procurement of representation formulae. The general study of the properties of (1.36) is still mostly open.

On the numerical side, [33] confirms that this model also exhibits time-periodic solutions in the delayed case with strong inhibition, like for (NNLIF), but also time-periodic solutions in the average excitatory case. A complete picture of the influence of each of the numerous parameters in (1.36) on the occurrence of self-sustained oscillations is currently lacking.

<footnote>$^4$Like mentioned before, up to a finite number of values of the parameters where it is up to multiplicity only.</footnote>
### 1.11.5 NNLIF with learning rule

In all NNLIF-type models we have described so far, the connectivity parameter is fixed throughout the time evolution. Following works on models for learning among a finite number of neurons [122, 121, 117, 118], attempts were made in [213] at describing learning in a large (mathematically infinite) network of neurons. The goal was to answer some questions inspired from [138]:

* can any pattern of neural activity be generated by a heterogeneous synaptic weights distribution?

* what is the equilibrium synaptic weights distribution for given learning rule and external inputs?

* can the network remember an external signal presented through a stationary weight distribution?

In order to investigate those questions, [213] proposes two NNLIF-like models. First, a model structured by voltage $v \in (-\infty, V_F]$ and synaptic weights $w \in \mathbb{R}$,

$$ 
\left\{ 
\begin{array}{l}
\displaystyle \frac{\partial \rho}{\partial t} + \frac{\partial}{\partial v} [(-v + J(w) + w\sigma(\bar{N}(t)))\rho] - a \frac{\partial^2 \rho}{\partial v^2} = \delta_{v=V_R} N(w, t), \\
\\
\displaystyle N(w, t) = -a \frac{\partial \rho}{\partial v}(V_F, w, t), \quad \bar{N}(t) = \int_{-\infty}^{+\infty} N(w, t) dw, \\
\\
\displaystyle \rho(V_F, w, t) = \rho(-\infty, w, t) = 0, \quad \rho(v, w, 0) = \rho^0(v, w) \ge 0, \\
\\
\displaystyle \int_{-\infty}^{V_F} \int_{-\infty}^{+\infty} \rho^0(v, w) dv dw = 1,
\end{array}
\right. 
\qquad (1.37) 
$$

where $a > 0$ is a constant, $\sigma \in C^2(\mathbb{R}_+, \mathbb{R}_+)$, $\sigma' \ge 0$, $\sup_{N \ge 0} \sigma(N) = \sigma_m < +\infty$. The neurons interact through the mean-firing rate $\bar{N}$ and the distribution along the connectivity parameter $w$. For a same average firing rate, neurons with large positive $w$ will receive strong excitation and neurons with a large negative $w$ will receive strong inhibition. The function $J$ is an external input. In this model, the synaptic distribution cannot change. The synaptic distribution can be represented by the function

$$ H(w) = \int_{-\infty}^{V_F} \rho(v, w, t) dv = \int_{-\infty}^{V_F} \rho^0(v, w) dv. $$

Note that

$$ \int_{-\infty}^{+\infty} H(w) dw = 1. $$

Second, they propose a kinetic model of the form

$$ 
\left\{ 
\begin{array}{l}
\displaystyle \frac{\partial \rho}{\partial t} + \frac{\partial}{\partial v} [(-v + J(w) + w\sigma(\bar{N}(t)))\rho] + \varepsilon \frac{\partial}{\partial w} [(\Phi - w)\rho] - a \frac{\partial^2 \rho}{\partial v^2} = \delta_{v=V_R} N(w, t), \\
\\
\displaystyle N(w, t) = -a \frac{\partial \rho}{\partial v}(V_F, w, t), \quad \bar{N}(t) = \int_{-\infty}^{+\infty} N(w, t) dw, \\
\\
\displaystyle \rho(V_F, w, t) = \rho(-\infty, w, t) = 0, \quad \rho(v, w, 0) = \rho^0(v, w) \ge 0, \\
\\
\displaystyle \int_{-\infty}^{V_F} \int_{-\infty}^{+\infty} \rho^0(v, w) dv dw = 1.
\end{array}
\right. 
\qquad (1.38) 
$$

This model allows for a change in $H$ over time. The scaling parameter $\varepsilon$ allows for a slow change in the weights distribution compared to the evolution of voltage in neurons. The learning rule is defined by the function $\Phi$. One possible choice is the celebrated Hebbian learning [134, 121]: the strength of the connectivity between two neurons increases when they have high activity simultaneously, which can be represented in the context of (1.38) by

$$ \Phi(w, t) = \bar{N}(t) N(w, t) K(w), $$

for some learning kernel $K(\cdot)$. The idea is that if the specific activity $N(w, t)$ of neurons with connectivity $w$ is high and at the same time the average activity represented by $\bar{N}(t)$ is high, then $\Phi$ is high, driving neurons towards stronger connectivities via a flux towards positive $w$.
Another, more general choice could be

$$ \Phi(w, t) = N(w, t) \int_{-\infty}^{+\infty} K(w, \omega) N(\omega, t) d\omega, $$

where $K(\cdot, \cdot)$ is a learning kernel, or other choices inspired from spike timing dependent plasticity (STDP) rules [76, 40, 173, 118].

**Theorem 1.47** (Perthame, Salort, Wainrib [213]). Assume $J \in L^\infty(\mathbb{R})$ and

$$ \int_{0}^{+\infty} w^2 H(w) dw < +\infty, $$

then there exists at least one stationary state of (1.37). If $H$ is supported on $(-\infty, 0]$, there is a unique stationary state.

From a stationary state $\rho_\infty$ with firing rate $N_\infty$, we can define the normalised output signal

$$ S(w) = \frac{N_\infty(w)}{N_\infty}. $$

Given a reasonable output signal $S(w)$, [213] proves that it is always possible to find a synaptic distribution $H(w)$ which allows the network to produce the output $S(w)$.

**Theorem 1.48** (Perthame, Salort, Wainrib [213]). Assume $J \in L^\infty(\mathbb{R}) \cap C^1(\mathbb{R})$ and $S \in L^1(\mathbb{R})$ satisfies $\int_{-\infty}^{+\infty} S(w)dw = 1$ and $\int_{-\infty}^{0} S(w)e^{-\gamma w}dw < +\infty$ for some $\gamma > \frac{\sigma_m}{a}(V_F - V_R)$. Then there exists a weigth distribution $H$ and a stationary state of (1.37) with normalised output signal $S(w)$. If $S$ is supported in $(-\infty, 0]$, then $H$ is unique.

A precise notion of discrimination property is also given and studied theoretically: from two input signals $J_1(w)$ and $J_2(w)$, we have under reasonable assumptions

$$ \int_{-\infty}^{+\infty} |N_1(w) - N_2(w)| dw \geq C \int_{-\infty}^{+\infty} L(w) N(w) |J_1(w) - J_2(w)| dw, $$

with $L > 0$ a well-chosen constant and $C > 0$ a constant depending on $H$, $J_1$, $J_2$ and other parameters.

For the kinetic variant (1.38), [213] proves that a given input $J$ can lead to many different stationary weight distributions $H$ after learning.

**Theorem 1.49** (Perthame, Salort, Wainrib [213]). Assume $J \in L^\infty(\mathbb{R}) \cap C^1(\mathbb{R})$, $wK(w) > 0$ for $w \neq 0$ and $0 < k_m \leq K(w) \leq K_M$ almost everywhere. Then there are infinitely many different stationary states of (1.38).

It is however possible to describe some aspects of these stationary states. In specific cases, [213] gives a way to organise and select them, by adding a Gaussian noise $\varepsilon\nu \frac{\partial^2 \rho}{\partial w^2}$ in the system in the fully inhibitory context. Then, they select a stationary state of (1.38) in the small noise limit $\nu \to 0$.

Numerical tests are then proposed by letting the system learn from an input $J_1$ in (1.38), then using the corresponding stationary weigth distribution $H$ as an initial condition in (1.37) subject to a different input $J_2$. The idea is to let a system learn the input $J_1$ on a long timescale and then, on a shorter timescale, to let it react to the other input $J_2$. Further numerics with a carefully crafted scheme have been proposed in [133].
### 1.11.6 Poissonian mean-field variant

In all the variants of the NNLIF model we have described above, the idea was to complexify the model in order to describe more realistic properties of neural networks. In two companion papers [252, 236], a simpler variant is proposed with the idea of doing more explicit computations and achieving a better characterisation of finite time blow-up. They assume that neurons are driven by Poisson-like noisy inputs and enter a refractory period after spiking, leading to what they name a delayed Poissonian mean-field dynamics (dPMF).

$$ 
\left\{ 
\begin{array}{l}
\displaystyle \frac{\partial \rho}{\partial t} = (\nu + bN(t)) \left( \frac{\partial \rho}{\partial v} + \frac{1}{2} \frac{\partial^2 \rho}{\partial v^2} \right) + N(t - \varepsilon_0) \delta_{v=v_R} \\[1.5em]
\displaystyle N(t) = \frac{1}{2} \frac{\partial \rho}{\partial v}(0, t), \\[1.5em]
\displaystyle \rho(0, t) = \rho(+\infty, t) = 0, \quad \rho(v, 0) = \rho^0(v) \ge 0, \quad \int_0^{+\infty} \rho^0(v) dv = 1,
\end{array}
\right. 
\qquad (\text{dPMF})
$$

with $b \in \mathbb{R}$ the connectivity parameter, $\nu > 0$ a constant drift and $\varepsilon_0 > 0$ the length of the refractory period. Note that they choose $0 = V_F < V_R$ and study the model in the non-negative half space $v \in [0, +\infty)$ for simplicity. For this system, the firing rate can be written as

$$ N(t) = \frac{\nu \partial_v \rho(0, t - \varepsilon_0)}{2 - b \partial_v \rho(0, t)} $$

which shows that as (NNLIF) with internal noise $a(N) = a_0 + a_1 N$, finite-time blow-up happens with a finite boundary derivative $\frac{1}{2} \partial_v \rho(V_F, T^*) = \frac{1}{b}$. Their method is to introduce the time change

$$ \tau = \Phi(t) = \nu + bF(t), \qquad F(t) = \int_{0}^{t} N(s)ds, $$

which allows to dilate the blow-up and continue the dynamics within. This idea is quite similar to the strategy of [91, 90] although developed with completely different methods and, up to our best knowledge, independently and simultaneously. After the time dilation, [91, 90] use more deterministic PDE tools while [252, 236] use more probabilistic methods. This idea was later used on a hyperbolic PDE in [46].

In the case of the model (dPMF), the time change allows to get precise information about the blow-up dynamics and to perform a vanishing refractory period limit $\varepsilon_0 \to 0$. A notable byproduct is a characterisation of the divergence speed for $N(t)$.

**Theorem 1.50** (Taillefumier, Whitman [252]). Under a full blow-up condition, blowing-up solutions to (dPMF) satisfy at the blow-up time $T^*$,

$$ N(t) \sim \frac{b}{\sqrt{2C(T^* - t)}}, $$

for some constant $C > 0$ which can be characterised.

The method also allows to characterise the mass of neurons participating into a full blow-up event. In the companion article [236], it is proved that there are global generalised solutions with infinitely many full blow-up events.

### 1.11.7 NNLIF-like models in finance modelling

Many researchers in financial mathematics [131, 162, 191, 192, 193, 129, 83] consider probabilistic NNLIF-type models in the form of McKean-Vlasov equations depending on singular hitting times; they have developed their own terminology and set of techniques. Although a careful description of their methods and results is beyond the scope of this review, we want to give a broad idea of the kind of
models they work with, provide some references, and briefly mention some results they have obtained that can be extended to the NNLIF model with enough technical work.

Default contagion in large financial networks, bank networks, or in a large portfolio can be modelled by the distance to default $X_t \in [0, +\infty)$, 0 being the default threshold. In the mean-field limit, the most simple probabilistic formulation [130] is

$$ egin{cases} X_t = X_0 + B_t - bL_t \\ \tau = \inf\{t \ge 0 \mid X_t \le 0\}, \\ L_t = \mathbb{P}(\tau \le t) \end{cases} $$

whose law is a solution to the PDE

$$ egin{cases} \frac{\partial \rho}{\partial t} - bN(t)\frac{\partial \rho}{\partial x} - \frac{1}{2}\frac{\partial^2 \rho}{\partial x^2} = 0 \\ N(t) = \frac{1}{2}\frac{\partial \rho}{\partial x}(0, t), \\ \rho(0, t) = 0, \quad \rho(x, 0) = \delta_{x=X_0}. \end{cases} $$

Many more realistic financial models are built on the same basis and mean-field limits have been proved. Contrarily to neurons, the financial agents, banks or portfolio entities are not resurrected to a reset value after defaulting. However, like the NNLIF model without delay, these finance models are prone to finite-time blow-up and a notion of *physical solution* can be constructed like in [82].

The size of the jumps in physical solutions is of a particular importance for finance applications as it is indicative of the severity of, for example, a financial crisis. It can be proved that the *physical solutions* have minimal jumps after blow-up [130, Theorem 1.2] and the article makes progresses towards a proof of uniqueness for the *physical solutions* [130, Theorem 1.8, discussion therein]. Recently, more methods were developed to extend solutions of similar stochastic equations in spite of the presence of blow-up [83, 162, 193]. This community also performed a careful study of how the decay of the solution near the boundary 0 has an impact on the regularity after blow-up events.

## 1.12 Numerical methods for NNLIF type systems

Several numerical methods have been proposed for approximating NNLIF equations and their variants. The first attempts by physicists and applied mathematicians have been applying Monte-Carlo methods based on the stochastic processes underlying these PDE models, see for instance [245, 220, 35, 224, 36, 34, 222, 196, 197]. The first numerical results by deterministic schemes of the NNLIF model (1.4)-(1.5) and their kinetic conductance-voltage counterpart in [27, 28] used explicit in time finite difference methods in voltage and conductance based on *Weighted Essentially Non-Oscillatory* (WENO) schemes [246] originating in conservation laws [247, 144]. They approximate the advection terms to high order in order to simulate efficiently and accurately large gradient regions in voltage or conductance that show up due to the blow-up of solutions for large connectivity. Second-order terms are discretized by centered finite differences while time discretization is performed via 3rd-order *Total Variation Diminishing Runge-Kutta* schemes. The resulting schemes have been very successful in capturing blow-up phenomena, kink formation and convergence to stationary states for both the NNLIF (1.4)-(1.5) without delay $d = 0$ [27] and the kinetic conductance-voltage model [28]. Moreover, they have been widely used to analyse the NNLIF (1.4)-(1.5) with delay in [32], the two excitatory-inhibitory population model (1.36) in [33], see Subsection 1.11.4, and the sequence of pseudo-equilibria describing the long-time behaviour of the NNLIF model (1.4)-(1.5) with large delay in [26], see Subsection 1.10.2. We refer to the corresponding previous subsections for details on the interpretation of the numerical results related to stability of steady states, periodic solutions and synchronization of the network.

A second family of deterministic numerical methods proposed in the literature are based on structure preserving methods, meaning schemes that preserve some properties of the evolution of the problem. In the case of the NNLIF model, the basic property to preserve at least in the linear case, $b = 0$, it is the
decay of the relative entropy to equilibrium. This family of methods based on finite difference/finite volume schemes for Fokker-Planck type equations can be found with the name Chang-Cooper method, see for instance [20], or Scharfetter-Gummel method, see for instance [145]. Both methods use the following rewriting of the Fokker-Planck equation (1.4) in terms of the Maxwellian

$$M(v, t) = e^{-\frac{(v-bN(t))^2}{2a(N(t))}}$$

as follows,

$$\frac{\partial p}{\partial t}(v,t) - a(N(t)) \frac{\partial}{\partial v} \left[ M(v,t) \frac{\partial}{\partial v} \left( \frac{p(v,t)}{M(v,t)} \right) \right] = N(t) \delta(v - V_R).$$

Then, the Chang-Cooper method performs a kind of $\theta$-finite difference approximation of $p/M$, see [20] for details, while the Scharfetter-Gummel method does a harmonic average, see [139] for details. Both methods have been used in [27, 28, 139, 133] for a variety of the computational neuroscience problems discussed in the previous sections, showing some of the conservation and dissipation properties of these schemes. In comparison with the WENO schemes, they are lower order, however they are able to keep known dissipation properties in the linear case $b = 0$ without network connectivity. Other numerical approaches have been introduced such as spectral methods [274], multiscale methods [94] mixing stochastic microscopic schemes with deterministic macroscopic schemes to analyse synchronization in NNLIF, and discontinuous-Galerkin methods [244, 243].

## 2 Fokker-Planck representation of a stochastic neural field and the modelling of grid cells

We present here a PDE framework recently developed for the study of the asymptotic behaviour and phase transitions in stochastic neural fields, with particular interest to its applications to patterning in noisy grid cells.

We first explain the main application goal of this family of models in the context of previous modelling attempts on grid cells. Then we describe the rigorous mean-field limit which gave rise to these systems. Finally, we introduce the mathematical setting and propose an overview of mathematical results on the solutions: stationary states and their bifurcations, existence and regularity of solutions, asymptotic behaviour.

### 2.1 Modelling the hexagonal firing pattern of grid cells

Grid cells, discovered in 2005 [127], have been described as the GPS in the navigation system for mammals. In fact, the hexagonal patterns they create play an important role in the integration of the moving paths of animals (including humans [86]) as demonstrated for mammalian brains in [189]. The understanding of their precise behaviour, see [234, 189] and the references therein, is still a matter of modern research. One of the main challenges highlighted in [234], was the effect of noise on grid cells and the robustness of these hexagonal patterns. The grid cell network is commonly described by deterministic continuous attractor network dynamics through a system of neural field models [101]. In order to better understand the formation of the hexagonal firing pattern and the role of grid cells in path integration, modelling studies have tried to represent in a simple framework how patterned activity arises through autonomous neuronal interactions and external inputs of speed and orientation.

Around 2009, mathematical models were designed to explain the patterns of grid cells, among which the works of Burak and Fiete, [22, 108]. In order to stay consistent with assumptions and notation, we do not present exactly their model but a slightly modified version. Consider neurons arranged in a two-dimensional sheet at locations $(x_1, \dots, x_N)$ on a square $Q = [0, L] \times [0, L]$ of length $L$. Each space location contains four neurons with orientation preference $\beta \in \{1, 2, 3, 4\}$, where these numbers represent North, West, South, East (all preference angles in $[0, 2\pi]$ are possible in reality; the model
uses only four directions for simplicity). The electrical activity (how much it fires per unit time) $s_i^\beta$ of a neuron at location $x_i$ with orientation preference $\beta$ evolves according to the differential equation

$$ \tau \frac{ds_i^\beta}{dt} = \Phi \left( \sum_{\beta'=1}^{4} \sum_{j=1}^{N} W_{i,j}^{\beta'} s_j^{\beta'} + B^\beta(t) \right) \quad (2.1) $$

where the connectivity parameter between locations $x_i$ and $x_j$ given the orientation preference $\beta$ is $W_{i,j}^\beta = W(x_i - x_j - r^\beta)$, computed according to the connectivity kernel $W$ (a function), periodically extended outside of the square domain $Q$, and the shift $r^\beta$ (a two-dimensional vector). The modulation function $\Phi$ can be thought of as a filter, like $\Phi(z) = \max(z, 0)$ (the so-called ReLU function), or a logistic / sigmoid function like $\phi(z) = (1 + e^{-x})^{-1}$. The characteristic response time for neurons is $\tau \approx 10$ms.

The input from the external world or other neuronal cells related to navigation, such as place or border cells, is contained in the terms $B^\beta$, which in turn depend on the velocity $v(t)$ and the orientation $\theta(t)$ of the animal. We consider here the simple form

$$ B^\beta(t) = B + \alpha v(t) \cos(\theta - \theta^\beta), $$

used in [22], where $B, \alpha \in [0, \infty)$ and $(\theta^1, \theta^2, \theta^3, \theta^4) = (\frac{\pi}{2}, \pi, \frac{3\pi}{2}, 2\pi)$. This setting allowed to reproduce the firing pattern of grid cells and to explain how they can be used by the brain for path integration; see [22, 72].

## 2.2 Rigorous mean-field limit towards a new model

Noise plays a role in many emerging brain phenomena and is more a feature than an accidental element of the way neural networks perform tasks; see [229] for an overview. In general, studies on understanding noisy neural field models have been lacking [13, Sec. 6] until very recently [149, 148, 154, 155, 171, 21, 260, 259, 14, 23, 103]. However, theoretical foundations are still sparse. In [42, 50], Brownian fluctuations were added to the attractor model (2.1), in order to investigate how robust the grid cell pattern is to noise and how it appears as a strong noise decreases to lower levels.

Consider $N$ space locations $x_1, \dots, x_N \in \mathbb{T}^d := \mathbb{R}^d/(L\mathbb{Z})^d$ of the neural cortex. The model is stated for a general dimension $d$ because it does not impact the results, but $d=2$ is taken for the application to grid cells. Note that $d=1$ yields interesting insights for the celebrated ring model [4]. The choice of a torus for $\mathbb{T}^d$ for the (cortical) spatial domain is motivated by recent insights from topological data analysis [120], which indicate that the electrical activity of a grid cell module can be represented by a toroidal topology. At each location, we now consider not just one neuron but a cortical column of $4M$ neurons ($M$ per orientation preference). The $k^{th}$ neuron at location $x_i$ with orientation preference $\beta \in \{1, 2, 3, 4\}$ has an activity level $s_{ik}^\beta(t)$ at time $t$, which evolves according to the coupled system

$$ \tau s_{ik}^\beta(t) = \tau s_{ik}(0) + \int_0^t \left[ -s_{ik}^\beta(\tau) + \Phi \left( \frac{1}{4NM} \sum_{\beta'=1}^4 \sum_{j=1}^N \sum_{l=1}^M W_{i,j}^{\beta'} s_{jl}^{\beta'}(\tau) + B^\beta(\tau) \right) \right] d\tau + \sqrt{2\sigma} B_{ik}^\beta(t) - \ell_{ik}^\beta(t), \quad (2.2) $$

$$ \ell_{ik}^\beta(t) = -|\ell_{ik}^\beta|(t), \quad |\ell_{ik}^\beta|(t) = \int_0^t 1_{\{s_{ik}^\beta(\tau)=0\}} d|\ell_{ik}^\beta|(\tau). \quad (2.3) $$

where the Brownian motions $(B_{ik}^\beta)_{\beta=1,...,4}$ can be correlated (see below).

For each $i, k$ and $\beta$, the reflection term $\ell_{ik}^\beta$ prevents the activity level $s_{ik}^\beta$ from taking negative values. It is equal for all time $t$ to the opposite of its total variation $\ell_{ik}^\beta(t) = -|\ell_{ik}^\beta|(t)$ which stays constant when $s_{ik}^\beta > 0$ and increases as $\int_0^t 1_{\{s_{ik}^\beta(\tau)=0\}} d|\ell_{ik}^\beta|(\tau)$ when $s_{ik}^\beta = 0$; the result is that $s_{ik}^\beta$ is pushed away from
the boundary value $s = 0$. Proving existence and uniqueness of this kind of term is a Skorokhod problem. Such reflection terms are rigorously constructed for the model (2.2)–(2.3) in [42], applying the strategy of Lions and Sznitman [164, 250].

Denote $\rho_{N,M}^{\beta}$ the empirical measure associated to the $4MN$ particles, that is

$$ \rho_{N,M}(t, dy, ds) = \frac{1}{NM} \sum_{j=1}^{N} \sum_{m=1}^{M} \delta_{(x_j, s_{jm}(t))} \quad \text{regarded as a measure on } \mathbb{T}^d \times \mathbb{R}^4. $$

where $s_{jm}(t) = (s_{jm}^1, s_{jm}^2, s_{jm}^3, s_{jm}^4) \in \mathbb{R}^4$.

We are then interested in the mean-field limit $M \to +\infty$, $N \to +\infty$. However, this is not a classical mean field problem with $4MN$ particles in $\mathbb{R}$ or $MN$ particles in $\mathbb{R}^4$. As noted in [42, 63] Its structure is more a collection of $N$ mean-field systems (the cortical columns) that interact with each other, each containing $M$ particles. For this modelling reason and for mathematical reasons related to measurability in $x$, the four-dimensional Brownian motions $B_{ik}$ and the initial conditions $s_{ik}^{\beta}(0)$ are not independent. There will be independence only between the $M$ neurons at the same location $x_j$ but there are correlations between initial conditions and noise in two different space locations $x_i$ and $x_j$. The authors of [42] show that it results in the convergence rate

$$ \frac{1}{\sqrt{M}} + \frac{1}{N^{\alpha}}, $$

for some positive $\alpha$ in the mean-field limit $M, N \to +\infty$ instead of the decay rate $\frac{1}{\sqrt{MN}}$ classicaly obtained with $MN$ particles.

Assuming the points $(x_1, \dots, x_N)$ are adequately chosen in the torus $\mathbb{T}^d$, and under other technical assumptions, it is shown in [42] that, up to technical modifications due to the measurability problems and the correlations, the system at each location $x \in \mathbb{T}^d$ converges to independent copies evolving according to the stochastic differential equation

$$ \begin{aligned} \tau s^{\beta}(x,t) &= s^{\beta}(x,0) + \int_{0}^{t} [-s^{\beta}(x,\tau) + \Phi^{\beta}(x,\tau)] d\tau + \sqrt{2\sigma} \mathcal{B}(x,t)^{\beta}(t) - \ell^{\beta}(x,t), \\ \Phi^{\beta}(x,t) &= \Phi \left( \frac{1}{4} \sum_{\beta'=1}^{4} \int_{\mathbb{T}^d} W^{\beta'}(x-y) \int_{\mathbb{R}^4} s^{\beta'}(x,y,t) \rho(t,y,ds) dy + B^{\beta}(t) \right), \\ \ell^{\beta}(x,t) &= -|\ell^{\beta}(x,\cdot)|(t), \quad |\ell^{\beta}(x,\cdot)|(t) = \int_{0}^{t} \mathbb{1}_{\{s^{\beta}(x,\tau)(\tau)=0\}} d|\ell^{\beta}(x,\cdot)|(\tau). \end{aligned} $$

where we denote $\rho(t, y, ds) = \text{Law}_{\mathbb{R}^4}(s(y, t))$, which induces the probability density $\rho(t, dx, ds)$ on $\mathbb{T}^d \times \mathbb{R}^4$ defined by

$$ \forall \phi \in C_b(\mathbb{T}^d \times \mathbb{R}^4), \quad \int_{\mathbb{T}^d \times \mathbb{R}^4} \phi(x, s) \rho(t, dx, ds) = \int_{\mathbb{T}^d} \int_{\mathbb{R}^4} \phi(x, s) \rho(t, x, ds) dx. $$

The marginal distributions $\rho^{\beta}(x, s, t) = \text{Law}_{\mathbb{R}}(s^{\beta}(x, t))$ are solution to the 4-population system of partial differential equations, for $\beta \in \{1, 2, 3, 4\}$,

$$ \tau \frac{\partial \rho^{\beta}}{\partial t} = -\frac{\partial}{\partial s} \left( [\Phi^{\beta}(x, t) - s] \rho^{\beta} \right) + \sigma \frac{\partial^2 \rho^{\beta}}{\partial s^2}, \quad (\text{GC}) $$

where $\Phi^{\beta}(x, t)$ is given by

$$ \Phi^{\beta}(x,t) = \Phi \left( \frac{1}{4} \sum_{\beta'=1}^{4} \int_{\mathbb{T}^d} W^{\beta'}(x-y) \int_{0}^{\infty} s \rho^{\beta'}(y, s, t) ds dy + B^{\beta}(t) \right), $$
and the system is supplemented with the nonlinear boundary condition

$$ \left. \left( \Phi^{\beta}(x, t) \rho^{\beta}(x, s, t) - \sigma \frac{\partial \rho^{\beta}}{\partial s}(x, s, t) \right) \right|_{s=0} = 0. $$

The next order in the expansion in terms of M and N in the mean field limit can be analysed by understanding the behaviour of a stochastic PDE for the fluctuations of the system; see [63]. See also [203] for a neural coding perspective on grid cell modelling.

## 2.3 Mathematical setting and definitions

Consider again the $d$-dimensional torus of length $L$, $T^d := \mathbb{R}^d/(L\mathbb{Z})^d$, $L > 0$. Writing proofs about the system (GC) being quite cumbersome, an efficient route is to first study a simpler model, without the orientation preference $\beta$, and then see what results can be extended to the general case. Thus, we consider the simpler system for one population given by

$$ \tau_c \frac{\partial \rho}{\partial t} = - \frac{\partial}{\partial s} \left[ \left( \Phi \left( \int_{T^d} W(x-y) \int_0^{+\infty} s \rho(y, s, t) ds dy + B \right) - s \right) \rho(x, s, t) \right] + \sigma \frac{\partial^2 \rho}{\partial s^2}, \quad (\text{SNF}) $$

with the nonlinear and non-local no-flux boundary condition

$$ \left. \left( \Phi \left( \int_{T^d} W(x-y) \int_0^{+\infty} s \rho(y, s, t) ds dy + B \right) \rho(x, s, t) - \sigma \frac{\partial \rho}{\partial s}(x, s, t) \right) \right|_{s=0} = 0, $$

and the mass normalisation condition

$$ \forall x \in T^d, \quad \int_0^{+\infty} \rho^0(x, s) ds = \frac{1}{L^d}. \qquad (2.4) $$

The system (SNF) can be thought of as a Stochastic Neural Field in the form of a partial differential equation. Other types of stochastic neural fields have been proposed and investigated with probabilistic techniques [149, 171, 16, 15, 154, 155, 61, 156]; see also books and reviews [13, 66, 68, 65] for discussions about the extensive literature on neural fields in general, and [23, 24, 67, 251, 240, 56] for efforts to build more realistic second-generation neural fields.

We will use the classical notation $u \ast v : x \mapsto \int_{T^d} u(x-y)v(y)dy$ for the convolution of $u$ and $v$ on the torus $T^d$. We also denote

$$ \Phi_{\bar{\rho}}(x,t) = \Phi(W * \bar{\rho}(x,t) + B), \quad \Phi'_{\bar{\rho}}(x,t) = \Phi'(W * \bar{\rho}(x,t) + B), \quad \Phi''_{\bar{\rho}}(x,t) = \Phi''(W * \bar{\rho}(x,t) + B), $$

with the mean activity level defined by

$$ \bar{\rho}(x, t) := \int_{0}^{+\infty} s \rho(x, s, t) ds. $$

The main equation (SNF) can then be rewritten as

$$ \tau_c \frac{\partial \rho}{\partial t} = - \frac{\partial}{\partial s} \left[ (\Phi_{\bar{\rho}} - s) \rho \right] + \sigma \frac{\partial^2 \rho}{\partial s^2}, \quad \sigma \frac{\partial \rho}{\partial s}(x, 0, t) = \Phi_{\bar{\rho}} \rho(x, 0, t). $$

The nonlinear non-local boundary condition is the principal theoretical difficulty: it destroys the gradient flow structure of the Fokker-Planck equation, complicates the search for stationary states and the construction of smooth solutions. However, note that weak solutions are already provided by the mean-field limit in [42].

For the spatial component $T^d$ of the system, we will use the Hilbert space $L^2(T^d)$ or the subspace $L_S^2(T^d)$ of coordinate-wise even $L^2$-functions:

$$ L_S^2(T^d) = \left\{ u \in L^2(T^d) \mid \forall i \in [1, d], \ u(x_1, \dots, -x_i, \dots, x_d) = u(x_1, \dots, x_i, \dots, x_d) \text{ a.e. in } T^d \right\}. $$

Similarly, for any $m \in \mathbb{N}$, the coordinate-wise even Sobolev space is $H_S^m(T^d) = H^m(T^d) \cap L_S^2(T^d)$.
**Assumption 2.1 (On $\Phi$, $W$ and $B$).** There exists $m \in \mathbb{N}$, such that $2m > d$ and $W \in H_S^m(\mathbb{T}^d)$. The function $\Phi$ is continuous on $\mathbb{R}$, $C^3$ on $\mathbb{R}^*$, increasing on $\mathbb{R}_+$ and $\Phi \equiv 0$ on $\mathbb{R}_-$. Moreover, $B > 0$ and $W_0 = \int_{\mathbb{T}^d} W(y)dy < 0$.

There are biological and mathematical reasons for these choices:

* Because $W \in H_S^m(\mathbb{T}^d)$, with $2m > d$, the smoothing property of the convolution and Sobolev embeddings imply that $W * h$ will be continuous for all $h \in L^2(\mathbb{T}^d)$. It also avoids pathological connectivity kernels with singularities.

* Consistent with the modelling purpose, the external signal $B$ is excitatory, the connectivity kernel $W$ is inhibitory on average and the modulation function $\Phi$ acts as a filter, transmitting only excitatory signals to the network.

For the evolution problem, we will use the following definitions.

**Definition 2.2.** (Fast decay at $s = +\infty$) Let $a \in \mathbb{R}$. We say that a function $\rho: \mathbb{T}^d \times [a, +\infty) \to \mathbb{R}$, is fast-decaying if for all polynomial functions $s \mapsto P(s)$, for all $x \in \mathbb{T}^d$, $\lim_{s \to +\infty} \rho(x, s)P(s) = 0$.

**Definition 2.3 (Compatible initial condition).** We say that the function $\rho^0 \in C^{0,1}(\mathbb{T}^d \times \mathbb{R}_+)$ is a compatible initial condition for (SNF) if $\rho^0$ is a non-negative probability density, if $\rho^0$, $\partial_s\rho^0$ are fast-decaying, and if

$$ \forall x \in \mathbb{T}^d, \quad \int_0^{+\infty} \rho^0(x, s)ds = \frac{1}{L^d}. $$

**Definition 2.4 (Classical solution).** We say that a non-negative function $\rho \in C^0(\mathbb{T}^d \times \mathbb{R}_+ \times [0, T_{\max}))$, $T_{\max} \in (0, +\infty]$ is a classical solution of (SNF) if

* it has the regularity $\rho \in C^{0,2,1}(\mathbb{T}^d \times (0, +\infty) \times (0, T_{\max})) \cap C^{0,1,1}(\mathbb{T}^d \times \mathbb{R}_+ \times (0, T_{\max}))$, and it satisfies the equation in the classical sense;

* for all $t \in [0, T_{\max})$, $\rho(\cdot, \cdot, t)$, $\partial_s\rho(\cdot, \cdot, t)$ are fast-decaying;

* for all $(x, t) \in \mathbb{T}^d \times [0, T_{\max})$, $\int_0^{+\infty} \rho(x, s, t)ds = \frac{1}{L^d}$.

## 2.4 Analytical study of the stationary problem

Numerical simulations in [50] show that the model (GC) reproduces the hexagonal firing pattern of grid cells for appropriate choice of parameters as seen in Figure 3. Although there is no rigorous proof of the existence of hexagonal patterns yet, partial answers were obtained rigorously via a bifurcation analysis in [53]. In the following, we present the theory for the system (SNF) for simplicity, but the technical extension for (GC) has been done in [50, 53].

### 2.4.1 Characterisation of stationary states

Integrating with respect to $s$ and using the boundary condition, the stationary states satisfy

$$ \sigma \partial_s \rho(x, s) = -(s - \Phi(W * \bar{\rho}(x) + B)) \rho(x, s), $$

which yields the form

$$ \rho(x, s) = \frac{1}{Z_{\rho}} \exp \left( - \frac{(s - \Phi(W * \bar{\rho}(x) + B))^2}{2\sigma} \right), \quad (2.5) $$

with, owing to (2.4), $Z_{\rho} = L^d \int_0^{+\infty} \exp \left( - \frac{(s - \Phi(W * \bar{\rho}(x) + B))^2}{2\sigma} \right) ds$. Stationary states must be zeros of

$$ \mathcal{G}(\rho, \sigma) = \rho - \frac{1}{Z_{\rho}} e^{-\frac{(s-\Phi\bar{\rho})^2}{2\sigma}}, \quad Z_{\rho} = L^d \int_0^{+\infty} e^{-\frac{(s-\Phi\bar{\rho})^2}{2\sigma}} ds. $$
Consequently, the average $x \mapsto \bar{\rho}(x)$ is a zero of the functional

$$ \bar{G}(\bar{\rho}, \sigma) = \bar{\rho} - \frac{1}{Z_{\rho}} \int_{0}^{+\infty} s e^{-\frac{(s-\Phi_{\bar{\rho}})^2}{2\sigma}} ds, \quad Z_{\rho} = L^d \int_{0}^{+\infty} e^{-\frac{(s-\Phi_{\bar{\rho}})^2}{2\sigma}} ds. \quad (2.6) $$

We can fully characterise the stationary states of (SNF) considering the zeros of the functional $\bar{G}(\bar{\rho}, \kappa)$: for any zero of the averaged functional, a stationary state can be defined through (2.5) and for any stationary state, its activity average is a solution to (2.6).

A magic ingredient for the study of this problem is the short calculation

$$ \begin{aligned} \rho(x, 0) &= - \int_{0}^{+\infty} \frac{\partial \rho}{\partial s}(x, s) ds = \frac{1}{\sigma} \int_{0}^{+\infty} (s - \Phi_{\bar{\rho}}(x)) \rho(x, s) ds \\ &= \frac{1}{\sigma} \left( \underbrace{\int_{0}^{+\infty} s \rho(x, s) ds}_{= \bar{\rho}(x)} - \Phi_{\bar{\rho}}(x) \underbrace{\int_{0}^{+\infty} \rho(x, s) ds}_{= \frac{1}{L^d}} \right), \end{aligned} $$

yielding a relation between the value at $s = 0$ of a stationary state, its firing rate and its mean value,

$$ \rho(x, 0) = \frac{1}{\sigma} \left( \bar{\rho}(x) - \frac{\Phi_{\bar{\rho}}(x)}{L^d} \right). \quad (2.7) $$

For space-homogeneous stationary states $\rho_{\infty}(s)$, the characterisation (2.6) simplifies to an easier scalar problem.

**Theorem 2.5** (Carrillo, Holden, Solem [50]). *Grant Assumption 2.1. Then, for all $\sigma > 0$, there exists a unique space-homogeneous stationary state $s \mapsto \rho_{\infty}^{\sigma}(s)$ of (SNF).*

The actual result is stated under slightly more general assumptions in [50]. Denoting

$$ \Phi_0 = \Phi(W * \bar{\rho}_{\infty} + B), \quad \Phi'_0 = \Phi'(W * \bar{\rho}_{\infty} + B), \quad \Phi''_0 = \Phi''(W * \bar{\rho}_{\infty} + B), $$

the equation (2.7) becomes the scalar identity

$$ \rho_{\infty}(0) = \frac{1}{\sigma} \left( \bar{\rho}_{\infty} - \frac{\Phi_{0}}{L^{d}} \right). $$

### 2.4.2 Noise-driven bifurcations

As we mentioned, noise is ubiquitous in the brain and plays an important role, which was the main reason for the development of this new PDE model. An interesting question is: How robust are the grid cell patterns to noise? If increasing the noise can flatten the activity and kill all stationary states but the space-homogeneous one, then while decreasing the noise bifurcations should occur.

The first study [50] on this PDE showed through numerical simulations that noise-driven bifurcations from the space-homogeneous stationary state do occur for the general four populations model (GC), and that the hexagonal grid cell pattern appears at some point via secondary bifurcations. They study a linearised version of the PDE and obtain a rigorous criterion for when the space-homogeneous stationary state can be guaranteed to be stable in the linear problem, which we will describe in more detail in Section 2.6.2. The bifurcation analysis was then made rigorous in [53] as follows, building upon methods in [49] for the McKean-Vlazov equation on the torus.

**Functional analysis setup** For simplicity in computation, the bifurcation parameter will be

$$ \kappa = \frac{1}{\sigma}. $$
<figure><img src="image_5.png" /><figcaption>Figure 3: Numerical simulation of the four population model for grid cells (GC) for different values of $\sigma$. Top Row is the repartition of the density for $s = 0.5$; middle row is the same for $s = 0$; bottom row represents in grey the path of a real rat and in red the prediction of the model of where in physical space one specific grid cell fires – the one at $x = (0, 0)$. Carrillo, Holden, Solem [50]</figcaption></figure>

Using the equivalence between the zeros of $\mathcal{G}$ and $\tilde{\mathcal{G}}$, we only need a functional setting for the space component. While working in the space $L_S^2(\mathbb{T}^d)$, we use the Hilbert basis $(\omega_k)_{k \in \mathbb{N}^d}$ defined by

$$ \omega_k(x) = \frac{\Theta(k)}{L^{\frac{d}{2}}} \prod_{i=1}^{d} \omega_{k_i}(x_i), \quad \omega_{k_i}(x) = \cos\left(\frac{2\pi k_i}{L} x_i\right), \quad \Theta(k) = \prod_{i=1}^{d} \sqrt{2 - \delta_{k_i, 0}}, \quad (2.8) $$

where $\delta_{i,j}$ is 0 if $i \neq j$ and 1 if $i = j$. The Fourier modes of $W$ are denoted by

$$ \tilde{W}(k) = \langle W | \omega_k \rangle = \int_{\mathbb{T}^d} W(x) \omega_k(x) dx, \quad k \in \mathbb{N}^d, $$

If $W \in L_S^2(\mathbb{T}^d)$, for all $g \in L_S^2(\mathbb{T}^d)$,

$$ W * g(x) = \int_{\mathbb{T}^d} W(x - y)g(y)dy = \frac{L^{\frac{d}{2}}}{\Theta(k)} \sum_{k \in \mathbb{N}^d} \tilde{W}(k)\tilde{g}(k)\omega_k(x). $$

In particular, for all $k \in \mathbb{N}^d$,

$$ W * \omega_k(x) = \frac{L^{\frac{d}{2}} \tilde{W}(k)}{\Theta(k)} \omega_k(x). $$

In the larger space $L^2(\mathbb{T}^d)$, we use the Hilbert basis $(\omega_k)_{k \in \mathbb{Z}^d}$ defined by (2.8),

$$ \omega_{k_i}(x) = \begin{cases} \cos\left(\frac{2\pi k_i}{L} x_i\right) & \text{if } k_i > 0, \\ 1 & \text{if } k_i = 0, \\ \sin\left(\frac{2\pi k_i}{L} x_i\right) & \text{if } k_i < 0. \end{cases} \quad \text{and} \quad \Theta(k) = \prod_{i=1}^{d} \sqrt{2 - \delta_{k_i, 0}}. $$
**Applying the Crandall-Rabinowitz theorem** The simplest way to prove the existence of noise-driven bifurcations is to apply the Crandall-Rabinowitz theorem to a modified version of the $\kappa$-rewritten functional

$$ \bar{G}(\bar{\rho}, \kappa) = \bar{\rho} - \frac{1}{Z_{\rho}} \int_{0}^{+\infty} se^{-\kappa \frac{(s-\Phi_{\bar{\rho}})^2}{2}} ds, \quad Z_{\rho} = L^d \int_{0}^{+\infty} e^{-\kappa \frac{(s-\Phi_{\bar{\rho}})^2}{2}} ds, $$

Denote $\rho_{\infty}^{\kappa}$ the space-homogeneous stationary state associated to the parameter $\kappa$, which has firing rate $\Phi_{0}^{\kappa}$, and define the modified functional

$$ \mathcal{H}: L_S^2(\mathbb{T}^d) \times (0, +\infty) \rightarrow L_S^2(\mathbb{T}^d) \\ (\bar{\rho}, \kappa) \mapsto \bar{G}(\bar{\rho}_{\infty}^{\kappa} + \bar{\rho}, \kappa). $$

Note that for all $\kappa > 0$, $\mathcal{H}(0, \kappa) = 0$. The Fréchet derivatives of $\mathcal{H}$ at $(0, \kappa)$ can be computed in terms of the ones of $\bar{G}$:

$$ D_{\bar{\rho}}\mathcal{H}(0, \kappa)[h] = D_{\bar{\rho}}\bar{G}(\bar{\rho}_{\infty}^{\kappa}, \kappa)[h], $$

$$ D_{\kappa}\mathcal{H}(0, \kappa) = D_{\kappa}\bar{G}(\bar{\rho}_{\infty}^{\kappa}, \kappa) + D_{\bar{\rho}}\bar{G}(\bar{\rho}_{\infty}^{\kappa}, \kappa) \left[ \frac{d\bar{\rho}_{\infty}}{d\kappa}(\kappa) \right], $$

$$ D_{\bar{\rho}\kappa}^2 \mathcal{H}(0, \kappa)[h] = D_{\bar{\rho}\kappa}^2 \bar{G}(\bar{\rho}_{\infty}^{\kappa}, \kappa)[h] + D_{\bar{\rho}\bar{\rho}}^2 \bar{G}(\bar{\rho}_{\infty}^{\kappa}, \kappa) \left[ \frac{d\bar{\rho}_{\infty}}{d\kappa}(\kappa), h \right], $$

for all $h \in L_S^2(\mathbb{T}^d)$. The road map for proving a bifurcation at parameter $\kappa^*$ consists in the following steps:

i) Check smoothness assumptions for $\mathcal{H}$ and that $D_x\mathcal{H}(0, \kappa^*)$ is a Fredholm operator of index 0.

ii) Check that $\dim(\ker(D_x\mathcal{H}(0, \kappa^*))) = 1$.

iii) Check a second order condition.

**Theorem 2.6** (Carrillo, Roux, Solem [53]). *Grant Assumption 2.1, and let $\kappa^* \in (\frac{2|W_0|^2}{L^{2d}\pi B^2}, +\infty)$. Assume that there exists a unique $k^* \in \mathbb{N}^d$ such that*

$$ \frac{\tilde{W}(k^*)}{\Theta(k^*)} = \frac{1}{L^{\frac{d}{2}}(\Phi_0^{k^*})' \left( \frac{1}{L^d} - L^d \bar{\rho}_{\infty}^{k^*} \left( \bar{\rho}_{\infty}^{k^*} - \frac{\Phi_0^{k^*}}{L^d} \right) \kappa^* \right)}. \quad (2.9) $$

Assume also that $(\Phi_0^{k^*})'' \ge 0$. Then $(\bar{\rho}_{\infty}^{k^*}, \kappa^*)$ is a bifurcation point of $\bar{G}(\bar{\rho}, \kappa) = 0$.

Moreover, in a neighbourhood $U \times (\kappa^* - \varepsilon, \kappa^* + \varepsilon)$ of $(\bar{\rho}_{\infty}^{k^*}, \kappa^*)$ in $L_S^2(\mathbb{T}^d) \times \mathbb{R}_+^*$, the average of the stationary state is either of the form $(\bar{\rho}_{\infty}^{k^*}, \kappa)$ or on the non-homogeneous solution curve

$$ \{ (\bar{\rho}_{\kappa(z)}, \kappa(z)) \mid z \in (-\delta, \delta), (\bar{\rho}_{\kappa(0)}, \kappa(0)) = (\bar{\rho}_{\infty}^{k^*}, \kappa^*), \delta > 0 \}, $$

with

$$ \bar{\rho}_{\kappa(z)}(x) = \bar{\rho}_{\infty}^{k(z)} + z\omega_{k^*}(x) + o(z), \quad x \in \mathbb{T}^d. $$

Notice that the new branch cannot be made up of space-homogeneous stationary states due to the presence of the vector $\omega_{k^*}$ in the asymptotic expression for small $z$. Moreover, the vector $\omega_{k^*}$ that spans the kernel of the functional derivative gives information about the shape of the new pattern close to the bifurcation point.
**Equivariant bifurcations** When there are symmetries in the connectivity kernel W, the uniqueness condition on the Fourier coefficient W(k*) may not be satisfied, but a bifurcation may still be found. Based on equivariant bifurcation theory [60, 85, 105], symmetries can be taken into account to give a more general classification.

**Theorem 2.7** (Carrillo, Roux, Solem [53]). Under the conditions of Assumption 2.1, consider $\mathcal{L}$ the trivial $d$-dimensional cubic lattice of the flat torus generated by the vectors $(L, 0, \dots, 0), \dots, (0, \dots, 0, L)$. Let $\Gamma$ be a subgroup of $E_d$. Assume $\mathcal{H}$ is $\Gamma$-equivariant on $L^2(\mathbb{T}^d)$ and that there exists a multi-index $k^* \in \mathbb{N}^d$ and $\kappa > 0$ such that (2.9) holds with $(\Phi_0^\kappa)'' \ge 0$. Denote $\mathcal{V} = \text{Ker}(D_{\bar{\rho}}\mathcal{H}(0, \kappa)) \subset L^2(\mathbb{T}^d)$. Then, for all isotropy subgroups $\Sigma$ of $\Gamma$, such that $\text{dim}(\text{Fix}_\mathcal{V}(\Sigma)) = 1$, there exists a unique branch of stationary states of (SNF) bifurcating from the trivial curve and whose average $\bar{\rho}$ in $s$ has the symmetry of $\Sigma$.

The application framework for grid cells is in dimension $d = 2$. The Euclidean group $E_2$ associated with the cubic lattice on the torus is the compact semidirect sum of the dihedral group $D_4$ and the compact group of translations on $\mathbb{T}^2$, $D_4 \dot{+} \mathbb{T}^2$. Note that the semidirect sum $\dot{+}$ is in fact a semidirect product, but noted differently for notational clarity in a geometry context. Assume W is radially symmetric. Then $\mathcal{H}$ is $\Gamma$-equivariant with $\Gamma = D_4 \dot{+} \mathbb{T}^2$. The patterns obtained from Theorem 2.7 can be classified using the procedure in [85, 84]:

* write the kernel $\mathcal{V} \subset L^2(\mathbb{T}^d)$ as a direct sum of $\Gamma$-irreducible subspaces $\mathcal{V} = \mathcal{V}_1 \oplus \dots \oplus \mathcal{V}_n$, which in turn implies, for any isotropy subgroup $\Sigma$ of $\Gamma$,

$$ \text{Fix}_{\mathcal{V}}(\Sigma) = \text{Fix}_{\mathcal{V}_1}(\Sigma) \oplus \dots \oplus \text{Fix}_{\mathcal{V}_n}(\Sigma). $$

* classify $\Gamma$-irreducible representations of $\mathcal{V}$ depending on its dimension.

In dimension 2, this work has already been done, for example in a general context [85, 84] or in the context of deterministic neural field models in [105, 266]. If W is radially symmetric and the terms $\tilde{W}(k)/\Theta(k)$ are unique up to permutations of $k$, then $\text{dim}(\mathcal{V}) = 4$ or $\text{dim}(\mathcal{V}) = 8$. The method allows to prove the existence of branches bifurcating from the spatially homogeneous stationary state with the following symmetries: stripes, also called rolls [50, Figure 4]; simple squares [84, Figure 1]; squares [84, Figure 2]; and anti-squares [84, Figure 3].

Numerical results indicate that the first bifurcation points are associated to the multi-indices $k = (4, 0), (0, 4), k = (4, 1), (1, 4)$, and $k = (3, 3)$. Although equivariant bifurcation theory does not allow us to prove the existence of stationary states with hexagonal symmetry, it is likely that the hexagonal firing pattern appears as a secondary bifurcation, since the indexes $(4, 1), (1, 4), (3, 3)$ are related to hexagonal patterns in many mathematical models exhibiting bifurcations [174]. It can also be noted that the simple addition of these modes,

$$ \omega(x_1, x_2) = \cos(6\pi x_1) \cos(6\pi x_2) + \cos(8\pi x_1) \cos(2\pi x_2) + \cos(2\pi x_1) \cos(8\pi x_2), $$

looks very close visually to the transient pattern in Figure 4 at $t = 1810ms$, and that

$$ \omega^{\text{hex}}(x_1, x_2) = \cos(6\pi x_1) \cos(6\pi x_2) + \cos(8\pi x_1) \cos(2\pi x_2) + \cos(2\pi x_1) \cos(8\pi x_2) \\ + \sin(6\pi x_1) \sin(6\pi x_2) - \sin(8\pi x_1) \sin(2\pi x_2) - \sin(2\pi x_1) \sin(8\pi x_2), $$

looks very close to the final hexagonal pattern in Figure 4 at $t = 2400ms$.

## 2.5 Global-in-time existence of classical solutions

In [54] global-in-time existence of smooth solutions is proved for both the one population system (SNF) and the general system (GC) under fairly general assumptions.
The proof consists in applying the same changes of variables than in the NNLIF existence proof. First, an affine change of variables to fix $\sigma = \tau_c = 1$. We use again the variables $y = e^t s$ and $\tau = \frac{1}{2}(e^{2t} - 1)$, and the function $\alpha(\tau) = (2\tau + 1)^{-\frac{1}{2}} = e^{-t}$, to define the new unknown

$$q(x, y, \tau) = \alpha(\tau)\rho(x, y\alpha(\tau), -\log(\alpha(\tau))) \quad \text{or equivalently} \quad \rho(x, s, t) = e^t q(x, e^t s, \frac{1}{2}(e^{2t} - 1)).$$

and we define $\beta(\tau) = B(-\log(\alpha(\tau)))$ and $\Psi(x, \tau) = \Phi(\alpha(\tau)W * \bar{q}(x, \tau) + \beta(\tau))\alpha(\tau)$. After the additional change of variable $z = y - \int_0^\tau \Psi(x, \eta)d\eta$, defining $u(x, z, \tau) = q(x, y, \tau)$ and

$$\gamma(x, \tau) = - \int_{0}^{\tau} \Psi(x, \eta) d\eta, \quad \bar{u} = \int_{\gamma(x, \tau)}^{+\infty} zu(z, \tau) dz,$$

we come to another nonlocal free-boundary Stefan-like problem,

$$\left\{
\begin{array}{ll}
\displaystyle \frac{\partial u}{\partial \tau}(x, z, \tau) = \frac{\partial^2 u}{\partial z^2}(x, z, \tau), & z \in (\gamma(x, \tau), +\infty), \tau \in \mathbb{R}_+, \\[2ex]
\displaystyle \gamma(x, \tau) = -\int_0^\tau \Psi(x, \eta) d\eta, & \tau \in \mathbb{R}_+, \\[2ex]
\displaystyle \frac{\partial u}{\partial y}(x, \gamma(x, \tau), \tau) = \Psi(x, \tau)u(x, \gamma(x, \tau), \tau), & \tau \in \mathbb{R}_+, \\[2ex]
\displaystyle \Psi(x, \tau) = \Phi(\alpha(\tau)W * [\bar{u}(x, \tau) - \gamma(x, \tau)] + \beta(\tau))\alpha(\tau), & \tau \in \mathbb{R}_+,
\end{array}
\right.$$

with $u(x, z, 0) = u^0(x, z)$ on $\mathbb{T}^d \times \mathbb{R}_+$ defined from $\rho^0$. Using again the heat kernel and Green identity

$$G(z, \tau, \xi, \eta) = \frac{1}{\sqrt{4\pi(\tau - \eta)}} e^{-\frac{(z-\xi)^2}{4(\tau-\eta)}}, \quad \frac{\partial}{\partial \xi} \left( G \frac{\partial u}{\partial \xi} - u \frac{\partial G}{\partial \xi} \right) - \frac{\partial}{\partial \eta} (Gu) = 0,$$

we can integrate for each $x \in \mathbb{T}^d$ the identity in the domain $(\gamma(x, t), +\infty) \times (0, \tau)$; we thereby obtain

$$u(x, z, \tau) = \int_{0}^{+\infty} G(z, \tau, \xi, 0)u^0(x, \xi)d\xi + \int_{0}^{\tau} \frac{\partial G}{\partial \xi}(z, \tau, \gamma(x, \eta), \eta)u(x, \gamma(x, \eta), \eta)d\eta. \quad (2.10)$$

Defining the boundary value term $v(x, \tau) = u(x, \gamma(x, \tau), \tau)$ and doing some work on (2.10), it is possible to obtain a closed integral system

$$\left\{
\begin{array}{ll}
v(x, \tau) = F_v[v, \gamma, \bar{u}](x, \tau) & \text{value on the free boundary,} \\
\gamma(x, \tau) = F_\gamma[v, \gamma, \bar{u}](x, \tau) & \text{position of the free boundary,} \\
\bar{u}(x, \tau) = F_{\bar{u}}[v, \gamma, \bar{u}](x, \tau) & \text{average in activity,}
\end{array}
\right.$$

where the activity $s$ has been eliminated. It allows to prove local-in-time existence of a solution to the original problem by application of the Banach fixed point theorem. A priori bounds and moment analysis lead to the conclusion that solutions are global-in-time for any Lipschitz continuous $\Phi$.

**Theorem 2.8** (Carrillo, Roux, Solem [54]). Assume $\Phi$ is locally Lipschitz, $B \in C^0(\mathbb{R}_+)$ and $W \in L^1(\mathbb{T}^d)$. For any compatible initial condition $\rho^0$, equation (SNF) admits a unique maximal solution in the sense of Definition 2.4. Moreover, if $\Phi$ is globally Lipschitz, the solution is global-in-time and

* if $\Phi$ non-negative, then

$$\forall t \in \mathbb{R}_+, \forall x \in \mathbb{T}^d, \quad \rho(x, 0, t) \le \frac{1}{L^d} \sqrt{\frac{2}{\pi\sigma}} \frac{1}{\sqrt{1 - e^{-2\tau_c t}}};$$

* if $\Phi$ is non-positive, then

$$\forall t \in \mathbb{R}_+, \forall x \in \mathbb{T}^d, \quad \bar{\rho}(x, t) \le \frac{1}{L^d} (1 + \sigma - \sigma e^{-2\tau_c t}) + e^{-2\tau_c t} \sup_{x \in \mathbb{T}^d} \int_0^{+\infty} s^2 \rho^0(x, s) ds.$$

The result is then extended to the four population Grid cells model (GC) with similar assumptions and definitions.
## 2.6 Asymptotic stability of stationary states

The relative entropy method [172] can be used to study the nonlinear stability of the stationary states of (GC) and (SNF). As above, the idea is to study it first for the simpler (SNF) and then to generalise to (GC), with the exception of linear stability, which was performed directly on (GC) in [50].

An important quantity for the study of asymptotic stability for these models is the variance associated to any stationary state $\rho_\infty$:

$$M_\infty(x) = \int_0^{+\infty} (s - L^d \bar{\rho}_\infty(x))^2 \rho_\infty(x, s) ds.$$

As proved in [53], we can obtain the following identity

$$\frac{M_\infty}{\sigma} = \frac{1}{L^d} g\left(\frac{\Phi_0}{\sqrt{2\sigma}}\right), \quad g(\eta) = 1 - \frac{2}{\sqrt{\pi}} \frac{\exp(-\eta^2)}{1 + \operatorname{erf}(\eta)} \left[ \frac{1}{\sqrt{\pi}} \frac{\exp(-\eta^2)}{1 + \operatorname{erf}(\eta)} + \eta \right]$$

where $g$ is a differentiable positive bounded function whose image lies in $(1 - \frac{2}{\pi}, 1)$.

### 2.6.1 Naive entropy method

The easiest way to get nonlinear stability is to do, like for the NNLIF model, a weakly-nonlinear analysis via a Lyapunov control on the dissipation of the relative entropy. The key step is to control the non-local and nonlinear effect of the firing rate.

**Theorem 2.9** (Carrillo, Roux, Solem [54]). Let $\Phi$ be Lipschitz, $W \in L^2(\mathbb{T}^d)$, and $\rho_\infty$ be a stationary state of (SNF). Assume that $\sigma$ satisfies

$$\|\Phi'\|_{\infty} \frac{\|W\|_{L^2(\mathbb{T}^d)}}{L^{\frac{d}{2}}} \left\| M_{\infty}^{\frac{1}{2}} \right\|_{\infty} < \frac{\sigma}{2} \tilde{\gamma}(\rho_{\infty})^{\frac{1}{2}}, \quad (2.11)$$

where $\tilde{\gamma}(\rho_\infty) = \inf_{x \in \mathbb{T}^d} \gamma(\rho_\infty(x))$, and $\gamma(\rho_\infty(x))$ is a Poincaré constant associated to $\rho_\infty(x)$. Then, the $L^1$ norm of the relative entropy

$$\int_{\mathbb{T}^d} \int_0^{+\infty} \left( \frac{\rho(x, s, t) - \rho_\infty(x, s)}{\rho_\infty(x, s)} \right)^2 \rho_\infty(x, s) ds dx, \quad (2.12)$$

decays exponentially fast whenever the compatible initial datum $\rho_0$ is close enough to the stationary state $\rho_\infty$ in relative entropy.

The idea of the proof is to write $h = \frac{\rho}{\rho_\infty}$ and choose $G(h) = \frac{1}{2}(h-1)^2$. Then the entropy dissipation can be bounded by

$$\frac{d}{dt} \int_{0}^{+\infty} \frac{1}{2}(h-1)^2 \rho_{\infty} ds \le \frac{1}{2\sigma} |\Phi_0 - \Phi_{\bar{\rho}}|^2 \int_{0}^{+\infty} (h-1)^2 \rho_{\infty} ds \\ - \frac{\sigma}{2} \int_{0}^{+\infty} \left(\frac{\partial h}{\partial s}\right)^2 \rho_{\infty} ds - (\Phi_0 - \Phi_{\bar{\rho}}) \int_{0}^{+\infty} \frac{\partial h}{\partial s} \rho_{\infty} ds.$$

The deviation of the firing rate from the stationary firing rate can be controlled by

$$|\Phi_0 - \Phi_{\bar{\rho}}| \le \|\Phi'\|_{\infty} \|W\|_{L^2(\mathbb{T}^d)} \sup_{x \in \mathbb{T}^d} M_{\infty}^{\frac{1}{2}}(x) \left( \int_{\mathbb{T}^d} \int_0^{+\infty} (h-1)^2 \rho_{\infty} ds dy \right)^{\frac{1}{2}}.$$

By applying a Poincaré inequality, we arrive at
$$ \frac{d}{dt} \int_{\mathbb{T}^d} \int_0^{+\infty} \frac{1}{2}(h-1)^2 \rho_\infty \mathrm{d}s \mathrm{d}y \le \frac{1}{2\sigma} C^2 \left( \int_{\mathbb{T}^d} \int_0^{+\infty} (h-1)^2 \rho_\infty \mathrm{d}s \mathrm{d}y \right)^2 + \left( C - \frac{\sigma}{2} \tilde{\gamma}(\rho_\infty)^{1/2} \right) \left( \int_{\mathbb{T}^d} \int_0^{+\infty} |\partial_s h|^2 \rho_\infty \mathrm{d}s \mathrm{d}x \right)^{1/2} \left( \int_{\mathbb{T}^d} \int_0^{+\infty} (h-1)^2 \rho_\infty \mathrm{d}s \mathrm{d}y \right)^{1/2}. $$

By assumption (2.11), $\tilde{C} := C - \frac{\sigma}{2} \tilde{\gamma}(\rho_\infty)^{1/2} < 0$. Thus, applying again the same Poincaré inequality,

$$ \frac{d}{dt} \int_{\mathbb{T}^d} \int_0^{+\infty} \frac{1}{2}(h-1)^2 \rho_\infty \mathrm{d}s \mathrm{d}y \le \left[ \frac{C^2}{2\sigma} \int_{\mathbb{T}^d} \int_0^{+\infty} (h-1)^2 \rho_\infty \mathrm{d}s \mathrm{d}y + \tilde{\gamma}(\rho_\infty)^{1/2} \tilde{C} \right] \int_{\mathbb{T}^d} \int_0^{+\infty} (h-1)^2 \rho_\infty \mathrm{d}s \mathrm{d}y. $$

The expression between the brackets will be negative as long as the initial relative entropy is small enough, and then there is exponential convergence in relative entropy.

As the condition (2.11) depends nonlinearly on $\sigma$, it is not easy to check theoretically. Results on Poincaré constants [190] and bounds on $M_\infty$ provide the clarified (less sharp) condition

$$ \frac{1}{L^d} \|\Phi'\|_{\infty} \|W\|_{L^2(\mathbb{T}^d)} < \frac{1}{2}. $$

The problem with this result is that, as it can be checked numerically, it applies only in cases where there is a unique stationary state (the space-homogeneous one). As only a norm of $W$ is involved, nothing can be concluded about nonlinear stability of the constant stationary state in terms of the shape of the connectivity kernel.

<figure><img src="image_6.png" /><figcaption>Figure 4: Numerical simulation of (SNF) starting from a small perturbation of a constant function, $t = 40, 220, 1500, 1810, 2190$, and $2400$ ms. Carrillo, Roux, Solem [53].</figcaption></figure>

## 2.6.2 Linear stability of the space-homogeneous stationary state

Stability of the space-homogeneous stationary state was first investigated in [50] by linearising the equation (GC). Let $\rho_\infty$ be a space-homogeneous stationary state. Denoting $\rho^\beta = \rho_\infty + u^\beta$, a linearised equation for the perturbation $(u^1, u^2, u^3, u^4)$ is,

$$ \tau \frac{\partial u^\beta}{\partial t} + \frac{\partial \rho_\infty}{\partial s} \frac{\Phi_0'}{4} \sum_{\beta'=1}^{4} W^{\beta'} * \bar{u}^{\beta'} + \frac{\partial}{\partial s} \left[ (\Phi_0 - s) u^\beta \right] - \frac{\partial^2 u}{\partial s^2} = 0. $$
Consider here an exponential basis for Fourier modes (not to be confused with the cosine and sine Fourier modes $\tilde{W}(k)$ above), for multi-indexes $k \in \mathbb{Z}^d$,

$$ \hat{W}(k) = \int_{\mathbb{T}^d} W(x)e^{-i2\pi k \cdot x} dx. $$

Assume that they satisfy $F(k) < \frac{\sigma}{M_{\infty}}$, with

$$ F(k) = \frac{\Phi'_{0}}{4} \hat{W}(k) \sum_{\beta=1}^{4} e^{-2i\pi k \cdot r^{\beta}}. $$

Notice that these quantities are real due to the symmetry by component of the connectivity $W$ and the shifts $r^{\beta}$. If the linear perturbation $u$ is taken in $L^2(\mathbb{T}^d \times \mathbb{R}_+)$ and decomposed in plane waves as

$$ u^{\beta}(x, s, t) = \rho_{\infty}(s) \sum_{k \in \mathbb{Z}^d} u_k^{\beta}(s, t) e^{-2i\pi k \cdot x}, $$

then, we can reduce the $L^2$-relative entropy (2.12) to estimate

$$ \int_{0}^{+\infty} \rho_{\infty}|U|^{2}ds \quad \text{with} \quad U = \sum_{\beta} \exp(-i\mathbf{k} \cdot \mathbf{r}^{\beta})u^{\beta} \quad (2.13) $$

for all $k \in \mathbb{Z}^d$. However, computing the evolution of this linearized $L^2$-entropy for one Fourier mode (2.13) over the linearized dynamics leads to a dissipation of the quantity

$$ \int_{0}^{+\infty} \rho_{\infty}|U|^{2}ds - \frac{F(\mathbf{k})}{\sigma} \left| \int_{0}^{\infty} \rho_{\infty}Uds \right|^{2} $$

This quantity can be rewritten as

$$ \int_{0}^{+\infty} \rho_{\infty}|V|^{2}ds + \left(1 - M_{\infty}\frac{F(\mathbf{k})}{\sigma}\right) M_{\infty}|c(t)|^{2} \quad (2.14) $$

by expanding $U$ into $U = V + c(t)(s - \bar{\rho}_{\infty})$. Here $\int_{0}^{+\infty} V \rho_{\infty} ds = 0$ and $c(t)$ is chosen such that $\int_{0}^{+\infty} V(s - \bar{\rho}_{\infty}) \rho_{\infty} ds = 0$ leading to

$$ c(t) = \frac{1}{M_{\infty}} \int_{0}^{\infty} \rho_{\infty} U ds, $$

Detailed computations in [50, Theorem 3.6] allowed to estimate both terms in the quantity (2.14) under the assumption (2.15) leading to the following result.

**Theorem 2.10** (Carrillo, Holden, Solem [50]). *Grant Assumption 2.1 and assume $d = 2$. The unique spatially-homogeneous stationary state $\rho_{\infty}$ of (GC) is linearly asymptotically stable in $L^2(\mathbb{T}^d \times \mathbb{R}_+)$ as long as for all multi-index $k \in \mathbb{Z}^d$,

$$ F(k) < \frac{\sigma}{M_{\infty}}. \quad (2.15) $$

By taking all shifts equal to zero, $B^{\beta} = B$ and all components $\rho_{\infty}^{\beta}$ equal, this implies a linear stability result for (SNF). It has then been proved in [53] that the loss of linear stability perfectly matches the first point of bifurcation on the space-homogeneous stationary states curve.
### 2.6.3 Nonlinear stability of the space-homogeneous stationary state

In order to go from linear stability to nonlinear stability, a possible method is to find nonlinear quantities which would, after linearisation, yield quantities similar to the linearly decaying (2.13). With this in mind, define

$$E(x, t) = \frac{1}{2} \int_{0}^{+\infty} \left( \frac{\rho - \rho_{\infty}}{\rho_{\infty}} \right)^2 \rho_{\infty} ds, \quad \text{and} \quad H(x, t) = \Phi^{\delta} \bar{\rho}^{\delta},$$

where $\Phi^{\delta} = \Phi_{\bar{\rho}} - \Phi_0, \bar{\rho}^{\delta} = \bar{\rho} - \bar{\rho}_{\infty}$. Linearising $\Phi^{\delta}$, we get

$$2E - \frac{H}{\sigma} \approx \int_{0}^{+\infty} \left( \frac{\rho - \rho_{\infty}}{\rho_{\infty}} \right)^2 \rho_{\infty} ds - \frac{\Phi_{0}'}{\sigma} W * (\bar{\rho} - \bar{\rho}_{\infty})(\bar{\rho} - \bar{\rho}_{\infty}).$$

Define then

$$Q = \int_{\mathbb{T}^d} 2E - \frac{H}{\sigma} dx, \quad \mathcal{E}(t) = \int_{\mathbb{T}^d} \int_0^{+\infty} \left( \frac{\rho(x, s, t) - \rho_{\infty}(x, s)}{\rho_{\infty}(x, s)} \right)^2 \rho_{\infty}(x, s) ds dx.$$

It is possible to show under the condition (2.17) that there are constants $\alpha, C_1, C_2, C_3 > 0$ such that

$$\frac{d}{dt} Q \le -\sigma C_1 \alpha \left( \alpha - C_2 Q^{\frac{1}{2}} \right) Q$$

and,

$$2\alpha\mathcal{E}(t) \le Q(t) \le C_3\mathcal{E}(t). \qquad (2.16)$$

The differential inequality allows to prove nonlinear decay of $Q$ to 0 at exponential speed, provided $Q(0)$ is small enough. The second inequality proves that the relative entropy $\mathcal{E}$ and the quantity $Q$ are equivalent. In the proof of (2.16), the condition (2.17) appears naturally as a requirement.

**Theorem 2.11** (Carrillo, Roux, Solem [54]). Let $\Phi \in C^2(\mathbb{R})$, $\nabla\Phi \in W^{1,\infty}(\mathbb{R})$, $W \in L_S^2(\mathbb{T}^d)$, and $\rho_{\infty}$ be a spatially homogeneous stationary state. Assume for a suitably small $\alpha \in (0, 1)$, that $\sigma > 0$ satisfies

$$\int_{\mathbb{T}^d} (1 - \alpha) g^2(x) - \frac{\Phi_g^\delta(x)}{\sigma} g(x) \, dx > 0, \quad g \in L^2(\mathbb{T}^d), \qquad (2.17)$$

where

$$\Phi_g^\delta(x) = \Phi(M_\infty W * g(x) + W_0 \bar{\rho}_\infty + B) - \Phi_0.$$

Then, if $\mathcal{E}(0)$ is small enough, $\mathcal{E}(t)$ converges to 0 exponentially fast.

In order to better understand the link between condition (2.17) and the linear stability condition, take for example the simplest case $\Phi(z) = z$ and send formally $\alpha$ to 0; then

$$\Phi_g^\delta(x) = M_\infty W * g(x) + W_0 \bar{\rho}_\infty + B - \Phi_0 = M_\infty W * g(x).$$

and condition (2.17) becomes

$$\int_{\mathbb{T}^d} g^2(x) > \frac{M_\infty}{\sigma} \int_{\mathbb{T}^d} \int_{\mathbb{T}^d} W(x-y)g(y)g(x) \, dy dx, \quad g \in L^2(\mathbb{T}^d),$$

If all Fourier modes of $W$ are negative (implying unconditional global linear stability), then for all $g \in L^2(\mathbb{T}^d)$,

$$0 > \int_{\mathbb{T}^d} \int_{\mathbb{T}^d} W(x - y)g(y)g(x) \, dy dx$$

and there is also local nonlinear stability. Condition (2.17) thus appears as a perturbation of the linear stability condition.

In [54], the asymptotic stability results for (SNF) are extended to the model (GC) under additional technical assumptions.
### 3 Rate models for decision making tasks

Understanding the mechanisms underlying decision-making has been one of the fundamental questions in behavioral neuroscience for the last decades. Optimizing tasking behaviour involves the constant making of decisions between alternative choices. Neurophysiological and theoretical neurosciences have started to reveal the neural mechanisms underlying decision-making. Identifying areas of the cortex of animals involved in particular decisions has been one of the cornerstones of neuroscience. The most basic studies focused on binary decision-making tasks, i.e., the subject is asked to make a choice between two alternatives according to an experimentally defined criterion based on sensory input information. The difficulty of making a decision can be manipulated by the level of uncertainty in the discrimination of the sensory input information influenced by regulating the signal to noise ratio.

In a wide range of such decision-making experiments, the behavioral response, i.e. performance and reaction times, can be properly described by a simple stochastic diffusion model or Fokker-Planck equation [225]. It is plausible to think that the underlying neuronal system perform decision-making by accumulation of evidences. Indeed, electrophysiological recordings from awake behaving monkeys performing decision-making have shown that trial-averaged spiking activity of some neurons shows a gradual, nearly linear increase shortly after stimulus onset until just before the response is made [230, 231, 242, 228, 141]. One example of these binary decision making tasks is the so-called random-dot task. The random-dot task is a two-choice visual decision-making task. In this task, the subject is confronted with a visual stimulus consisting of a field of randomly moving dots. On any one trial, a fixed fraction of the dots, determined by the experimentalist, moves coherently in one of two directions. The subject must discriminate in which direction the majority of dots are moving. Electrophysiological recordings from awake behaving monkeys performing this task have shown that trial-averaged spiking activity of neurons in the lateral intraparietal (LIP) cortex reflects the accumulation of information mentioned above [242, 228, 141]. The difference in the firing rates of cells with opposing preferred directions is a likely measure of the available visual evidence for the direction of coherent motion of the moving dot stimulus, further indicating that the area LIP is involved in this decision-making task.

A computational model of the cortical circuit underlying LIP cell activity has been proposed and studied numerically [268, 272, 167]. This model was able to explain qualitatively the experimentally observed trial-averaged spiking activity. The model consists of two populations of spiking neurons, within which interactions are mediated by excitatory synapses and between which interactions occur principally through an intermediary population of inhibitory interneurons. Sensory input may bias competition in favour of one of the populations, potentially resulting in a gradually developing decision in which neurons in the chosen population exhibit increased activity while activity in the other population is inhibited. Here, activity is measured as the global firing rate of each of the populations of spiking neurons. When the activity of one of the two populations exceeds a pre-defined threshold, then a behavioural response is generated.

There are two possible mechanisms by which the decision can be made which correspond to two different dynamical working points, see Figure 5. In one mechanism, the appearance of the stimulus destabilises the so-called spontaneous state, in which all neurons show low activity. In the other mechanism, even after stimulus presentation, the network can still sustain low activity, but random fluctuations eventually drive the collective activity of the network to an activated state, in which a fraction of the neurons fires at a high rate. Let us remark that the behavioural performance of the model (i.e. the fraction of trials in which the 'correct' direction is chosen) and corresponding reaction time to make a decision match well the data from both monkeys and humans performing the moving-dot task [235]. Therefore, it may be possible to establish a direct link between the neuronal and behavioural correlates of decision making.

In order to perform the explicit link between a neuronal circuit processing decision-making and a diffusion model able to describe the behavioural level, we consider a mean-field approximation consisting of two competing rate models. The precise model uses a Wilson-Cowan [270] type system, describing
<figure><img src="image_7.png" /><figcaption>Figure 5: Dynamics in the unbiased ($\lambda_i = 0$) case of (3.1) in the deterministic system $\beta = 0$ (straight lines) and the stochastic system $\beta > 0$ (wiggled line). Straight lines emphasize the location of the equilibrium points. The black bars are the equilibrium points of the deterministic differential system. The single particle realization of the stochastic system starts its dynamics at point (5, 5), moves almost straight forward towards a slow-manifold nearby the spontaneous state, and then slowly oscillates towards one of the decision states. Carrillo, Cordier, Mancini [44].</figcaption></figure>

the time evolution of two population of neurons firing rates $\nu_i$, $i = 1, 2$:

$$ \tau \frac{d\nu_i(t)}{dt} = -\nu_i(t) + \phi \left( \lambda_i + \sum_{j=1,2} w_{ij} \nu_j(t) \right) + \beta \mathcal{B}_i(t), \quad i = 1, 2, \qquad (3.1) $$

where $\tau$ is the typical time relaxation and $\mathcal{B}_i(t)$, $i = 1, 2$, represent normalized independent Brownian motions. In (3.1) the function $\phi(x)$ has a sigmoidal shape determining the response function of the neuron population to a mean excitation $x$ given by $x_i(t) = \lambda_i + \sum_j w_{ij}\nu_j$, $i = 1, 2$, in each population:

$$ \phi(x) = \frac{\nu_c}{1 + \exp(-\alpha(x/\nu_c - 1))}, $$

where $\lambda_i$ are the external stimuli applied to each neuron population and $w_{ij}$ are the connection coefficients. The parameter $\nu_c$ represents both the maximal activity rate of the population and the frequency input needed to drive the population to half of its maximal activity.

Following [226, 158, 79], we assume that neurons within a specific population are likely to correlate their activity, and to interact via strong recurrent excitation with a dimensionless weight $w_+ > 1$ greater than a reference baseline value established to 1. Analogously, neurons in two different populations are likely to have anti-correlated activity expressed by an excitatory weight lower than the baseline, $w_- < 1$. Furthermore, we assume that there is global feedback inhibition, as a result of which all neurons are mutually coupled to all other neurons in an inhibitory fashion (we will denote this inhibitory weight by $w_I$). As a result, the synaptic connection coefficients $w_{ij}$, representing the interaction between population $i$ and $j$, are the elements of a $2 \times 2$ symmetric matrix $W$ given by

$$ W = \begin{bmatrix} w_+ - w_I & w_- - w_I \\ w_- - w_I & w_+ - w_I \end{bmatrix}, $$
The typical synaptic values considered are such that $w_{-} < w_{I} < w_{+}$ leading to cross-inhibition and self-excitation.

Standard stochastic calculus methods show that the probability density $p = p(t, \nu)$ of finding the neurons of both populations firing at averaged rates $\nu = (\nu_{1}, \nu_{2})$ at $t > 0$, satisfies a Fokker-Planck equation of the form

$$ \partial_t p + \nabla \cdot ([-\nu + \Phi(\Lambda + W \cdot \nu)] p) - \frac{\beta^2}{2} \Delta p = 0 \quad (3.2) $$

where $\nu \in \Omega = [0, \nu_{m}] \times [0, \nu_{m}]$, $\Lambda = (\lambda_{1}, \lambda_{2})$, $\Phi(x_{1}, x_{2}) = (\phi(x_{1}), \phi(x_{2}))$, $\nabla = (\partial_{\nu_{1}}, \partial_{\nu_{2}})$ and $\Delta = \Delta_{\nu}$. We choose to complete equation (3.2) by the following no-flux boundary conditions:

$$ \left( [-\nu + \Phi(\Lambda + W \cdot \nu)] p - \frac{\beta^2}{2} \nabla p \right) \cdot n = 0 $$

where $n$ is the outward normal to the domain $\Omega$. Physically, these boundary conditions mean that neurons cannot spike with arbitrarily large firing rates and thus there is a typical maximal value of the averaged firing rate $\nu_{m}$ and that the solution to (3.2) is a probability density function, i.e.,

$$ \int_{\Omega} p(t, \nu) d\nu = 1. $$

In order to simplify notations, we will from now on consider the vector field $F = (F_{1}, F_{2})$, representing the flux in the Fokker-Planck equation:

$$ F = -\nu + \Phi(\Lambda + W \cdot \nu) = \begin{pmatrix} -\nu_{1} + \phi(\lambda_{1} + w_{11}\nu_{1} + w_{12}\nu_{2}) \\ -\nu_{2} + \phi(\lambda_{2} + w_{21}\nu_{1} + w_{22}\nu_{2}) \end{pmatrix} $$

then, equation (3.2) and boundary conditions (3.3) read:

$$ \partial_t p + \nabla \cdot \left( F p - \frac{\beta^2}{2} \nabla p \right) = 0 \quad (3.4) $$

$$ \left( F p - \frac{\beta^2}{2} \nabla p \right) \cdot n = 0 \quad (3.5) $$

The nonlinear nature of the SDE system (3.1), however, hinders analytical progress. For this reason, the main analysis of such noise driven probabilistic decision-making systems remains based on numerical investigations, which are time consuming because of the need for sufficiently many trials to generate statistically meaningful data. This is an advantage of the mean-field Fokker-Planck equation due to the existence of efficient numerical solvers. Moreover, the Fokker-Planck framework (3.4)-(3.5) is amenable to nonlinear analysis.

## 3.1 Analysis of the rate model

The analysis of (3.4)-(3.5) was performed in [44] where we refer for further details. Let us elaborate on the most important structural property of this model. Assume that $u_{1}$ and $u_{2}$ are smooth solutions of (3.4)-(3.5) with $u_{1} > 0$ in $\bar{\Omega}$. Given any convex smooth function $H = H(\omega)$, where $\omega = u_{2}/u_{1}$, we have

$$ \frac{d}{dt} [u_1 H(\omega)] = \frac{\beta^2}{2} \left( \Delta (u_1 H(\omega)) - u_1 H''(\omega) |\nabla \omega|^2 \right) - \nabla \cdot [F u_1 H(\omega)] . \quad (3.6) $$

In order to check this, let us develop the left hand side of (3.6), using (3.4) to obtain

$$ \begin{aligned} \frac{d}{dt} [u_1 H(\omega)] = & - \nabla \cdot [F u_1 H(\omega)] + u_1 F \cdot \nabla H(\omega) \\ & + \frac{u_2}{u_1} H'(\omega) \nabla \cdot (u_1 F) - H'(\omega) \nabla \cdot (u_2 F) \end{aligned} $$
$$+ \frac{\beta^2}{2} \left[ H(\omega) - \frac{u_2}{u_1} H'(\omega) \right] \Delta u_1 + \frac{\beta^2}{2} H'(\omega) \Delta u_2.$$

We separate now the computation in two parts: the one concerning first order derivatives (I) and the one concerning second order derivatives (II). We start by computing (I). Leaving the first term unchanged and developing the following ones, we get:

$$(I) = - \nabla \cdot [F u_1 H(\omega)] + H'(\omega) F \cdot \left( \nabla u_2 - \frac{u_2}{u_1} \nabla u_1 \right) \\ + \frac{u_2}{u_1} H'(\omega) \nabla u_1 \cdot F + H'(\omega) u_2 \nabla \cdot F \\ - H'(\omega) \nabla u_2 \cdot F - H'(\omega) u_2 \nabla \cdot F = - \nabla \cdot [F u_1 H(\omega)].$$

Concerning the second part (II), we can compute $\Delta (u_1 H(\omega))$ to obtain the relation

$$\Delta (u_1 H (\omega)) - u_1 |\nabla \omega|^2 H'' (\omega) = \left[ H (\omega) - \frac{u_2}{u_1} H' (\omega) \right] \Delta u_1 + H' (\omega) \Delta u_2$$

Multiplying by $\frac{\beta^2}{2}$ and collecting (I) and (II) completes the proof of (3.6). Let us now define the relative entropy between two solutions as

$$\mathcal{H}(u_2|u_1) = \int_{\Omega} u_1 H(\omega) \, d\nu.$$

For any $u_1$ and $u_2$ strong solutions of (3.4)-(3.5), and any $H$ smooth convex function then, the relative entropy is decreasing in time and

$$\frac{d}{dt} \mathcal{H}(u_2|u_1) = -\frac{\beta^2}{2} \mathcal{D}(u_2|u_1) = -\frac{\beta^2}{2} \int_{\Omega} u_1 H''(\omega) \left| \nabla \left( \frac{u_2}{u_1} \right) \right|^2 d\nu \le 0. \quad (3.7)$$

This is simply obtained by integrating (3.6) over the domain $\Omega$ and showing that the boundary terms disappear due to (3.5).

Under an additional assumption on the flux $F$ one can prove that the stationary problem has a unique classical solution. More precisely, assuming that

$$F \in C^1(\bar{\Omega}, \mathbb{R}^2) \text{ with } F \cdot n < 0 \text{ on } \partial\Omega, \quad (3.8)$$

there exists an unique probability density function $u_\infty \in H^4(\Omega)$, $u_\infty(\nu) > 0$ in $\bar{\Omega}$ satisfying

$$\left\{ \begin{array}{ll} \nabla \cdot \left( F p - \frac{\beta^2}{2} \nabla p \right) = 0 & \text{in } \Omega, \\ \\ \left( F u - \frac{\beta^2}{2} \nabla u \right) \cdot n = 0 & \text{on } \partial \Omega. \end{array} \right.$$

The consequences of the existence of this family of Liapunov functionals given in (3.7) for (3.4)-(3.5) have already been explored for several equations in [178, 179] and it is usually referred as the general relative entropy (GRE) inequalities. The same conclusions apply here as in those cases apply here and we just list them here refering to [44] for more details.

**Corollary 3.1** (Carrillo, Cordier, Mancini [44]). Given $F$ satisfying (3.8) and any solution $u$ with normalized initial data $u_0$ to (3.4)-(3.5), then the following properties hold:

i) Nonnegativity: If $u_0$ is nonnegative, then the solution $u$ of problem (3.4)-(3.5) is nonnegative.

ii) Contraction principle:

$$\int_{\Omega} |u(t, \nu)| \, d\nu \le \int_{\Omega} |u_0(\nu)| \, d\nu.$$
iii) $L^p$ bounds, $1 < p < \infty$:

$$ \int_{\Omega} u_{\infty}(\nu) \left| \frac{u(t, \nu)}{u_{\infty}(\nu)} \right|^p d\nu \le \int_{\Omega} u_{\infty}(\nu) \left| \frac{u_0(\nu)}{u_{\infty}(\nu)} \right|^p d\nu. $$

iv) Pointwise estimates:

$$ \inf_{\nu \in \Omega} \frac{u_0(\nu)}{u_{\infty}(\nu)} \le \frac{u(t, \nu)}{u_{\infty}(\nu)} \le \sup_{\nu \in \Omega} \frac{u_0(\nu)}{u_{\infty}(\nu)}. $$

This corollary is a consequence of the GRE inequality (3.7) with $H(\omega) = \omega^-$ (the negative part of $\omega$), $H(s) = |s|$, $H(s) = |s|^p$, and $H(s) = (s-k)_+^2$ respectively by approximation from smooth convex functions. Moreover, the GRE inequality gives the convergence of the solution $u(t)$ to the stationary state $u_\infty$.

**Corollary 3.2 (Long time asymptotic).** Given $F$ satisfying (3.8) and any solution $u$ with normalized initial data $u_0$ to (3.4)-(3.5), then

$$ \lim_{t \to \infty} \int_{\Omega} |u(t, \nu) - u_{\infty}(\nu)|^2 d\nu = 0. $$

Convergence rates have been identified in terms of Poincaré type inequalities in [51].

## 3.2 Reduction to a diffusion model

One dimensional Fokker-Planck effective reductions near bifurcation points were first obtained in [235]. They proposed to locally approximate the evolution on the slow manifold at the bifurcation point by Taylor expanding the nonlinearities with respect to the slow variable based on the different time scales. In this way, they obtained an effective potential of degree 4 for the one dimensional Fokker-Planck dynamics on the slow manifold.

In [45, 43] a different strategy was used instead by not performing any Taylor expansion and rather approximate the slow manifold in terms of the nonlinearities defining the evolution using the slow-fast character of the dynamical system, see [6]. Once, the slow manifold is approximated then one should project the Fokker-Planck dynamics on it to obtain a one dimensional Fokker-Planck reduction possibly valid beyond the local character of the expansion in [235]. In this way, an approximated full reduced potential was obtained together with an explicit formula for the stationary state distribution on the slow manifold. This effective reduced potential covers the other two stable equilibrium points, and not only the central equilibrium point changing its stability character as in [235].

We illustrate this strategy in the particular example studied in [43] of a neuronal population model with two pools with self-excitation and cross-inhibition as in [79] related to (3.1). The firing rates $\nu_1$ and $\nu_2$ of the neuronal networks are determined by the stochastic dynamical system

$$ \begin{cases} d\nu_1 = [-\nu_1 + \phi(\lambda_1 + w_+\nu_1 + w_I\nu_2)] dt + \beta dB_t^1 \\ d\nu_2 = [-\nu_2 + \phi(\lambda_2 + w_I\nu_1 + w_+\nu_2)] dt + \beta dB_t^2 \end{cases}, \quad (3.9) $$

with $\nu_1, \nu_2 \ge 0$. Here the applied stimuli are, $\lambda_1 = 33$, $\lambda_2 = \lambda_1 - \Delta\lambda$, with the bias $\Delta\lambda \in [0, 10^{-3}]$, the inhibitory connectivity coefficient is $w_I = 1.9$, the standard deviation of the Brownian motion is $\beta = 3 \cdot 10^{-3}$, and the excitatory connectivity coefficient $w_+$ is the bifurcation parameter, its range of values will be discussed later on. Moreover, the sigmoid (response) function $\phi(z)$ is given by

$$ \phi(z) = \frac{\nu_c}{1 + \exp(-bz + \alpha)}, $$

with $\nu_c = 15$, $b = 0.25$, and $\alpha = 11.1$.
<figure><img src="image_8.png" /><figcaption>Figure 6: Bifurcation diagrams with respect to $w_+$ for $\Delta\lambda = 0$ and $10^{-3}$. The spontaneous state is the black line. There are other two stable equilibria: the blue and light-blue lines. The zoom shows the behavior at the bifurcation point $w_+ = 2.5695$. Carrillo, Cordier, Deco, Mancini [43].</figcaption></figure>

Let us first consider the underlying dynamical system (3.9) without noise. In terms of $w_+$, it presents a subcritical Hopf bifurcation whose bifurcation diagram is shown in Figure 6. For values of the excitatory coefficient around $w_+ = 1.4$ the system goes from a single asymptotically stable equilibria to a situation in which there are three stable (continuous lines) and two unstable equilibria (hashed lines). The second local bifurcation, where the central asymptotic equilibria disappears, happens at around $w_+ = 2.5695$ for $\Delta\lambda = 10^{-3}$.

<figure><img src="image_9.png" /><figcaption>Figure 7: Change of variable from the phase space $(\nu_1, \nu_2)$ to $(x, y)$. Carrillo, Cordier, Deco, Mancini [43].</figcaption></figure>

An approximation of the slow manifold joining the spontaneous $S_0$ and the decision states $S_1$ and $S_2$ is needed (Figure 7). This curve is found by introducing the linearization of the dynamical system around the spontaneous equilibria $S_0$ as a new set of variables defined by $X = P^{-1}(\nu - S_0)$, with $X = (x, y)^T$, or equivalently $\nu = S_0 + PX$ with $\nu = (\nu_1, \nu_2)^T$. Here, $P$ is the matrix diagonalizing the jacobian of the dynamical system (3.1) at the spontaneous state $S_0$, see [45, 43] for more details.

This change of variables, sketched in Figure 7, is a natural way to introduce slow $y$ and fast $x$ variables in the system determined by the eigenvalues of the linearization at the spontaneous state $S_0$. By rewriting the dynamical system (3.9) above in these new variables together with the new time variable $\tau = \varepsilon t$, we have an evolution governed by

$$ \begin{cases} \varepsilon \frac{dx}{d\tau} = f(x, y) \\ \frac{dy}{dt} = g(x, y) \end{cases}, $$

where $\varepsilon$ is a small parameter. Assuming that the scaling ratio is zero, $\varepsilon = 0$, we find the implicit relation $f(x, y) = 0$ that determines the approximated slow manifold $x^*(y)$. Inverting this identity to find the
approximated curves restrict the set of possible values of the bifurcation parameter in order to keep positive firing rates, in the particular case chosen above, one needs $w_+ \ge 1.9$ approximately.

Given the approximation of the slow manifold $x^*(y)$, we can restrict the dynamics in (3.9) to a single effective Langevin equation. This equation is determined by a potential $G(y)$ obtained from the evaluation of the dynamics over the slow manifold approximation, leading to

$$dy = g(x^*(y), y) dt + \beta_y dW_t^2,$$

with $\beta_y$ properly obtained in terms of $\beta$ via the change of variables, see [45] for details. The effective potential is determined by the relation $\partial_y G(y) = -g(x^*(y), y)$.

<figure><img src="image_10.png" /><figcaption>Figure 8: Potential associated to the 1D reduction. Left: $\Delta\lambda = 10^{-3}$, right: zoom for $\Delta\lambda = 10^{-3}$ around the spontaneous state. Carrillo, Cordier, Deco, Mancini [43].</figcaption></figure>

Returning to the example (3.9), the previous described reduction strategy yields an approximated potential as in Figure 8 for the case of the subcritical Hopf bifurcation in Figure 6. This is in contrast with the approach in [235] leading to a potential of degree 4 approximating an interval around the central equilibrium point, see [235, Figure 2]. Our strategy is valid as long as the approximated slow manifold lies fully in the relevant biological range of positive rates leading to a restricted range of the bifurcation parameter $w_+$.

The 1D effective computed potential $G(y)$ for the biased case ($\Delta\lambda = 10^{-3}$) with respect to the slow variable $y$ is plotted in Figure 9 for $w_+ = 2.0, \dots, 2.9$ (left), together with a zoom at the spontaneous state (right). We note that increasing $w_+$ beyond the second bifurcation point, the spontaneous state pass from local minimum to local maximum.

<figure><img src="image_11.png" /><figcaption>Figure 9: Full 1D potential along the $y$ variable, for $w_+ = 2.0, \dots, 2.9$ and $\Delta\lambda = 10^{-3}$ (left) and a zoom around the spontaneous state (right). Carrillo, Cordier, Deco, Mancini [43].</figcaption></figure>

These reduced Fokker-Planck equations carry information about reaction times and performance of the decision tasks since one can compute probabilities for these events in terms of the stationary distributions. We refer for more details to [45, 43].