## Part II

# Hyperbolic equations of transport type

## 4 The Time Elapsed model for neural networks

Although the introduction of the model (NNLIF) was motivated by observations of self-sustained periodic activity in inhibitory networks with sparse firing, the Time Elapsed model (TE) described below, on the contrary, was elaborated from observations of periodic activity in excitatory networks with intermittent intense spiking.

Multiple neuron models have exhibited self-sustained oscillations [264, 215, 216, 153], leading to think that the core mechanism relies upon a few key aspects of neural dynamics. The Time Elapsed model focuses on one key point: the probability for a neuron to spike in a time window $[t, t+dt)$ can be modelled according to the time $s$ elapsed since the last action potential, the intrinsic refractory period of the neuron and the synaptic integration of the recent history of network activity.

We will first explain the constitutive ideas of Time Elapsed neuron models and propose a set of mathematical assumptions used throughout the section. Then, we will discuss the first attempts, through direct calculations and characteristics, at constructing solutions and obtaining convergence towards a stationary state. Many of these results are obtained in simplified cases where formulae, relative entropies and decay rates can be obtained via direct computations. Then, we present more abstract methods adapted from semigroup theory or the Doeblin-Harris machinery, which allow to prove results about existence of solutions and long-time asymptotics in more general frameworks. We then focus on the theoretical and numerical insights into the emergence and stability of periodic solutions. Finally, we discuss the rigorous derivation of this family of equations and present many variants with additional features.

### 4.1 Neurons structured by time elapsed since the last discharge

Summarising [215, 216], the article [205] introduces and studies rigorously a hyperbolic partial differential equation of transport type modelling an infinite neural network structured by the time $s$ elapsed since the last time a neuron has emitted an action potential. If $n^0(s)$ is the initial probability density for finding in the network at time $t=0$ a neuron which did not fire since a time $s$, then this density is transported at time 1 and neurons are continuously removed upon firing:

$$ \frac{\partial n}{\partial t}(s,t) + \frac{\partial n}{\partial s}(s,t) + p(s,t)n(s,t) = 0, $$

where $p(s,t)$ is at time $t$ the firing rate of neurons which did not fire since a time $s$. The firing function $p(s,t) = p(s, X(t))$ can depend both on the time $s$ elapsed since the last discharge and the past history of network activity $X(t)$. This model is a nonlinear version of renewal equations [125, 208, 116]. The key idea for the Time Elapsed model is that when $s$ is small for a neuron, i.e. when this neuron just fired, it is not able to fire again due to the refractory period. This refractory period has an absolute component (physical impossibility to fire) and a relative component dependent on the network activity. The Time Elapsed model aims at understanding how self-sustained oscillations can arise from the interplay between refractoriness of neurons and network activity.

The instantaneous firing rate of the network, i.e. the flux of neurons firing at time $t$ in the network, is given by

$$ N(t) = \int_{0}^{+\infty} p(s,t)n(s,t)ds. $$

Then those neurons that just fired must be reset at time elapsed $s=0$ through the boundary condition $n(0,t) = N(t)$. Network activity is defined through a memory kernel represented by the measure $\mu_X$,

$$ X(t) = \int_{0}^{t} N(t-\tau)d\mu_{X}(\tau). $$
The transmission can be delayed (smoothly if $\mu_X(\tau) = J\alpha(\tau)d\tau$ or singularly with $\mu_X = J\delta_{\tau=d}$, $d > 0$, which is $X(t) = JN(t-d)$), or instantaneous ($\mu_X = J\delta_{\tau=0}$, which is $X(t) = JN(t)$). In both cases, the strength of interconnections is represented by a parameter $J \ge 0$.

## 4.2 Mathematical setting and the choice of parameters

Let us summarise the Time Elapsed model into one PDE system.

$$ 
\left\{
\begin{array}{ll}
\displaystyle \frac{\partial n}{\partial t}(s,t) + \frac{\partial n}{\partial s}(s,t) + p(s, X(t)) n(s,t) = 0, & s, t > 0, \\
\\
\displaystyle N(t) = \int_{0}^{+\infty} p(s, X(t)) n(s,t) ds, \quad X(t) = \int_{0}^{+\infty} N(t-\tau) \mu_X(d\tau), & t \ge 0, \\
\\
n(0,t) = N(t), \quad n(s,0) = n^0(s) \ge 0, \quad \int_{0}^{+\infty} n^0(s) ds = 1, & s, t \ge 0.
\end{array}
\right.
\quad (\text{TE})
$$

Assuming enough regularity, direct integration yields the conservation of mass

$$ \forall t \ge 0, \quad \int_{0}^{+\infty} n(s,t)ds = \int_{0}^{+\infty} n^0(s)ds = 1. \quad (4.1) $$

For the time elapsed neuron model, this feature is at the core of many methods, as the use of characteristics lines in the free transport part provides handy representation formulae for $N(t)$ and $n(s, t)$.

The discharge function $(s, X) \mapsto p(s, X)$ represents the firing rate for neurons with age $s$, given that the feedback from activity in the network is $X$. In the excitatory case, it should increase with activity and with age (note that the inhibitory case, although considered later, has also been defined and studied rigorously [255]). We expect that, because of the refractory period, neurons will have zero probability to fire between $s = 0$ and some $X$-dependent age $\sigma(x)$. This firing rate should also be bounded. These properties can be summarised as

$$ \frac{\partial p}{\partial s} \ge 0, \quad \frac{\partial p}{\partial X} \ge 0, \quad p \equiv 0 \text{ on } [0, \sigma(X)], \quad p > 0 \text{ on } (\sigma(X), +\infty), \quad \|p\|_{\infty} < +\infty. $$

In order to perform explicit computations and take full advantage of the simple characteristics in the transport part of (TE), many studies make the following choice for $p(s, X)$.

**Assumption 4.1 (Simplified choice for $p(s, X)$).** The discharge function has the shape

$$ p(s, X) = 1_{s>\sigma(X)}, $$

where $\sigma: \mathbb{R}_+ \to \mathbb{R}_+$ is a Lipschitz function such that $\sigma'(x) \le 0$ almost everywhere and

$$ \forall X \in \mathbb{R}_+, \quad 0 < \sigma^- \le \sigma(X) \le \sigma^+ = \sigma(0) < 1. $$

In contrast, it is possible to choose a smooth function for the discharge function $p$.

**Assumption 4.2 (Smooth firing function $p(s, X)$).** The discharge function satisfies $p \in W^{2,\infty}(\mathbb{R}^2)$ and

$$ \frac{\partial p}{\partial s} \ge 0, \quad \frac{\partial p}{\partial X} \ge 0, \quad 0 < p^- = \lim_{s \to +\infty} p(s, 0) \le \lim_{s, X \to +\infty} p(s, X) = p^+ < +\infty. $$

Note that Assumption 4.2 excludes Assumption 4.1 because of the higher regularity it requires. Although existence and long time behaviour can be studied under more general assumptions, many studies choose well-behaved initial conditions in order to get global-in-time bounds for the solution and simplify the computations.

**Assumption 4.3 (Agreeable initial condition).** The initial condition satisfies $n^0 \in L^\infty(\mathbb{R}_+) \cap L^1(\mathbb{R}_+)$ and

$$ 0 \le n^0(\cdot) \le 1, \quad \int_{0}^{+\infty} n^0(s)ds = 1. $$
The impact of firing neurons on the dynamics is represented through the quantity $X(t)$, itself defined through the integration of the firing rate $N(t)$ under a measure $\mu_X$. The choice of this measure plays a critical role in the modelling goals and the behaviour of the solutions. The simplest choice is the non-delayed Time Elapsed neuron model where $\mu_X(d\tau) = J\delta_{\tau=0}$ for some connectivity strength $J \ge 0$. This translates to

$$X(t) = JN(t).$$

In the same fashion, one can choose a constant delay for the feedback between firing neurons and the dynamics, with $\mu_X(d\tau) = J\delta_{\tau=d}$ for some $d > 0$, yielding

$$X(t) = JN(t - d).$$

In larger networks where the transmission delay and synaptic integration are not evenly distributed, a random delay with distribution $\alpha$ and average connectivity strength $J$ can be considered.

**Assumption 4.4 (Smoothly delayed feedback).** The measure $\mu_X$ can be written as $\mu_X = J\alpha(\tau)d\tau$ with $J \ge 0$, $\alpha \in L^q(\mathbb{R}_+) \cap L^1(\mathbb{R}_+)$ for some $q \in (1, +\infty)$ and

$$\int_0^{+\infty} \alpha(\tau)d\tau = 1.$$

For some results we will require stonger decay on the delay feedback.

**Assumption 4.5 (Exponential delay).** The measure $\mu_X$ can be written as $\mu_X = J\alpha(\tau)d\tau$ with $J \ge 0$, $\alpha \in L^1(\mathbb{R}_+)$ and there exists $\lambda > 0$ such that

$$\int_0^{+\infty} \alpha(\tau)d\tau = 1, \quad \int_0^{+\infty} e^{\lambda\tau}(\alpha(\tau) + |\alpha'(\tau)|)d\tau < +\infty.$$

Solutions are to be understood in the sense of characteristics or in the sense of distributions, as is typical for transport equations. If the distributional solution has a density, for example if it belongs in $L^1$ for all times, then it will be a weak-solution in the following sense:

$$\int_0^{+\infty} \phi(s, T)n(s, T)ds = \int_0^{+\infty} \phi(s, 0)n^0(s)ds + \int_0^T \int_0^{+\infty} \left(\frac{\partial \phi}{\partial t}(s, t) + \frac{\partial \phi}{\partial s}(s, t)\right) n(s, t)dsdt - \int_0^T \int_0^{+\infty} \phi(s, t)p(s, X(t))n(s, t)dsdt + \int_0^T N(t)\phi(0, t)dt, \\ N(t) = \int_0^{+\infty} p(s, X(t))n(s, t)ds, \quad X(t) = \int_0^T N(t - \tau)\mu_X(d\tau),$$

for any test function $\phi \in C^1(\mathbb{R}_+ \times \mathbb{R}_+)$ with appropriate decay at $+\infty$.

Last, if we want to talk about solutions to (TE) in the sense of measures, we have to introduce the following notions.

Let $(E, \mathcal{A})$ be a measurable space. Denote $\mathcal{M}(E)$ the set of finite measures on $E$, $\mathcal{M}_+(E)$ the set of positive finite measures on $E$. Recall that the total variation norm $\|\cdot\|_{\text{TV}}$ on $\mathcal{M}(E)$ is defined by

$$\|\mu\|_{\text{TV}} = \int_E \mu_+ + \int_E \mu_-,$$

where $\mu_+, \mu_-$ is the Hahn-Jordan decomposition of the measure $\mu$ into its positive and negative parts. Then, $(\mathcal{M}(\mathbb{R}_+), \|\cdot\|_{\text{TV}})$ is a Banach space and the topological dual of $C_{\to 0}(\mathbb{R}_+)$. The weak-* topology on $\mathcal{M}(\mathbb{R}_+)$ can be metrised on any tight subset of $\mathcal{M}(\mathbb{R}_+)$ with bounded total variation by the bounded Lipschitz norm

$$\|\mu\|_{\text{BL}} = \sup_{\psi \in W^{1, +\infty}} \int_0^{+\infty} \psi(x)\mu(dx).$$

Then we denote $C([0, T], \mathcal{M}(\mathbb{R}_+))$ the set of functions on $[0, T)$ with values in $\mathcal{M}(\mathbb{R}_+)$ which are continuous for the bounded Lipschitz norm. We denote the translation semigroup generated on $(\mathcal{M}(\mathbb{R}_+), \|\cdot\|_{\text{BL}})$ by the operator $-\partial_s$, which means by abuse of notation $T_t n^0 = n(\cdot - t)$.
**Definition 4.6.** Consider the case $\mu_X = \delta_{\tau=0}$, i.e. $X(t) = N(t)$. A couple $(n, N)$ is said to be a measure solution to (TE) on $[0, T)$ with initial condition $n^0 \in \mathcal{M}_+(\mathbb{R}_+)$, $N^0 \in \mathbb{R}_+$ if $n \in C([0, T), \mathcal{M}_+(\mathbb{R}_+))$, $N \in C([0, T), \mathbb{R}_+)$, $n(0) = n^0$, $N(0) = N^0$, and for all $t \in [0, T)$,

$$n(t) = T_t n^0 + \int_0^t T_{t-\tau} (p(\cdot, N(\tau))n(\cdot, \tau)) (s) d\tau + \int_0^t T_{t-\tau} (N(\tau)\delta_{s=0}) d\tau.$$

$$\text{and } N(t) = \int_0^{+\infty} p(s, N(\tau))n(t)(ds).$$

In the definition, we write $n(\cdot, t)$ by abuse of notation even if $n(t) \in \mathcal{M}_+(\mathbb{R}_+)$ does not always have a density. The time integrals have to be understood as Bochner integrals. It can be readily checked that the crucial mass conservation (4.1) holds for measure solutions too, which implies by positivity $\|n(t)\|_{\text{TV}} = \|n^0\|_{\text{TV}}$. Note also that

$$\int_0^t T_{t-\tau} (N(\tau)\delta_{s=0}) (s) d\tau = \int_0^t N(\tau)\delta_{s=t-\tau} d\tau = N(t-s)1_{[0, +\infty)}(t-s),$$

and that (4.6) can be turned into the easier equivalent form(see arguments in [39]): for all $\phi \in C_c^\infty(\mathbb{R}_+)$,

$$\frac{d}{dt} \int_0^{+\infty} \phi(s)n(t)(ds) = \int_0^{+\infty} \frac{\partial \phi}{\partial s}(s)n(t)(ds) - \int_0^{+\infty} \phi(s)p(s, X(t))n(t)(ds) + N(t)\phi(0).$$

In the context of measure solutions, we introduce another assumption on the firing function $p$.

**Assumption 4.7 (Firing function for measure solutions).** The firing function $p$ belongs in $W^{1,\infty}(\mathbb{R}_+ \times \mathbb{R}_+)$, $p \ge 0$ and there exists $L > 0$ such that

$$\forall s, X_1, X_2 \ge 0, \quad |p(s, X_1) - p(s, X_2)| \le L|X_1 - X_2|.$$

**Assumption 4.8 (Preparation for decay in measures).** The firing function $p \in W^{1, +\infty}(\mathbb{R}_+ \times \mathbb{R}_+)$ satisfies, for some constants $\sigma, p^-, p^+ > 0$,

$$\forall s, X \ge 0, \quad \frac{\partial p}{\partial s}(s, X) \ge 0, \quad p^{-}1_{[\sigma, +\infty)} \le p(s, X) \le p^{+}.$$

Before proceeding with the description of the numerous analytical results on (TE), let us emphasise that two equivalent ways of making a weak-nonlinearity assumption are used in the literature: one can choose a Lipschitz condition on $p$ like $\|\sigma'\|_\infty$ small enough in Assumption 4.1 or $L$ small enough in Assumption 4.7, and then $X$ is not scaled as in $X(t) = N(t)$ or $X(t) = \int_0^t N(t-\tau)\alpha(\tau)d\tau$, $\int_0^{+\infty} \alpha(\tau)d\tau = 1$; or one could choose $J$ small enough with $X(t) = JN(t)$ or $J$ in Assumption 4.4 or Assumption 4.5. This switch between a Lipschitz condition on the discharge function $p$ and a smallness condition on the connectivity parameter $J$ allows to simplify computations depending on the method. We try, for each reference described, to stay compatible with the global framework proposed above for the Time Elapsed model, which may require some mental gymnastic on the part of the reader.

## 4.3 Global existence and well-posedness

### 4.3.1 Uniformly bounded solutions

In [205], global-in-time solutions are constructed via a fixed point argument. Let us choose the smoothly delayed assumption 4.4 and fix some time $T > 0$. Given some external input $Y(t)$ in place of the nonlinear $X(t)$, there is a unique solution $n$ with firing rate $N$ to the linear problem

$$\left\{
\begin{array}{ll}
\displaystyle \frac{\partial n}{\partial t}(s, t) + \frac{\partial n}{\partial s}(s, t) + p(s, Y(t)) n(s, t) = 0, & s, t > 0, \\
\\
\displaystyle N(t) = \int_{0}^{+\infty} p(s, Y(t)) n(s, t) ds, & t \ge 0, \\
\\
\displaystyle n(0, t) = N(t), \quad n(s, 0) = n^0(s) \ge 0, \quad \int_{0}^{+\infty} n^0(s) ds = 1, \quad s, t \ge 0.
\end{array}
\right.
\qquad (4.2)$$
Given such a solution, we can define the functional

$$ \begin{cases} F: L^{\infty}([0, T]) \rightarrow L^{\infty}([0, T]) \\ Y: \mapsto J \int_{0}^{+\infty} \alpha(\tau) N(t-\tau) d\tau. \end{cases} $$

Consider two external inputs $Y_1, Y_2 \in L^\infty([0, T])$, for $T$ small enough but independent from $n^0$. Because the differential operator is speed one transport, it is possible to exploit the characteristics: for all $s \ge t$, the initial condition is propagated along characteristics while being dampened by the loss of firing neurons,

$$ n(s,t) = n^0(s-t) e^{-\int_0^t p(\tau+s-t, Y(\tau)) d\tau}, \quad s \ge t. \quad (4.3) $$

On the contrary, when $s \le t$, the neuron have already fired and reset and come from the boundary condition at $s = 0$, which can be computed as

$$ n(s,t) = N(t-s) e^{-\int_0^s p(\tau, Y(\tau+t-s)) d\tau}, \quad s \le t. \quad (4.4) $$

Combining (4.3) and (4.4) gives the representation formula

$$ N(t) = \int_{0}^{t} p(s, Y(t)) N(t-s) e^{-\int_{0}^{s} p(\tau, Y(\tau+t-s)) d\tau} ds \\ + \int_{t}^{+\infty} p(s, Y(t)) n^{0}(s-t) e^{-\int_{0}^{t} p(\tau+s-t, Y(\tau)) d\tau} ds. $$

Applying this formula for both $Y_1$ and $Y_2$, it is then possible to decompose the difference $F(Y_1) - F(Y_2)$ into many terms, each of those one can control by the norm $\|Y_1 - Y_2\|_\infty$. It is possible to choose $T$ small enough independently of $n^0$ such that

$$ \|F(Y_1) - F(Y_2)\|_{\infty} \le \frac{1}{2} \|Y_1 - Y_2\|_{\infty}. $$

Hence, by applying the Banach fixed point theorem, there is a unique fixed-point for $F$, which yields a unique local-in-time solution of (TE). Since $T$ is independent from $n^0$, the procedure can be iterated to extend the local solution into a global-in-time solutions. Uniform bounds can be derived along the way.

**Theorem 4.9** (Pakdaman, Perthame, Salort [205]). *Grant Assumptions 4.1 and 4.3 and 4.4 (simple p, good n⁰, smooth delay). There exists a unique solution $n \in C(\mathbb{R}_+, L^1(\mathbb{R}_+))$ to (TE) and for all $s, t \ge 0$,*

$$ 0 \le n(s,t) \le 1, \quad 0 \le N(t) \le 1, \quad J(1-\sigma^+) \le X(t) \le J. $$

**Remark 4.10.** The authors of [205] claim that, if $p$ is taken to be Lipschitz continuous, then the same proof works when the smooth delay Assumption 4.4 is replaced by a fixed delay assumption $\mu = J\delta_{\tau=d}$, $d \ge 0$, $J \ge 0$. The case $d = 0$ yields the non-delayed Time Elapsed model where $X(t) = JN(t)$.

### 4.3.2 Solutions in the sense of measures

In [39], measures solution in the sense of definition 4.6 are constructed for well-behaved firing functions $p$ and initial conditions with small enough total variation norm, in the special case $X(t) = N(t)$.

Following [38, 205], a fixed point argument is applied in an appropriate space. Fix $C > 0$ and consider

$$ \mathcal{X} = \{n \in C([0, T], \mathcal{M}_+(\mathbb{R}_+)) \mid n(0) = n^0, \forall t \in [0, T], \|n(t)\|_{\text{TV}} \le C\}, $$

endowed with the norm $\|n\|_{\mathcal{X}} = \sup_{t \in [0, T]} \|n(t)\|_{\text{TV}}$. Beware that the continuity in $C([0, T], \mathcal{M}_+(\mathbb{R}_+))$ is not with respect to the total variation norm $\|\cdot\|_{\text{TV}}$ but according to the bounded Lipschitz norm $\|\cdot\|_{\text{BL}}$. The fixed point is sought for the operator $\Psi: \mathcal{X} \to \mathcal{X}$ defined as

$$ \Psi[n](t) = T_t n^0 + \int_0^t T_{t-\tau} (p(\cdot, N(\tau)) n(\cdot, \tau)) d\tau + \int_0^t T_{t-\tau} (N(\tau) \delta_{s=0}) (s) d\tau, $$
where

$$N(t) = \int_{0}^{+\infty} p(s, N(t))n(t)(ds).$$

An important step is to check that $N$ is uniquely defined, allowing the operator $\Psi$ to be well-defined as well.

**Lemma 4.11** (Cañizo, Yoldaş [39]). Let $n \in M_+(\mathbb{R}_+)$. Grant Assumption 4.7 with $L$ satisfying

$$L \le \frac{1}{4 \|n\|_{\text{TV}}}. \qquad (4.5)$$

Then there exists a unique solution $N \ge 0$ to the fixed-point problem $N = \int_0^{+\infty} p(s, N)n(ds)$. Given $n_1, n_2 \in M_+(\mathbb{R}_+)$ satisfying (4.5), the associated $N_1, N_2$ satisfy

$$|N_1 - N_2| \le \frac{\|p\|_{\infty}}{1 - L \|n_1\|_{\text{TV}}} \|n_1 - n_2\|_{\text{TV}}.$$

Using this lemma, the mapping $\Psi$ is well-defined and can be proved to be a contraction on $\mathcal{X}$ for $T$ small enough, yielding an existence result in the space of finite measures.

**Theorem 4.12** (Cañizo, Yoldaş [39]). Let $n^0 \in \mathcal{P}(\mathbb{R}_+)$. Grant Assumption 4.7 with

$$L \le \frac{1}{4 \|n^0\|_{\text{TV}}}. \qquad (4.6)$$

Consider the non-delayed case $\mu_X = \delta_{\tau=0}$, i.e. $X(t) = N(t)$. There exists a unique global-in-time measure solution $(n, N)$ to (TE) in the sense of Definition 4.6. If $n^{0,1}, n^{0,2} \in \mathcal{P}(\mathbb{R}_+)$ satisfying (4.6), then the associated solutions $n_1$ and $n_2$ satisfy

$$\forall t \ge 0, \quad \|n_1(t) - n_2(t)\|_{\text{TV}} \le e^{4\|p\|_{\infty}t} \|n^{0,1} - n^{0,2}\|_{\text{TV}}.$$

We enforce $\int n^0 = 1$ in this theorem for the sake of consistency with our formulation (TE) but this assumption is unnecessary in [39].

The Lipschitz condition on $p$ could be translated onto a connectivity parameter $J$ taking $X(t) = JN(t)$. Note that this smallness condition on $L$ is important, as will be seen in Section 4.7 where solutions with discontinuous firing rate are constructed using the non-uniqueness of $N(t)$ for large $L$ (translated to large $\|\sigma'\|_{\infty}$ in the context of Assumption 4.1). See also the Lipschitz assumption on $\sigma$ in Theorem 4.20.

### 4.3.3 Improved existence proof including the strongly inhibitory case

In the article [241], the existence results in [205] and [39] are revisited and improved. They proceed by considering a cutoff and solutions to an approximate problem, then going back to the full (TE).

In the non-delayed case, they obtain:

**Theorem 4.13** (Sepulveda, Torres, Villada [241]). Assume $\mu_X = \delta_{\tau=0}$, i.e. $X(t) = N(t)$, $p \in W^{1, +\infty}(\mathbb{R}_+ \times \mathbb{R}_+)$ and let $\gamma = \sup_{s,X \ge 0} \partial_X p(s, X)$. For all initial condition $n^0 \in L^1(\mathbb{R}_+) \cap \mathcal{P}(\mathbb{R}_+)$ such that $\gamma \|n^0\| < 1$, there exists a unique solution $n \in C(\mathbb{R}_+, L^1(\mathbb{R}_+))$, $N \in C(\mathbb{R}_+)$ to (TE).

It is claimed in [241] that the regularity can be relaxed to include step functions such as Assumption 4.1 or the case $p(s, X) = \phi(X)1_{s>\sigma}$ used in [255]. A key advantage of this result is that, because there is no absolute value or norm in the definition of $\gamma$, the condition $\gamma = \sup_{s,X \ge 0} \partial_X p(s, X)$ is always satisfied when $\gamma \le 0$, which included the strongly inhibitory framework that is studied for example in [255]. Similarly, for the model (NNLIF), solutions are always global-in-time in the inhibitory case as described in Section 1.

In the distributed delayed case, they obtain a more general existence result:
**Theorem 4.14** (Sepulveda, Torres, Villada [241]). Assume $\mu_X = \alpha(\tau)d\tau$, $\alpha \in L^1(\mathbb{R}_+)$, and $p \in W^{1,+\infty}(\mathbb{R}_+ \times \mathbb{R}_+)$. For all initial condition $n^0 \in L^1(\mathbb{R}_+) \cap \mathcal{P}(\mathbb{R}_+)$, there exists a unique solution $n \in C(\mathbb{R}_+, L^1(\mathbb{R}_+))$, $N, X \in C(\mathbb{R}_+)$ to (TE).

Note that the delayed assumption is slightly more general than Assumption 4.4 used in Theorem 4.9 as the delayed kernel is not assumed to be in some $L^q$, $q > 1$.

## 4.4 Study of the stationary equation

Under either Assumption 4.4 or for fixed delay $\mu_X = J\delta_{\tau=d}$, $J, d \ge 0$, a stationary firing rate implies that $X_\infty = JN_\infty$. Then stationary states are solutions of

$$ \frac{dn_{\infty}}{ds}(s) = -p(s, X_{\infty})n_{\infty}(s), \quad X_{\infty} = J \int_{0}^{+\infty} p(s, X_{\infty})n_{\infty}(s)ds, \quad \int_{0}^{+\infty} n_{\infty}(s)ds = 1, $$

which has solutions of the form

$$ n_{\infty}(s) = \frac{e^{-\int_{0}^{s} p(\eta, X_{\infty})d\eta}}{\int_{0}^{+\infty} e^{-\int_{0}^{s} p(\eta, X_{\infty})d\eta} ds}. $$

Under Assumption 4.1 for $p$, it simplifies down to the full characterisation

$$ n_{\infty}(s) = \frac{e^{-(s-\sigma(X_{\infty}))+}}{1 + \sigma(X_{\infty})}, \quad X_{\infty} = J(1 + \sigma(X_{\infty})). \qquad (4.7) $$

The number of stationary states thus depends on the number of solutions $X_\infty$ in (4.7); there can be several of them for some choices of $\sigma$. For small enough $J$, a Banach fixed point argument yields uniqueness. For larger $J$, the criterion is more involved.

**Theorem 4.15** (Pakdaman, Perthame, Salort [205]). *Grant Assumption 4.1 and assume either 4.4 or* $\mu_X = J\delta_{\tau=d}$, $J, d \ge 0$.

* If $J$ is small enough, there is a unique stationary state of (TE).

* There exists a unique stationary state of (TE) if $J$ is large enough and

$$ \limsup_{J \to +\infty} \sup_{X > \frac{J}{1+\sigma^+}} \frac{J|\sigma'(X)|}{(1+\sigma(x))^2} < 1. $$

If we choose a smoother function $p$, going for Assumption 4.2, we can also obtain existence of stationary states and uniqueness for small connectivities.

**Theorem 4.16** (Mischler, Weng [187]). *Grant Assumption 4.1 and assume either 4.5 or $\mu_X = J\delta_{\tau=0}$, $J \ge 0$. Then, there exists at least one stationary state ($n_\infty, N_\infty, X_\infty$) of (TE) such that for some constant $C > 0$,*

$$ \forall s \ge 0, \quad 0 \le n_{\infty}(s) \le C e^{-\frac{p}{2}s}, \quad \left| \frac{\partial n_{\infty}(s)}{\partial s} \right| \le C e^{-\frac{p}{2}s}. $$

There exists $J_0 > 0$ small enough such that for all $J \in [0, J_0)$, there is a unique stationary state to (TE).

This result has been extended to the strongly nonlinear case in [184]. In the context of measures solutions, stationary states can be found in the sense of Definition 4.6.

**Theorem 4.17** (Cañizo, Yoldaş [39]). *Grant Assumptions 4.7 and 4.8 and assume $\mu_X = \delta_{\tau=0}$, i.e.* $X(t) = N(t)$. If

$$ L < \frac{1}{\left(\frac{\sigma^2}{2} + \frac{\sigma}{p^-} + \frac{1}{(p^-)^2}\right) (p^+)^2}, $$

then there exists a unique stationary state of (TE) in the space of probability measures $\mathcal{P}(\mathbb{R}_+)$.

Further insights into existence and uniqueness of measure stationary states can be obtained via the Doeblin-Harris method, as described later in Subsection 4.6.3
## 4.5 Asymptotic relaxation to a stationary state for the step function case

We first discuss here the long-time asymptotic results under Assumption 4.1, when the firing function $p(s, X) = 1_{s>\sigma(X)}$ is a step function, or when we choose a firing function $p(s, X) = \phi(X)1_{s>\sigma}$ which separates into a network influence $\phi$ and a constant refractory period $1_{s>\sigma}$.

### 4.5.1 Convergence in relative entropy

**Linear case:** let us first showcase the method on the simpler linear case, which corresponds, under Assumption 4.1, to $\sigma(X) \equiv \sigma \in (0, 1)$ constant, and we write for clarity $\bar{p}(s) = p(s, X) = 1_{s>\sigma}$. Then, there is a unique stationary state to (TE) which can be written as

$$n_{\infty}(s) = \frac{e^{-\int_{0}^{s} \bar{p}(\tau)d\tau}}{\int_{0}^{+\infty} e^{-\int_{0}^{s} \bar{p}(\tau)d\tau}} ds.$$

Note that this function $n_{\infty}$ is constant on $[0, \sigma]$ and then is exponentially decaying.

Taking inspiration from the general methods in [179, 178], we want to find some function $\Psi$ and some $\mu > 0$ such that

$$\frac{d}{dt} \int_{0}^{+\infty} \Psi(s) |n(s, t) - n_{\infty}(s)| ds \leq -\mu \int_{0}^{+\infty} \Psi(s) |n(s, t) - n_{\infty}(s)| ds. \quad (4.8)$$

First, we can notice that by linearity, $\tilde{n} = n - n_{\infty}$ still satisfies (TE), and then it is also the case for $|\tilde{n}|$. By conservation of mass, $\int_{0}^{+\infty} \tilde{n}(s, t) ds = 0$ for all $t \geq 0$, which allows the computation, for any $\nu \in (0, 1)$,

$$|\tilde{n}(0, t)| = |\tilde{N}(t)| = \left| \int_{0}^{+\infty} \bar{p}(s) \tilde{n}(s) ds \right| = \left| \int_{0}^{+\infty} (\bar{p}(s) - \nu) \tilde{n}(s) ds \right|.$$

Using integration by part (with a density argument on $n^0$ with compact support if necessary),

$$\frac{d}{dt} \int_{0}^{+\infty} \Psi(s) |n(s, t) - n_{\infty}(s)| ds = \int_{0}^{+\infty} \left( \frac{d\Psi}{ds}(s) - \bar{p}(s)\Psi(s) \right) |n(s, t) - n_{\infty}(s)| ds \\ - \Psi(0) \left| \int_{0}^{+\infty} (\bar{p}(s) - \nu) \tilde{n}(s) ds \right|$$

Hence, it is possible to prove (4.8) by choosing $\Psi$ as a solution of

$$\frac{d\Psi}{ds}(s) - (\bar{p}(s) - \mu)\Psi(s) = -\Psi(0)|\bar{p}(s) - \nu|. \quad (4.9)$$

It is possible to prove that there is a choice of $\mu, \nu$ and $\Psi(0) = 1$ such that there exists a non-negative solution $\Psi$ to this problem.

**Theorem 4.18** (Pakdaman, Perthame, Salort [205]). *Grant Assumptions 4.1 and 4.3 with $\sigma(X) \equiv \sigma \in (0, 1)$ constant. There exists $\mu, B > 0$ and a non-negative function $\Psi: \mathbb{R}_+ \to \mathbb{R}_+$ such that $\Psi(0) = 1$, $\Psi(\cdot) \geq B$, $\Psi$ grows as $e^{\mu s}$ as $s \to +\infty$ and if $\int_{0}^{+\infty} \Psi(s) n^0(s) ds + \infty$, then for all $t \geq 0$,*

$$\int_{0}^{+\infty} \Psi(s) |n(s, t) - n_{\infty}(s)| ds \leq e^{-\mu t} \int_{0}^{+\infty} \Psi(s) |n^{0}(s) - n_{\infty}(s)| ds.$$

**Weakly-nonlinear and strongly nonlinear case:** For $J \neq 0$ small enough (and $\sigma$ not constant), it is still possible to construct the function $\Psi$ as a solution of (4.9), which can be interpreted as spectral gap property.

Denoting $f(t) = \int_{0}^{+\infty} \Psi(s) |n(s, t) - n_{\infty}(s)| ds$, it is possible to write a differential inequality of the form

$$\frac{df}{dt}(t) \leq -\mu f(t) + C \|\sigma'\|_{L^{\infty}([J(1-\sigma^{+}), J])} |X(t) - X_{\infty}|.$$
The difficulty is to control the difference $|X - X_{\infty}|$ which depends on past values of $N(t)$. However, according to Theorem 4.9, we have the control

$$J(1 - \sigma^{+}) \le X(t) \le J, \quad J(1 - \sigma^{+}) \le X_{\infty} \le J.$$

Hence $\frac{df}{dt} \le -\mu f + C \|\sigma'\|_{L^{\infty}([J(1-\sigma^{+}),J])} J$. Under either a smallness assumption on $J$ or taking $J$ large and assuming that $J \|\sigma'\|_{L^{\infty}([J(1-\sigma^{+}),J])}$ is small enough when $J$ is large enough, we can control the size of $f$ and feed this control into an iteration procedure for further refined estimation of $|X - X_{\infty}|$.

**Theorem 4.19 (Pakdaman, Perthame, Salort [205]).** *Grand Assumptions 4.1, 4.3 and 4.4. If either $J$ is small enough or if $J$ is large enough and the quantity*

$$\limsup_{J \to +\infty} \sup_{X \in [J(1-\sigma^{+}), J]} J|\sigma'(X)|$$

is finite and small enough, there exists $\mu, B > 0$ and a non-negative function $\Psi: \mathbb{R}_{+} \to \mathbb{R}_{+}$ such that $\Psi(0) = 1$, $\Psi(\cdot) \ge B$, $\Psi$ grows as $e^{\mu s}$ as $s \to +\infty$. For any initial condition such that $\int_{0}^{+\infty} \Psi(s)|n^{0}(s) - n_{\infty}(s)|ds < +\infty$, the solution of (TE) satisfies

$$\lim_{t \to +\infty} \int_{0}^{+\infty} \Psi(s)|n(s, t) - n_{\infty}(s)|ds = 0,$$

where $n_{\infty}$ is the unique stationary state.

Without more information on $\alpha$ in Assumption 4.4, it is not possible to conclude about the rate of convergence to the stationary states. However, in the simpler non-delayed case $\mu = J\delta_{\tau=0}$, i.e. $X(t) = JN(t)$, [205] proves under similar hypotheses that if $J$ is small enough, for some positive $\beta$,

$$\int_{0}^{+\infty} \Psi(s)|n(s, t) - n_{\infty}(s)|ds \le e^{-\beta t} \int_{0}^{+\infty} \Psi(s)|n^{0}(s) - n_{\infty}(s)|ds.$$

### 4.5.2 Convergence of the firing rate

Because of the simple characteristic flow of the transport operator $\partial_t + \partial_s$, it is also possible to work directly at the level of the firing rate $N(t)$. If we know that $N(t)$ converges to a stationary firing rate $N_{\infty}$, then we can conclude that the density $n$ will also become stationary as a solution to a linear equation with external force $X_{\infty}$. Following the approach of [206] we choose the simplified non-delayed case $\mu = J\delta_{\tau=0}$ with $J = 1$, which is equivalent to $X(t) = N(t)$. Note that the general case $J \ne 1$ can be retrieved by a suitable modification of $\sigma$. The connectivity parameter is thus transmitted onto the derivative of $\sigma$, which we will assume bounded:

$$0 \le -\sigma'(X) \le \|\sigma'\|_{\infty}.$$

**Linear case:** Consider again the linear case $\sigma(X) = \sigma \in (0, 1)$ constant, under Assumption 4.1. Note first that then the unique stationary state has the firing rate

$$N_{\infty} = \frac{1}{1 + \sigma}.$$

A fundamental idea is that for $t \ge s$, the initial condition is gone through the free speed one transport at point $(s, t)$ and if $s \le \sigma$ there is no firing yet. Hence

$$n(s, t) = N(t - s), \quad t \ge s, s \le \sigma.$$

Then the conservation of mass can be rewritten as

$$1 = \int_{0}^{+\infty} n(s, t)ds = \int_{0}^{\sigma} n(s, t)ds + \int_{\sigma}^{+\infty} 1_{s \ge \sigma}n(s, t)ds,$$
which yields

$$N(t) + \int_{t-\sigma}^{t} N(s)ds = 1, \quad t \ge \sigma. \qquad (4.10)$$

Because the stationary firing rate satisfies this identity and because, by Theorem 4.9, $0 \le N, N_\infty \le 1$, a contraction can be obtained by iteration:

$$N(t) - N_{\infty} = \int_{t-\sigma}^{t} (N_{\infty} - N(s))ds \implies \forall t \ge k\sigma, \quad |N(t) - N_{\infty}| \le \sigma^k.$$

Looking further at (4.10), we can also see that $N(t)$ oscillates around $N_\infty$ while exponentially converging towards it.

**Weakly-nonlinear case:** The same procedure can be applied when $\sigma$ is not constant but varies slowly $\|\sigma'\|_\infty < 1$. The idea is to obtain and use again a formula like (4.10). First, note that $N(t)$ is defined by

$$N(t) = \int_{\sigma(N(t))}^{+\infty} n(s,t)ds \qquad (4.11)$$

In the case $\|\sigma'\|_\infty < 1$, there is a unique solution to this problem. As we will see later, it is not the case in other parameter regimes and non-uniqueness of $N$ can induce oscillations. From this fixed point formula we have

$$\forall t \in \mathbb{R}_+, \quad N(t) + \int_{0}^{\sigma(N(t))} n(s,t)ds = 1. \qquad (4.12)$$

Using (4.12) and the equation for $n$, it is possible to prove that $N$ is a Lipschitz function. It then has a almost everywhere derivative that satisfies

$$N'(t)[1 + \sigma'(N(t))n(\sigma(N(t)), t)] = -N(t) + n(\sigma(N(t)), t) \quad \text{a.e. in } [\sigma^+, +\infty).$$

One can then obtain $N'(t)\sigma'(N(t)) \le \|\sigma'\|_\infty$, which leads to the nonlinear version of (4.10) after some tedious developments involving characteristics:

$$N(t) + \int_{t-\sigma(N(t))}^{t} N(s)ds = 1. \qquad (4.13)$$

Application of the Banach fixed point theorem to (4.13) proves that there is a unique stationary firing rate $N_\infty$. The argument above for the linear case can then be adapted to get the contraction

$$|N(t) - N_{\infty}| \le \left( \frac{\sigma^{+}}{1 - \|\sigma'\|_{\infty} N_{\infty}} \right)^{k}.$$

**Theorem 4.20 (Pakdaman, Perthame, Salort [206]).** *Grant Assumptions 4.1 and 4.3. Assume that $\mu = \delta_0$, i.e. $X(t) = N(t)$, and $\|\sigma'\|_\infty < 1$ and*

$$\sigma(x) < \sigma^{+} < 1 - \|\sigma'\|_{\infty} N_{\infty}.$$

*Then there exists a unique stationary firing rate $N_\infty$ and $\nu > 0$ such that the solution of (TE) satisfies* $\frac{d\sigma(N(t))}{dt} < 1$ and

$$\forall t \ge 2\sigma, \quad |N(t) - N_{\infty}| \le e^{-\nu t}.$$
### 4.5.3 Reduction to a delay differential equation

Consider the non-delayed case $X(t) = N(t)$, and take now a firing function with the form $p(s, X) = \varphi(X)1_{s>\sigma}$, where $\sigma > 0$ is a constant and $\varphi$ a positive bounded function $0 < p^- \le \varphi \le p^+$ (then, by definition of $N(t)$, $0 \le N(t) \le p^+$). The constant $\sigma > 0$ is a fixed refractory period. When $\varphi$ is non-decreasing, the network is excitatory as more activity implies more discharges after the refractory time is elapsed. On the contrary, if $\varphi$ is non-increasing, the network influence is inhibitory: more activity in the network diminishes the quantity of neurons undergoing discharges.

Following [255], let us define the auxiliary function $\Psi: \mathbb{R}_+ \to \mathbb{R}_+$ and its derivative by

$$ \psi(X) = \frac{X}{\varphi(X)}, \quad \psi'(X) = \frac{\varphi(X) - X\varphi'(X)}{\varphi^2(X)}. $$

As we will see, complex dynamics can occur only when $\psi'$ changes sign on $[0, p^+]$.

* In the inhibitory case $\varphi' \le 0$, we always have $\psi'(X) > 0$ for $X \in [0, p^+]$.

* In the excitatory case $\varphi' \ge 0$, we can distinguish between a weakly-excitatory case if $\psi'(X) > 0$ and a strongly-excitatory case when it changes sign on $[0, p^+]$. A sufficient condition for being in the weakly-excitatory case is for example $0 \le \varphi'(X) \le \frac{p^-}{p^+}$.

The main idea to exploit the properties of $\psi$ in the study of (TE) is to reduce the full system to a delayed differential equation on the firing rate $N(t)$ that involves $\psi'$. This can be done using mass conservation and the methods of characteristics.

**Lemma 4.21** (Cáceres, Perthame, Salort, Torres [255]). Assume $\mu_X = \delta_{\tau=0}$, i.e. $X(t) = N(t)$, $p(s, X) = \varphi(X)1_{s>\sigma}$, where $\varphi \in C^1(\mathbb{R}_+)$, $\sigma > 0$ and $0 < p^- \le \varphi \le p^+$. Then, given a solution to (TE) with initial condition $n^0$,

* for $0 < t < \sigma$,

$$ \int_{0}^{t} N(\tau)d\tau + \int_{0}^{\sigma-t} n^{0}(s)ds + \psi(N(t)) = 1, $$

and if $N$ is differentiable on $(0, \sigma)$,

$$ \psi'(N(t))N'(t) = -N(t) + n^{0}(\sigma - t); $$

* for $\sigma < t$,

$$ \int_{t-\sigma}^{t} N(\tau)d\tau + \psi(N(t)) = 1, \quad (4.14) $$

and if $N$ is differentiable on $(\sigma, +\infty)$,

$$ \psi'(N(t))N'(t) = -N(t) + N(t - \sigma). \quad (4.15) $$

Then, [255] proves that from a solution to the integral equation (4.14) and with compatibility conditions, it is possible to construct a solution to the full system (TE).

**Theorem 4.22** (Cáceres, Perthame, Salort, Torres [255]). Assume $\mu_X = \delta_{\tau=0}$, i.e. $X(t) = N(t)$, $p(s, X) = \varphi(X)1_{s>\sigma}$, where $\varphi \in C^1(\mathbb{R}_+)$, $\sigma > 0$ and $0 < p^- \le \varphi \le p^+$. Let $N \in L^\infty(\mathbb{R}_+)$ such that $\psi \circ N \in C([\sigma, +\infty)) \cap C^1((0, \sigma))$. Then if

* for all $t \in (0, \sigma)$, $\tilde{n}(t) := N(\sigma - t) + \frac{d\psi(N)}{dt}(\sigma - t) \ge 0$

* the function $N$ is solution to (4.14) on $[\sigma, +\infty)$.

for any initial condition $n^0$ of mass 1 such that for all $s \in (0, \sigma)$, $n^0(s) = \tilde{n}(s)$, there exists a unique solution to (TE) with firing rate $N$.
This theorem is proved by inputing $N$ in the linear system (4.2) and checking that the solution is also a solution to the nonlinear problem (TE).

From the delay differential equation formulation, we can conclude that in both the inhibitory and the weakly-excitatory case, the firing rate $N(t)$ oscillates an infinite number of times. Indeed, if we reason by contradiction and assume that the firing rate is, for example, increasing over an interval of length $\sigma$, then by (4.15),

$$ \psi'(N(t))N'(t) = -N(t) + N(t - \sigma) < 0. $$

Using that $\psi'(N(t)) > 0$, we get $N'(t) < 0$, which is a contradiction with $N$ increasing. This way, [255] proves that $\psi' > 0$ on $[0, p^+]$ implies the existence of an infinite quantity of local maxima and minima.

Then, it can be proved that the amplitude of oscillations is actually decreasing: define $I_k = [k\sigma, (k+1)\sigma]$, then

$$ \max_{t \in I_{k-1}} N(t) \ge \max_{t \in I_k} N(t), \quad \min_{t \in I_{k-1}} N(t) \le \min_{t \in I_k} N(t). $$

Last, defining the sequence $N_n(t) = N(t+n)$, by compactness there exists a subsequence converging to a limit $N_\infty(t)$ which defines a $C^1$ bounded function of time. Using $\psi' > 0$ and Lemma 4.21, $N_\infty$ is proved to be a constant function and in fact the firing rate of the only possible stationary state. Then, it it proved in [255] that $N(t) \to N_\infty$ and $n(\cdot, t) \to n_\infty$.

**Theorem 4.23** (Cáceres, Perthame, Salort, Torres [255]). Assume $\mu_X = \delta_{\tau=0}$, i.e. $X(t) = N(t)$, $p(s, X) = \varphi(X)1_{s>\sigma}$, where $\varphi \in C^1(\mathbb{R}_+)$, $\sigma > 0$ and $0 < p^- \le \phi \le p^+$. Assume that for all $X \in [0, p^+]$, $\psi'(X) > 0$. Then, there exists a unique stationary state $(n_\infty, N_\infty)$ to (TE) and for any solution to (TE), $N(t)$ converges to $N_\infty$ when $t \to +\infty$.

As we have seen, the method relies on the sign of $\psi'$. We will see later in 4.7 that when the system is not inhibitory nor weakly-excitatory, the oscillations can still happen but are not necessarily dampened towards the stationary state, as explicit periodic solutions can be constructed.

However, even in the strongly excitatory case, converge to a (now not necessarily unique) stationary state can happen: [255] proves that under technical conditions on a function $N \in C([0, \sigma])$, a solution $(n, N)$ that prolongs $N$ on $[0, +\infty)$ can be constructed using Theorem 4.22 and the extended $N$ monotonically converges towards a stationary state $N_\infty$. Numerical simulations in [255] show that it happens when the solution remains in an area where the sign of $\psi'$ is constant. For example, they choose

$$ \varphi(X) = \frac{1}{1 + e^{-9X + 3.5}}, \quad \sigma = \frac{1}{2}. $$

In this case, $\psi$ changes sign and there are three stationary states $N_\infty^1, N_\infty^2, N_\infty^3$. Even more surprising, all of them can attract solutions from the same initial density. Choose $n^0(s) = \frac{1}{2}e^{-(s-1)^+}$. Then, both the theoretical and the numerical solution are not uniquely defined, because the equation

$$ N(0) = \varphi(N(0)) \int_{\sigma}^{+\infty} n^0(s)ds $$

has three different solutions $N_1(0), N_2(0)$ and $N_3(0)$. Each of them leads to a different numerical solution $n_1, n_2, n_3$ converging to $N_\infty^1, N_\infty^2, N_\infty^3$ respectively. Two of them converge while oscillating and the third does so monotonically. For the same $\phi$ and $\sigma$, some other initial conditions exhibit jumps on the firing rate. Those jumps preserve the continuity of $\psi(N(t))$.

Contrary to the case of Assumption 4.1 described before, the speed of convergence is unknown in Theorem 4.23; [255] conjectures that the speed is exponential.

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

In the first works [205, 206] on the model (TE), asymptotic convergence is studied in the case where the firing function $(s, X) \mapsto p(s, X)$ has the simple form $p(s, X) = 1_{s \ge \sigma(X)}$. This choice allows for
explicit computations, predictable characteristics paths and an explicit construction for relative entropy decay. In later works [187, 39], other methods are introduced which allow to prove weakly-nonlinear stability of a stationary state for other families of firing functions p.

### 4.6.1 Abstract semigroup approach for the weakly-nonlinear problem

In [187], the strategy of perturbation of semigroup initiated in [182] and previously used in a neuroscience context in [183] was applied to the Time Elapsed neuron model. It allows to prove a general weakly-nonlinear stability result.

**Theorem 4.24** (Mischler, Weng [187]). *Grant Assumption 4.2. Assume either $\mu_X = J\delta_{\tau=0}$, $J \ge 0$, i.e. $X(t) = JN(t)$ or Assumption 4.5. Then, there exists $J_0 > 0$, $\nu > 0$, $C \ge 1$, $\eta > 0$ such that for all $J \in [0, J_0)$, for all initial condition $n^0 \in L^1(\mathbb{R}_+) \cap \mathcal{P}(\mathbb{R}_+)$ and satisfying $\|n^0 - n_\infty\|_{L^1} < \eta$,*

$$ \forall t \ge 0, \quad \|n(\cdot, t) - n_\infty(\cdot)\|_{L^1} \le Ce^{-\nu t}. $$

The method proceeds as follows.

**Case without delay.** Consider the case $\mu_X = J\delta_0$, which yields $X(t) = JN(t)$. Let us write the stationary state perturbation

$$ n(s, t) = n_{\infty}(s) + m(s, t), \quad N(t) = N_{\infty} + M(t). $$

Then the linearised system for the pair $(m, M)$ is

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial m}{\partial t} + \frac{\partial m}{\partial s} + p(s, JN_{\infty})m + J \frac{\partial p}{\partial X}(s, JN_{\infty})n_{\infty}M = 0, \\ M(t) = \int_{0}^{+\infty} \left( p(s, JN_{\infty})m(s, t) + J \frac{\partial p}{\partial X}(s, JN_{\infty})n_{\infty}(s)M(t) \right) ds, \\ m(0, t) = M(t), \quad m(s, 0) = m^0(s). \end{array} \right. \quad (4.16) $$

Define

$$ \kappa = J \int_{0}^{+\infty} \frac{\partial p}{\partial X}(s, JN_{\infty})n_{\infty}(s)ds. $$

For $J$ small enough, $\kappa < 1$, and solving for $M(t)$ in the second equation of (4.16) yields

$$ M(t) = \frac{1}{1-\kappa} \int_{0}^{+\infty} p(s, JN_{\infty})m(s,t)ds := \mathcal{M}_J[m(\cdot, t)]. $$

Defining this linear operator $\mathcal{M}_J$ allows to rewrite the linearised system with a linear operator acting on $m$:

$$ \frac{\partial m}{\partial t} = L_J m := -\frac{\partial m}{\partial s} - p(s, JN_{\infty})m - J \frac{\partial p}{\partial X}(s, JN_{\infty})n_{\infty} \mathcal{M}_J[m]. $$

A semigroup $S_{L_J}(t)$ with generator $L_J$ acting on $L^1(\mathbb{R}_+)$ is then classically associated to (4.16) and it has the domain

$$ D(L_J) = \{m \in W^{1,1}(\mathbb{R}_+) \mid m(0) = \mathcal{M}_J(m)\}. $$

The linearised system can also be rewritten as

$$ \frac{\partial m}{\partial t} = \Lambda_J m := -\frac{\partial m}{\partial s} - p(s, JN_{\infty})m - J \frac{\partial p}{\partial X}(s, JN_{\infty})n_{\infty} \mathcal{M}_J[m] + \delta_{s=0} \mathcal{M}_J[m], $$

allowing to define a semigroup $S_{\Lambda_J}(t)$ with generator $\Lambda_J$ acting on the space of bounded Radon measures

$$ \chi = M^1(\mathbb{R}_+) = \{m \in C_{\to 0}(\mathbb{R})' \mid \text{supp}(m) \subset \mathbb{R}_+\} $$
endowed with the weak-* topology $\sigma(M^1, C_{\to 0})$. It can be checked by duality that $S_{\Lambda_J|L^1(\mathbb{R}_+)} = L_J$.
A core idea of [187] is to split the operator $\Lambda_J$ as

$$ \Lambda_J = \mathcal{A}_J + \mathcal{B}_J, \quad \mathcal{A}_J m := \nu_J \mathcal{M}_J[m], \quad \mathcal{B}_J m := -\frac{\partial m}{\partial s} - p(s, JN_\infty)m, $$

where $\gamma_J = \delta_{s=0} + J \partial_X p(s, JN_\infty)n_\infty$. The authors of [187] then propose to adapt the strategy in [185, 181] relying on Weyl's theorem and the Spectral Mapping theorem. However, an obstacle is that $\Lambda_J$ does not generate a strongly continuous semigroup on $\mathcal{X}$, but only a weakly-* continuous semigroup. Among diverse ways to circonvent the issue, [187] chooses to rely on the splitting and the Banach structure of $L^1(\mathbb{R}_+) \subset \mathcal{X}$.

**Definition 4.25** (Hypodissipativity [124]). A closed operator $L$ on a Banach space $(E, \|\cdot\|)$ with dense domain $D(L)$ is $\beta$-hypodissipative if there exists a norm $|||\cdot|||$ on $E$ equivalent to $\|\cdot\|$ such that

$$ \forall x \in D(L), \exists \phi \in F_{|||\cdot|||}(x), \quad \mathcal{R}e \langle \phi | (L - \beta) x \rangle \le 0 $$

where

$$ F_{|||\cdot|||}(x) = \{\phi \in E' | \langle \phi | x \rangle = |||x|||^2 = |||\phi|||_{E'}\}. $$

An operator $L: D(L) \to E$ being hypodissipative is equivalent to the growth estimate

$$ \|S_L(t)\|_{\mathcal{B}(E)} \le Me^{\beta t}, $$

on the semigroup $S_L$ associated to $L$, where $\mathcal{B}(E)$ is the space of linear and bounded operators on $E$.

**Lemma 4.26** (Mischler, Weng [187]). *Grant Assumption 4.2.* Denote $\beta_0 = -\frac{p_-}{2}$. For all $J \ge 0$,

* $A_J \in \mathcal{B}(W^{-1,1}(\mathbb{R}_+), \mathcal{X})$.

* $\mathcal{B}_J$ is $\beta_0$-hypodissipative in both $L^1(\mathbb{R}_+)$ and $\mathcal{X}$.

* The family of operators $S_{\mathcal{B}_J} * A_J S_{\mathcal{B}_J}$ satisfies

$$ \forall \beta \in (\beta_0, +\infty), \quad \|(S_{\mathcal{B}_J} * A_J S_{\mathcal{B}_J})(t)\|_{\mathcal{B}(\mathcal{X}, \mathcal{Y})} \le C_\beta e^{\beta t}, $$

where $C_\beta > 0$, $\mathcal{Y} = BV(\mathbb{R}_+) \cap L_1^1(\mathbb{R}_+)$, $BV$ is the space of functions with bounded variation and $L_1^1$ the Lebesgue space with weight $s \mapsto \langle s \rangle$.

Applying this lemma, adapting Weyl's theorem and the Spectral Mapping theorem [185, 181] and taking advantage of the splitting introduced above, the article [187] characterises the spectrum $\Sigma(\Lambda_J)$ of $\Lambda_J$ on $\mathcal{X}$.

**Theorem 4.27** (Mischler, Weng [187]). *Grant Assumption 4.2.* Denote $\beta_0 = -\frac{p_-}{2}$. For all $J \ge 0$, the operator $\Lambda_J$ is the generator of a weekly-* continuous semigroup $S_{\Lambda_J}$ acting on $\mathcal{X}$ endowed with the weak-* topology. There exist a finite rank projector $\Pi_J$ which commutes with $S_{\Lambda_J}$, $j \in \mathbb{N}$ and complex numbers

$$ \xi_1, \dots, \xi_j \in \Delta_{\beta_0} := \{z \in \mathbb{C} \mid \mathrm{Re}z > \beta_0\}, $$

such that $\Sigma(\Lambda_J|\Pi_J\mathcal{X}) \cap \Delta_{\beta_0} = \{\xi_1, \dots, \xi_j\}$ and for any $\beta \in (\beta_0, +\infty)$, there exists $C_\beta$ such that the remainder semigroup satisfies

$$ \forall t \ge 0, \quad \|S_{\Lambda_J}(I_{\mathcal{X}} - \Pi_J)\|_{\mathcal{B}(\mathcal{X})} \le C_\beta e^{\beta t}. $$

In the linear case $J = 0$, the semigroup is positive and [187] proves the following version of the Krein-Rutman theorem: there exists $\beta < 0$, $C > 0$ such that $\Sigma(\Lambda_0) \cap \Delta_\beta = \{0\}$ and for any $m^0 \in \mathcal{X}$ of average 0,

$$ \forall t \ge 0, \quad \|S_{\Lambda_0}(t)m^0\|_{\mathcal{X}} \le C e^{\beta t} \|m^0\|_{\mathcal{X}}. $$

Then, it is possible to use the stability theory for semigroups [147, 182, 186, 261] to apply a perturbation argument and prove stability for the weakly connected linearised system.
**Theorem 4.28** (Mischler, Weng [187]). *Grant Assumption 4.2. There exists $J_0$, $\beta < 0$, $C > 0$ such that for all $J \in [0, J_0)$, $\Sigma(\Lambda_J) \cap \Delta_\beta = \{0\}$ and for any $m^0 \in \mathcal{X}$ of average 0,*

$$ \forall t \ge 0, \quad \|S_{\Lambda_0}(t)m^0\|_{\mathcal{X}} \le C e^{\beta_0 t} \|m^0\|_{\mathcal{X}}. $$

For the weakly-nonlinear case, consider the functional

$$ \Phi[m, N] = \int_{0}^{+\infty} p(s, JN) - N. $$

Using the $W_1$ Wasserstein distance on the set $\mathcal{P}(\mathbb{R}_+)$ of probability measures, [187] proves (similar to the key uniqueness result in [206] under Assumption 4.1) that for $J$ small enough, there is a function $\phi_J: \mathcal{P}(\mathbb{R}_+) \to \mathbb{R}$, Lipschitz continuous for the weak topology of probability measures, such that $N = \phi_J(m)$ is the unique solution to $\Phi(m, N) = 0$, and thus the unique stationary firing rate associated to the perturbation $m$. Then, using this result, the decay for the linearised system and a Grönwall estimate, it is possible to prove that the perturbation $m$ decays exponentially fast in $L^1(\mathbb{R}_+)$ for the nonlinear problem, hence proving Theorem 4.24 in the non-delayed case.

**Case with delay.** Consider the case of a delayed interaction with an exponentially decaying regular time kernel as in Assumption 4.5. The strategy described above for the non-delayed case can be adapted using a modified functional framework. First, to write the linearised system as a time-autonomous problem, introduce a local time $\tau$ and, given the perturbation $M(t)$ of the firing rate $N(t) = N_\infty + M(t)$, transport the value in time through

$$ \frac{\partial v}{\partial t} + \frac{\partial v}{\partial \tau} = 0, \quad v(0, t) = M(t), \quad v(\tau, 0) = 0. $$

The solution is $v(\tau, t) = M(t - \tau)1_{0 \le \tau \le t}$. The variation $x$ of $X = X_\infty + Y(t)$ is characterised by

$$ Y(t) = \mathcal{D}[v] := J \int_{0}^{+\infty} v(\tau, t)\alpha(\tau)d\tau $$

and we have, given the perturbation $n(s, t) = n_\infty(s) + m(s, t)$, the linearised equation

$$ M(t) = \int_{0}^{+\infty} p(s, X_{\infty})m(s, t)ds + Y(t)J \int_{0}^{+\infty} \frac{\partial p}{\partial X}(s, X_{\infty})n_{\infty}(s)ds, $$

which can be rewritten in operator form

$$ M(t) = \mathcal{O}_J[m, v] := \mathcal{N}_J[m] + \kappa \mathcal{D}[v], \quad \kappa = \int_0^{+\infty} \frac{\partial p}{\partial X}(s, X_\infty)n_\infty(s)ds. $$

This leads to rewriting the linearised system in the delayed case in the form

$$ \frac{\partial}{\partial t} \begin{pmatrix} m \\ v \end{pmatrix} = \mathcal{L}_J \begin{pmatrix} m \\ v \end{pmatrix}, \quad \mathcal{L}_J \begin{pmatrix} m \\ v \end{pmatrix} = \begin{pmatrix} -\partial_s m - p(s, X_\infty)m - \partial_X p(s, X_\infty)n_\infty \mathcal{D}[v] \\ -\partial_\tau v \end{pmatrix}, $$

where the domain of $\mathcal{L}_J$ is

$$ \mathcal{D}(\mathcal{L}_J) = \{(m, v) \in W^{1,1}(\mathbb{R}_+) \times W^{1,1}(\mathbb{R}_+, \omega) \mid m(0) = v(0) = \mathcal{O}_J[m, v]\}, \quad \omega(s) = e^{-\lambda s}, $$

with $\lambda > 0$ from Assumption 4.5. The semigroup $S_{\mathcal{L}_J}$ acts on $L^1(\mathbb{R}_+) \times L^1(\mathbb{R}_+, \omega)$. As in the non-delayed case, let us introduce the larger space $\mathcal{X} = M^1(\mathbb{R}_+) \times M^1(\mathbb{R}_+, \omega)$, and the extended operator $\Lambda_J$ defined as

$$ \Lambda_J \begin{pmatrix} m \\ v \end{pmatrix} = \begin{pmatrix} -\partial_s m - p(s, X_\infty)m - \partial_X p(s, X_\infty)n_\infty \mathcal{D}[v] + \delta_{s=0} \mathcal{O}_J[m, v] \\ -\partial_\tau v + \delta_{\tau=0} \mathcal{O}_J[m, v] \end{pmatrix}. $$

It is again possible to decompose it in the form $\Lambda_J = \mathcal{A}_J + \mathcal{B}_J$ and to obtain a semigroup satisfying $S_{\Lambda_J|L^1(\mathbb{R}_+) \times L^1(\mathbb{R}_+, \omega)} = S_{\mathcal{L}_J}$. From there, it is possible to adapt the arguments in the non-delayed case and to prove the other part of Theorem 4.24. Let us mention the main intermediary results about the operator $\Lambda_J$ and its spectrum, as they are useful in themselves.
**Lemma 4.29** (Mischler, Weng [187]). *Grant Assumptions 4.2 and 4.5. Let $\beta_0 = \max(-\frac{p_-}{2}, -\lambda) < 0$.*

* $A_J \in \mathcal{B}(W^{-1,1}(\mathbb{R}_+) \times W^{-1,1}(\mathbb{R}_+, \omega), \mathcal{X})$.

* $B_J$ is $\beta_0$-hypodissipative in both $L^1(\mathbb{R}_+) \times L^1(\mathbb{R}_+, \omega)$ and $\mathcal{X}$.

* The family of operators $S_{B_J} * A_J S_{B_J}$ satisfies

$$ \forall \beta \in (\beta_0, +\infty), \quad \|(S_{B_J} * A_J S_{B_J})(t)\|_{\mathcal{B}(\mathcal{X}, \mathcal{Y})} \le C_\beta e^{\beta t}, $$

$$ \text{where } C_\beta > 0, \mathcal{Y} = BV(\mathbb{R}_+) \cap L_1^1(\mathbb{R}_+) \times BV(\mathbb{R}_+, \omega) \cap L_1^1(\mathbb{R}_+, \omega). $$

Using this lemma, it is possible to obtain the same spectral result as in the non-delayed case.

**Theorem 4.30** (Mischler, Weng [187]). *Grant Assumptions 4.2 and 4.5. The conclusion of Theorem 4.27 holds with $\beta_0 = \max(-\frac{p_-}{2}, -\lambda) < 0$.*

### 4.6.2 Back to non-smooth firing functions and the strong connectivity case

In [184], the method of [187, 269] and ideas from [205, 212] are combined with refined spectral theory [185], leading to convergence results for more general firing functions $p$ and allowing to consider the strong connectivity regime. It namely extends the strong connectivity convergence result of [205] described above in Theorem 4.19 and extends the weakly-nonlinear convergence results of Theorem 4.24 to less regular firing functions like those satisfying Assumption 4.1.

Denote $L_q^1(\mathbb{R}_+)$ the space of measurable functions $f$ such that $\|f\|_{L_q^1(\mathbb{R}_+)} := \int_0^{+\infty} (1+s^q)|f|ds < +\infty$ and $L_w^1$ the space $L^1(\mathbb{R}_+)$ endowed with the weak topology $\sigma(L^1, L^\infty)$.

Let us define a generalised weak nonlinearity condition.

**Assumption 4.31.** There exists $L > 0$ small enough such that for all $x_0 > 0$ there exists $J_0 > 0$ such that

$$ \forall J \in (0, J_0), \forall X_1, X_2 \in [0, \frac{x_0}{J}), \quad \int_0^{+\infty} |p(s, X_1) - p(s, X_2)|ds \le \frac{L}{J}|X_1 - X_2|. $$

Similarly we have the strong nonlinearity condition

**Assumption 4.32.** There exists $L > 0$ small enough such that for all $x_\infty > 0$ there exists $J_\infty > 0$ such that

$$ \forall J \in (J_\infty, +\infty), \forall X_1, X_2 \in (\frac{x_\infty}{J}, +\infty), \quad \int_0^{+\infty} |p(s, X_1) - p(s, X_2)|ds \le \frac{L}{J}|X_1 - X_2|. $$

Note that these assumptions cover many of the firing functions in Assumption 4.1 used in [205, 206].

**Theorem 4.33** (Mischler, Quiñinao, Weng [184]). *Assume $p \ge 0$, $\partial_s p \ge 0$, $\partial_X p \ge 0$,*

$$ 0 < \lim_{s \to +\infty} p(s, 0) \le \lim_{s, X \to +\infty} p(s, X) < +\infty, $$

and for all $s \ge 0$, $\int_0^s p(w, \cdot)dw \in C(\mathbb{R}_+)$. Then for all $q > 0$, for all $n^0 \in L_q^1(\mathbb{R}_+) \cap L^\infty(\mathbb{R}_+) \cap \mathcal{P}(\mathbb{R}_+)$, if one of the following three conditions hold:

1. $\mu_X = J\alpha$ with $\alpha \in L^1(\mathbb{R}_+)$ satisfying $\int_0^{+\infty} \alpha(\tau)e^{\delta\tau}d\tau < +\infty$;

2. $\mu_X = J\delta_{s=0}$ with $J > 0$ small enough and Assumption 4.31 holds;

3. $\mu_X = J\delta_{s=0}$ with $J > 0$ large enough, Assumption 4.32 holds and $\int_0^{+\infty} p(s, 0)n^0(s)ds > 0$;
there exists a global-in-time weak solution $n \in C(\mathbb{R}_+, L_w^1(\mathbb{R}_+)) \cap L^\infty(\mathbb{R}_+, L_q^1(\mathbb{R}_+)) \cap L^\infty(\mathbb{R}_+^2)$, $N, X \in C(\mathbb{R}_+)$ to (TE). When case 2. or 3. is assumed, this solution is unique.

The reason the solution is not proved to be unique in case 1. is that a Schauder fixed point theorem is used. Cases 2. and 3. yield uniqueness as they rely upon the Banach fixed point theorem.

The proof of existence of stationary states ($n_\infty, N_\infty, X_\infty$) with $n_\infty \in W^{1,+\infty}(\mathbb{R}_+)$ in [187] is also extended to more general firing functions $p$ in [184]. When $J$ is small enough or large enough, they prove that there exists a unique stationary state.

**Theorem 4.34** (Mischler, Quiñinao, Weng [184]). *Grant the assumptions of Theorem 4.33. There exists $J_0 > 0$ small enough and $J_\infty > 0$ large enough, $\nu > 0$ and $C \ge 1$ such that if $J \in (0, J_0) \cup (J_\infty, +\infty)$, the solution to (TE) satisfies*

$$ \forall t \ge 0, \quad \|n(\cdot, t) - n_\infty\|_{L^1} \le Ce^{-\nu t}. $$

### 4.6.3 Convergence using the Doeblin-Harris method

In [39], a general convergence result is achieved in the space of probability measures using the Doeblin-Harris machinery. After the seminal [132], these methods have been revisited between 2000 and 2010, namely in [116, 250, 128], and they have since been applied to a variety of linear and nonlinear partial differential equations; see [273] for a review. We start by reproducing here the introduction provided by [39].

**Definition 4.35** (Stochastic operators and semigroups). Let $(E, \mathcal{A})$ be a measurable space. We call a stochastic operator any linear operator $S : \mathcal{M}(E) \to \mathcal{M}(E)$ such that

$$ \forall \mu \in \mathcal{M}(E), \quad \int_E S\mu = \int_E \mu, \quad \mu \ge 0 \implies S\mu \ge 0. $$

A stochastic semigroup is a semigroup $(S_t)_{t \ge 0}$ of operators $S_t : \mathcal{M}(E) \to \mathcal{M}(E)$ such that for all $t \ge 0$, $S_t$ is a stochastic operator. An equilibrium or stationary state of a stochastic semigroup is a finite measure $\mu_\infty \in \mathcal{M}(E)$ such that for all $t \ge 0$, $S_t\mu_\infty = \mu_\infty$.

Equivalently, stochastic operators on $E$ preserve the space of probability measures $\mathcal{P}(E)$.

**Definition 4.36** (Doeblin condition [116]). Let $(E, \mathcal{A})$ be a measurable space. We say that a stochastic operator $S$ satisfies the Doeblin condition if there exists $\alpha \in (0, 1)$ and a probability measure $\nu \in \mathcal{P}(E)$ such that,

$$ \forall \mu \in \mathcal{P}(E), \quad S\mu \ge \alpha\nu. $$

**Theorem 4.37** (Doeblin's theorem [116]). Let $(E, \mathcal{A})$ be a measurable space and $(S_t)_{t \ge 0}$ a stochastic semigroup on $\mathcal{M}(E)$. If there exists $t_0 > 0$ such that $S_{t_0}$ satisfies the Doeblin condition with constant $\alpha \in (0, 1)$, then the semigroup has a unique equilibrium $n_\infty \in \mathcal{P}(E)$, and

$$ \forall n^0 \in \mathcal{P}(E), \quad \|S_t(n^0 - n_\infty)\|_{\text{TV}} \le \frac{1}{1-\alpha} e^{-\lambda t} \|n^0 - n_\infty\|_{\text{TV}} $$

$$ \text{where } \lambda = -\frac{1}{t_0} \log(1 - \alpha) > 0. $$

In order to apply Doeblin's theorem, [39] focuses first on the linear problem $p(s, X) = \bar{p}(s)$. Then, as seen in Subsections 4.3.1 and 4.5.2 for a simpler $\bar{p}$, the method of characteristics provides explicit estimates on the solution and allows to construct solutions without any smallness assumption. This extends to the framework of measure solutions: in the linear case, there exists a global-in-time measure solution to (TE) for any probability $n^0 \in \mathcal{P}(\mathbb{R}_+)$ and any $L$ in Assumption 4.7. Using the linear characteristics and performing lower bounds on them, the Doeblin condition can be proved in the form

$$ n(2\sigma) \ge p^{-}\sigma e^{-2p^{+}\sigma} 1_{(0,\sigma)} $$
**Lemma 4.38** (Cañizo, Yoldaş [39]). *Grant Assumptions 4.7 and 4.8. Assume for all s, X ≥ 0, p(s, X) = p̄(s). Then there exists a stochastic semigroup (S_t)_{t≥0}, such that for all n⁰ ∈ P(R⁺), t ↦ S_t n⁰ is the unique measure solution to (TE). Moreover, (S_t)_{t≥0} satisfies the Doeblin condition with t₀ = 2σ and α = p⁻σe⁻²p⁺σ, ν = 1/σ¹(0,σ).*

Then, Doeblin's theorem allows to straightforwardly conclude that the linear Time Elapsed equation has a spectral gap in the space of probability measures.

For the weakly-nonlinear case, the idea is to rewrite and decompose the evolution problem as

$$ \frac{\partial n}{\partial t} = \mathcal{L}_N n = \mathcal{L}_{N_\infty} n + h $$

with the operator

$$ \mathcal{L}_N n = -\frac{\partial n}{\partial s} - p(\cdot, N)n + \delta_{s=0} \int_0^{+\infty} p(u, N)n(du). $$

and the perturbation

$$ h = [p(s, N_\infty) - p(s, N)]n + \delta_{s=0} \int_0^{+\infty} [p(s, N) - p(s, N_\infty)]n(du) $$

It can be proved that under a smallness condition on L (condition (4.18) below),

$$ \|h(t)\|_{\text{TV}} \le \frac{2p^+L}{1-L} \|n(t) - n_\infty\|_{\text{TV}}, \quad \int_0^{+\infty} h(t)(ds) = 0 \quad (4.17) $$

The second step is to use the Duhamel formula

$$ n(t) = S_t n^0 + \int_0^t S_{t-\tau} h(\cdot, \tau) d\tau, $$

which yields after subtracting the stationary state and taking the total variation norm

$$ \|n(t) - n_\infty\|_{\text{TV}} \le \|S_t(n^0 - n_\infty)\|_{\text{TV}} + \left\| \int_0^t S_{t-\tau} h(\cdot, \tau) d\tau \right\|_{\text{TV}}. $$

Then, using the linear decay result for the first term in the right-hand side, (4.17) and Grönwall's inequality for the second part of the right-hand side, one can obtain the following result.

**Theorem 4.39** (Cañizo, Yoldaş [39]). *Let n⁰ ∈ P(R⁺). Grant Assumptions 4.7 and 4.8 and μ_X = δ_τ=0 with*

$$ L \le \frac{1}{4 \|n^0\|_{\text{TV}}}. \quad (4.18) $$

Then, there exists a unique measure stationary solution (n_∞, N_∞) in the sense of Definition 4.6 and there exists constants C ≥ 1, λ > 0 independent of n⁰ such that the solution to (TE) satisfies

$$ \forall t \ge 0, \quad \|n(t) - n_\infty\|_{\text{TV}} \le C e^{-\lambda t} \|n^0 - n_\infty\|_{\text{TV}}. $$

The constants C and λ in this result are constructive; in [39], explicit formulae are provided.

# 4.7 Self-sustained oscillations

## 4.7.1 Explicit construction in the step function case

Theorem 4.20 can be considered sharp in some sense, as there are analytic examples of periodic solutions when its assumptions are not satisfied.

Consider again the non-delayed case X(t) = N(t) and the step firing function p(s, X) = 1_{s>σ(X)} as in Assumption 4.1. Given the firing rate N(t), the system (TE) becomes linear. Hence, a way to
construct a periodic solution is to find a suitable T-periodic function $t \mapsto N(t)$ such that the solution to the linear hyperbolic boundary-value problem

$$ \frac{\partial n}{\partial t} + \frac{\partial n}{\partial s} + p(s, N(t)) = 0, \quad n(0, t) = N(t), \quad n(s, 0) = n(s, T), $$

is a solution to the full system (TE). Because there is no access to $n$, the mass conservation has to be translated as a property on $N$ like in (4.13). More precisely, one can prove:

**Theorem 4.40** (Pakdaman, Perthame, Salort [206]). *Grant Assumption 4.1 and assume $\mu_X = \delta_0$, i.e. $X(t) = N(t)$. Let $N \in L^\infty(\mathbb{R})$ a periodic function such that, in the sense of distributions,*

$$ \frac{d}{dt}(\sigma \circ N) \le 1, \quad N(t) + \int_0^{\sigma(N(t))} N(t-s)ds = 1. \quad (4.19) $$

Then there exists a periodic solution of (TE) with firing rate $N$.

Proving existence of periodic solutions to (TE) is then reduced to finding periodic functions satisfying (4.19). Such functions can have discontinuity points. Then, $\frac{d}{dt}(\sigma \circ N)$ in the sense of distributions means that $\sigma(N)$ has a negative jump, which is a positive discontinuity jump for $N$ given that $\sigma$ is non-increasing by Assumption 4.1. The reason for these discontinuities is the possibility for non-uniqueness of $N$ in the fixed point formulation (4.11) when $\|\sigma'\|_\infty$ exceeds one. Those positive jumps are a feature of the system and a core mechanism in non-delayed periodic solutions to (TE), as can be seen on Figure 10 and the explicit analytical examples below.

<figure><img src="image_12.png" /><figcaption>Figure 10: Numerical simulation of the firing rate $N(t)$ of a solution to (TE) with $\mu_X = J\delta_{\tau=0}$, $J = 20$, $\sigma$ as in (4.22) with $\alpha = 3$, $p(s, X) = 1_{s \ge \sigma(X)}$, $n^0(s) = e^{-s}$, original simulation by us, numerical scheme from [206].</figcaption></figure>

The article [206] proposes two kinds of explicit construction using Theorem 4.40.

First, choose a number $n$ of discontinuities and an increasing sequence $0 = \alpha_0 < \alpha_1 < \dots < \alpha_{n+1}$. Then, define the quantities

$$ N_n^+ = \frac{1}{1 + \sum_{j=1}^{n-1} (e^{\alpha_{j+1} - \alpha_j} - 1) + \alpha_{n+1} - \alpha_n}, \quad N_i^+ = e^{\alpha_{i+1} - \alpha_i} N_n^+, \quad 0 \le i \le n-1, $$

and the periodic function

$$ N(t) = \begin{cases} N_i^+ e^{\alpha_i t}, & t \in (\alpha_i, \alpha_{i+1}), \quad 0 \le i \le n-1 \\ N_n^+, & t \in (\alpha_n, \alpha_{n+1}), \\ N(t - k\alpha_{n+1}), & t \in [k\alpha_{n+1}, (k+1)\alpha_{n+1}), \quad k \in \mathbb{N}, k > 0. \end{cases} \quad (4.20) $$
If we define

$$ \underline{N} = N_n^+, \quad \bar{N} = \sup_{0 \le i \le n} N_i^+, \quad \sigma_0 = \alpha_{n+1} - \alpha_1 + \log \left( \frac{N_0^+}{\bar{N}} \right), $$

and the Lipschitz function

$$ \sigma(X) = \begin{cases} \sigma_0 - \log(\underline{N}) + \log(\bar{N}), & X \in [0, \underline{N}], \\ \sigma_0 - \log(x) + \log(\bar{N}), & X \in [\underline{N}, \bar{N}], \\ \sigma_0, & X \in [\bar{N}, +\infty]. \end{cases} \quad (4.21) $$

there exists a periodic solution to (TE) with firing rate N as defined in (4.20)

The authors of [206] construct a second kind of periodic solution. Choose a number $\alpha > 0$ and set

$$ \underline{N} = \frac{1}{2e^{\alpha} - 1}, \quad \bar{N} = \frac{e^{\alpha}}{2e^{\alpha} - 1}. $$

Similar to (4.21) we use the Lipschitz function

$$ \sigma(X) = \begin{cases} 2\alpha, & X \in [0, \underline{N}], \\ \alpha - \log(x) + \log(\bar{N}), & X \in [\underline{N}, \bar{N}], \\ \alpha, & X \in [\bar{N}, +\infty]. \end{cases} \quad (4.22) $$

Then, choose any $q \in (0, \alpha)$. There is a solution $\delta = \delta(q) \in [q + \alpha, 2\alpha]$ to the problem

$$ e^{\alpha - \delta} (-2\alpha + 1 + e^{\alpha} + \delta) + e^{-\alpha} (\delta - q - \alpha + e^q - 1) = 2 $$

and there exists a solution of (TE) with firing rate the $2\alpha$-periodic function defined on $[0, 2\alpha]$ by

$$ N(t) = \begin{cases} \bar{N}e^{-t}, & t \in [0, \alpha], \\ e^q \underline{N}e^{-t+\alpha}, & t \in (\alpha, q + \alpha], \\ \underline{N}, & t \in (q + \alpha, \delta(q)], \\ \bar{N}e^{-t+\alpha}(t - 2\alpha + e^{\alpha}), & t \in (\delta(q), 2\alpha]. \end{cases} $$

Note the particular case $q = \alpha$, $\delta(q) = 2\alpha$.

These numerous analytical examples and further numerical insights [205, 206] show that for fixed parameters there can be a great variety of periodic solutions. In the delayed cases, the firing rate may be periodic without exhibiting discontinuities. Characterising the stability profile of all those co-existing periodic solutions is an open question.

### 4.7.2 Construction from a delayed differential equation

Let us now come back to the framework of Subsection 4.5.3, which is the non-delayed case $\mu_X = \delta_{\tau=0}$ ($X(t) = N(t)$), and the firing function $p(s, X) = \varphi(X)1_{s>\sigma}$, $\sigma > 0$, $0 < p^- \le \phi \le p^+$. Recall the auxiliary functions

$$ \psi(X) = \frac{X}{\varphi(X)}, \quad \psi'(X) = \frac{\varphi(X) - X\varphi'(X)}{\varphi^2(X)}. $$

We have seen already that when $\psi' > 0$ on $[0, p^+]$ (inhinitory case or weakly-excitatory case) there is only convergence to the stationary state.

However, in the strongly excitatory case, which is when $\psi'$ changes sign, explicit examples of periodic solutions are described in [255]. The main idea is to use the reconstruction theorem 4.22 to obtain a periodic solution to (TE) from a well-chosen function N on $[0, \sigma]$.

The first kind of solution [255] finds is piecewise constant $\sigma$-periodic solutions.
**Theorem 4.41** (Cáceres, Perthame, Salort, Torres [255]). Assume $\mu_X = \delta_{\tau=0}$, i.e. $X(t) = N(t)$, $p(s, X) = \phi(X)1_{s>\sigma}$, where $\phi \in C^1(\mathbb{R}_+)$, $\sigma > 0$ and $0 < p^- \le \phi \le p^+$. Assume $\psi'$ changes sign on $[0, p^+]$. Let $N_1, N_2 \in (0, p^+]$ such that $N_1 \ne N_2$ and $\psi(N_1) = \psi(N_2)$. If there exists $\beta \in (0, \sigma)$ such that

$$ \beta N_1 + (\sigma - \beta)N_2 + \psi(N_1) = 1, \quad (4.23) $$

then there exists a $\sigma$-periodic solution to (TE) with firing rate defined on $[0, \sigma)$ by

$$ N(t) = \begin{cases} N_1, & t \in [0, \beta), \\ N_2, & t \in [\beta, \sigma). \end{cases} $$

Note that the condition (4.23) is equivalent to satisfying (4.14) for this piece-wise-constant stationary state. When such a $\beta$ exists, there is necessarily at least one stationary state

$$ N_{\infty} \in [\min(N_1, N_2), \max(N_1, N_2)]. $$

A similar construction can be made with more than one jump in $[0, \sigma)$, provided there are values $N_1, \dots, N_k$ such that $\psi(N_1) = \dots = \psi(N_k)$. In the special case where there exists an interval $[a, b] \subset [0, p^+]$ where $\phi$ is linear, i.e. $\phi(X) = CX$ on $[a, b]$, and if

$$ a\sigma < 1 - \frac{1}{C} < b\sigma, $$

then the same method allows to construct a smoother $\sigma$-periodic solution with firing rate

$$ N(t) = \frac{1}{\sigma} \left( 1 - \frac{1}{C} \right) + \beta u(t). $$

where $u(t)$ is any $\sigma$-periodic function such that $\int_0^\sigma u(t)dt = 0$ and the constant $\beta > 0$ is chosen small enough depending upon $u$ to keep $N(t)$ in $[a, b]$.

A second, more involved construction is proposed in [255], which allows to obtain piece-wise smooth $2\sigma$-periodic solutions to (TE). Similar to what we have described for the delayed NNLIF model in Section 1.7, it is possible to use arguments from [126] in the study of the delay equation

$$ \frac{d\psi \circ N}{dt}(t) = -N(t) + N(t - \sigma) \quad (4.24) $$

that we find in Lemma 4.21. Assume that $\psi$ is strictly convex around a local minimum $\underline{N}$. Then we can choose $\underline{N}^- < \underline{N} < \underline{N}^+$ such that $\psi(\underline{N}^-) = \psi(\underline{N}^+)$, $\psi$ decreasing on $[\underline{N}^-, \underline{N}]$ and increasing on $[\underline{N}, \underline{N}^+]$. Then the idea of [255] is to seek a periodic solution of (4.24) satisfying

* $N \in (\underline{N}^-, \underline{N})$ and $\underline{N}' < 0$ on $(0, \sigma)$; $N \in (\underline{N}, \underline{N}^+)$ and $\underline{N}' < 0$ on $(\sigma, 2\sigma)$;

* $N(0) = \underline{N}^+$, $N((2\sigma)^-) = \underline{N}^-$;

* $\psi(\underline{N}(\sigma^-)) = \psi(\underline{N}(\sigma^+))$.

Then, the solution on $[\sigma, 2\sigma)$ is given by the solution on $[0, \sigma)$ as an initial condition in the DDE (4.24). Such a procedure can be reduced to finding a fixed-point of a well-chosen functional $T: C \to C$ defined over the Banach space

$$ C = \{N \in C[0, \sigma] | N(0) = \underline{N}^+, N \text{ non-increasing, } \psi \circ N \ge \psi(\underline{N}^+) - (\underline{N}^+ - \underline{N}^-)\sigma\}. $$

Given $N \in C$, the corresponding extension on $[\sigma, 2\sigma)$ with the DDE (4.24) is chosen as a solution $M$ to the backward ODE problem

$$ \frac{d\psi \circ M}{dt}(t) = -M(t) + N(t - \sigma), \quad M(2\sigma) = \underline{N}^{-}. $$
Then, this defines a new portion $L$ of the solution on $[2\sigma, 3\sigma]$ as a solution to the forward ODE

$$ \frac{d\psi \circ L}{dt}(t) = -L(t) + M(t - \sigma), \quad L(2\sigma) = N^+. $$

The natural choice to get a periodic solution of (4.24) as a fixed point of $T$ is then to define

$$ T[N](t) = L(t + 2\sigma). $$

When $\sigma$ is small enough, it can be proved that $T$ is a contraction on $C$. Provided additional technical assumptions, a solution to the full system (TE) can be reconstructed from this periodic solution of (4.24) via Theorem 4.22.

**Theorem 4.42** (Cáceres, Perthame, Salort, Torres [255]). Assume $\mu_X = \delta_{\tau=0}$, i.e. $X(t) = N(t)$, $p(s, X) = \phi(X)1_{s>\sigma}$, where $\phi \in C^1(\mathbb{R}_+)$, $\sigma > 0$ and $0 < p^- \le \phi \le p^+$. Assume $\psi'$ changes sign on $[0, p^+]$ and there exist $N, \varepsilon > 0$ such that $\psi$ is strictly convex on $(N - \varepsilon, N + \varepsilon)$ and $N$ is a local minimum of $\psi$. Then,

* If $\sigma$ is small enough, there exists a $2\sigma$-periodic solution $N$ of (4.24) such that $\psi \circ N \in W^{1,\infty}$, $N$ is decreasing on $(0, \sigma)$ and on $(\sigma, 2\sigma)$ with a discontinuity at $t = \sigma$.

* If in addition $\psi(N) < 1$ and $\psi(N \pm \varepsilon) > 1$, then there exists a solution of (TE) with firing rate $N$.

<figure><img src="image_13.png" /><figcaption>Figure 11: Construction of a periodic firing rate yielding a solution to (TE) from a local minimum of $\psi$ as in Theorem 4.42. The arrows represent the path of $N(t)$ over a cycle of length $2\sigma$. The horizontal lines are the jump discontinuities at $t = \sigma$ and $t = 2\sigma$; note that they preserve the continuity of $\psi(N)$. Cáceres, Perthame, Salort, Torres [255].</figcaption></figure>

Note that given the fixed point operator is defined through a succession of ODEs, depending on the regularity of $\psi$, more regularity can be proved on the solution $N$ on $(0, \sigma)$ and $(\sigma, 2\sigma)$, meaning that this solution can be in fact smooth by part around the discontinuity point $t = \sigma$. The same procedure allows to construct solutions from a local maximum of $\psi$, but this time the solution will be increasing by part.

Numerical simulations in [255] confirm that different types of periodic solutions can be attractive. The solution to (TE) converges to a periodic solution with constant by part firing rate for the choice

$$ \varphi(X) = \max(\min(1.6X, 1), 0.25), \quad \sigma = 1, \quad n^0(s) = e^{-s}. $$
With the choice

$$ \varphi(X) = \frac{10X^2}{1 + X^2}, \quad \sigma = 1, \quad n^0(s) = e^{-(s-1)} 1_{s>1}. $$

the solution converges towards a periodic solution with a jump discontinuity and two non-constant paths in each period.

Note that only periodic solutions of period a multiple of $\sigma$ have been constructed so far. Given that multiple periodic solutions are possible for a same choice of $\psi$ and $\sigma$, [255] conjectures that the stable solutions are the one with few jump discontinuities. The speed of convergence towards the stable periodic solutions is unknown.

### 4.7.3 Further numerical evidence with theoretical guarantees on the scheme

In the article [241], numerical schemes are carefully constructed for the model (TE) in both the delayed and non-delayed cases. They propose first order upwind schemes based on a finite volume framework. These schemes are proved to be stable, consistent and convergent for small enough compactly supported initial conditions in the space $BV(\mathbb{R}_+)$ and $W^{1,+\infty}$ firing functions $p$.

The improved and theoretically validated numerical simulations in [241] confirm the previous observations of stable periodic solutions [206, 255] and exhibit some other cases. For example, periodic solutions are obtained in the delayed strongly inhibitory case in the framework of [255] with singular delay $d = \frac{1}{2}$, with

$$ p(s, X(t)) = \varphi(X)1_{s>\sigma}, \quad X(t) = N\left(t - \frac{1}{2}\right), \quad \phi(N) = e^{-9N}, \quad \sigma = \frac{1}{2}, \quad n^0(s) = \frac{1}{2}e^{-(s-1)^+}. $$

These oscillations are also observed in delayed settings. Moreover, new delayed excitatory oscillations are shown in [241] with detailed information about the numerical density $n(s, t)$.

## 4.8 Rigorous derivation from microscopic models

In [59], a rigorous connection is made between point processes used by neuroscientists to fit spike data (Poisson, Wold, Hawkes) and the macroscopic Time Elapsed model for neural networks.

Consider a locally finite point process $(T_z)_{z \in \mathbb{Z}}$ on $(\mathbb{R}, \mathcal{B}(\mathbb{R}))$, i.e. a random set of points such that there is almost surely a finite number of points in any bounded interval. For a measurable set $A \in \mathcal{B}(\mathbb{R})$, $N_A \in \mathbb{N} \cup \{+\infty\}$ denotes the number of points in $A$. Then, $t \mapsto N_t := N_{(0,t]}$, defined on $\mathbb{R}_+$, is the counting process associated to $N$. The natural filtration of $N$ is the family $(\mathcal{F}_t)_{t \in \mathbb{R}}$ of sigma-algebras $\mathcal{F}_t = \sigma(N \cap (-\infty, t])$ and the predictable filtration is $(\mathcal{F}_{t-}^N)_{t \in \mathbb{R}}$ defined by $\mathcal{F}_{t-}^N = \sigma(N \cap (-\infty, t))$. The former contains all the information given by the point process at times $t$ and the latter only the information strictly before $t$, thus does not contain information about whether or not $t$ is a point. In most cases, a stochastic intensity $\lambda(t, \mathcal{F}_{t-}^N)$ can be associated to the point process $(N_t)_{t \ge 0}$. Then, $\lambda(t, \mathcal{F}_{t-}^N)dt$ represents the probability to have a new point in $[t, t + dt)$ given the past points.

Those point processes here correspond to spike times. Then, $T_{N_t}$ is the last discharge for the neuron that occurred at or before time $t$. The time since the last discharge is then $S_t = t - T_{N_t} \ge 0$, and for all $t \ge 0$, $S_t$ is $\mathcal{F}_{t-}^N$ measurable. However, as the conditional intensity cannot use knowledge exactly at $t$, and thus cannot use $S_t$, we also need the variable $S_{t-} = t - T_{N_{t-}}$ which is $\mathcal{F}_{t-}^N$ measurable; $T_{N_{t-}}$ is the time of the last spike emitted strictly before $t$. Given this framework, the firing function $p(s, X(t))$ and the stochastic intensity $\lambda(t, \mathcal{F}_{t-}^N)$ play the same role: they give the firing rate of neurons with time $s$ elapsed since the last discharge given the full information about strictly past activity. The article [59] then proposes three different point processes to model the time since the last discharge.

* Poisson processes: if the Poisson process is homogeneous, which is $\lambda(t, \mathcal{F}_{t-}^N) = \lambda$ independent of $t$ and of the past activity, it corresponds to $p(s, X(t)) = \lambda$ as

$$ \mathbb{P}(S_{t-} > s) = \mathbb{P}(N_{[t-s,t)} = 0) = e^{-\lambda s}, $$
which corresponds to the stationary distribution $n_∞(s) = λe^{-λs}$ in (TE) with $p(s, X(t)) = λ$. The Poisson process can then be chosen inhomogeneous by letting $λ(t)$ depend on $t$. If we want to take into account the refractory period, the intensity has to depend on age $λ(t, \mathcal{F}_{t-}^N) = f(S_{t-})$. It is then equivalent to the setting $p(s, X(t)) = \bar{p}(s)$, which is the linear case of (TE). This renewal process enjoys i.i.d. inter-spike intervals (ISI).

* Wold processes: consider $A_t^1$ the time between the last discharge and the penultimate discharge. Then $A_t^1 = T_{N_{t-}} - T_{N_{t-}-1}$ is a Wold process characterised by $λ(t, \mathcal{F}_{t-}^N) = f(S_{t-}, A_t^1)$ [159]. This would correspond to $p(s, X(t)) = p(s, a)$ where a is the time elapsed since the penultimate discharge. This extension of (TE) has later been investigated at the PDE level in [256] (see also Subsection 4.9.4 below). One could define in the same way the intervals $A_t^k$ between the last $k^{th}$ and $(k+1)^{th}$ discharges, which was also to be investigated at the PDE level in [88] (see Subsection 4.9.5)

* Hawkes processes: let us consider now a univariate Hawkes process

$$λ(t, \mathcal{F}_{t-}^N) = μ + \int_{-∞}^{t-} h(t - τ)N(dτ),$$

with $μ$ the spontaneous rate and $h$ an interaction kernel supported in $\mathbb{R}_+$ such that $\int_0^{+∞} h(τ)dτ < 1$. Under mild assumptions the second part can be rewritten as $\int_0^{t-} h(t - τ)N(dτ)$ and represents the influence of the network activity $X(t)$ in $p(s, X(t))$. Although we expect that $p(s, X(t)) = μ + X(t)$ in this framework, [59] proves that it does not hold and provides more specific rigorous information about the link between these Hawkes processes and the model (TE).

Rigorous links between these point processes and the model (TE) are proved in [59] after technical considerations beyond the scope of this review. Their main tool is Ogata's thinning algorithm [163, 200].

In [66], the approach of [59] is generalised and a rigorous link between generalised Hawkes processes and the Time Elapsed model is proved.

Another microscopic approach has been proposed for the Time Elapsed model: in [218] a jump process interacting through a set of random global activity variables is proved to enjoy propagation of chaos and to converge toward (TE). The convergence rate is quantified by use of a coupling method.

## 4.9 Variants of the Time Elapsed model

### 4.9.1 Inhomogeneous distribution of the refractory period

In many cases, the refractory period can be divided into an absolute refractory period (physiological impossibility to fire) and a relative refractory period (the firing threshold is increased), ranging on average from 0.5 to 1 ms for the former and around 10ms for the later. Different types of nerve cells can have different absolute and relative refractory periods [115].

A possible extension of the Time Elapsed neuron model is then to add an other structuring variable $r \in \mathbb{R}$, and to choose a firing function $p(r, s, X)$ where for each $r \in \mathbb{R}$ the refractory period is different. The article [146] proposed the model

$$\left\{
\begin{array}{l}
\displaystyle \frac{\partial n}{\partial t}(r, s, t) + \frac{\partial n}{\partial s}(r, s, t) + p(r, s, X(t)) n(r, s, t) = 0, \\
\\
\displaystyle N(r, t) = \int_{0}^{+\infty} p(r, s, X(t)) n(r, s, t) ds, \quad X(t) = J \int_{-\infty}^{+\infty} N(r, t) dr, \\
\\
\displaystyle n(r, 0, t) = N(r, t), \quad n(r, s, 0) = n^0(r, s) \ge 0, \quad \int_{-\infty}^{+\infty} \int_{0}^{+\infty} n^0(r, s) ds = 1.
\end{array}
\right.
\qquad (4.25)$$

For the firing function, they make the choice

$$p(r, s, X) = \begin{cases} a(r) & s \le \sigma(r, X), \\ b(r) & s > \sigma(r, X), \end{cases} \quad a(r) < b(r), \qquad (4.26)$$
with

$$0 < a_0 = \inf_{r \in \mathbb{R}} a(r), \quad b_0 = \sup_{r \in \mathbb{R}} b(r) < +\infty, \quad \sigma \in L^\infty(\mathbb{R}, C_b^1(\mathbb{R}_+)), \quad \partial_X \sigma \le 0. \quad (4.27)$$

The quantity $\sigma(r, X)$ represents the refractory period for neurons of type $r$ under global network activity $X$. The firing rate during the refractory period is taken to be a small value $a(r)$ and not 0, in part to account for stochastic variability in firing and in part for purely mathematical reasons.

From the initial condition $n^0$, the distribution in $r$ is $\rho(r) = \int_0^{+\infty} n^0(r, s)ds$; it is preserved over time by the equation. Let us then assume

$$0 \le n^0 \le b(r)\rho(r). \quad (4.28)$$

The article [146] claims that standard arguments and the previous work [205] allow to construct a distributional solution $n \in C(\mathbb{R}_+, L^1(\mathbb{R}))$ such that $0 \le n(r, s, t) \le b(r)\rho(r)$. Stationary states are then of the form

$$n_{\infty} = \frac{\rho(r)}{\frac{1}{a(r)} - \left(\frac{1}{a(r)} - \frac{1}{b(r)}\right) e^{-a(r)\sigma(r, X_{\infty})}}.$$

There exists at least one stationary network activity $X_\infty$ defining a stationary state $(n_\infty, N_\infty, X_\infty)$, and [146] proves that if the average connectivity strength $J$ is small enough, the stationary network activity, and thus the stationary state, is unique.

Then, similar to the standard model (TE), it is possible to prove weakly-nonlinear stability of the (then unique) stationary state.

**Theorem 4.43 (Kang, Perthame, Salort [146]).** Assume the conditions (4.26), (4.27) and (4.28). If there exists $\beta > 0$ such that

$$\forall r \in \mathbb{R}, \quad a(r) - \frac{2Jb(r) \left\| \int_{-\infty}^{+\infty} \rho(r)b(r)(b(r) - a(r))\partial_X \sigma(r, \cdot) dr \right\|_{\infty}}{1 - J \left\| \int_{-\infty}^{+\infty} b(r)(b(r) - a(r))\partial_X \sigma(r, \cdot) dr \right\|_{\infty}} \ge \beta > 0, \quad (4.29)$$

then for all $t \ge 0$,

$$\int_{-\infty}^{+\infty} \int_{0}^{+\infty} |n(r, s, t) - n_{\infty}(r, s)| ds dr \le e^{-\beta} \int_{-\infty}^{+\infty} \int_{0}^{+\infty} |n^0(r, s) - n_{\infty}(r, s)| ds dr.$$

Note that (4.29) will be satisfied for $J$ small enough given (4.27).

## 4.9.2 Adaptation and fatigue of neurons modelled as a fragmentation kernel

An important aspect of neuronal dynamics which is overlooked in the models (TE) and (4.25) is that the refractory period of neurons can depend on their past activity. Following previous literature on fragmentation-type equations [93, 210, 177, 37, 208], the article [207] proposes to incorporate a fragmentation kernel into the Time Elapsed model.

$$\left\{
\begin{array}{l}
\displaystyle \frac{\partial n}{\partial t}(s, t) + \frac{\partial n}{\partial s}(s, t) + p(s, N(t)) n(s, t) = \int_{0}^{+\infty} K(s, a)p(a, N(t))n(a, t)da, \\[2ex]
\displaystyle N(t) = \int_{0}^{+\infty} p(s, N(t))n(s, t)ds, \\[2ex]
\displaystyle n(0, t) = 0, \quad n(s, 0) = n^0(s) \ge 0, \quad \int_{0}^{+\infty} n^0(s)ds = 1.
\end{array}
\right.
\quad (4.30)$$

The kernel $K(s, a)$, written as a function by abuse of notation, can be any positive measure with total mass 1 on $\mathbb{R}_+ \times \mathbb{R}_+$. Its values represent the probability density to reset to the state $s$ a neuron which fired while being in a state $a$. The variable $s$ then becomes a more general state variable not necessarily related to actual time elapsed since the last discharge. Note that the $a$-independent choice $K = \delta_{s=0}$ corresponds<fnref n="5" /> to system (TE) with $X(t) = N(t)$. An other simple choice would be $K = \delta_{s=\psi(a)}$, which

<footnote n="5">⁵In such case, the condition $n(0, t) = 0$ has to be overlooked or replaced by $n(0, t) = N(t)$.</footnote>
means that upon firing a neuron is deterministically reset to a state determined as a function $\psi(a)$ of its previous state $a$. Taking $(s, a) \mapsto K(s, a)$ as a smooth function instead adds randomness to which state neurons are reset to.

For the function $p$, [207] chooses a regularisation of Assumption 4.1.

**Assumption 4.44 (Structural assumptions on the firing function).** Assume there exists $p_M > 0$ such that for all $N$,

$$0 \le p(s, N) \le p_M, \quad \left\| \frac{\partial p}{\partial s} \right\|_{\infty} + \left\| \frac{\partial^2 p}{\partial s^2} \right\|_{\infty} + \left\| \frac{\partial^2 p}{\partial s \partial N} \right\|_{\infty} + \int_{0}^{+\infty} \left| \frac{\partial p}{\partial s}(s, N) \right| ds < +\infty,$$

and a bounded measurable function $\sigma : (0, +\infty) \to (0, \sigma^+]$ such that $p(s, N) \equiv p_M$ for $s \ge \sigma(N)$ and $p(s, N) < p_M$ for $s < \sigma(N)$.

Regarding the fragmentation kernel $K$, the natural properties to enforce are to not reset neurons at a larger time elapsed than just before they fire, and to preserve the mass of the equation (redistribute the firing neurons without creation nor annihilation). Also, following [161], an important assumption is that neurons firing a short time after the previous discharge suffer from fatigue and are reset at lower value $s$, but neurons which didn't fire since a long time are less fatigued and are reset at larger values of $s$, thus enjoying a shorter refractory period. A last assumption is that the change of state is significant when neurons fire. All these assumptions are summarised in mathematical form in Assumption 4.45 hereafter.

**Assumption 4.45 (Structural assumptions on the fragmentation kernel).** Assume that the measure $K$ satisfies, by abuse of notation

$$K(s, a) \ge 0, \quad K(s, a) = 0, s > a, \quad \int_{0}^{a} K(s, a) ds = 1.$$

and, there exists $\theta \in [0, 1)$ such that

$$-\frac{\partial}{\partial a} \int_{0}^{s} K(w, a) dw \ge 0, \quad \frac{\partial}{\partial a} \int_{0}^{s} s K(s, a) ds \le \theta < 1.$$

As we have seen for the standard model (TE), it is not to be taken for granted that the firing rate $N(t)$ is uniquely defined and multiple solutions for $N$ can be the driving force for the arising of periodic solutions. Smallness assumptions can however enforce uniqueness both for the time-dependent firing rate $N(t)$ and the stationary firing rate $N_\infty$. Since (4.30) does not contain a connectivity parameter, a weak-nonlinearity assumption has to be expressed as a condition on the variations of $p$. If there exists $\eta$ small enough such that

$$\left\| \frac{\partial p}{\partial N} \right\|_{\infty} < \eta, \qquad (4.31)$$

then the firing rate $N(t)$ is uniquely defined and [207] proves that there exists a unique stationary state $(n_\infty, N_\infty)$. We can then define the functions and quantities $p_\infty(s) = p(s, N_\infty)$, $\sigma_\infty = \sigma(N_\infty)$, and

$$B_\infty = e^{\int_0^{\sigma_\infty} p_\infty(s) ds} \left[ \sigma_\infty \int_0^{\sigma_\infty} |p'_\infty(s)| ds + \theta \int_0^{\sigma_\infty} p_\infty(s) ds \right] > 0.$$

Consider then the following assumption on the stationary state

$$B_\infty < 1, \quad -\frac{e^{\int_0^{\sigma_\infty} p_\infty(s) ds}}{1 - B_\infty} \int_0^{\sigma_\infty} \frac{\partial}{\partial a} \int_0^s K(w, a) dw da - \int_{\sigma_\infty}^{+\infty} \frac{\partial}{\partial a} \int_0^s K(w, a) dw da < p_M, \quad (4.32)$$

which has to be understood as a smallness assumption on $\sigma$, and thus a smallness assumption on $p$. Following the approach of [161], [207] introduces the quantities

$$m(s, t) = n(s, t) - n_\infty(s), \quad M(s, t) = \int_0^s m(w, t) dw, \quad J(s, t) = \frac{\partial M}{\partial t}(s, t).$$

and applies a similar method to prove convergence in the linear case.
**Theorem 4.46** (Pakdaman, Perthame, Salort [207]). *Grant Assumptions 4.44 and 4.45. Assume (4.32),* $p(s, N) = p(s)$ and $\sigma(N) = \sigma$ are independent of $N$ and there exists $\lambda > 0$

$$ \sup_{u \le \sigma} p(u) \left( \int_{0}^{u} K(s, u) e^{\lambda(u-s)} ds - 1 \right) < \lambda $$

There exists $b > 0$ and a function $P(s) > b$ such that if

$$ \int_{0}^{+\infty} P(s) |M(s, 0)| ds + \int_{0}^{+\infty} P(s) |J(s, 0)| ds < +\infty, $$

there exists $C > 0$ and $\nu > 0$ such that

$$ \int_{0}^{+\infty} P(s) |n(s, t) - n_{\infty}(s)| ds \le C e^{-\nu t} \int_{0}^{+\infty} P(s) \left( |M(s, 0)| + |J(s, 0)| + \mathbb{1}_{s \le \sigma} |m(s, 0)| \right) ds. $$

This linear result then allows to prove a weakly-nonlinear result.

**Theorem 4.47** (Pakdaman, Perthame, Salort [207]). *Grant Assumptions 4.44 and 4.45. Assume (4.31) with $\eta$ small enough and (4.32). There exists $b > 0$ and a function $P(s) > b$ such that if*

$$ \int_{0}^{+\infty} P(s) |M(s, 0)| ds + \int_{0}^{+\infty} P(s) |J(s, 0)| ds < +\infty, $$

there exists $C > 0$ and $\nu > 0$ such that

$$ \begin{aligned} & \int_{0}^{+\infty} |n(s, t) - n_{\infty}(s)| ds + \int_{0}^{+\infty} P(s) |M(s, t)| ds + |N(t) - N_{\infty}| + |N'(t)| \\ & \le C e^{-\nu t} \left( 1 + \int_{0}^{+\infty} P(s) |M(s, 0)| ds \right). \end{aligned} $$

The method of [39] presented in Subsections 4.3.2 and 4.6.3 for the standard model (TE) also allows to prove existence and weakly-nonlinear stability of the stationary state in the space of probability measures. Besides Assumptions 4.7 and 4.8 already introduced, we need, similar to Assumption 4.45, to ensure that the mass is preserved over time, that the kernel is a positive measure and that neurons are reset in a physical way. We now require a stronger condition, which is that for any state $u$ before firing, a least some neurons are reset in a fixed neighbourhood of $s = 0$. We sum up these requirements as

**Assumption 4.48** (Fragmentation kernel for measure solutions). *For all $u \ge 0$, $K(\cdot, u)$ is a probability measure supported on $[0, u]$. There exists $\varepsilon > 0$, $\sigma > 0$ and $0 < \delta < \sigma$ such that*

$$ \forall u \ge \sigma, \quad K(\cdot, u) \ge \varepsilon \mathbb{1}_{[0, \delta]}. $$

Using this assumption, [39] extends the existence of measure solution to (4.30) with a similar fixed-point method.

**Theorem 4.49** (Cañizo, Yoldaş [39]). Let $n^0 \in \mathcal{P}(\mathbb{R}_+)$. Grant Assumptions 4.7 and 4.48 with

$$ L \le \frac{1}{4 \|n^0\|_{\text{TV}}}. \qquad (4.33) $$

There exists a unique global-in-time measure solution $(n, N)$ to (4.30). If $n^{0,1}, n^{0,2} \in \mathcal{M}_+(\mathbb{R}_+)$ satisfying (4.33), then the associated solutions $n_1$ and $n_2$ satisfy

$$ \forall t \ge 0, \quad \|n_1(t) - n_2(t)\|_{\text{TV}} \le e^{4\|p\|_{\infty}t} \|n^{0,1} - n^{0,2}\|_{\text{TV}}. $$

Then, the linear and weakly-nonlinear proofs involving Doeblin's theorem can also be extended with similar arguments, leading to
**Theorem 4.50** (Cañizo, Yoldaş [39]). *Grant Assumptions 4.7, 4.8 and 4.48 with same constant σ and L small enough. Then, there exists a unique stationary state $(n_∞, N_∞)$ to (4.30), there exist constants C ≥ 1, λ > 0 such that for all probability density n⁰, if L ≤ (1/(4||n⁰||ₜₗ)), then*

$$ \forall t \ge 0, \quad \|n(t) - n_{\infty}\|_{\mathrm{TV}} \le C e^{-\lambda t} \|n^0 - n_{\infty}\|_{\mathrm{TV}}. $$

In [207], the strongly-nonlinear case is studied numerically with the choice $K(s, u) = \delta_{s=\frac{u}{2}}$, where a neuron is reset to half the time elapsed before its discharge. They observe that this fragmentation kernel has a smoothing effects on the periodic solutions constructed and observed for the standard model (which is $K(s, u) = \delta_{s=0}$): there is a unique continuous periodic firing rate attracting all solutions. Periodic solutions with bursting pattern can be observed for some parameters. Note that the choice $K(s, u) = \delta_{s=\frac{u}{2}}$ is incompatible with the kernel assumption 4.48. To the best of our knowledge, the numerical exploration of periodic solutions in the framework of [39] has not been performed yet.

### 4.9.3 Adaptation and fatigue modelled as an additional variable

In [110], another point of view is proposed for the modelling of adaptation and fatigue of neurons. Instead of resetting the neurons through a fragmentation kernel, they keep the physical meaning of the time elapsed since the last discharge s and add a separate leaky memory variable m ≥ 0. It allows to take into account with more flexibility phenomena like short-term plasticity and spike-frequency adaptation [152, 254, 262].

In [239], the following model is obtained rigorously through a mean-field limit of age and memory dependent Hawkes processes.

$$ \left\{ \begin{array}{l} \frac{\partial n}{\partial t} + \frac{\partial n}{\partial s} + \frac{\partial}{\partial m} (-\lambda m n) + p(s, m, X(t)) n = 0, \\ n(0, m, t) = \mathbb{1}_{m \ge \gamma(0)} |(\gamma^{-1})'(m)| \int_{0}^{+\infty} p(s, \gamma^{-1}(m), X(t)) n(s, \gamma^{-1}(m), t) ds, \\ X(t) = J \int_{0}^{t} \int_{0}^{+\infty} \int_{0}^{+\infty} h(t - \tau, s, m) p(s, m, X(\tau)) n(s, m, \tau) dm ds d\tau, \\ n(s, m, 0) = n^0(s, m) \ge 0, \quad \int_{0}^{+\infty} \int_{0}^{+\infty} n^0(s, m) dm ds = 1. \end{array} \right. \quad (4.34) $$

As in the previous Time Elapsed models we presented, p is the firing function, $X(t)$ the network activity at time t and J the connectivity strength. The network activity is computing according to a memory kernel h and the history of all the neurons spiking between time 0 and t; this variant of the Time Elapsed model is then delayed. If $h(\tau, s, m) = h(\tau)$, the network activity simplifies down to the familiar

$$ X(t) = \int_{0}^{t} h(t - \tau) \mathcal{N}(\tau) d\tau, \quad \mathcal{N}(t) = \int_{0}^{+\infty} \int_{0}^{+\infty} p(s, m, X(\tau)) n(s, m, \tau) dm ds. $$

However, more general biologically relevant forms are possible for h, for example:

* Spike frequency adaptation: Take $\gamma(m) = m + \Gamma$ for some constant $\Gamma > 0$, $h(t, s, m) = \tilde{h}(t)$ and $p(s, m, t) = \tilde{p}(\eta(s) - m + X)$ with $\tilde{p}$ non-decreasing and $\eta$ a bounded refractory Kernel such that $\eta(+\infty) = 0$. This can represent the fact that neurons spiking at high frequency will undergo spike frequency adaptation and be unable to spike for some time [5]. With these parameters, the model exhibits bursting at the population level as can be seen on Figure 12.

* Short term synaptic depression: when neurons fire a lot in recent past, it can induce short term depression at the level of the synapse [263], it can be modelled by restricting $m \in (0, 1)$ and choosing $h(t, s, m) = \tilde{h}(t)(1-m)$, $\gamma(m) = 1 - v + vm$ for some fixed $v \in (0, 1)$. Neurons with a leaky memory m closed to 1 will have very little impact on the network activity $X(t)$.
<figure><img src="image_14.png" /><figcaption>Figure 12: Self sustained bursts for large $J$ and relaxation to stationary state for small $J$, seen as the network activity $X(t)$ in model (4.34) or the stochastic spikes in the underlying particle system [239]. Fonte Sánchez, Schmutz [110].</figcaption></figure>

The resetting of neurons depends on an increasing $C^1$-diffeomorphism $\gamma: \mathbb{R}_+ \to (0, +\infty)$ called the jump mapping. A neuron firing with state $(s, m)$ is reset to $(0, \gamma(m))$, and then the boundary condition at $(0, m)$ is written in term of $\gamma^{-1}$.

It can be formally checked and rigorously proved that system (4.34) preserves the mass. Define the weight $\omega(s, m) = m + 1$, consider the weighted space $L_+^1(\omega)$ and the set of hypotheses

**Assumption 4.51** (For well-posedness). There exists $L$ such that for all $(s_1, m_1, X_1), (s_2, m_2, X_2) \in \mathbb{R} + \times (0, +\infty) \times \mathbb{R}_+$,

$$|p(s_1, m_1, X_1) - p(s_2, m_2, X_2)| \le L(|s_1 - s_2| + |m_1 + m_2| + |X_1 - X_2|).$$

There exists a positive bounded function $\Gamma$ such that for all $m > 0$, $\gamma(m) = m + \Gamma(m)$.

Then, similar to the strategy in Subsection 4.3.1 two applications of the Banach fixed point theorem to an appropriate functional yields existence of solutions in the weighted space $L_+^1(\omega)$.

**Theorem 4.52** (Fonte Sánchez, Schmutz [110]). Grant Assumption 4.51. There exists constants $\alpha$ and $b \in \mathbb{R}$ such that for any probability density $n^0 \in L_+^1(\omega)$, there exists a unique weak-solution $n \in C(\mathbb{R}_+, L_+^1(\omega))$ to (4.34) and

$$\forall t \ge 0, \quad \|n(\cdot, t)\|_{L_+^1(\omega)} \le \|n^0\|_{L_+^1(\omega)} e^{-\alpha t} + \frac{b}{\alpha}(1 - e^{-\alpha t}).$$

Note that well-posedness is also proved in $L^1$ not weighted in [110], but the weighted version is required for the Doeblin-Harris method.

**Assumption 4.53** (On firing and resetting). There exists $\sigma > 0, p_0 > 0$ such that if $s > \sigma, p(s, \cdot, \cdot) \ge p_0$. There exists $C_\gamma \in (0, 1]$ such that $C_\gamma \le \gamma' \le 1$. The function $\bar{h}: (s, m) \mapsto \int_0^{+\infty} h(t, s, m)dt$ is bounded.

**Theorem 4.54** (Fonte Sánchez, Schmutz [110]). Grant Assumptions 4.51 and 4.53. There exist at least one stationary state $(n_\infty, X_\infty)$ to (4.34). If $|J|$ is small enough, the stationary state is unique.
Then, similar to [39] but adding the weighted framework and using Harris' theorem instead of Doeblin theorem, [110] proves that when $J$ is small enough, there is global convergence towards the unique stationary state.

**Theorem 4.55** (Fonte Sánchez, Schmutz, [110]). *Grant Assumptions 4.51 and 4.53. Assume there exists $C_1, c_2 > 0$ such that $h(t, \cdot, \cdot) \le C_1 e^{-c_2t}$. For any $A > 0$, if $|J|$ is small enough with respect to $A$, then there exists $C \ge 1$, $\nu > 0$ such that for all initial condition satisfying $\|n^0\|_{L_+^1(\omega)} \le A$, the solution to (4.34) satisfies*

$$ \|n(\cdot, t) - n_{\infty}\|_{L_{+}^{1}(\omega)} + |X(t) - X_{\infty}| \le C e^{-\nu t} \left( 1 + \|n^0 - n_{\infty}\|_{L_{+}^{1}(\omega)} \right). $$

A similar result can be proved when instead of $\gamma(m) = m + \Gamma(m)$ we assume that $\gamma$ is bounded.

### 4.9.4 Time elapsed since the last two discharges

Another direction for extending the standard Time Elapsed model is to incorporate another variable $a$ representing, for a neuron with time $s$ elapsed since the last discharge, the time elapsed since the penultimate discharge. Such a model is proposed and explored in [256]. Multiple renewal equations had previously been explored in [114]. Define

$$ \Omega = \{(s, a) \in \mathbb{R}^2 | 0 \le s \le a\}. $$

The multiple time elapsed model reads

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial n}{\partial t}(s, a, t) + \frac{\partial n}{\partial s}(s, a, t) + \frac{\partial n}{\partial a}(s, a, t) + p(s, a, X(t)) n(s, a, t) = 0, \\[1.5em] \displaystyle N(a, t) = \int_a^{+\infty} p(s, a, X(t)) n(a, u, t) du, \quad X(t) = \int_0^{+\infty} N(a, t) da, \\[1.5em] \displaystyle n(0, a, t) = N(a, t), \quad n(s, a, 0) = n^0(s, a) \ge 0, \quad \int_0^{+\infty} \int_s^{+\infty} n^0(s, a) da ds = 1. \end{array} \right. \qquad (4.35) $$

As for the standard Time Elapsed model, we have conservation of mass for weak solutions

$$ \int_{0}^{+\infty} \int_{s}^{+\infty} n^{0}(s, a) \mathrm{d}a \mathrm{d}s = \int_{0}^{+\infty} \int_{s}^{+\infty} n(s, a, t) \mathrm{d}a \mathrm{d}s = 1. $$

Note that if $p(s, a, X) = p(s, X)$ is independent of $a$, then the quantity $m(s, t) = \int_{s}^{+\infty} n(s, a, t) da$ solves (TE). The model (4.35) can hence exhibit all the behaviours of the standard Time Elapsed model, like periodic solutions and the existence of multiple stationary states.

By studying the linear problem and then applying a fixed point argument, [256] proves the well-posedness of this model in the weakly-nonlinear regime.

**Theorem 4.56** (Perthame, Salort, Torres, [256]). *Assume that $n^0 \in L^1(\Omega)$ is a probability density and that $p \in W^{1,\infty}(\Omega \times \mathbb{R})$. If $\|\partial_X p\|_{\infty} < 1$, there exists a unique solution of (4.35) satisfying*

$$ n \in C(\mathbb{R}_+, L^1(\Omega)) \cap L^\infty(\Omega \times \mathbb{R}_+), \quad N \in C(\mathbb{R}_+, L^1(\mathbb{R}_+)) \cap L^\infty(\mathbb{R}_+^2), \quad X \in C(\mathbb{R}_+) \cap L^\infty(\mathbb{R}_+). $$

Then, for the linear system $p(s, a, X) = p(s, a)$, they prove the following entropy dissipation: for all convex function $H \ge 0$ with $H(0) = 0$, for all solution $n$ to the linear system,

$$ \frac{d}{dt} \int_{0}^{+\infty} \int_{s}^{+\infty} H \left( \frac{n(s, a, t)}{n_{\infty}(s, a)} \right) n_{\infty}(s, a) \mathrm{d}a \mathrm{d}s = -D_{H}[n](t) \le 0, $$

where $(n_\infty, N_\infty, X_\infty)$ is the unique linear stationary state and

$$ D_H[n](t) = \int_{0}^{+\infty} \int_{s}^{+\infty} p(s, a) H \left( \frac{n(s, a, t)}{n_{\infty}(s, a)} \right) n_{\infty}(s, a) \mathrm{d}a \mathrm{d}s - \int_{s}^{+\infty} H \left( \frac{N(a, t)}{N_{\infty}(a)} \right) N_{\infty}(a) \mathrm{d}a. $$
However, a requirement for obtaining this entropy dissipation is that the stationary state satisfies $n_\infty > 0$, $N_\infty > 0$ which is not true in general (it can be proved to be false for example with $p(s, a, X) = 1_{s>1}$). Moreover, it is not easy to prove an appropriate Poincaré-like inequality to bound $D_H[n]$. For these reasons, [256] turns to the Doeblin theory used in [39] and described in Subsection 4.6.3. In the linear case $p(s, a, X) = \bar{p}(s, a)$, they prove that there exists a unique stationary state $(n_\infty, N_\infty, X_\infty)$ of (4.35) and the solution n satisfies

$$ \|n(t) - n_{\infty}\|_{L^1} \le \frac{1}{1-\beta} e^{-\lambda t} \|n^0 - n_{\infty}\|_{L^1}, $$

with $\beta = \frac{1}{2}p_0^2\sigma^2e^{-3p_\infty\sigma}$ and $\lambda = -\frac{1}{3\sigma}\log(1-\beta) > 0$. In the weak-connectivity regime $\partial_X p$ small, there is also a unique stationary state for the nonlinear problem. The linear result can be extended to the weakly-nonlinear case.

**Theorem 4.57** (Perthame, Salort, Torres, [256]). Assume that $p$ is smooth, $\|\partial_X p\|_\infty$ is small enough and that there exists $p_0 > 0$, $p_\infty > 0$, $\sigma > 0$ such that $p_0 1_{a>s\ge\sigma} \le p \le p_\infty$. Then, there exists a unique stationary state $(n_\infty, N_\infty, X_\infty)$ of (4.35) and constants $C \ge 1$, $\lambda > 0$ such that if $n^0 \in L^1(\Omega)$ is a probability density, the solution n to (4.35) satisfies

$$ \forall t \ge 0, \quad \|n(t) - n_{\infty}\|_{L^1} \le C e^{-\lambda t} \|n^0 - n_{\infty}\|_{L^1}. $$

Numerical exploration in [256] shows that (4.35) also displays solutions with jump discontinuities in the strongly-nonlinear regime. It was observed that convergence to a periodic solution in the case $p(s, a, X) = p(s, X)$ is replaced by steady asymptotics when the dependency in $a$ is added. Complex dynamics elicited by the $a$ variable has yet to be observed and studied. One can hope to find a reduction to a 2-delays equation similar to the findings of [255] for (TE).

### 4.9.5 Convergence towards an infinite renewal model

The idea of taking into account both the time since the last and the penultimate discharges is pushed further in [88]: they consider for each neuron $N$ variables $s_1, \dots, s_N$ representing the time elapsed since the last $N$ discharges. Then, the firing function can depend upon more information about the past history of the neurons. With the aim of sending $N$ towards $+\infty$, a particular structure is proposed, which gives less and less weight to later discharges

$$ p_N(s_1, \dots, s_N) = \sum_{i=1}^{N} \phi_i(s_1, \dots, s_N), \quad \phi_i > 0, \quad \sum_{i=1}^{+\infty} \|\phi_i\|_{\infty} \le a_+ < +\infty. \quad (4.36) $$

Note that, for the sake of mathematical tractability, $p_N$ does not depend on the network activity $X(t)$. The model for $N$ time elapsed variables can be written as

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial n_N}{\partial t} + \sum_{i=1}^{N} \frac{\partial n_N}{\partial s_i} + p_N(s_1, \dots, s_N)n_N = 0. \\[1.5ex] \displaystyle n_N(0, s_2, \dots, s_N, t) = \int_0^{+\infty} p_N(s_2, \dots, s_N, u)n_N(s_2, \dots, s_N, u, t)du, \\[1.5ex] \displaystyle n_N(s_1, \dots, s_N, 0) = n^0(s_1, \dots, s_N) \ge 0, \quad \int_{C_N} n_N(s_1, \dots, s_N)ds_1 \dots ds_N = 1. \end{array} \right. \quad (4.37) $$

where we define the domains

$$ C_N = \{0 \le s_1 \le \dots \le s_N\}, \quad C_\infty = \{0 \le s_1 \le \dots \le s_N \le \dots\}. $$

As we have mentioned in Subsection 4.8, there is a connection between (4.37) and Wold processes. However, there is currently, to our knowledge, no framework for Wold processes with an infinite number of variables.
In order to improve readability, especially when $N \to +\infty$, let us use $[s]_N$ to denote $(s_1, \dots, s_N) \in C_N$ and $[s]_\infty = (s_1, \dots, s_N, \dots) \in C_\infty$. Consider also the shift operator $\tau$ defined by

$$ \tau(s_1, s_2, \dots, s_N) = (0, s_1, s_2, \dots, s_{N-1}), \quad \tau(s_1, s_2, \dots, s_N, \dots) = (0, s_1, s_2, \dots, s_N, \dots). $$

A weak-solution to (4.37) can be defined with the shift operator: $n \in C(\mathbb{R}_+, L^1(C_N))$ if for all $\psi \in C_b^1([0, T] \times C_N)$,

$$ \int_0^T \int_{C_N} n_N([s]_N, t) \left[ -\frac{\partial \psi}{\partial t} - \nabla \cdot \psi + (\psi([s]_N, t) - \psi(\tau([s]_N, t)))p_N([s]_N) \right] d[s]_N dt \\ = - \int_{C_N} \psi(\tau([s]_N, T)) n_N([s]_N, T) d[s]_N + \int_{C_N} \psi(\tau([s]_N, 0)) n_N([s]_N, 0) d[s]_N $$

As for the model with $N = 2$ in [256] and other Time Elapsed models [39, 110], the Doeblin-Harris method allows to prove convergence towards a unique stationary state.

**Theorem 4.58** (Dou, Perthame, Qi, Salort, Zhou [88]). Assume $0 < a_- \le p_N \le a_+$, with constants $a_- > 0, a_+ > 0$. Then, for all appropriate initial condition, for all $[s]_N \in C_N$, $n_N([s]_N, t) \ge a_- e^{-s_N a_+}$. There exist a unique stationary state $n_{N,\infty}$ and constants $c_N > 0, \lambda_N > 0$ such that

$$ \forall t \ge 0, \quad \|n_N(t) - n_{N,\infty}\|_{L^1(C_N)} \le c_n e^{-\lambda_N t} \|n_N^0 - n_{N,\infty}\|_{L^1(C_N)}. $$

Moreover, the following tightness estimate holds for the stationary state:

$$ \int_{C_N} \sigma_N([s]_N) n_{N,\infty} d[s]_N \le \frac{2}{a_-}, \quad \sigma_N([s]_N) := \sum_{i=1}^N \frac{s_i}{2^i}. $$

The constants $c_N, \lambda_N$ are explicitly constructed in [88]. Unfortunately, they converge towards 0 when $N$ goes to $+\infty$. Thus, it is not possible to use such decay estimates for the limit $N \to +\infty$. Another difficulty is to write a rigorous formulation for the limit of (4.37). Formally, it would be an infinite renewal equation with variables $[s]_\infty \in C_\infty$ in the form

$$ \frac{\partial n_\infty}{\partial t} + \sum_{i=1}^{+\infty} \frac{\partial n_\infty}{\partial s_i} + p_\infty([s]_\infty) n_\infty = 0. $$

But because there is no "last" variable, it is not easy to define a solution, or simply to prescribe the boundary condition at $s_1 = 0$. Therefore, [88] proposes to study the $K$-s of the solution to this limiting infinite renewal equation, similar to BBGKY hierarchies in classical kinetic theory [57].

First, let us write the marginals for the solution $n_N$ to the problem with a finite number $N$ of variables $s_1, \dots, s_N$. For $K \le N$, denote $n_N^{(K)}$ the first $K$ marginals of $n_N$:

$$ n_N^{(K)}([s]_K, t) = \int_{s_K \le s_{K+1} \le \dots \le s_N} n_N([s]_N, t) ds_{K+1} \dots ds_N. $$

Note that $n_N^{(K)} = \int_{s_K}^{+\infty} n_N^{(K+1)} ds_{K+1}$.

**Definition 4.59** (Consistent sequence). A sequence of measures $(\mu_N)_{N \ge 1}$ is said to be consistent if for all $K < N$, the marginals satisfy $\mu_K = \mu_N^{(K)}$. A measure $\mu_\infty$ on $C_\infty$ is consistent if for all $K \ge 1$, $\mu_\infty^{(K)} = \mu_\infty^{(K+1)} ds_{K+1}$.

Given the coupling term

$$ E_N^{(K)}([s]_K, t) = \sum_{i=K+1}^{N} \int_{s_{K+1}}^{+\infty} \dots \int_{s_N=0}^{+\infty} \phi_i([s]_i) n_N^{(i)}([s]_i, t) d[s]_{K+1,i}, $$
the marginals $n_N^{(K)}$ satisfy

$$ \frac{\partial n_N^{(K)}}{\partial t} + \sum_{i=1}^{K} \frac{\partial n_N^{(K)}}{\partial s_i} + p_K([s]_K) n_N^{(K)} + E_N^{(K)} = 0, $$

$$ n_N^{(K)}(0, [s]_{2,K}, t) = \int_0^{+\infty} \left[ p_K([s]_{2,K}, u) n_N^{(K)}([s]_{2,K}, u, t) + E_N^{(K)}([s]_{2,K}, u, t) \right] du $$

It leads to the infinitely coupled hierarchy system

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial n_{\infty}^{(K)}}{\partial t} + \sum_{i=1}^{K} \frac{\partial n_{\infty}^{(K)}}{\partial s_i} + p_K([s]_K) n_{\infty}^{(K)} + E_{\infty}^{(K)} = 0. \\[1.5em] \displaystyle n_{\infty}^{(K)}(0, s_2, \dots, s_K, t) = \int_{0}^{+\infty} \left[ p_K([s]_{2,K}, u) n_{\infty}^{(K)}([s]_{2,K}, u, t) + E_{\infty}^{(K)}([s]_{2,K}, u, t) \right] du, \\[1.5em] \displaystyle n_{\infty}^{(K)}([s]_K, 0) = n_{\infty}^{0,(K)}([s]_K). \end{array} \right. \quad (4.38) $$

The coupling system is obtained by passing to the limit in $E_N^{(K)}$, which yields

$$ E_{\infty}^{(K)}([s]_K, t) = \sum_{i=K+1}^{+\infty} \int_{s_{K+1}=0}^{+\infty} \dots \int_{s_i=0}^{+\infty} \phi_i([s]_i) n_{\infty}^{(i)}([s]_i, t) d[s]_{K+1,i}. $$

**Theorem 4.60** (Dou, Perthame, Qi, Salort, Zhou [88]). Assume (4.36). Consider a sequence of initial conditions $n_N^0$ of mass 1 on $C_N$ such that there is a consistent family $(n_{\infty}^{0,(K)})_K$ such that

$$ \lim_{N \to +\infty} \sup_K \left\| n_N^{0,(K)} - n_{\infty}^{0,(K)} \right\| = 0. $$

* For all $T > 0$ and $K \in \mathbb{N}$, $(n_N^{(K)})_{N \in \mathbb{N}}$ is a Cauchy sequence in $C([0, T], L^1(C_K))$ and thus has a consistent limit $n_{\infty}^{(K)} \in C([0, T], L^1(C_K))$.

* $E_N^{(K)}$ converges towards $E_{\infty}^{(K)}$ in $C([0, T], L^1(C_K))$ as $N \to +\infty$.

* $(n_{\infty}^{(K)})_{K \in \mathbb{N}}$ is the unique consistent weak solution of the hierarchy (4.38).

Such a solution to the hierarchy (4.38) is obtained as a local-in-time limit on $[0, T]$ of the models with finite number of variables. Under additional smallness assumptions, [88] also proves a uniform-in-time, using the Doeblin method. It allows to approximate the stationary solution of the infinite renewal equation by the stationary state $n_{N,\infty}$ for large $N$:

$$ \lim_{N \to +\infty} \left\| n_{N,\infty} - n_{\infty,\infty}^{(N)} \right\|_{L^1(C_N)} = 0. $$

Under smallness assumption on the functions $\phi_i$ in (4.36), they obtain uniqueness of the stationary state $n_{\infty,\infty}$ of the hierarchy.

Last, [88] proposes a further interpretation of the limiting infinite renewal equation by introducing solutions as measures on $C_{\infty}$ which are weak limits of the solutions for finite $N$. The time-dependent infinite-dimensional measure $n_{\infty}(t) \in \mathcal{P}(C_{\infty})$ can be built using the Kolmogorov extension theorem. This notion of solution is proved to be equivalent to the hierarchy (4.38). It allows in particular to use the Monge-Kantorovich distance to write and prove a long-term convergence result for the infinite system.
### 4.9.6 Spatial extension and evolution of the connectivity

The standard Time Elapsed model does not take into account spatial variability of the connectivity $J$ nor evolution of the connectivity in time. Following previous spatial extensions in Fitzhugh-Nagumo models [73, 74] (see also Section 8) and foreshadowing spatial extension in Fokker-Planck models (see Section 2), [257] proposes a spatial extension of the Time Elapsed model where neurons at locations $x \in \Omega$ and $y \in \Omega$ interact through the time-dependent connectivity $w(x, y, t)$.

$$ egin{cases} 
\frac{\partial n}{\partial t}(x, s, t) + \frac{\partial n}{\partial s}(x, s, t) + p(s, X(x, t)) n(x, s, t) = 0, \\ 
N(x, t) = \int_{0}^{+\infty} p(s, X(t)) n(x, s, t) ds, \\ 
X(t) = \int_{\Omega} w(x, y, t) N(y, t) dy + I(x, t), \\ 
\frac{\partial w}{\partial t}(x, y, t) = -w(x, y, t) + \gamma G[N(x, t), N(y, t)], \quad w(x, y, 0) = w^0(x, y). \\ 
n(0, t) = N(t), \quad n(s, 0) = n^0(s) \ge 0, \quad \int_{0}^{+\infty} n^0(s) ds = 1. 
\end{cases} 
\quad (4.39) $$

Neurons at location $x$ are stimulated by a network activity $X(x, t)$ which is determined as an integral on the spatial domain $\Omega$ of the firing rate $N(y, t)$ with weight given by the connectivity function $w(x, y, t)$, plus an external input $I(x, t)$ which is a given function. The connectivity kernel $w$ evolves according to its own equation, depending on the firing rates at locations $x$ and $y$. The smooth function $G: \mathbb{R}^2 \to \mathbb{R}$ is the learning rule for the connectivity, and $\gamma > 0$ is the learning rate. We can take for example a representation of Hebbian learning like $G[N_1, N_2] = N_1 N_2$. The following assumption is made for simplicity

$$ \|G\|_{\infty} + \|\nabla G\|_{\infty} \le 1. \quad (4.40) $$

There is no loss of generality as $\gamma$ can be scaled to enforce it. With a fixed point argument, [257] proves that the spatially extended problem is well-posed in the weak-interaction regime for small enough initial condition.

**Theorem 4.61** (Salort, Torres, [257]). Assume (4.40), $p \in W^{1,+\infty}(\mathbb{R}_+ \times \mathbb{R})$ and that there exists $p_0, p_\infty, \bar{\sigma}$ such that $p_0 1_{s>\bar{\sigma}} \le p \le p_\infty$. For all initial condition $n^0 \in C_b(\Omega, L^1(\mathbb{R}_+))$, $w^0 \in C_b(\Omega \times \Omega)$ with $n^0$ of total mass 1, if

$$ |\Omega| \left\| \int_{0}^{+\infty} n^0(\cdot, s) ds \right\|_{\infty} \left\| \frac{\partial p}{\partial X} \right\|_{\infty} \max (\|w^0\|_{\infty}, \gamma) < 1, $$

then the system (4.39) has a unique solution such that $n \in C_b(\Omega \times \mathbb{R}_+, L^1(\mathbb{R}_+))$ is non-negative, $N \in C_b(\Omega \times \mathbb{R}_+)$ and $w \in C_b(\Omega \times \Omega \times \mathbb{R}_+)$.

The result can also be proved for a larger class of firing functions $p$. A fixed point argument also allows [257] to prove that, given an initial spatial distribution of neurons $g(x) = \int_{0}^{+\infty} n^0(\cdot, s) ds$, there exists a unique continuous stationary state in the weak-connectivity and small learning rate regime (it can be related to the homogeneous brain assumption for the grid cells model in Section 2). Note that since the connectivity kernel is evolving in time, a kernel converging towards a small function in $L^\infty$ norm means that the nonlinearity through $X(x, t)$ will be small too, hence the uniqueness of the $n_\infty$ stationary density.

Following [179, 178, 208, 205, 257] proves the convergence towards the unique stationary state in the weak-connectivity regime via a relative entropy method.

**Theorem 4.62** (Salort, Torres, [257]). Assume (4.40), $p \in W^{1,+\infty}(\mathbb{R}_+ \times \mathbb{R})$ and that there exists $p_0, p_\infty, \bar{\sigma}$ such that $p_0 1_{s>\bar{\sigma}} \le p \le p_\infty$. For all initial condition $n^0 \in C_b(\Omega, L^1(\mathbb{R}_+))$, $w^0 \in C_b(\Omega \times \Omega)$
with $n^0$ of total mass 1, if $\gamma$, $\|\partial_X p\|_\infty$ are small enough with respect to $n^0$, $w^0$, then there exists $C \ge 1$, $\nu > 0$ such that

$$ \forall t \ge 0, \quad \|n(\cdot, \cdot, t) - n_\infty\|_{L^1} + \|w(\cdot, \cdot, t) - w_\infty\|_{L^1} \le Ce^{-\nu t} (\|n^0 - n_\infty\|_{L^\infty} + \|w^0 - w_\infty\|_{L^1}), $$
and $X \to X_\infty, N \to N_\infty$ both exponentially fast.

Then, [257] builds on the Doeblin theory in [39] and improves the relative entropy result, proving under similar conditions that

$$ \|n(\cdot, \cdot, t) - n_\infty\|_{L_x^\infty L_s^1} + \|w(\cdot, \cdot, t) - w_\infty\|_\infty \le Ce^{-\nu t} (\|n^0 - n_\infty\|_{L_x^\infty L_s^1} + \|w^0 - w_\infty\|_\infty). $$

The article [257] also provides theoretical insights into the effect of large inputs $I(x, t)$ and slow learning dynamics compared to the evolution of network activity. For the later, they consider a slow timescale represented by a parameter $\varepsilon > 0$ and they change the first equation of the system with

$$ \varepsilon\partial_t n + \partial_s n + p(s, X(x, t))n = 0. $$

Then, [257] rigorously proves that when $\max(\|w^0\|, \gamma)$ is small enough, the solution $n_\varepsilon$ converges in $L^1$ towards a solution to the limiting problem

$$ \partial_s n + p(s, X(x, t))n = 0, $$

and that this limiting problem is well-posed in the weakly-nonlinear regime. They also prove asymptotic convergence towards a unique stationary state for the weakly-nonlinear limiting problem.

The theoretical analysis of the strongly-nonlinear regime is yet unexplored for this model, but [257] provides numerical simulations showing that a non-constant external input $I(x)$ can lead to asymptotic convergence towards spatial patterns. When $I$ is constant, convergence towards spatially uniform stationary states are observed even for strong nonlinearities. The numerical experiments are performed with a Hebbian learning rule as explicitly above.

### 4.9.7 Master equation and link to the NNLIF model

A series of papers [99, 98, 100] proposed a connection between the NNLIF model described throughout Section 1 and the Time Elapsed model.

First, in the linear case, consider the auxiliary equation

$$ \frac{\partial q}{\partial s} + \frac{\partial}{\partial v}[-vq] - a \frac{\partial q}{\partial v^2} = 0, \quad q(V_F, s) = q(-\infty, s) = 0, \quad q(v, 0) = \delta_{v=V_R}, \quad (4.41) $$

posed for voltage $v \in (-\infty, V_F]$ and time since the last discharge $s \in [0, +\infty)$. The solution $q$ represents, when a neuron fires and is reset at voltage $v = V_R$, how its probability distribution in voltage evolves until the next spike. In the limit $s \to +\infty$, the neuron has more and more probability to have fired a second time, as represented by the decay of $q(v, s)$ towards 0 in age. The quantity

$$ S(s) = \int_{-\infty}^{V_F} q(v, s)dv $$

is the survival probability at age $s$ since the last discharge. It has been used in [143] for parameter estimation of the linear NNLIF parameters.

**Theorem 4.63** (Dumont, Henri, Tarniceriu, [99, 98]). Let $\rho$ a solution to (NNLIF) with $b = 0$ from an initial condition $\rho^0$, $n$ a solution to (TE) with $\mu_X = 0$ from an initial condition $n^0$ and $q$ a solution of (4.41). If

$$ \forall v \in (-\infty, V_F], \quad \rho^0(v) = \int_0^{+\infty} \frac{q(s, v)}{\int_{-\infty}^{V_F} q(v, a)da} n^0(s)ds, $$
then

$$ \forall (v, s, t) \in (-\infty, V_F] \times \mathbb{R}_+ \times \mathbb{R}_+, \quad \rho(v, t) = \int_0^{+\infty} \frac{q(s, v)}{\int_{-\infty}^{V_F} q(v, a) da} n(s, t) ds. $$

Moreover, for all $t$ such that $0 < s < t$,

$$ n(s, t) = \frac{\partial}{\partial s} \int_{-\infty}^{V_F} q(v, s) \rho(v, t-s) dv. $$

and for $s \ge t$,

$$ \int_{t}^{+\infty} n(s, t) ds = \int_{-\infty}^{V_F} q(v, s) \rho^0(v) dv. $$

This correspondence can be interpreted as a Bayes rule: the kernel in the integral transform represents the probability density for a neuron at electric potential $v$ given that it did not spike at time $s$ since the last discharge.

In [100], a unifying approach is proposed in the nonlinear case $b \ne 0$, $\mu_X \ne 0$. The idea is to consider the stochastic system

$$ dV_t = (-v + bN(t))dt + \sigma dB_t, \quad dS_t = dt, $$

with $N(t)$ the firing rate of the network. When a neuron has voltage $V_{t^-} = V_F$, it is reset at $V_{t^+} = V_R$, $S_{t^+} = 0$. The two-dimensional kinetic Fokker-Planck equation for the density, that we also call the master equation is

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial M}{\partial t} + \frac{\partial M}{\partial s} + \frac{\partial}{\partial v} [(-v + bN(t))M] - a \frac{\partial M}{\partial v^2} = 0, \\[1.5ex] \displaystyle N(t) = -a \int_0^{+\infty} \frac{\partial M}{\partial v}(V_F, s, t) ds, \\[1.5ex] \displaystyle M(V_F, s, t) = M(-\infty, s, t) = M(v, +\infty, t) = 0, \quad M(v, 0, t) = N(t) \delta_{v=V_R}, \\[1.5ex] \displaystyle M(v, s, 0) = M^0(v, s) \ge 0, \quad \int_0^{+\infty} \int_{-\infty}^{V_F} M^0(v, s) dv ds = 1. \end{array} \right. $$

where $a = \frac{\sigma^2}{2}$. Integrating the density $M$ with respect to age $s$ yields a solution $\rho$ to (NNLIF):

$$ \rho(v, t) = \int_{0}^{+\infty} M(v, s, t) ds. $$

A solution $n$ to (TE) with $p(s, t) = \frac{-a \partial_v q(V_F, s, t)}{\int_{-\infty}^{V_F} q(v, s, t) dv}$ can also be recovered by computing

$$ M(v, s, t) = \frac{q(v, s, t)}{\int_{-\infty}^{V_F} q(v, s, t)} n(s, t), $$

where $q$ solves the first passage time problem

$$ \frac{\partial q}{\partial t} + \frac{\partial q}{\partial s} + (-v + bN(t)) \frac{\partial q}{\partial v} - a \frac{\partial q}{\partial v^2} = 0, \quad q(v, 0, t) = \delta_{v=V_R}, \quad q(V_F, s, t) = 0. $$

Both these fact are rigorously proved in the linear case $b = \mu_X = 0$ in [100]. Formal integration shows they hold in the nonlinear case $b \ne 0$, but then, the strong coupling through $N(t)$ between the equations for $q$ and $M$ has to be taken into account. Also, the reduction to (TE) is not standard as the firing function does not depends straightforwardly on the network activity.
# 5 Jump based Leaky Integrate & Fire model

An issue with many of the models we have presented so far is that they rely on a mean-field assumption of an infinite number of all-to-all connected neurons. For example, in the NNLIF model in Section 1, the K neurons in the particle system make jumps of size $\frac{b}{K}$ upon receiving a spike. In some contexts, it can be useful to keep jumps at a non-vanishing size during the mean-field limit. The jump-based Leaky Integrate & Fire model is an example of this idea.

We will discuss how it has been obtained by physicists, then propose a mathematical setting, discuss well-posedness theory in both function and measure spaces, and talk about the occurrence of finite time blow-up in strongly excitatory networks. Last, we will explain how this model can be connected back to the NNLIF model through heuristic arguments in the vanishing jumps limit.

## 5.1 Macroscopic description of deterministic IF neurons

Like for the derivation of the model (NNLIF) in Section 1, let us start from a (now renormalised) Laplacian Integrate and Fire model

$$ \frac{dV}{dt}(t) = -V(t) + h \sum_{j \in \mathbb{N}} \delta_{t=t_j}, \quad V(t) \ge V_F \implies V(t^+) = V_R, $$

where a representative neuron of voltage $V(t)$ receives spikes from other neurons at times $t_j$. The quantity $h$ represents an average jump in voltage upon receiving a spike.

Instead of assuming Poisson discharges and replacing the incoming spikes by an Ornstein-Uhlenbeck process of same mean and variance like in the NNLIF derivation, another popular choice [202, 151, 122] to obtain a macroscopic description consistent with behaviours in particle simulations [204] is to formally consider a large number of $N$ replicas of the population of neurons and to take an average over all the replicas. If for the $k^{th}$ replica, there are $n_k(v, t)dv$ neurons in a voltage element $[v, v+dv)$, then the macroscopic probability density $\rho$ is constructed as

$$ n(v, t) = \lim_{N \to +\infty} \frac{1}{N} \sum_{N \to +\infty} n_k(v, t), \quad \rho(v, t) = \frac{n(v, t)}{\int_{-\infty}^{V_F} n(w, t) dw}. $$

The initial conditions for the replicas are randomly drawn but there is no source of noise in the time evolution of each of them, leading to a hyperbolic equation with no intrinsic noise source. Note an important distinction between the model (NNLIF) and this approach: because we do not take the mean-field limit for a number $N$ of all-to-all connected neurons going to $+\infty$, we do not need to introduce a mean-field scaling $\frac{h}{N}$ in the process. The macroscopic model obtained with the replica ansatz then has finite-sized, non-vanishing jumps when excitatory spikes are received at the macroscopic level, leading to a jump-based structure in the PDE.

After constructing each quantity involved through averages over the replicas, [202] obtains the evolution equation

$$ \frac{\partial \rho}{\partial t}(v, t) + \frac{\partial}{\partial v}[-v\rho(v, t)] + \sigma(t)[\rho(v, t) - \rho(v - h, t)] = r(t)\delta_{v_R}, $$

where the function $\sigma(t)$ is the deterministic flux of excitatory spikes reaching the population and $r(t)$ is the quantity of neurons spiking and being reset at time $t$. Then, $\sigma(t)$ is itself a sum of the external stimulation $\sigma_0 \ge 0$ from outside the network and the internal excitatory stimulation $J_E r(t)$ coming from neurons in the network spiking, with $J_E \ge 0$ the excitatory connectivity strength.

The quantity $\rho(v, t) - \rho(v - h, t)$ represents the neurons in a window $[v - h, v]$ which are excited by the spikes. It can be better understood in divergence form

$$ \sigma(t) [\rho(v, t) - \rho(v - h, t)] = \frac{\partial}{\partial v} \left( \sigma(t) \int_{v-h}^{v} \rho(w, t) dw \right). $$
In particular, the flux $r(t)$ of neurons spiking at time $t$ corresponds to the part of the network making the $h$-sized jump from $V_F - h$ to $V_F$:

$$r(t) = \sigma(t) \int_{V_F-h}^{V_F} \rho(w,t) dw.$$

Recurrent inhibition can be taken into account by adding the term $I(t)[\rho(v+h) - \rho(v)]$, to the right-hand side of the equation, where $I(t) = J_I r(t)$ is the inhibitory current and $J_I \ge 0$ is the inhibitory connectivity parameter. A comparison between the numerical simulation of the PDE (jLIF) and a Monte-Carlo simulation of a particle system can be found in [96].

Instead of the instantaneous transmission $\sigma(t) = \sigma_0 + J_E r(t)$, $I(t) = J_I r(t)$ a delayed transmission can be considered like for the model (TE):

$$\sigma(t) = \sigma_0 + J_E \int_0^t r(t-\tau)\alpha(\tau)d\tau, \quad I(t) = J_I \int_0^t r(t-\tau)\alpha(\tau)d\tau, \quad \int_0^{+\infty} \alpha(\tau)d\tau = 1.$$

Since there is no noise (no diffusion) an the deterministic field $-v$ does not send the mass in negative voltage space, the deterministic case $J_E \ge 0$, $J_I = 0$ can be equivalently posed on $[0, V_F]$ by choosing an initial condition supported in $[0, V_F]$.

## 5.2 Mathematical setting

Following several references on the problem [96, 97, 95], we make the simplifying choice $V_F = 1$ and sum up the jump-based Leaky Integrate and Fire model as the closed system

$$\left\{
\begin{array}{l}
\displaystyle \frac{\partial \rho}{\partial t} + \frac{\partial}{\partial v} [-v\rho] + \sigma(t) [\rho(v,t) - \rho(v-h,t)] = I(t) [\rho(v+h,t) - \rho(v,t)] + r(t)\delta_{v=v_R} \\
\\
\displaystyle r(t) = \sigma(t) \int_{1-h}^{1} \rho(v,t)dv, \\
\\
\displaystyle \sigma(t) = \sigma_0(t) + J_E \int_{0}^{t} r(t-\tau)\alpha(\tau)d\mu_{\sigma}, \quad I(t) = J_I \int_{0}^{t} r(t-\tau)\alpha(\tau)d\mu_{\sigma} \\
\\
\displaystyle \rho(1,t) = 0, \quad \rho(v,0) = \rho^0(v) \ge 0, \quad \int_{-\infty}^{1} \rho^0(v,t)dv = 1.
\end{array}
\right.
\qquad (\text{jLIF})$$

The abstract measure$^6$ $\mu_\sigma$ allows to represent both instantaneous transmission with $\mu_\sigma = \delta_{\tau=0}$ and smoothly delayed transmission with $\mu_\sigma = \alpha(\tau)d\tau$. Note that when there is instantaneous transmission, $\sigma(t)$ solves the closed equation

$$\sigma(t) = \frac{\sigma_0}{1 - J_E \int_{1-h}^{1} \rho(v,t) dv}.$$

In most cases, solutions are to be understood in the weak sense: for all suitable $C^1$ test function $v \mapsto \phi(v)$, we require in the sense of distributions

$$\frac{d}{dt} \int_{-\infty}^{1} \phi(v) \rho(v,t) dv + \int_{-\infty}^{1} v \frac{d\phi}{dv}(v) \rho(v,t) dv + \sigma(t) \int_{-\infty}^{1} [\phi(v) - \phi(v+h)] \rho(v,t) dv \\ + \sigma(t) \int_{1-h}^{1} \rho(v,t) \phi(v+h) dv = I(t) \int_{-\infty}^{1} [\phi(v-h) - \phi(v)] \rho(v,t) dv + \phi(V_R) r(t). \quad (5.1)$$

Note that this formulation is obtained by extending $\rho$ and $\phi$ by 0 past $V_F = 1$ in the computations. Taking $\phi(v) = 1$ in (5.1) yields mass conservation:

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

<footnote>$^6$Beware that contrary to our convention for (TE) we do not include the connectivity strength $J$ into the transmission measure; in both cases we chose the most suitable notation to encompass the largest proportion of references on the topic.</footnote>
When working with measure solutions, we use the same formalism as for (TE) in Section 4. We denote $M(E)$ the set of finite measures on a measurable space $E$, $M_+(E)$ the set of positive finite measures, $P(E)$ the set of probability measures and $\|\cdot\|_{TV}$ the total variation norm.

**Definition 5.1** (measure solutions). Let $T > 0$, we say that a family $(\rho_t)_{t \ge 0}$ of $P([0, 1])$ is a solution to the excitatory ($J_I = 0, J_E \ge 0$) system (jLIF) on $[0, T)$ with initial datum $\rho_0 = \rho^0$ if

* $t \mapsto \sigma(t) := \frac{\sigma_0(t)}{1 - J_E \rho_t([1-h, 1])}$ is positive and locally integrable;

* $t \mapsto \rho_t$ is weak-* continuous on $[0, T]$;

* for all $\phi \in C^1([0, 1])$, for all $t \in [0, T)$,

$$ \int_0^1 \phi(v) d\rho_t = \int_0^1 \phi(v) d\rho_0 + \int_0^t \int_0^1 \left( -v \frac{d\phi}{dv}(v) + \sigma(s) \left[ \phi(v+h) \mathbb{1}_{[0,1-h]}(v) - \phi(v) + \phi(V_R) \mathbb{1}_{[1-h,1]}(v) \right] \right) d\rho_s ds. $$

Define the operators $A_\sigma$ and $B$ acting on $C^1([0, 1])$ by

$$ A_{\sigma}\phi(v) = -v \frac{d\phi}{dv}(v) + \sigma B\phi(v), \quad B\phi(v) = \phi(v+h)\mathbb{1}_{[0,1-h]}(v) - \phi(v) + \phi(V_R)\mathbb{1}_{[1-h,1]}(v). $$

Then, denoting $\rho_t\phi = \int_0^1 \phi(v)d\rho_t$, the last part of the definition of a measure solution reduces to

$$ \rho_t \phi = \rho_0 \phi + \int_0^t \rho_s A_{\sigma(s)} \phi ds. $$

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

### 5.3.1 Case with distributed delay

Consider the excitatory case $J_E \ge 0$ and $J_I = 0$. Let us first define a capped problem using a function $D(x) = \min(M, \max(0, x))$ for some constant $M \ge 0$:

$$ \frac{\partial \rho}{\partial t}(v, t) + \frac{\partial}{\partial v} [-v\rho(v, t)] + D(\sigma(t)) [\rho(v, t) - \rho(v - h, t)] = D(\sigma(t)) \int_{1-h}^{1} \rho(w, t) dw \delta_{v_R}. \quad (5.2) $$

Then, a solution to (5.2) is a solution to (jLIF) as long as $0 \le \sigma(t) \le M$.

Like for (TE), the method to prove well-posedness is to first consider the linear problem with $J = 0$. The solution can be constructed via a fixed-point argument in the Banach space $X = C([0, T], L^1(0, 1))$ endowed with its natural Banach norm on the mapping $\Lambda : q \to \rho$ where $\rho$ is the solution to

$$ \frac{\partial \rho}{\partial t}(v, t) + \frac{\partial}{\partial v} [-v\rho(v, t)] + D(\sigma(t)) [\rho(v, t) - q(v - h, t)] = D(\sigma(t)) \int_{1-h}^{1} q(w, t) dw \delta_{v_R}, $$

with same initial and boundary condition as (jLIF). For small enough $T$, $\Lambda$ is a contraction mapping and the argument can be iterated to obtain a global-in-time solution to the linear problem.

Then, the solution to the (capped) nonlinear problem can be obtained with a second fixed point argument in $C([0, T])$ applied to the mapping

$$ F : \sigma \mapsto \sigma_0 + J \int_0^t G[\sigma](t-\tau)\alpha(\tau)d\tau, \quad G[\sigma](t) = D(\sigma) \int_{1-h}^1 \rho(w,t)dw, $$

with $\rho$ solution to the linear problem with external input $\sigma$.

Last, by deriving local in time bounds of the form

$$ r(t) \le \|\sigma_0\|_{\infty} e^{J\|\alpha\|_{\infty} t}, $$

the cap $D$ can be removed on any finite time interval $[0, T]$ by letting $M$ be large enough.
**Theorem 5.2** (Dumont, Henry [96]). Assume $\mu_\sigma = \alpha(\tau)d\tau$, with $\alpha \in L^\infty(\mathbb{R}_+) \cap L^1(\mathbb{R}_+)$ and $\int_0^{+\infty} \alpha(\tau)d\tau = 1$, $\sigma_0 \in C(\mathbb{R}_+) \ge 0$, $J_I = 0$ and $J_E \ge 0$ (excitatory case). For all $\rho^0 \in L_+^1(0,1)$ of unit mass, there exists a unique solution $\rho \in C(\mathbb{R}_+, L_+^1(0,1))$ to (jLIF), and it satisfies

$$r(t) \le \|\sigma_0\|_{\infty} e^{J\|\alpha\|_{\infty}t}.$$

If $J_E < 1$, it satisfies the improved bound

$$r(t) \le \frac{\|\sigma_0\|_{\infty}}{1 - J_E}.$$

In the inhibitory case, the same strategy applies, but now the equation has to be considered on $(-\infty, 1]$ instead of just $[0, 1]$.

**Theorem 5.3** (Dumont, Henry [96]). Assume $\mu_\sigma = \alpha(\tau)d\tau$, with $\alpha \in L^\infty(\mathbb{R}_+) \cap L^1(\mathbb{R}_+)$ and $\int_0^{+\infty} \alpha(\tau)d\tau = 1$, $\sigma_0 \in C(\mathbb{R}_+) \ge 0$, $J_I \ge 0$ and $J_E = 0$ (inhibitory case). For all $\rho^0 \in L_+^1(-\infty, 1)$ of unit mass, there exists a unique solution $\rho \in C(\mathbb{R}_+, L_+^1(-\infty, 1))$ to (jLIF).

### 5.3.2 Instantaneous transmission

As was already noticed within a less rigorous formalism in [202], when $\sigma(t) = \sigma_0 + J\sigma(t)$, finite time blow-up can happen. It is easier to see in the reformulation

$$\sigma(t) = \frac{\sigma_0}{1 - J \int_{1-h}^{1} \rho(w, t) dw},$$

as when $J \int_{1-h}^{1} \rho(w, t) dw$ approaches 1, $\sigma(t)$ ceases to exist in the classical sense.

However, since $\rho$ is a probability density, as long as $J < 1$, we should have $J \int_{1-h}^{1} \rho(w, t) dw < 1$ and the firing rate $\sigma(t)$ should be well-defined.

**Theorem 5.4** (Dumont, Henry [96]). Assume $\mu_\sigma = \delta_{\tau=0}$, $\sigma_0 \in C(\mathbb{R}_+) \ge 0$, $J_I = 0$ and $J_E \in [0, 1)$ (weakly-excitatory case). For all $\rho^0 \in L_+^1(0, 1)$ of unit mass, there exists a unique solution $\rho \in C(\mathbb{R}_+, L_+^1(0, 1))$ to (jLIF), and it satisfies

$$r(t) \le \frac{\|\sigma_0\|_{\infty}}{1 - J_E}.$$

The proof follows the same pattern as in the delayed case, with some shortcuts. However, for a general excitatory connectivity $J_I = 0$, $J_E \ge 0$, as the solution may blow-up (see later), it is only possible to obtain a local-in-time solution. It can be sought in the space of admissible functions

$$\mathcal{X} = \left\{ \rho \in L_{+}^{1}(0, 1) \mid \int_{0}^{1} \rho(v) dv = 1, \quad J_{E} \int_{1-h}^{1} \rho(v) dv < 1 \right\}.$$

Then, it is possible to apply another fixed point procedure to a version of (jLIF) with a mollified impact of $\sigma$ and obtain the following result.

**Theorem 5.5** (Dumont, Henry [97]). Assume $\sigma_0 \in C_b(\mathbb{R}_+, \mathbb{R}_+)$, $J_I = 0$, $J_E \ge 0$ (excitatory case). For all $n^0 \in L_+^1(0, 1)$ of unit mass, there exists a unique maximal solution $\rho \in C([0, T^*], L_X^1(0, 1))$ to (jLIF).

The inhibitory case also proceeds similarly, but there is no issue with finite time blow-up even when $J_I$ is large, allowing to achieve unconditional global-in-time existence.

**Theorem 5.6** (Dumont, Henry [96]). Assume $\mu_\sigma = \delta_{\tau=0}$, $\sigma_0 \in C(\mathbb{R}_+) \ge 0$, and $J_I \ge 0$, $J_E = 0$ (inhibitory case). For all $\rho^0 \in L_+^1(-\infty, 1)$ of unit mass, there exists a unique solution $\rho \in C(\mathbb{R}_+, L_+^1(-\infty, 1))$ to (jLIF).
### 5.3.3 Solutions in the sense of measures

Consider the case $\mu_\sigma$ (instantaneous transmission) and $J_I = 0, J_E \in [0, 1)$ (weakly-excitatory). Following the duality method applied in [116] on a toy example, [96] constructs global-in-time solutions to (jLIF). See also the construction of measure solutions to the system (TE) in Subsection 4.3.2 which lies in the same toolbox and was performed a year earlier in [39]. If $\psi(s, t, v)$ satisfies

$$ \frac{\partial \psi}{\partial t} = v \frac{\partial \psi}{\partial v} + \frac{\sigma_0}{1 - J_E \mu_0 \psi(0, s, \cdot)} \left[ \psi(s, t, v) - \psi(s, t, v + h) \mathbb{1}_{[0,1-h)}(v) - \psi(s, t, V_R) \mathbb{1}_{[1-h,1]}(v) \right], $$

with terminal condition $\psi(t, t, v) = \mathbb{1}_{[1-h,1]}(v)$, then $\rho_t([1-h, 1]) = \rho_0\psi(0, t, \cdot)$. Solving the dual equation allows to obtain a Duhamel formula and to desingularise the problem, proceeding similarly to the construction of the linear semigroup below.

However, because of the discontinuity introduced by terms like $\mathbb{1}_{[1-h,1]}$, and because of the need to solve for $\psi$ in the space of continuous functions in order to get to $M([0, 1])$ by duality, the existence proof involves intermediate mollified $\psi_n$ equations and a lot of technicalities that we omit in our account.

**Theorem 5.7** (Dumont, Gabriel [95]). Assume $\mu_\sigma = \delta_{\tau=0}$ (instantaneous transmission), $J_I = 0, J_E < 1$ (weakly-excitatory case). Then for all $\rho_0 \in P([0, 1])$, there exists a unique global-in-time measure solution to (jLIF) in the sense of Definition 5.1.

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

In [97], it is proved that finite time blow-up does occur for the excitatory (jLIF). Unlike in the model (NNLIF), blow-up does not arise in weakly-connected networks. However, the universal blow-up property for high enough connectivity still holds.

**Theorem 5.8** (Dumont, Henry [97]). If $\mu_\sigma = \delta_{\tau=0}$ (instantaneous transmission) and

$$ J_I = 0, \quad J_E \ge \frac{1 - V_R}{h} + 1, \quad h\sigma_0 > 1, $$

then for all initial condition $\rho^0 \in L_+^1(0, 1)$ with unit mass, the solution $\rho$ to (jLIF) cannot be global-in-time. The maximal time of existence $T^*$ satisfies

$$ T^* \le \frac{1}{h\sigma_0 - 1}, \quad \limsup_{t \to T^*} r(t) = +\infty. $$

**Remark 5.9.** Under those conditions, there is in particular no stationary state to (jLIF).

The proof proceeds similar to the one of Theorem 1.12, by considering the weak formulation with test function $\phi(v) = e^{\mu v}$, for a constant $\mu > 0$ large enough (not to be confused with the measure $\mu_\sigma$). The exponential moment $M_\mu(t) = \int_0^1 \rho(v, t)e^{\mu v}dv$ then satisfies, after some computations and using the hypotheses,

$$ \frac{dM_\mu}{dt}(t) \ge K_\mu M_\mu(t) + \mu r(t) \left( JhM_\mu(t) - \frac{e^{1+h} - e^{V_R}}{\mu} \right). $$

for some constant $K_\mu > 0$. We have

$$ \frac{dM_\mu}{dt}(t) \ge \mu r(t) \left( JhM_\mu(t) - \frac{e^{1+h} - e^{V_R}}{\mu} \right). $$

So if $M_\mu(0) > \frac{e^{1+h} - e^{V_R}}{\mu}$, then $M_\mu$ is non-decreasing and

$$ \frac{dM_\mu}{dt}(t) \ge K_\mu M_\mu(t), $$

which leads to a contradiction with the natural bound $M_\mu(t) \le e^{\mu} \int_0^1 \rho(v, t)dv = e^{\mu}$. The problem is reduced to prove that for any initial condition there exists a time $t^*$ such that $M_\mu(0) > \frac{e^{1+h} - e^{V_R}}{\mu}$, which happens due to the strength of the external drift $\sigma_0$ in the assumptions.
## 5.5 Long time behaviour

In [95], the long time asymptotic of solutions in the linear and weakly-connected excitatory regimes are investigated for (jLIF). Like in [39], the Doeblin method allows to do so in the space of measures and to obtain quantitative global asymptotic stability of a unique invariant measure. Since we already described some parts of this procedure, we will not enter all technical details; see Subsection 4.6.3.

### 5.5.1 Study of the linear asymptotic

Starting with the spectral gap in the linear case, [95] builds on the method of [116] based on the dual equation

$$ \frac{\partial f}{\partial t} + v \frac{\partial f}{\partial v} + \sigma_0 f = \sigma_0 [f(v+h,t)1_{[0,1-h)}(v) + f(V_R,t)1_{[1-h,1]}(v)], $$

which they prove to be well-posed in the space of continuous functions. Then, it is possible to define a positive conservative semigroup $(M_t)_{t \ge 0}$ acting on $C([0, 1])$ such that $M_t f^0 = f(\cdot, t)$. It allows in turn to define by duality a semigroup $(M_t)_{t \ge 0}$ acting on $M([0, 1]) = C([0, 1])'$, which is also positive and conservative (because of the duality construction, this semigroup acts left). It preserves the set of probabilities $P([0, 1])$ and is a contraction in total variation norm.

After some technicalities involving regularisation, [96] proves that this semigroup provides solutions to the linear problem. They also prove that it satisfies Doeblin's condition, which leads, applying Doeblin theorem, to the following result.

**Theorem 5.10** (Dumont, Gabriel [95]). Assume $\mu_\sigma = \delta_{\tau=0}$ (instantaneous transmission), $\sigma_0(t) \equiv \sigma_0 > 0$ is constant and $J_I = J_E = 0$ (linear case). Then the equation (jLIF) generates a positive contraction left-acting weak-* continuous semigroup $(M_t)_{t \ge 0}$ on $M([0, 1])$. The set of probability measures $P([0, 1])$ is invariant under $(M_t)_{t \ge 0}$. It admits a unique invariant probability measure $\rho_\infty$ and for all $\rho^0 \in P([0, 1])$,

$$ \|\rho^0 M_t - \rho_\infty\|_{TV} \le e^{-\nu(t-t_0)} \|\rho^0 - \rho_\infty\|_{TV}, $$

where

$$ t_0 = \log \left( \frac{4}{h} \right) > 0, \quad \nu = - \frac{\log \left( 1 - \frac{\sigma_0}{2} \left( \frac{h}{4} \right)^{\sigma_0} \right)}{\log \left( \frac{4}{h} \right)} > 0. $$

### 5.5.2 Nonlinear stability and existence of stationary states

Existence of a nonlinear stationary state in the excitatory case with instantaneous transmission is equivalent to finding a constant $\sigma_\infty$ such that

$$ \sigma_{\infty} = \frac{\sigma_0}{1 - J_E \rho^{\sigma_{\infty}}([1-h, h])}, $$

where $\rho^{\sigma_\infty}$ is the invariant measure to the linear problem with external input $\sigma_\infty$. This done in [95] by defining the mappings $F, G$ by

$$ F(\sigma) = \rho^{\sigma}([1-h, h]), \quad G(\sigma) = \frac{1}{J_E} \left(1 - \frac{\sigma_0}{\sigma}\right), $$

and by looking for a solution to $F(\sigma_\infty) = G(\sigma_\infty)$.

First, one can prove that for fixed constant external $\sigma$, the invariant measure is absolutely continuous with respect to Lebesgue measure and enjoys some regularity;

$$ v \rho_{\infty}^{\sigma} \in W^{1,\infty}([0, V_R) \cup (V_R, 1]), \quad \forall v \in (0, 1), \ 0 < \rho_{\infty}^{\sigma}(v) \le \min \left( \frac{\sigma v^{\sigma-1}}{h^{\sigma}}, \frac{\sigma}{v} \right). $$
Note that then the nonlinear stationary state also enjoys such regularity. From this regularity and the decomposition of the measure formulation with $A_\sigma$ and $B$, [95] proves asymptotic properties on $F$:

$$ \lim_{\sigma \to 0} F(\sigma) = 0, \quad \lim_{\sigma \to +\infty} F(\sigma) = \frac{1}{1 + \left\lfloor \frac{1-V_R}{h} \right\rfloor}. $$

Combining these results, using the Duhamel formula and the contraction of the linear semigroup $(M_t)_{t \ge 0}$, and working some more yields

**Theorem 5.11** (Dumont, Gabriel [95]). Assume $\mu_\sigma = \delta_{\tau=0}$ (instantaneous transmission), $\sigma_0(t) \equiv \sigma_0 > 0$ is constant and $J_I = 0$, $J_E \ge 0$ (excitatory case).

* There exist at least two stationary states to (jLIF) under the condition

$$ J_E > 1 + \left\lfloor \frac{1 - V_R}{h} \right\rfloor, \quad \sigma_0 < \frac{1 - h}{4J_E}. $$

* There exists at least one stationary state to (jLIF) under the condition

$$ J_E < 1 + \left\lfloor \frac{1 - V_R}{h} \right\rfloor. $$

* If

$$ J_E < (5 - 2\sqrt{6}) \left(\frac{h}{4}\right)^{\sigma_0+1}, $$

there exists a unique stationary state $\rho_\infty$, and there are constants $t_0 \ge 0, \nu > 0$ such that for all initial condition $\rho_0 \in \mathcal{P}([0, 1])$ the solution to (jLIF) satisfies,

$$ \|\rho_t - \rho_\infty\|_{\text{TV}} \le e^{-\nu(t-t_0)} \|\rho_0 - \rho_\infty\|_{\text{TV}}. $$

## 5.6 Heuristic link to the NNLIF model

As noted in the physics literature [202] and explained in detail in [96], the singular structure of the jump mechanism in (jLIF) can be approximated by a Taylor expansion when the jumps of size $h$ are small enough, leading to a diffusive model which happens to be the model (NNLIF) from Section 1. Consider the excitatory case $J_I = 0, J_E \ge 0$. Denote $b_0 + bN(t) = h\sigma(t)$ for some $b_0 > 0, b \ge 0$. If $h$ is small enough, we can write the second order expansion

$$ \begin{aligned} \sigma(t)[\rho(v-h,t) - \rho(v,t)] &\approx \sigma(t)h \frac{\partial \rho}{\partial v} - \sigma(t) \frac{h^2}{2} \frac{\partial^2 \rho}{\partial v^2} \\ &\approx (b_0 + bN(t)) \frac{\partial \rho}{\partial v} - (b_0 + bN(t)) \frac{h}{2} \frac{\partial^2 \rho}{\partial v^2}, \end{aligned} $$

and, using that $\rho(1, t) = 0$

$$ \begin{aligned} r(t) = \sigma(t) \int_{1-h}^{1} \rho(v,t) dv &\approx \sigma(t) h \rho(1,t) - \sigma(t) \frac{h^2}{2} \frac{\partial \rho}{\partial v}(1,t) \\ &\approx -(b_0 + bN(t)) \frac{h}{2} \frac{\partial \rho}{\partial v}(1,t). \end{aligned} $$

Note then that for instantaneous transmission $\mu_\sigma = \delta_{\tau=0}$,

$$ b_0 + bN(t) = h\sigma(t) = h(\sigma_0 + J_E r(t)) = h\sigma_0 + hJ_E r(t), $$
which allows to identify $b_0 = \sigma_0 h$, $N(t) = r(t)$, $b = hJ_E$. Plugging these approximations into (jLIF), denoting

$$a(N(t)) = \frac{h}{2}(b_0 + bN(t)),$$

and rescaling the voltage $v$ so as to make a term $b_0$ disappear in the $-v$, eventually yields

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

with

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

It is exactly the model (NNLIF) in the case of intrinsic noise. Note that the right-hand side in the condition $J_E < (5 - 2\sqrt{6}) (\frac{h}{4})^{\sigma_0+1}$ for unconditional asymptotic stability in Theorem 5.11 goes to 0 as $h \to 0$, consistent with the fact that the limiting NNLIF model exhibits blow-up for any positive connectivity.

# 6 Some other hyperbolic models

We provide here a short description of some other hyperbolic PDEs used in mathematical neuroscience.

## 6.1 Leaky Integrate and Fire with random firing

Another hyperbolic model that was rigorously obtained as a mean-field limit [113, 78] is the Leaky Integrate and Fire with random firing. At the particle level, the neurons of potential $V_t^1, \dots, V_t^N$ evolve according to the system of SDEs

$$V_t^i = V_0^i + \int_0^t (-V_s^i + j_0)ds + \frac{J}{N} \sum_{j \neq i} \int_0^t \int_0^{+\infty} 1_{z \le \phi_{\varepsilon}(V_s^j)} N^j(ds, dz) - \int_0^t \int_0^{+\infty} V_{s-}^i 1_{z \le \phi_{\varepsilon}(V_s^i)} N^i(ds, dz),$$

where the $N^i(ds, dz)$ are independent Poisson measures with intensity $dsdz$, $v \mapsto \phi_{\varepsilon}(v)$ is a firing function, $J > 0$ a connectivity parameter, $j_0 > 0$ an external excitatory input. In the limit $N \to +\infty$, the system converges towards

$$V_t = V_0^i + \int_0^t (-V_s + j_0)ds + J \int_0^t \mathbb{E}[\phi_{\varepsilon}(X_s)] ds - \int_0^t \int_0^{+\infty} V_{s-} 1_{z \le \phi_{\varepsilon}(V_{s-}^i)} N(ds, dz), \quad (6.1)$$

whose law is a solution of

$$\left\{
\begin{array}{l}
\displaystyle \frac{\partial \rho}{\partial t} + \frac{\partial}{\partial v} [(-v + j_0 + Jr(t))\rho] + \phi_{\varepsilon}(v)\rho = 0, \\
r(t) = \displaystyle \int_0^{+\infty} \phi_{\varepsilon}(v)\rho(v, t)dv, \\
\rho(0, t) = \displaystyle \frac{r(t)}{j_0 + Jr(t)}, \quad \int_0^{+\infty} \rho^0(v)dv = 1, \quad \rho(v, 0) = \rho^0(v) \ge 0,
\end{array}
\right.
\qquad (6.2)$$

The model (6.2) shares similarities both with the LIF with jump (jLIF) and the NNLIF model with random firing (1.33).

In [70], the model is studied at the level of the stochastic mean-field equation (6.1). They prove that for small enough $J$ there is a unique invariant measure and that the solutions converge towards it when $t \to +\infty$. In [71], it is rigorously proved that the model (6.2) undergoes a Hopf bifurcation in terms of the parameter $\varepsilon$ for well-chosen $\phi_{\varepsilon}$, which proves the existence of time-periodic solutions (see Figure
<figure><img src="image_15.png" /><figcaption>Figure 13: Periodic solution in the LIF model with random firing. Left is the firing rate $N(t)$ in the PDE (6.2), right represents spikes from the mean-field particle system (6.1). Cormier, Tanré, Veltz [71].</figcaption></figure>

13). For both studies [70, 71], a crucial tool is a Volterra equation related to the system. Further results about long-time behaviour and bistability were recently proved in [69].

Note that (6.2) was also studied from the PDE point of view in the linear case in [212] as the noiseless limit of (7.13) below in Subsection 7.7.1. Conditions are given for linear instability of the stationary state in the nonlinear case, which can be related to the rigorously obtained periodic solutions in [71].

## 6.2 Continuity equation for pulse coupled oscillators

A simple model for interacting neurons are the so-called pulse coupled oscillators first introduced in the context of modelling cardiac cells [214]. It consists in a system of ordinary differential equations (ODE) coupled by singular interaction. Consider $V^1(t), \dots, V^n(t)$ the potential at time $t$ of $n$ neurons. Each of them evolves according to

$$ \frac{dV^i}{dt}(t) = f(V^i(t)), $$

and when $V^i(t^-) \ge V_F$, $V^i(t^+) = V_R = 0$ and the other neurons receive a kick $V^j(t^+) = V^j(t^-) + \frac{J}{n}$, $j \ne i$, $J \ge 0$. This is similar to the construction of (NNLIF) in Section 1, but in a fully deterministic framework. This deterministic particle system has been investigated in depth, with the milestone [180] establishing simple conditions for synchronisation or relaxation to the stationary state. A common strategy is to rescale the system from voltage variable $v \in [0, V_F]$ to phase variable $\phi \in [0, \Phi_F]$ with a constant vector field and a modified jump function:

$$ \frac{d\phi^i}{dt}(t) = 1, \quad \phi^i(t^-) \ge V_F \implies \phi^i(t^+) = 0, \phi^j(t^+) = \phi^j(t^-) + \frac{1}{n}K(\phi^j), \quad j \ne i. $$

A formal mean-field limit [1, 157] yields the continuity equation

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial \rho}{\partial t} + \frac{\partial}{\partial \phi} [(1 + K(\phi)N(t))\rho] = 0, \\[1.5em] \displaystyle N(t) = (1 + K(\Phi_F)N(t))\rho(\Phi_F, t), \\[1.5em] \displaystyle (1 + K(\Phi_F)N(t))\rho(\Phi_F, t) = (1 + K(0)N(t))\rho(0, t), \\[1.5em] \displaystyle \int_0^{\Phi_F} \rho^0(\phi) d\phi = 1, \quad \rho(\phi, 0) = \rho^0(\phi) \ge 0. \end{array} \right. \qquad (6.3) $$
Notice that as in (jLIF), we have for the firing rate N the closed expression

$$N(t) = \frac{\rho(t, \Phi_F)}{1 - K(\Phi_F)\rho(t, \Phi_F)}.$$

which indicates the model can blow-up in finite time.

This nonlinear mean-field problem has scarcely been investigated from the mathematical perspective. Criteria for the existence of stationary states and partial results about the synchronisation-relaxation dichotomy have been provided first in [176]. In particular, they show, similar to the findings for the particle system [180], that synchronisation occurs for increasing K and relaxation to a stationary state for decreasing K and adequate initial conditions. However, their results require many technical assumptions. Another approach has been proposed in [46]: consider the repartition function $F(\phi, t) = \int_0^\phi \rho(\psi, t)d\psi$ and its left-continuous pseudo-inverse

$$Q(\eta) = \inf\{\phi \ge 0 : F(\phi) \ge \eta\}, \quad \eta \in [0, 1].$$

The function Q is also known as the quantile function of ρ. Consider the time dilation dτ = N(t)dt; then, Q(η, τ) solves

$$\left\{
\begin{array}{ll}
\partial_{\tau} Q + \partial_{\eta} Q = \frac{1}{N(\tau)} + K(Q), & \tau > 0, \eta \in (0, 1), \\
Q(\tau, 0) = 0, & \tau > 0, \\
\frac{1}{N(\tau)} = \partial_{\eta} Q(\Phi_F, \tau) - K(\Phi_F), & \tau > 0, \\
Q(\eta, 0) = Q^0(\eta), & \eta \in [0, 1].
\end{array}
\right.
\qquad (6.4)$$

Note that this kind of time dilation has been used in previous [91, 252, 236] and subsequent [90] works. The simpler drift and the decoupled nonlinearities 1/N and K(Q) in (6.4) have allowed [46] to propose generalised versions of the results in [176] and to pave the road for continuation of solutions to (6.3) after blow-up.

# Part III

## Kinetic Fokker-Planck equations

### 7 The kinetic Voltage-Conductance model

We discuss here some kinetic voltage-conductance descriptions of neural networks based on the Morris-Lecar model for individual neurons. The degenerate diffusion, nonlinearity, and peculiar boundary conditions in these models make them difficult to study analytically, and many questions are still widely open.

We start by describing the historical construction of kinetic voltage-conductance models at the network level. Then we propose a mathematical framework, and we discuss how the already difficult linear case can be tackled with methods from kinetic theory or the use of probabilistic tools like Doeblin-Harris techniques. We then turn to the nonlinear case and discuss both a priori estimates relying on direct integration and moment estimates, and a recent, more abstract approach which allows to prove the existence of solutions for weakly-connected networks. Finally, we address the issues around numerical simulations of this model and present some asymptotic limits and variants.

#### 7.1 Fluctuations at the level of the conductance

The models we have discussed so far are overlooking an important aspect of the dynamics of individual neurons: the evolution of the membrane potential v of individual cell depends on underlying variables
with faster timescales. In large and mostly inhibitory networks observed over long periods of time, the voltage or activity descriptions provide good approximations, which is why the delayed inhibitory NNLIF model in Section 1 or the activity-based grid cell model in Section 2 are good at capturing macroscopic behaviours like self-sustained oscillations and the hexagonal firing pattern of grid cells.

However, as we have seen, those models are less precise for the modelling of excitatory ensembles of neurons. For example, in the primary visual cortex V1, so-called simple neurons can respond in a linear fashion to frequency in visual stimulation, when others, called complex, respond nonlinearly, see [140]. It is believed that this simple or complex nature is not an individual characteristic of those cells but rather an emerging property from the network they are embedded in [58]. The behaviour of the complex type of V1 neuron has been obtained in a model with a strong excitation; The simple type has been found to be better represented by inhibitory models; see [253].

More generally, accurately describing the fluctuation-driven dynamics of highly nonlinear networks of neurons with slow-postsynaptic responses requires to work at smaller scales for space and time. Around 2005, in the pioneering works [36, 35, 221], based upon previous literature about deterministic and stochastic representation of coarse-grained neural networks (see e.g. [112, 150, 198, 199, 202, 201, 19]), Cai, McLaughlin, Rangan, Shelley and Tao developed an accurate kinetic description for excitatory and excitatory-inhibitory ensembles of neurons.

Let us focus on the fully excitatory case. They start with a coarse-grained patch of N all-to-all connected excitatory neurons described by their voltage V and conductance G: assume each neuron labelled by j ∈ {1, ..., N} evolves according to the Integrate & Fire equations

$$ 
\begin{aligned}
\tau \frac{dV^j}{dt} &= -(V^j - V_L) - G^j(V^j - V_E), \\
\sigma_E \frac{dG^j}{dt} &= -G^j + S_{\text{ext}} \sum_{l=1}^{+\infty} \delta(t - T_j^l) + \frac{S_E}{N} \sum_{k=1, k \neq j}^{N} \sum_{l=1}^{+\infty} \delta(t - t_k^l),
\end{aligned}
$$

where $V_E$ is the excitatory reversal potential, $\tau$ is the leakage time constant, $\sigma_E$ the time constant of the synaptic receptors, $S_{\text{ext}}$ the strength of connections from outside the network, $S_E$ the strength of synapses inside the network, $T_j^l$ is the time of the $l^{th}$ exterior spike coming on neuron j and $t_k^l$ is the time of the $l^{th}$ spike coming from the neuron k onto the neuron j. We can compare with the simpler Laplacian model (1.1) that was used in [19] to derive the model (NNLIF): here the incoming current is replaced by a more complex dynamics involving the conductance of the cell membrane.

When neurons reach the voltage firing threshold $V_F$, which is $V(t^-) = V_F$, they are reset to the reset potential $V_R$, which is $(V(t^+), G(t^+)) = (V_R, G(t^-))$. From each individual neuron k, the output train statistics are not Poisson in general. However, the inputs received from the network by each neuron j is a spike train summed over many output spike trains from the other neurons. If we assume that these incoming spikes are asymptotically independent from each other, the total input from the network can be assumed to be Poisson [62] or [159, pp. 549–606]. Applying this simplification and considering the discrete probability density

$$ \rho^N(v, g, t) = \mathbb{E} \left[ \frac{1}{N} \sum_{j=1}^{N} \delta(v = V^j(t)) \delta(g = G^j(t)) \right], $$

where $\mathbb{E}$ is the expectation with respect to all possible sets of N independent realisations of the exterior Poisson spikes ($T_j^l$)$_{j \in \{1, ..., N\}, l \in \mathbb{N}}$, [36, 35] derive by a moment closure approach a (2+1)-dimensional kinetic Fokker-Planck equation for the density $\rho(v, g, t)$ in the limit $N \to +\infty$. Up to rescaling and assuming that $0 = V_R = V_L < V_F < V_E$, it writes

$$ \frac{\partial \rho}{\partial t} + \frac{\partial}{\partial v} [J_v(v, g)\rho(v, g, t)] + \frac{\partial}{\partial g} [J_g(g, \mathcal{N}(t), t)\rho(v, g, t)] - \frac{a(\mathcal{N}(t), t)}{\sigma_E} \frac{\partial^2 \rho}{\partial g^2} = 0, \quad (7.1) $$

where the fluxes are defined by

$$ J_v(V_F, t) = -g_L v + g(V_E - v), \quad J_g(g, \mathcal{N}, t) = \frac{1}{\sigma_E} (S_{\text{ext}}\nu(t) + S_E\mathcal{N} - g), $$
the firing rate $N(t)$ and the intensity of the synaptic noise $a(N(t), t)$ are characterised by

$$ \mathcal{N}(t) = \int_{0}^{+\infty} N(g, t)dg, \quad N(g, t) = J_v(V_F, g)\rho(V_F, g, t), \\ a(\mathcal{N}(t), t) = \frac{1}{2\sigma_E} \left( S_{\text{ext}}^2\nu(t) + \frac{S_E}{N_E}\mathcal{N}(t) \right). $$

and the system is supplemented with initial and boundary constraints

$$ J_g(0, \mathcal{N}(t), t)\rho(v, 0, t) - \frac{a(\mathcal{N}(t), t)}{\sigma_E} \frac{\partial \rho}{\partial g}(v, 0, t) = 0 \quad \text{and} \quad \text{BC (7.3)} \qquad (7.2) $$

and

$$ \rho(v, g, 0) = \rho^0(v, g) \ge 0, \quad \int_0^{V_F} \int_0^{+\infty} \rho^0(v, g)dgdv = 1. $$

Let us comment on the boundary conditions (7.2). The first part is necessary for the conductance to stay positive: we impose a dynamic no-flux condition on $g = 0$. We could also impose this condition at $g = +\infty$, but in order to avoid technical complications, we will simply assume fast decay at $+\infty$ for $\rho$ and $\partial_v\rho$. The second boundary condition (7.3), that we have not described yet, is more tricky. There is no diffusion in $v$, so the fluxes entering or exiting through the boundary at $v = 0$ or $v = V_F$ depend only on $J_v\rho$. The model prescribes that nothing comes in the domain from $\{V_F\} \times (0, +\infty)$ and anything that goes out through this boundary line must come in through the line $\{0\} \times (0, +\infty)$. Looking at the drift $J_v$, we see that it vanishes on a curve$^7$

$$ g = \frac{g_L v}{V_E - v} $$

and is always non-negative when

$$ g \ge g_F := \frac{g_L V_F}{V_E - V_F}. $$

Below this critical conductance ($g \le g_F$), the fluxes are coming inside the domain from both sides, so they are set to 0 (if nothing comes out at $v = V_F$, nothing should come in at $v = 0$). Hence, we impose the following flux conditions in the $v$ direction:

$$ \left\{ \begin{array}{ll} \rho(v, g, t) = \rho(V_F, g, t) = 0, & g \in (0, g_F], \\ J_v(0, g)\rho(0, g, t) = J_v(V_F, g)\rho(V_F, g, t) = N(g, t), & g \in (g_F, +\infty). \end{array} \right. \qquad (7.3) $$

Given the profusion of parameters, let us also recapitulate about their meaning.

* The exterior environment (other networks, visual stimuli, etc) acts through the given function $v(t) > 0$, that we will later take as constant, amplified by the synaptic strength $S_{\text{ext}} > 0$.

* The internal parameters are the average synaptic strength $S_E > 0$, the normalisation factor $N_E$, the leak conductance $g_L$ and the decay rate $\sigma_E$ for the conductance.

* The neurons are driven towards the leak potential $V_L = 0$ with conductance $g_L$ and towards the excitatory reversal potential $V_E > V_F$ with conductance $g$; when a neuron reaches the potential $V_F$ on its way to $V_E$, it fires (see the relaxed model (7.22) below for a smoother approach).

Although we will not discuss it further, note that [35, Section 5] proposes a model for excitatory-inhibitory networks. It consists in introducing two separate conductance variables $g_E$ and $g_I$, two excitatory and inhibitory reversal potentials $V_I < V_E$ and a $v$-flux defined as

$$ J_v(v, g_E, g_I) = g_L(V_L - v) + g_E(V_E - v) + g_I(V_I - v). $$

<footnote>$^7$This will be an issue for mathematical analysis later.</footnote>
## 7.2 Notation for mathematical analysis

Consider, performing renormalisations and taking constant external inputs, the fluxes and noise term

$$J_v(v, g) = -g_L v + g(V_E - v), \quad J_g(g, \mathcal{N}) = \frac{s_0 + S_E \mathcal{N}(t) - g}{\sigma_E}, \quad a(\mathcal{N}) = \frac{a_0 + a_1 S_E^2 \mathcal{N}}{\sigma_E},$$

with $g_L, V_E, s_0, a_0, \sigma_E > 0$ positive parameters and $a_1, S_E > 0$ non-negative connectivity coefficients. Recall the crucial g threshold

$$g_F := \frac{g_L V_F}{V_E - V_F},$$

where the boundary condition switches. We can sum up the Voltage-Conductance kinetic model as

$$\left\{
\begin{array}{l}
\displaystyle \frac{\partial \rho}{\partial t} + \frac{\partial}{\partial v} [(-g_L v + g(V_E - v))\rho] + \frac{\partial}{\partial g} \left[ \frac{s_0 + S_E \mathcal{N}(t) - g}{\sigma_E} \rho \right] - \frac{a(\mathcal{N}(t))}{\sigma_E} \frac{\partial^2 \rho}{\partial g^2} = 0, \\
\\
\displaystyle \mathcal{N}(t) = \int_0^{+\infty} N(g, t) dg, \quad N(g, t) = J_v(V_F, g) \rho(V_F, g, t), \\
\\
\displaystyle \rho(0, g, t) = \rho(V_F, g, t) = 0, \quad g \le g_F, \quad J_v(0, g) \rho(0, g, t) = J_v(V_F, g) \rho(V_F, g, t), \quad \text{(V-C)} \\
\\
\displaystyle J_g(0, \mathcal{N}(t)) \rho(v, 0, t) - \frac{a(\mathcal{N}(t))}{\sigma_E} \frac{\partial \rho}{\partial g}(v, 0, t) = 0, \\
\\
\displaystyle \rho(v, g, 0) = \rho^0(v, g) \ge 0, \quad \int_0^{V_F} \int_0^{+\infty} \rho^0(v, g) dg dv = 1.
\end{array}
\right.$$

In most of the results about (V-C), the solutions have to be understood in the distributional sense, or in the following weak sense first introduced in [87]: consider the test function space

$$\Phi = \left\{ \phi \in C_c^\infty([0, V_L] \times \mathbb{R}_+ \times \mathbb{R}_+) \mid \phi(0, g, t) = \phi(V_F, g, t), \frac{\partial g}{\partial g}(v, 0, t) = 0 \right\}.$$

then a weak-solution $\rho \in L^\infty(\mathbb{R}_+, L_+^1([0, V_L] \times \mathbb{R}_+))$ with firing rate $\mathcal{N} \in L_{\text{loc}}^1(R_+)$ satisfies

$$\int_0^{V_F} \int_{-\infty}^{+\infty} \phi(v, g, T) \rho(v, g, T) dv dg = \int_0^{V_F} \int_{-\infty}^{+\infty} \phi(v, g, 0) \rho(v, g, 0) dv dg \\ + \int_0^T \int_0^{V_F} \int_{-\infty}^{+\infty} \left[ \frac{\partial \phi}{\partial t} + J_v \frac{\partial \phi}{\partial v} + J_g(g, \mathcal{N}(t)) \frac{\partial \phi}{\partial g} + \frac{a(\mathcal{N}(t))}{\sigma_E} \frac{\partial \phi}{\partial g^2} \right] \rho(v, g, t) dv dg dt.$$

Note that the test functions can have any values $\phi(v, 0, t)$, which is crucial to ensure the boundary condition $J_g(0, \mathcal{N}(t))\rho(v, 0, t) - \frac{a(\mathcal{N}(t))}{\sigma_E}\partial_g\rho(v, 0, t) = 0$.

Most of the computations in articles about (V-C), for example in [211, 212, 146, 89], are done in the sense of distributions but would work in function spaces with this notion of weak solutions. It relieves a large part of the notational burden to just integrate and differentiate everything a priori and then prove a posteriori higher $L^q$ or Sobolev regularity from the propagation of finiteness in norms.

## 7.3 The linear Voltage-Conductance model

Consider the simplified case of a linear equation (V-C) by taking $S_E = 0$; this case has been studied first in [211], later improved in [89] and recently generalised to weak interaction in [111]. Then, there is no interaction between neurons. Yet, because of the degenerate structure of the differential operator, the vanishing v-flux and the complex boundary conditions, the mathematical analysis of solutions requires some care.

Indeed, consider a stationary state $(\rho_\infty, \mathcal{N}_\infty)$ of the linear problem: contrary to the stationary states of (NNLIF), (SNF) or (TE), there is no semi-implicit formulation or parabolic regularity to ensure they are more than $L^1$ functions. They have to be constructed in a more abstract way, akin to the stationary solutions of (jLIF). As we will explain, it is possible to obtain more regularity for linear stationary states by carefully performed iteration methods [211, 89].
**Theorem 7.1** (Dou, Perthame, Salort, Zhou [89]). Assume $S_E = 0$. There exists a unique stationary state $\rho_\infty \in L^\infty([0, V_F] \times \mathbb{R}_+)$ to (V-C).

Although we include it in the theorem for clarity, existence was supposed a priori in [211, 89] and proved in the linear case first in [87] and then in [111] which extends it to the weakly-nonlinear case.

### 7.3.1 A priori properties of linear stationary states

Denote

$$a = \frac{a_0}{\sigma_E}, \quad M(g) = e^{-\frac{(g-s_0)^2}{2a}}, \quad Z = \int_0^{+\infty} M(g)dg.$$

The stationary equation for (V-C) can be rewritten

$$\frac{\partial}{\partial v}[J_v(v, g) \rho_\infty(v, g)] - \frac{a}{\sigma_E} \frac{\partial}{\partial g} \left[ M(g) \frac{\partial}{\partial g} \left( \frac{\rho_\infty(v, g)}{M(g)} \right) \right] = 0.$$

By direct integration and simple computations,

$$\int_0^{V_F} \rho_\infty(v, g)dv = \frac{M(g)}{Z}.$$

Then, by multiplying the equation by $\log \rho_\infty$ and integrating, it is possible to prove:

**Lemma 7.2** (Perthame, Salort [211]). Assume $S_E = 0$. There exist positive constants $K_0, K_1$ depending on $N_\infty$ such that

$$\int_0^{V_F} \int_0^{+\infty} \frac{1}{\rho(v, g)} \left( \frac{\partial \rho}{\partial g}(v, g) \right)^2 dg dv \le K_0, \quad \int_0^{V_F} \int_0^{+\infty} e^{\frac{g^2}{8a}} \left| \frac{\partial \rho}{\partial g}(v, g) \right| dg dv \le K_1.$$

From there, the strategy of [211] was to prove that $J_v \rho_\infty$ was lying in the Besov space

$$B_\infty^{\frac{1}{2},1}((0, V_R) \times \mathbb{R}_+) = \left\{ f \in L^1 \left| \sup_{h \le 1} h^{-\frac{1}{2}} \|f(\cdot + h) - f(\cdot)\|_{L^1} \right. \right\}.$$

Then, by Besov injection, they could get $J_v \rho_\infty \in L^q$ for all $1 \le q \le \frac{4}{3}$. Unfortunately, because $J_v$ vanishes on a curve in the domain, they could not get better than $\rho_\infty \in L^q$ for $1 \le q < \frac{8}{7}$.

A significant improvement was proposed in [89]. They bypass the Besov spaces by writing

$$\rho_\infty(v, g) = - \int_g^{+\infty} \frac{\partial \rho_\infty}{\partial g} dg \le \int_0^{+\infty} \left| \frac{\partial \rho_\infty}{\partial g}(v, g') \right| dg',$$

and performing the direct bound

$$\int_0^{V_F} \sup_{g \in \mathbb{R}_+} \rho_\infty(v, g) dv \le \int_0^{V_F} \int_0^{+\infty} \left| \frac{\partial \rho_\infty}{\partial g}(v, g') \right| dg' \le K_1.$$

After more computations and using Lemma 7.2, they build a set of useful estimates.

**Lemma 7.3** (Dou, Perthame, Salort, Zhou [89]). Assume $S_E = 0$. There is a constant $K^*$ such that

$$\|\rho_\infty\|_{L^1([0, V_F]; L^\infty(\mathbb{R}_+))} \le K_1, \quad \|J_v^2 \rho_\infty\|_{L^\infty([0, V_F]; L^1(\mathbb{R}_+))} \le K^*, \quad \|J_v \rho_\infty\|_{L^2} \le K_1 K^*$$

and for all $q \in (0, \frac{4}{3})$,

$$\int_0^{V_F} \int_0^{+\infty} (g + g_L)^2 \rho(v, g)^q dg dv < +\infty.$$
The direct sup<sub>g</sub> bound is possible due to the smoother parabolic g structure of the equation, embedded in the constant K<sub>1</sub> of Lemma 7.2. It is not possible to do the same in v, which is why [89] has to derive the weaker sup<sub>v</sub> bound

$$ \| J_v^2 \rho \|_{L^\infty([0, V_F]; L^1(\mathbb{R}_+))} = \sup_{v \in [0, V_F]} \int_0^{+\infty} J_v(v, g)^2 \rho_\infty(v, g) dg \le K^*, $$

before getting to the L<sup>2</sup> bound on J<sub>v</sub> ρ<sub>∞</sub>. Note that this constitutes a slight improvement from J<sub>v</sub>ρ<sub>∞</sub> ∈ L<sup>4</sup><sub>3</sub><sup>-</sup> in [211] to J<sub>v</sub>ρ<sub>∞</sub> ∈ L<sup>2</sup> in [89]. After carefully overcoming the curve of zeros of J<sub>v</sub>, it allows to go from ρ<sub>∞</sub> ∈ L<sup>8</sup><sub>7</sub><sup>-</sup> to (g + g<sub>L</sub>)ρ<sub>∞</sub> ∈ L<sup>4</sup><sub>3</sub><sup>-</sup>.

This does not sound impressive, but the difference allows to make more progress. Using the Alikakos iteration method [2] starting from L<sup>4</sup><sub>3</sub><sup>-</sup> and using specific weights, [89] manages to get to ρ<sub>∞</sub> ∈ L<sup>∞</sup>. Uniqueness of linear stationary states had already been obtained at hand in [211].

### 7.3.2 Asymptotic stability in the linear case via an entropy method

Let H ∈ C<sup>2</sup> a convex function. The entropy dissipation for solutions of (V-C) in the linear case S<sub>E</sub> = 0 writes

$$ \frac{d}{dt} \int_0^{V_F} \int_0^{+\infty} H\left(\frac{\rho}{\rho_\infty}\right) \rho_\infty dv dg \le - \int_0^{V_F} \int_0^{+\infty} H''\left(\frac{\rho}{\rho_\infty}\right) \left[\frac{\partial}{\partial g}\left(\frac{\rho}{\rho_\infty}\right)\right]^2 \rho_\infty dv dg \le 0. $$

Provided ρ<sup>0</sup><sub>∞</sub> ∈ L<sup>∞</sup>, following the method of [209, Chapter 8] and using the regularity in Theorem 7.1, one can deduce from the entropy dissipation a bound of the form

$$ C_{-}^{0} \le \rho(v, g, t) \le C_{+}^{0}, \qquad (7.4) $$

with C<sub>-</sub><sup>0</sup>, C<sub>+</sub><sup>0</sup> > 0 depending on ρ<sup>0</sup><sub>∞</sub>. The dissipation not being coercive, it is not possible to get exponential decay in relative entropy. Instead, [89] uses a compactness argument by defining

$$ p_k(v, g, t) = \rho(v, g, t + k), \quad h_k(v, g, t) = \frac{\rho(v, g, t + k)}{\rho_\infty(v, g, t + k)}. $$

By (7.4), h<sub>k</sub> is bounded in L<sup>∞</sup>. Thus, there is a subsequence h<sub>k(n)</sub> converging weakly towards some h<sub>∞</sub> when n → +∞.

With the choice H(x) = (x - 1)<sup>2</sup>, the entropy dissipation allows first to prove that h<sub>∞</sub>(v, g, t) does not depend on g, and then that ∂<sub>t</sub>h<sub>∞</sub> + ∂<sub>v</sub>h<sub>∞</sub> = 0, implying by mass conservation that h<sub>∞</sub> ≡ 1. In turn, this allows to prove time compactness for h<sub>k</sub> and p<sub>k</sub>, and to conclude using Aubin-Lions lemma and the structure of the equation that p<sub>k</sub> converges towards ρ<sub>∞</sub> almost everywhere and in L<sup>q</sup><sub>loc</sub>, 1 ≤ q < +∞.

The last step is then to use these converging sequences to prove the decay of the relative entropy.

**Theorem 7.4** (Dou, Perthame, Salort, Zhou [89]). Assume we are in the linear case S<sub>E</sub> = 0. Assume there exists C > 0 such that ρ<sup>0</sup> ≤ Cρ<sub>∞</sub>. Then, for all time t ≥ 0, ρ(·, ·, t) ≤ Cρ<sub>∞</sub> and

$$ \lim_{t \to +\infty} \int_{0}^{V_F} \int_{0}^{+\infty} \left( \frac{\rho(v, g, t) - \rho_{\infty}(v, g)}{\rho_{\infty}(v, g)} \right)^2 \rho_{\infty}(v, g) dg dv = 0. $$

Note however that in the linear case, the entropy method does allow to prove exponential convergence to a stationary state if we first integrate the voltage variable. Denote

$$ M_{\infty}(g) = \frac{M(g)}{Z}, \quad \varphi(g, t) = \int_{0}^{V_F} \rho(v, g, t) dv. $$

Then, [211] computes,

$$ \frac{d}{dt} \int_0^{+\infty} \left( \frac{\varphi(g, t) - M_\infty(g)}{M_\infty(g)} \right)^2 M_\infty(g) dg \le - \frac{a}{\sigma_E} \int_0^{+\infty} \left[ \frac{\partial}{\partial g} \left( \frac{\varphi(g, t)}{M_\infty(g)} \right) \right]^2 M_\infty(g) dg, $$

which by a Poincaré-like inequality and standard computations yields for some ν > 0,

$$ \int_{0}^{+\infty} \left( \frac{\varphi(g, t) - M(g)}{M_{\infty}(g)} \right)^2 M_{\infty}(g) dg \le e^{-\nu t} \int_{0}^{+\infty} \left( \frac{\varphi(g, 0) - M(g)}{M_{\infty}(g)} \right)^2 M_{\infty}(g) dg. $$
### 7.3.3 Probabilistic reformulation and exponential ergodicity

There are two issues in the approach from [211, 89] described so far. First, the knowledge of the properties of the stationary state is a priori and this is not enough to prove existence. Second, the convergence to a prospective stationary state is obtained by compactness, and there is thus no rate of convergence. To circumvent the difficulty of the degenerate diffusion in (V-C), [87] proposes to use a probabilistic reformulation (see also the similar approach in [135] described below in Subsection 7.8.2) in order to take advantage of the ergodicity of the microscopic trajectories.

Consider again $S_E = 0$ and fix for the sake of simplicity $a_0 = \sigma_E = s_0 = V_F = 1$. Denote $(V_t, G_t)$ the solution of the system of stochastic differential equations

$$ egin{cases} 
dV_t = J_v(V_t, G_t)dt, \\ 
dG_t = -(G_t - 1)dt + \sqrt{2}dB_t + dL_t, \\ 
V_{t-} = 1 \implies V_t = 0, 
\end{cases} 
\quad (7.5) $$

with $(B_t)_{t \ge 0}$ a standard Brownian motion and $L_t$ the local time preventing the conductance $G_t$ from becoming negative (see [227, Chapter 6] or the discussion on a similar issue in Section 2).

Since we are in the linear case, there is no feedback on $G_t$ from the behaviour of $V_t$. Hence, as we can see in (7.5), the equation for $G_t$ is autonomous and constitutes a standard reflected Ornstein-Uhlenbeck SDE, whose solution exist and enjoy many good properties. Thus, following [82], [87] considers the sequence of jump times of $V_t$,

$$ \tau_k = \inf\{t > \tau_{k-1} \mid \lim_{s \to t^-} V_s = 1\}, $$

and proves that $\mathbb{P}(\lim_{k \to +\infty} \tau_k) = 1$, which implies that there exists a trajectory-wise almost surely unique solution to (7.5).

Then, [87] proceeds to prove that the law of a solution to (7.5) is a weak solution to the linear (V-C), a key difficulty being the boundary behaviour, they prove separately that

$$ \mathbb{P}(V_t = 0) = \mathbb{P}(V_{t-} = 1) = \mathbb{P}(G_t = 0) = 0. $$

Given the hypoellipticity of the linear operator defining the PDE and the stochastic process, there is enough interior regularity to ensure that the law of $(V_t, G_t)$ is absolutely continuous with respect to Lebesgue measure, which allows to build a link to weak solutions of the PDE using Itô's formula for semi-martingales with jumps [217, Theorem 2.33].

Last, they prove that the v-boundary condition is satisfied for weak solutions in the sense that for test functions $\psi \in C_c^\infty(\mathbb{R}_+^2)$,

$$ \lim_{v \to 0^+} \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} J_v(v, g) \rho(v, g, t) \psi(g, t) dg dt = \lim_{v \to 1^-} \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} J_v(v, g) \rho(v, g, t) \psi(g, t) dg dt. $$

To prove the exponential convergence of the density associated to $(V_t, G_t)$, [87] uses the Doeblin-Harris method. See our presentation of Doeblin's theorem in the context of the Time Elapsed neuron model in Subsection 4.6.3. For (7.5), Doeblin theorem is not enough. Hence, like for example [110] for the system (4.34), they use Harris theorem [128], which requires Doeblin's minoration condition (Definition 4.36) on a compact subset only, but as a trade-off requires on top a Lyapunov control on the stochastic semigroup $(S_t)_{t \ge 0}$:

$$ S_t \omega \le \alpha_1 \omega + \alpha_2, $$

for some weight function $\omega(v, g) \ge 0$ and constants $\alpha_1 \in (0, 1)$ and $\alpha_2 \ge 0$. The idea is that the Lyapunov function sends back the mass towards the compact subset where the Doeblin condition is applied, which is exploiting the ergodicity of the process and the good properties at a location it visits by ergodicity. Given both Doeblin's condition on a compact subset and the Lyapunov control with weight $\omega$, Harris theorem says that the stochastic semigroup admits a unique invariant measure and that for some
$\beta \ge 0$, for all initial probability densities, $S_t \rho^0$ converges exponentially fast to the invariant measure in the weighted total variation norm

$$ \|\mu\|_{\text{TV}(\beta\omega)} = \int_{[0,1) \times \mathbb{R}_+} (1 + \beta\omega) d|\mu|. $$

The domain being bounded in the $v$ direction, [87] constructs the Lyapunov weight $\omega(v, g) = \omega(g)$ as a function of $g$ only and proves that $\omega(g) = (g-1)^2$ satisfies the Lyapunov condition. For Doeblin's condition, in some area $g \le R$, they exploit the interaction between the deterministic dynamical system, sending the particles towards the curve where $J_v(v, g) = 0$, and the Brownian noise making them move away. A careful construction yields a measure $\nu$ satisfying the minoration condition.

**Theorem 7.5** (Dou, Kong, Xu, Zhou [87]). Let $\omega(g) = (g-1)^2$. Assume the linear case $S_E = 0$. There exists a unique invariant measure $\rho_\infty$ and constants $C \ge 1, \lambda > 0, \beta > 0$ such that for all $\rho^0 \in \mathcal{P}([0, 1) \times \mathbb{R}_+)$, there exists a unique solution to (7.5) associated for all $t > 0$ to a density $\rho_t$ which is a weak solution to (V-C) and satisfies

$$ \forall t \ge 0, \quad \|\rho_t - \rho_\infty\|_{\text{TV}(\beta\omega)} \le C e^{-\lambda t} \|\rho_0 - \rho_\infty\|_{\text{TV}(\beta\omega)}. $$

## 7.4 A priori estimates on solutions in the nonlinear case

### 7.4.1 Existence of nonlinear stationary state

Given the existence of a linear stationary state, it is possible to investigate the existence of nonlinear stationary state by a fixed point on the stationary firing rate $\mathcal{N}_\infty$. Consider the functional $\Psi: \mathbb{R}_+ \to \mathbb{R}_+$ defined by

$$ \mathcal{N} \mapsto \Psi(\mathcal{N}) = V_E \int_0^{+\infty} \rho_\infty^{\mathcal{N}}(0, g)dg, $$

where $\rho_\infty^{\mathcal{N}}$ is the unique stationary state from the linear system with $\mathcal{N}$ as a constant external input. It is continuous and $\Psi(0) > 0$. Hence, to prove the existence of a fixed point (and hence of a nonlinear stationary firing rate yielding a stationary state) it is enough to find $\mathcal{N}_1$ such that $\Psi(\mathcal{N}_1) < \mathcal{N}_1$. To do so, [211] proves the bounds

$$ -V_F + (V_E - V_F) \left( g_\infty + \sqrt{\frac{a_\infty}{2\pi}} e^{-\frac{g_\infty^2}{2a_\infty}} \right) \le \Psi(\mathcal{N}) \le V_E \left( g_\infty + \sqrt{\frac{2a_\infty}{\pi}} e^{-\frac{g_\infty^2}{2a_\infty}} \right), $$

with $g_\infty = s_0 + S_E \mathcal{N}$, $a_\infty = a_0 + S_E^2 \mathcal{N}$. Analysing the explicit $\mathcal{N}$ dependency in the bounds allows to prove

**Theorem 7.6** (Perthame, Salort [211]). If $\frac{V_E}{V_F} S_E < 1$, there exists at least one stationary state to (V-C). However, if

$$ \frac{V_E - V_F}{V_F} S_E > 1, \quad \text{and} \quad (V_E - V_F)s_0 > V_F^2, $$

then there is no stationary state to (V-C).

The $L^\infty$ regularity from Theorem 7.1 translates to these nonlinear stationary states, as they are linear stationary states with a modified $\tilde{s}_0 = s_0 + S_E \mathcal{N}_\infty$.

### 7.4.2 Propagation of moments and global bounds on the firing rate

For any $k \in \mathbb{N}$, define the moments

$$ \Psi(t) = \int_0^{V_F} \int_0^{+\infty} v \rho(v, g, t) dg dv, $$
$$H_k(t) = \int_0^{V_F} \int_0^{+\infty} g^k \rho(v, g, t) dg dv \quad \text{and} \quad f(t) = \int_0^{V_F} \rho(v, 0, t) dv.$$

Then, providing these quantities exist, multiplying (7.1) by the appropriate weight and integrating, yields the moment equations

$$\sigma_E \frac{dH_1}{dt}(t) = -H_1(t) + s_0 + S_E \mathcal{N}(t) + a(t)f(t),$$

$$\sigma_E \frac{dH_2}{dt}(t) = -2H_2(t) + 2[s_0 + S_E \mathcal{N}(t)] H_1(t) + 2a(t),$$

$$\sigma_E \frac{dH_k}{dt}(t) = -kH_k(t) + k[s_0 + S_E \mathcal{N}(t)]h_{k-1}(t) + k(k-1)a(t)h_{k-2}(t), \quad k \ge 3$$

$$\frac{d\Psi}{dt}(t) = -g_L \Psi(t) - \int_0^{V_F} \int_0^{+\infty} gv\rho(v, g, t)dgdv + V_E H_1(t) - V_F \mathcal{N}(t).$$

Note that since $\rho$ is a probability density, $H^1(t)^2 \le H_2(t)$. Define

$$\omega = \frac{1}{\sigma_E} \left( \frac{V_E}{V_F} [S_E + a_1 S_E^2] - 1 \right).$$

From this system of differential equations for the moments and using the signs, [211] proves that if the initial moments are finite, it is propagated in time for any solution of the nonlinear problem (V-C).

**Theorem 7.7** (Perthame, Salort [211]). Let $\rho$ be a solution to (V-C). Assume that $H_k(0) < +\infty$ for all $k \in \mathbb{N}$. Then there exists $C > 0$ such that

$$H_1(t) \le C \max(1, e^{\omega t}), \quad \int_0^t \mathcal{N}(s)ds \le C \max(1, e^{\omega t}),$$

and for all $T > 0$ there exists $\tilde{C}(T) > 0$ such that for all $t \in [0, T]$,

$$H_k(t) \le \tilde{C}(T), \quad \int_0^T \int_0^{+\infty} (1+g)^{k-1} N(g, t) dg dt \le \tilde{C}(T).$$

This allows to prove that the firing rate is unbounded in the high connectivity regime. Using the moment equations, we can write

$$\frac{d\Psi}{dt}(t) \ge (V_E - V_F)H^1(t) - g_L\Psi(t) - V_F\mathcal{N}, \quad \frac{dH^1}{dt}(t) \ge \frac{1}{\sigma_E}(s_0 + S_E\mathcal{N}(t) - H^1(t)).$$

Multiplying the first inequality by $\frac{S_E}{V_F \sigma_E}$ and merging the constants into $\bar{C}, \zeta > 0$, it leads to

$$\frac{d}{dt} \left( \frac{S_E}{V_F \sigma_E} \Psi + H^1 \right) \ge \zeta \left( \frac{S_E}{V_F \sigma_E} \Psi + H^1 \right) - \bar{C}.$$

For large enough $H^1(0)$ it implies a uniform lower bound on $\frac{S_E}{V_F \sigma_E} \Psi + H^1$, which in turn yields exponential growth of $\mathcal{N}$ by re-injecting the lower bound in the differential inequality for $\Psi$.

**Theorem 7.8** (Perthame, Salort [211]). Assume that

$$\zeta := \frac{1}{\sigma} \left( \frac{S_E(V_E - V_F)}{V_F} - 1 \right) > 0.$$

Then, if $H_1(0)$ is large enough, there exists $C > 0$ such that for any solution to (V-C),

$$\int_0^t \mathcal{N}(s)ds \ge Ce^{\zeta t}.$$
Given that several over neuroscience models like (NNLIF), (jLIF) or (6.3) exhibit finite time blow-up at the level of the firing rate, that reciprocally some of them are proved to have global-in-time solutions as long as the firing rate stays in $L_{\text{loc}}^{\infty}$ [48, 31, 46], a natural question is whether or not the firing rate $N(t)$ of (V-C) is locally uniformly bounded. The article [211] provides the following partial answer.

**Theorem 7.9** (Perthame, Salort [211]). Let $q \ge 2$, and $\varepsilon > 0$. If

$$ \int_0^{V_F} \int_0^{+\infty} (1+g)^{2q+\varepsilon-1} \rho^0(v,g)^q dv dg < +\infty, \quad \int_0^{V_F} \int_0^{+\infty} (1+g)^2 \rho^0(v,g) dv dg < +\infty, $$

then for all $T \in \mathbb{R}_+$, any solution to (V-C) satisfies

$$ \int_0^T \mathcal{N}(t)^q dt < +\infty. $$

This is a $L_{\text{loc}}^q$ bound for $q \ge 2$ and not a $L_{\text{loc}}^{\infty}$ bound, but it should be enough to prove that blow-up cannot happen, provided a suitable well-posedness framework. Note that [233] has proved that an $L^q$ bound, $q > 2$, on the firing rate $N(t)$ is enough$^8$ to ensure global existence of solutions to (NNLIF).

## 7.5 Existence theory in the weak interaction regime

As we have seen, [87] has provided a framework for existence of weak solutions to (V-C) in the linear case by relying upon a probabilistic setting for the corresponding system of SDEs. Until very recently, there was no result about existence of time-dependent solutions for the nonlinear case. In [111], a partial answer is proposed: weak solutions are constructed in the case of weak interactions (small enough $S_E$).

### 7.5.1 Ultracontractivity estimates

One way to handle the structure of (V-C) is by using the hypocoercivity structure of the kinetic Fokker-Planck equation [267, 7]. However, the article [111] goes another way: following [41, 238], they prove ultracontractivity estimates and rely upon the Doeblin-Harris method (see Subsections 4.6.3 and 7.3.3), in order to construct a framework for well-posedness.

Consider, for external source terms $\alpha, \beta \in L^{\infty}(0, +\infty)$, the linear kinetic Voltage-Conductance equation

$$ \frac{\partial \rho}{\partial t} = \mathcal{L}_{\alpha, \beta} \rho, \quad \mathcal{L}_{\alpha, \beta} \rho = \frac{\partial}{\partial v} [J_v(v, g)\rho] + \frac{\partial}{\partial g} [(\beta(t) - g)\rho] + \alpha(t) \frac{\partial^2 \rho}{\partial g^2}, \quad (7.6) $$

with boundary conditions as in (V-C). Taking either $\omega(v, g) = (1+g^2)^{\frac{k}{2}}$, $k \in \mathbb{N}$, $k > 1$ or $\omega(v, g) = e^{\alpha g}$, $\alpha > 0$ as a weigth function, [111] formally integrates the equation and proves that there is another equivalent weigth $\tilde{\omega}$ and a constant $\kappa \ge 0$ depending only on $\omega, \alpha, \beta$ satisfying the a priori estimate:

$$ \|\rho(\cdot, \cdot, t)\|_{L^p(\tilde{\omega})} \le e^{\kappa t} \|\rho^0\|_{L^p(\tilde{\omega})}. \quad (7.7) $$

Then, using the penalisation of the $v$-boundary $\gamma(v) = \min(v, V_F - v)$ and building step by step from the proof of (7.7), they prove that for the exponential weigth $\omega(v, g) = e^{\alpha g}$, $\alpha > 0$, solutions to the linear equation (7.6) satisfy

$$ \int_0^T \phi^2 \int_0^{V_F} \int_0^{+\infty} \left[ \frac{1}{\delta^{\frac{1}{27}}} + \sqrt{1+g^2} \right] \rho \omega \, dv \, dg \, dt \le C \int_0^T (\phi \phi'_+ + \phi^2) \int_0^{V_F} \int_0^{+\infty} \rho^2 \omega^2 \, dt \, dv \, dg, \quad (7.8) $$

with $C \ge 1$ a constant and any test function $\phi \in C_c^1((0, T))$.

<footnote>$^8$The proof in [233] is done for $q = 3$ for the sake of clarity, but the argument works for any $q > 2$, and does not work for $q = 2$ or lower.</footnote>
From there, [111] considers $\bar{\rho} = \rho\phi(t)\chi(v)\omega(g)$, $\phi \in C_c^1((0, T))$, $\chi \in C_c^1((0, V_F))$ and performs a regularisation argument [267, Appendix] in order to get a bound of $\int_0^T \|\bar{\rho}\|_{L^2}^2 dt$ in function of $\int_0^T \|\rho\|_{L^2(\omega)}^2 dt$ and the regularisation functions. Combining with (7.8), they conclude that there exists $r > 2$ such that for any $\phi \in C^1((0, T)) \cap C_{\to 0}((0, T))$, $T \in (0, 1)$,

$$ \int_0^T \phi^2 \|\rho(\cdot, \cdot, t)\|_{L^r(\omega)} dt \le \tilde{C} \|(\phi + |\phi'|)\rho\omega\|_{L^2([0, V_F] \times \mathbb{R}_+ \times \mathbb{R}_+)} $$

Using Nash's argument [195], it leads to the bound

$$ \forall t \in (0, 1), \quad \|\rho(\cdot, \cdot, t)\|_{L^2(\omega)} \le \frac{1}{t^{\nu_1}} \|\rho^0\|_{L^1(\omega)}, $$

for some $\nu_1 > 0$. Last, they build from this result using the Nash-Moser iteration procedure and estimate (7.7) in order to get a full ultracontractivity estimate for the linear system (7.6).

**Theorem 7.10** (Fonte Sanchez, Mischler [111]). Assume $a_1 = 1$. Let $\omega(v, g) = e^{\alpha g}$, $\alpha > 0$ and $\rho^0 \in L^2(\omega)$. Then there exists $C \ge 1$, $\nu > 0$ and $\kappa \in [0, C(1 + \|\alpha\|_{\infty} + \|\beta\|_{\infty})]$ such that any solution $\rho$ to (7.6) satisfies the ultracontractivity estimate

$$ \|\rho(\cdot, \cdot, t)\|_{L^{\infty}(\omega)} \le C \frac{e^{\kappa(t-s)}}{(t-s)^{\nu}} \|\rho(\cdot, \cdot, s)\|_{L^1(\omega)}. $$

### 7.5.2 Well-posedness in a weighted $L^\infty$ framework

Adapting classical methods from kinetic theory, [111] proves that the linear system (7.6) with external inputs $\alpha, \beta$ admits for all $T > 0$ a unique weak-solution $\rho \in C([0, T], L^2(\omega)) \cap L^2((0, T) \times (0, V_F); H^1(0, +\infty))$ and that this solution admits a trace $\gamma(\rho) \in L^2(\Gamma_{v,t}, \omega^2(1+g^2)^{-\frac{1}{2}}J_v^2)$ on the $(v, t)$-boundary $\Gamma_{v,t} = (0, T) \times (\{0\} \cup \{V_F\}) \times \mathbb{R}_+$. Using a kinetic Harnack inequality in the spirit of [123], those solutions to the linear problem are also proved to enjoy a higher $C^\eta$ regularity in the strict interior of the time-voltage-conductance domain, for some $\eta > 0$.

All this work allows to prove that for small enough $S_E$ with respect to $\rho^0$, the mapping that takes an external firing rate $\mathcal{N}_1$ defining $\alpha, \beta$ and returns the firing rate $\mathcal{N}_2$ defined from the trace on the $(v, t)$-boundary of $\rho$ solution to (7.6) enjoys enough regularity to apply to Schauder-Tychonov fixed-point theorem. This provides a solution to the nonlinear problem.

**Theorem 7.11** (Fonte Sanchez, Mischler [111]). Assume $a_1 = 1$. Let $\omega(v, g) = (1+g^2)^{\frac{k}{2}}$, $k \in \mathbb{N}$, $k > 1$ or $\omega(v, g) = e^{\alpha g}$, $\alpha > 0$. There exists $\eta(\omega) > 0$ such that for all initial condition $\rho^0 \in L^\infty(\omega)$, if $(S_E + S_E^2) \|\rho^0\|_{L^\infty(\omega)} < \eta$, there exists at least one solution

$$ \rho \in L^\infty(\mathbb{R}_+, L^\infty(\omega)) \cap C(\mathbb{R}_+, L^2(\omega)) $$

to (V-C) which conserves mass and positivity.

Last, using a reformulation of the Doeblin-Harris method [128] in a Banach lattice, and taking advantage of the ultracontractivity estimate in Theorem 7.10, they prove the best result known so far on the existence and regularity of nonlinear stationary states.

**Theorem 7.12** (Fonte Sanchez, Mischler [111]). Assume $a_1 = 1$. Let $\omega_0(v, g) = (1+g^2)^{\frac{k}{2}}$, $k \in \mathbb{N}$, $k > 1$. There exists $\eta(\omega_0) > 0$ such that if $S_E < \eta(\omega_0) \|\omega_0^{-1}\|_{L^1}$, for either $\omega = \omega_0$ or $\omega(v, g) = e^{\alpha g}$, $\alpha > 0$, there exists a unique stationary state

$$ \rho_{\infty} \in L^{\infty}(\omega) \cap C((0, V_F) \times (0, +\infty)). $$

It is linearly exponentially asymptotically stable in $L^p(\omega)$ for all $p \in [1, +\infty]$.

Like in many other results obtained from Doeblin-Harris approaches, the constants are all constructive in the linear exponential contraction estimates.
## 7.6 Numerical simulations of the Voltage-Conductance equation

Numerical methods have also been introduced in order to explore the different features of the voltage-conductance kinetic equations (V-C). In [28] a numerical method based on Chang-Cooper schemes as in [20], see subsection 1.12, was developed and validated against Monte Carlo simulations. These numerical methods were used to showcase the appearance of bi-stability in this model showing the appearance of two steady states with low and high activity respectively with the appearance of a phase transition leading to complicated dynamics. They also observe the existence of periodic solutions or damped oscillations in the firing rate depending on the values of $S_E$. It is an open problem to show the existence of periodic solutions in this model even if numerically observed by multiple authors using different numerical approaches. We refer to [28] for more details on the numerical experiments and particular values in the parameter regime leading to periodic solutions. Finally, the authors look at the validity limits of moment closures of the kinetic-conductance model introduced in [223]. This moment-closure system reads as

$$ 
\begin{aligned}
\partial_t \rho_v(t, v) &= \partial_v \left\{ \left[ \left( \frac{v - V_R}{\tau} \right) + \mu_1(t, v) \left( \frac{v - V_E}{\tau} \right) \right] \rho_v(t, v) \right\} \\
\partial_t \mu_1(t, v) &= -\frac{1}{\sigma_E} [\mu_1(t, v) - s_0 - S_E N(t)] + \frac{a(N(t))}{\sigma_E \rho_v(t, v)} \partial_v \left[ \left( \frac{v - V_E}{\tau} \right) \rho_v(t, v) \right] \\
&\quad + \left[ \left( \frac{v - V_R}{\tau} \right) + \mu_1(t, v) \left( \frac{v - V_E}{\tau} \right) \right] \partial_v \mu_1(t, v).
\end{aligned}
$$

where $\rho_v(t, v)$ and $\rho_v(t, v)\mu_1(t, v)$ are the first two moments in $g$ of the kinetic distribution $\rho(t, v, g)$. This closure assumption was numerically verified to be valid in steady states away from the boundaries in the voltage variable.

## 7.7 The random firing voltage-conductance model and one variable reductions

Because of the difficulty of the complex v-boundary condition and the singular firing at $V_F$, a strategy to simplify the study of the model (V-C) is to consider a softer firing mechanism [212, 146]. Instead of letting the neurons fire as they cross the threshold potential $V_F$, they can be allowed to go into the space $[V_F, V_E)$, with again $V_E$ the excitatory reversal potential where the flux $J_v$ is always inward in the domain. The neurons with voltage $V_F \le v < V_E$ will then fire randomly with high probability $\phi_F(v)$, which translates into a removal term $\phi_F(v)\rho(v, g, t)$ in the deterministic framework. The firing rate is then the quantity of neurons over all the $(v, g)$ domain which are removed from the equation:

$$ \mathcal{N}(t) = \int_{0}^{+\infty} N(g, t)dg = \int_{0}^{+\infty} \int_{0}^{V_E} \Phi_F(v)\rho(v, g, t)dvdg. $$

Note that this approach has also been implemented into the model NNLIF [29, 90] (see also Subsection 1.11.2). The model writes

$$ 
\left\{
\begin{array}{l}
\displaystyle \frac{\partial \rho}{\partial t} + \frac{\partial}{\partial v} \left[ \left( -g_L v + g(V_E - v) \right) \rho \right] + \frac{\partial}{\partial g} \left[ \frac{s_0 + S_E \mathcal{N}(t) - g}{\sigma_E} \rho \right] - \frac{a}{\sigma_E} \frac{\partial^2 \rho}{\partial v^2} + \phi_F \rho = 0, \\
\displaystyle \mathcal{N}(t) = \int_0^{+\infty} N(g, t)dg, \quad N(g, t) = \int_0^{V_E} \phi_F(v)\rho(v, g, t)dv, \\
\displaystyle J_v(0, g)\rho(0, g, t) = N(g, t), \quad \rho(V_E, g, t) = 0, \quad J_g(0, \mathcal{N}(t))\rho(v, 0, t) = \frac{a}{\sigma_E} \frac{\partial \rho}{\partial g}(v, 0, t), \\
\displaystyle \rho(v, g, 0) = \rho^0(v, g) \ge 0, \quad \int_0^{V_E} \int_0^{+\infty} \rho^0(v, g)dgdv = 1,
\end{array}
\right.
\quad (7.11) $$

where we simplify $a(\mathcal{N}) \equiv a > 0$ and still denote

$$ J_v(v, g) = -g_L v + g(V_E - v), \quad J_g(g, \mathcal{N}) = \frac{s_0 + S_E \mathcal{N}(t) - g}{\sigma_E}. $$
Note that we do not have anymore the complex switch at $g = g_F$ on the boundary $v = V_F$: it is replaced by a homogeneous Dirichlet boundary condition at $v = V_E$. The firing function is chosen as $\phi_F \in W^{1,+\infty}(0, V_E)$ and such that $\phi_F \ge 0$, $\phi'_F \ge 0$.

### 7.7.1 Fast conductance limit

In [212], the Voltage-Conductance model with random firing (7.11) is used to derive rigorously a voltage-only model through the limit of fast conductance, which corresponds in our setting to the limit of vanishing conductance time constant $\sigma_E \to 0$.

Let us denote $\rho^\epsilon$ a solution obtained for $\sigma_E = \epsilon > 0$. The main idea is that when $\epsilon$ goes to 0, in the formal limit $\rho^\epsilon$ converges to a solution of

$$ \frac{\partial}{\partial g} \left[ (s_0 + S_E \mathcal{N}(t) - g) \rho \right] - a \frac{\partial \rho}{\partial g^2} = 0, \quad (s_0 + S_E \mathcal{N}(t)) \rho(v, 0, t) - a \frac{\partial \rho}{\partial g}(v, 0, t) = 0. $$

The solution can then be written in terms of a v-dependent part and a normalised Maxwellian

$$ \rho(v, g, t) = n(v, t)M(g, \mathcal{N}(t)), \quad M(v, \mathcal{N}(t)) = \frac{e^{-\frac{(s_0 + S_E \mathcal{N}(t) - g)^2}{2a}}}{\int_0^{+\infty} e^{-\frac{(s_0 + S_E \mathcal{N}(t) - g)^2}{2a}} dg}. $$

The quantity $n(v, t)$ should then be the limit of

$$ n^{\epsilon}(v, t) = \int_{0}^{+\infty} \rho^{\epsilon}(v, g, t) dg. $$

Integrating (7.11), we find that this quantity $n^\epsilon$ is solution to

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial n^{\epsilon}}{\partial t} + \frac{\partial}{\partial v} \left[ \left( -g_L v + (V_E - v) \int_0^{+\infty} g \frac{\rho^{\epsilon}(v, g, t)}{n^{\epsilon}(v, t)} dg \right) n^{\epsilon}(v, t) \right] + \phi_F n^{\epsilon}(v, t) = 0, \\[1.5em] \displaystyle \mathcal{N}(t) = \int_0^{V_E} \phi_F(v) n^{\epsilon}(v, t) dv, \\[1.5em] \displaystyle V_E \int_0^{+\infty} g \frac{\rho^{\epsilon}(v, g, t)}{n^{\epsilon}(v, t)} dg n^{\epsilon}(0, t) = \mathcal{N}(t), \quad n^{\epsilon}(V_E, t) = 0. \end{array} \right. \qquad (7.12) $$

The idea of [212] for performing a rigorous slow-fast limit in (7.12) is to rewrite the v-drift as

$$ -g_L v + (V_E - v) \int_0^{+\infty} g \frac{\rho^{\epsilon}(v, g, t)}{n^{\epsilon}(v, t)} dg = G_{\epsilon}(v, t)(V_{\epsilon}(v, t) - v), $$

with

$$ G_{\epsilon} = \int_{0}^{+\infty} (g_{L} + g) \frac{\rho^{\epsilon}(v, g, t)}{n^{\epsilon}(v, t)} dg dg, \quad V_{\epsilon} = \frac{1}{G_{\epsilon}(v, t)} \int_{0}^{+\infty} g V_{E} \frac{\rho^{\epsilon}(v, g, t)}{n^{\epsilon}(v, t)} dg. $$

The rewriting $G_\epsilon(V_\epsilon - v)$ can then be seen as a conductance acting on the difference between the current potential and an equilibrium potential $g_{eq}(t)(v_{eq}(t) - v)$, and the idea is to take advantage of the inequalities $0 < V_\epsilon \le V_E$ and $G_\epsilon \ge g_L$.

The quotient $\rho^\epsilon/n$ should then converge towards $M(v, \mathcal{N}(t))$. In the limit $\epsilon \to 0$, [212] proves that

$$ \mathcal{N}_{\epsilon} \to N(t), \quad G_{\epsilon} \to G := g_L + \int_{0}^{+\infty} g M(g, N(t)) dg, \quad V_{\epsilon} \to \frac{1}{G} V_E \int_{0}^{+\infty} g M(g, N(t)) dg. $$