# Well-posedness and numerical analysis of an elapsed time model with strongly
coupled neural networks

Mauricio Sepúlveda, Nicolas Torres  and Luis Miguel Villada  CI2MA and DIM,
Universidad de Concepción, Concepción-Chile. E-mail: mauricio@ing-
mat.udec.clCorresponding author. Departamento de Matemática Aplicada,
Universidad de Granada-España. E-mail<EMAIL_ADDRESS>de
Matemática, Universidad del Bío-Bío and CI2MA-Universidad de Concepción,
Concepción-Chile. E-mail<EMAIL_ADDRESS>
###### Abstract

The elapsed time equation is an age-structured model that describes dynamics
of interconnected spiking neurons through the elapsed time since the last
discharge, leading to many interesting questions on the evolution of the
system from a mathematical and biological point of view. In this work, we
first deal with the case when transmission after a spike is instantaneous and
the case when there exists a distributed delay that depends on previous
history of the system, which is a more realistic assumption. Then we study the
well-posedness and the numerical analysis of the elapsed time models. For
existence and uniqueness we improve the previous works by relaxing some
hypothesis on the non-linearity, including the strongly excitatory case, while
for the numerical analysis we prove that the approximation given by the
explicit upwind scheme converges to the solution of the non-linear problem. We
also show some numerical simulations to compare the behavior of the system in
the case of instantaneous transmission with the case of distributed delay
under different parameters, leading to solutions with different asymptotic
profiles.

2010 Mathematics Subject Classification. 35A35, 35F20, 35R09, 65M06.

Keywords and phrases. Structured equations; Mathematical neuroscience; Delay
differential equations, Well-posedness, Numerical analysis; Periodic
solutions.

##  1 Introduction

Structured equations have been widely studied in the modeling of biological
systems. In particular in the context of neuroscience age-structured models
are an interesting approach to modeling the dynamics of interconnected spiking
neurons. The study of the precise mechanisms of brain processes, leading to
synchronous regular or irregular activities, has been a challenge for
mathematician and biologists with many interesting models including discrete
systems, differential equations and stochastic processes (for a reference see
for example [1]).

One of these equations is the well-known elapsed time model, where a given
population of neurons is described by the time elapsed since their last
discharge. In this network neurons are submitted to random discharges, which
is related to changes in the membrane potential, stimulating other neurons to
spike. One the main motivations of this model is to predict the brain activity
through the previous history of spikes and the key element that determines the
evolution of the system is the way neurons interact within the network,
leading to different possible behaviors of the neural activity and patterns
formation. This equation is a mean-field limit of a microscopic model that
establishes a bridge of the dynamics of a single neuron with a population-
based approach, whose aspects has been investigated in [2, 3, 4, 5, 6, 7].

The elapsed time model has been studied by many authors. The pioneer works on
this model were studied by Pakdaman et al. [8, 9, 10] and some important
results on exponential convergence to the equilibrium for weak non-linearities
were proved in [11, 12, 13] through different techniques such as the entropy
method, semi-group theory and spectral arguments. Results on strong non-
linearities have been studied in [14], where existence of periodic solution
with jump discontinuities was established. Moreover, different extensions of
the elapsed time model have been studied by incorporating new elements such as
the fragmentation equation [10], spatial dependence with connectivity kernel
in [15], a multiple-renewal equation in [16] and a leaky memory variable in
[17].

The classical elapsed time model with instantaneous transmission, which we
will call throughout this article as ITM, is given by the following non-linear
age-structured equation

$\textrm{(ITM)}\begin{cases}\partial_{t}n+\partial_{s}n+p(s,N(t))n=0&t>0,\>s>0\\\
N(t)=n(t,s=0)=\int_{0}^{+\infty}p(s,N(t))n(t,s)\mathrm{d}s&t>0\\\
n(0,s)=n^{0}(s)\geq 0&s\geq 0,\\\
\int_{0}^{\infty}n^{0}(s)\mathrm{d}s=1,\end{cases}$ (1)

where $n(t,s)$ is the probability density of finding a neuron at time $t$,
whose elapsed time since last discharge is $s\geq 0$ and the function $N(t)$
represents the flux of discharging neurons.

For this equation, we assume that when a neuron spikes, its interactions with
other neurons are instantaneous so that we assume for simplicity in this model
that the total activity of the network is simply given by the value of $N(t)$.
The crucial non-linearity is given by the function
$p\colon[0,\infty)\times[0,\infty)\mapsto[0,\infty)$, which is called as the
hazard rate and it describes the susceptibility of neurons to spike. We assume
that $p$ depends on the elapsed time $s$ and the activity $N$ and without loss
of generality we consider that $p\in
W^{1,\infty}(\mathbb{R}^{+}\times\mathbb{R}^{+})$, though the regularity is
not a crucial assumption as we will see in the numerical examples. We stick to
these names of each term of the system, though they may differ in the
literature.

In this setting, neurons discharge at the rate given by $p$ and then the
elapsed time is immediately reset to zero, as it is stated by the integral
boundary condition of $n$ at $s=0$. Following the terminology of age-
structured equations, the elapsed time corresponds to ”age” of neuron. When a
neuron discharge, it is considered to ”die” where $p(s,N)n$ is the
corresponding ”death” term. After a neuron spikes, it instantaneously re-enter
the cycle and it is considered to be ”reborn” with the boundary condition at
$s=0$ representing the ”birth” term.

Moreover, we assume that the rate $p$ is increasing with respect to the age
$s$, which means that neurons are more prone to spike when the elapsed time
since last discharge is large. According to the dependence of the rate $p$ on
the total activity different regimes are possible. When $p$ is increasing with
respect to $N$ we say that the network is excitatory, which is means under a
high activity neurons are more susceptible to discharge. Similarly, when $p$
is decreasing we say that the network is inhibitory and we have the opposite
effect on the network. Moreover, if the following conditions holds

$\|\partial_{N}p\|_{\infty}<1,$ (2)

we say that the network is under a weakly interconnected regime, which means
that the non-linearity is weak.

For the initial data $n^{0}\in L^{1}(\mathbb{R}^{+})$ we assume that is a
probability density and we formally have the following mass-conservation
property

$\int_{0}^{\infty}n(t,s)\,\mathrm{d}s=\int n^{0}(s)\,\mathrm{d}s,\qquad\forall
t\geq 0,$ (3)

which will be crucial in the analysis of Equation (1). Throughout this article
we consider solutions in the weak sense but for simplicity we simply refer to
them as solutions.

######  Remark 1.1.

Since we look for solution $n\in\mathcal{C}([0,T],L^{1}(\mathbb{R}^{+}))$ for
some $T>0$, observe that $N(0)$ formally satisfies the following equation

$N(0)=\int_{0}^{\infty}p(s,N(0))n^{0}(s)\,\mathrm{d}s.$ (4)

In the inhibitory case and the weak interconnections regime this equation has
a unique solution for $N(0)$, while for the excitatory case we may have
multiple solutions and thus the solution of Equation (1) is not unique.
Moreover, we remark that in general $N(0)\neq n^{0}(0)$, which might imply
that $n$ has discontinuities along the line $\\{(t,s)\in\mathbb{R}^{2}\colon
t=s\\}$.

Concerning stationary solutions of the non-linear System (1), these are given
by solutions the problem

$\begin{cases}\partial_{s}n+p(s,N)n=0&s>0,\vspace{0.15cm}\\\
N=n(s=0)\coloneqq\int_{0}^{\infty}p(s,N)n(s)\,\mathrm{d}s,&\vspace{0.15cm}\textstyle\\\
\int_{0}^{\infty}n(s)\,\mathrm{d}s=1,\quad n(s)\geq 0.\end{cases}$ (5)

If the activity $N$ is given, we can determine the stationary density through
the formula

$n(s)=Ne^{-\int_{0}^{s}p(u,N)\,\mathrm{d}u}.$ (6)

Thus by integrating with respect to $s$, we get that $(n,N)$ corresponds to a
stationary solution of System (1) if the activity satisfies the fixed point
equation

$N=F(N)\coloneqq\left(\int_{0}^{\infty}e^{-\int_{0}^{s}p(u,N)\,\mathrm{d}u}\,\mathrm{d}s\right)^{-1}.$
(7)

So that depending on the rate $p$, we have a unique solution in the inhibitory
case and in the weak interconnections regime. For the excitatory case we may
have multiple steady-states.

We also remark that a prototypical form of the function $p$ is given by

$p(s,N)=\varphi(N)\chi_{\\{s>\sigma\\}},$ (8)

which represents a hazard rate with an absolute refractory period $\sigma>0$
so that for an age $s<\sigma$ neurons are not susceptible to discharge. For
$s>\sigma$ neurons are able to discharge and the density $n$ decays
exponentially according to $\varphi(N)$. The function $\varphi$ is assumed to
be smooth and satisfies the following bounds

$p_{0}\leq\varphi(N)\leq p_{1}\qquad\forall N\geq 0,$

for some constants $p_{0},p_{1}>0$. This special case has been studied which
has been studied in [14, 8], where they proved convergence for the inhibitory
case and constructed periodic solutions for the excitatory case. Another
possible example is to consider a variable refractory period depending on the
total activity

$p(s,N)=\chi_{\\{s>\sigma(N)\\}}.$ (9)

This type of hazard rate has been studied in [9], where existence of periodic
solutions was studied as well.

From a biophysical point of view, when neurons spike it is reasonable to
consider a delay in the transmission to other neurons. In order to take into
account this effect, Pakdaman et al. [8] considered a modification of the
elapsed time model incorporating a distributed delay, which corresponds to the
following variant of Equation (1) that we will call as DDM

$\textrm{(DDM)}\begin{cases}\partial_{t}n+\partial_{s}n+p(s,X(t))n=0&t>0,\>s>0\\\
N(t)=n(t,s=0)=\int_{0}^{+\infty}p(s,X(t))n(t,s)\mathrm{d}s&t>0\\\
X(t)=\int_{0}^{t}\alpha(t-\tau)N(\tau)\,\mathrm{d}\tau&t>0\\\
n(0,s)=n^{0}(s)\geq 0&s\geq 0,\\\
\int_{0}^{\infty}n^{0}(s)\mathrm{d}s=1.&\end{cases}$ (10)

The kernel $\alpha\in L^{1}(\mathbb{R}^{+})$ with $\alpha\geq 0$, corresponds
to the distributed delay and for simplicity we may assume that $\alpha$ is
smooth and uniformly bounded with
$\int_{0}^{\infty}\alpha(\tau)\mathrm{d}\tau=1$, but the theoretical results
are still valid if we only assume that $\alpha$ is integrable. For this model
$N(t)$ is the flux of discharging neurons and $X(t)$ is the total activity,
which depends on the values taken by $N$ in the past (i.e. in the interval
$[0,t]$) through the convolution with $\alpha$. Unlike ITM, the rate $p$
depends on the total activity $X(t)$ instead of the discharging flux.
Properties like the mass-conservation remain valid for this modified model. We
remark that under the condition
$\int_{0}^{\infty}\alpha(\tau)\mathrm{d}\tau=1$, we consider as steady-states
of DDM equation (10) the same as those of ITM equation.

In particular when $\alpha(t)$ approaches in the sense of distributions to the
Dirac’s mass $\delta(t)$, then we formally get $X(t)=N(t)$ and thus we recover
the ITM equation (1). An important example studied in [8] is the exponential
delay given by $\alpha(t)=\frac{1}{\lambda}e^{-t/\lambda}$ so that $X(t)$
satisfies the following differential equation

$\begin{cases}\lambda X^{\prime}(t)+X(t)=N(t),\\\ X(0)=0.\end{cases}$ (11)

giving a simple way to compute numerical solution of this system.

Similarly, when $\alpha(t)$ approaches in the sense of distributions to the
Dirac’s mass $\delta(t-d)$, then we formally get $X(t)=N(t-d)$ and we recover
a version of the classical elapsed time equation with a single discrete delay
$d$

$\begin{cases}\partial_{t}n+\partial_{s}n+p(s,N(t-d))n=0&t>0,\>s>0\\\
N(t)=n(t,s=0)=\int_{0}^{+\infty}p(s,N(t-d))n(t,s)\mathrm{d}s&t>0\\\
n(0,s)=n^{0}(s)\geq 0&s\geq 0,\\\
\int_{0}^{\infty}n^{0}(s)\mathrm{d}s=1.\end{cases}$ (12)

Given the flexibility of the DDM model through the kernel $\alpha$, in this
article we will focus on the ITM and DDM equations, but the results and
techniques are presented in this work also valid for Equation (12).

This article is devoted to two aspects of the elapsed time model: well-
posedness and numerical analysis. On one hand, we aim to give a
straightforward proof on well-posedness for both ITM and DDM equations,
improving the proofs given in [8, 11]. In the work of Pakdaman et al. [8] they
proved well-posedness when the rate $p$ is of the form given by (8) or (9)
under some asymptotic conditions in the growth of $p$ with respect to
discharging flux $N$ (or the total activity $X$ in the DDM equation), while in
the work of Cañizo et al. [11] they mainly focused in the weak
interconnections regime for the ITM equation (1). In this context, we give a
proof for a wider class of hazard rates $p$ with some simple and general
assumptions.

On the other hand, we aim to make a numerical analysis of these non-linear
models by proving the convergence of an explicit upwind scheme of first-order.
Previous works on the numerical analysis of age-structured equations includes
the works [18, 19] and further generalizations to solutions in the space of
positive regular measures $\mathcal{M}^{+}(\mathbb{R}^{+})$ has been
investigated by [20, 21, 22] through the particle method. This method consists
that measure solutions are approximated by a sum Dirac’s masses and then
transported according the structured equation, so that by compactness through
tightness of measures they proved that the approximation actually convergences
to a solution of the non-linear equation. Following the spirit of these ideas,
we study the evolution of the finite-volume approximation and then by an
estimate of the bounded variation norm, we get the necessary compactness to
conclude the convergence of the numerical method. This BV-estimate to prove
correctness of the scheme is a novelty that was missing in the literature.

The article is organized as follows. In Section 2, we study the ITM equation
by giving a proof on well-posedness in the inhibitory case (including the weak
interconnections case as well) and explaining how the arguments can be
extended for the general excitatory case. Then we proceed to explain the
scheme to solve numerically ITM equation 1 and prove the necessary estimates
that ensure the convergence of the numerical method. In Section 3, we make the
analogous analysis for the DDM equation by adapting the ideas of the arguments
applied to the ITM equation. Finally in Section 4, we present numerical
simulations to compare both ITM and DDM equations under different choice of
parameters, including the inhibitory and excitatory regime. In particular, we
consider different types of the delay kernel $\alpha$ in order to observe the
limit cases when $\alpha(t)$ is approaching a Dirac’s mass and the possible
asymptotic behavior thereof. This extends the numerical simulations made by
Pakdaman et al. [8], where they consider mainly the exponential kernel in the
DDM equation (10).

##  2 Instantaneous Transmission Model (ITM)

###  2.1 Well-posedness of the ITM

In this subsection we prove that the solution of Equation (1) is well-posed in
the inhibitory case and the weak interconnections regime. We improve the ideas
of Pakdaman et al. [8] by giving consider more general forms for the rate $p$
and we improve the result of the weak interconnection regime of [11] by
extending the existence and uniqueness of a solution when we drop the absolute
value in the Condition (2). The main idea of the proof is to propose the
appropriate fixed point problem that eventually leads to a solution of ITM
equation 1 through the contraction principle.

######  Theorem 2.1.

Consider a non-negative $n^{0}\in L^{1}(\mathbb{R}^{+})$. Assume that $p\in
W^{1,\infty}(\mathbb{R}^{+}\times\mathbb{R}^{+})$ and let
$\gamma\coloneqq\sup_{s,N}\partial_{N}p(s,N)$ with $\gamma\|n^{0}\|_{1}<1$,
then Equation (1) has a unique solution
$n\in\mathcal{C}\left([0,\infty),L^{1}(\mathbb{R}^{+})\right)$ and
$N\in\mathcal{C}[0,\infty)$. Moreover $n$ satisfies the mass-conservation
property (3).

######  Remark 2.1.

We remark that the regularity of the rate $p$ is not fundamental for the proof
and Theorem 2.1 is still valid for a wider class of functions such as

$p(s,N)=\varphi(N)\chi_{\\{s>\sigma(N)\\}},$

with $\varphi$ and $\sigma$ Lipschitz bounded functions, as they are studied
in [8, 14]. So under a similar conditions for the inhibitory and weakly
excitatory regimes, we can apply the arguments used in Theorem 2.1 to get
well-posedness in these cases.

Moreover, from Remark 1.1 we know that in the excitatory case multiple
solutions may arise [14] and the proof of the theorem can be replicated to
prove the existence of solutions with
$n\in\mathcal{C}\left([0,T],L^{1}(\mathbb{R}^{+})\right)$ for some $T>0$.
Indeed, from the proof Lemma 2.2 we can apply the implicit function theorem as
long as the following invertibility condition holds

$\Psi(N(t),n(t,\cdot)):=1-\int_{0}^{\infty}\partial_{N}p(s,N(t))n(t,s)\mathrm{d}s\neq
0\qquad\forall t\in[0,T],$ (13)

where $\Psi:\mathbb{R}^{+}\times L^{1}(\mathbb{R}^{+})\longmapsto\mathbb{R}$,
so that we obtain existence of a solution (or possible branches of solutions)
of Equation (1) defined locally in time by applying the arguments in the proof
of Theorem (2.1). If $\Psi(N(t^{*}),n(t^{*},\cdot))=0$ for some $t^{*}>0$,
then the continuity of solutions is not ensured and jump discontinuities might
arise. We explore this aspect in the section of numerical simulations.

For the proof we need the following lemmas, which will be the key idea
throughout this article. We start with following result on the linear case.

######  Lemma 2.1.

Assume that $n^{0}\in L^{1}(\mathbb{R}^{+})$ and $p\in
L^{\infty}(\mathbb{R}^{+}\times\mathbb{R}^{+})$. Then for a given
$N\in\mathcal{C}[0,T]$, the linear equation

$\begin{cases}\partial_{t}n+\partial_{s}n+p(s,N(t))n=0&t>0,\>s>0,\\\
n(t,s=0)=\int_{0}^{+\infty}p(s,N(t))n(t,s)\mathrm{d}s&t>0,\\\
n(0,s)=n^{0}(s)\geq 0&s\geq 0,\end{cases}$ (14)

has a unique weak solution $n\in\mathcal{C}([0,T],L^{1}(\mathbb{R}^{+}))$.
Moreover $n$ is non-negative and verifies the mass conservation property (3)
for $t\in[0,T]$.

###### Proof.

For a proof see the linear theory developed in [23, 8, 11]. ∎

By using the implicit function theorem, we prove the following lemma that is
the keystone in the proof of our theorem.

######  Lemma 2.2.

Consider a non-negative function $n^{0}\in L^{1}(\mathbb{R}^{+})$. Let
$\gamma\coloneqq\sup_{s,N}\partial_{N}p(s,N)$ and assume that
$\gamma\|n^{0}\|_{1}<1$. Then there exists a unique solution for $N$ of the
equation

$N=F(N)\coloneqq\int_{0}^{\infty}p(s,N)n^{0}(s)\,\mathrm{d}s,$ (15)

that we call $N\coloneqq\psi(n^{0})$, where the map $\psi\colon
L^{1}(\mathbb{R}^{+})\mapsto\mathbb{R}$ satisfies the following estimate

$|\psi(n^{1})-\psi(n^{2})|\leq\frac{\|p\|_{\infty}}{1-\gamma\|n^{0}\|_{1}}\int_{0}^{\infty}|n^{1}-n^{2}|(s)\,\mathrm{d}s$
(16)

for non-negative integrable functions $n^{1},n^{2}$ with
$\|n^{1}\|_{1}=\|n^{2}\|_{1}=\|n^{0}\|_{1}$.

###### Proof.

Observe that $F$ is a continuous and bounded with respect to $N$. Indeed, for
all $N$ we have

$0\leq F(N)\leq\|p\|_{\infty}\|n^{0}\|_{1}$

and hence there exists $N\in[0,\|p\|_{\infty}\|n^{0}\|_{1}]$ such that
$N=F(N)$. Moreover the function $g(N)=N-F(N)$ is strictly increasing. Indeed,

$g^{\prime}(N)=1-F^{\prime}(N)=1-\int_{0}^{\infty}\partial_{N}p(s,N)n^{0}(s)\,\mathrm{d}s\geq
1-\gamma\|n^{0}\|_{1}>0,$

and therefore the $N=F(N)$ has a unique solution that we call
$\overline{N}=\psi(\overline{n})$. Consider the set $U$ defined by

$U\coloneqq\\{n\in L^{1}(\mathbb{R}^{+})\colon\gamma\|n\|_{1}<1\\}$

And consider the map $G\colon X\times\mathbb{R}\mapsto\mathbb{R}$ defined by

$G(n,N)=N-\int_{0}^{\infty}p(s,N)n^{0}(s)\,\mathrm{d}s.$

Observe that for $(n,N)\in U\times\mathbb{R}$ with $G(n,N)=0$ we get
$\partial_{N}G(n,N)>0$. By the implicit function theorem, we notice that
$\psi$ it is a differentiable map on $U$. Moreover $D\psi\colon
L^{1}(\mathbb{R})\mapsto\mathbb{R}$ at a point $(n,N)\in U\times[0,\infty)$ is
given by

$D\psi[h]=\frac{\int_{0}^{\infty}p(s,N)h(s)\,\mathrm{d}s}{1-\int_{0}^{\infty}\partial_{N}p(s,N)n^{0}(s)\,\mathrm{d}s},$

and we have the following estimate in the operator norm at the point
$(n^{0},N^{0})$

$\|D\psi\|\leq\frac{\|p\|_{\infty}}{1-\gamma\|n^{0}\|_{1}},$

thus for $n^{1},n^{2}$ with the same norm $n^{0}$, the inequality (16) readily
follows. ∎

With this lemma, we continue the proof of the Theorem 2.1.

###### Proof.

Consider $T>0$ and a fixed non-negative $N\in\mathcal{C}[0,T]$. For this
function $N$ we define $n[N]\in\mathcal{C}([0,T],L^{1}(\mathbb{R}^{+})$ as the
weak solution of the linear equation (14), which can be expressed trough the
method of characteristics

$n(t,s)=n^{0}(s-t)e^{-\int_{0}^{t}p(s-t+t^{\prime},N(t^{\prime}))\mathrm{d}t^{\prime}}\chi_{\\{s>t\\}}+N(t-s)e^{-\int_{0}^{s}p(s^{\prime},N(t-s+s^{\prime}))\mathrm{d}s^{\prime}}\chi_{\\{t>s\\}}.$

From the mass conservation property we know that
$\|n[N](t)\|_{1}=\|n^{0}\|_{1}$.

On the other hand, for a given $n\in\mathcal{C}([0,T],L^{1}(\mathbb{R}^{+}))$,
let
$\psi\colon\mathcal{C}([0,T],L^{1}(\mathbb{R}^{+}))\mapsto\mathcal{C}[0,T]$
the unique solution $X\in\mathcal{C}[0,T]$ of the equation

$X=\int_{0}^{\infty}p(s,X)n(t,s)\,\mathrm{d}s.$

Under this setting we get a solution of the non-linear equation (1) if only if
we find $N\in\mathcal{C}[0,T]$ such that is a solution of the equation

$N=H(N)\coloneqq\psi\left(n[N]\right),$

where $H\colon\mathcal{C}[0,T]\mapsto\mathcal{C}[0,T]$.

We assert that for $T$ small enough the map $H$ is a contraction. For non-
negative functions $N_{1},N_{2}\in\mathcal{C}[0,T]$, we get from the method of
characteristics we get the following inequality

$\int_{0}^{\infty}|n[N_{1}]-n[N_{2}]|(t,s)\,\mathrm{d}s\leq
A_{1}+A_{2}+A_{3},$

where $A_{1},A_{2},A_{3}$ are given by

$\displaystyle A_{1}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}n^{0}(s)\left|e^{-\int_{0}^{t}p(s+t^{\prime},N_{1}(t^{\prime}))\mathrm{d}t^{\prime}}-e^{-\int_{0}^{t}p(s+t^{\prime},N_{2}(t^{\prime}))\mathrm{d}t^{\prime}}\right|\,\mathrm{d}s$
$\displaystyle A_{2}$ $\displaystyle=$
$\displaystyle\int_{0}^{t}|N_{1}-N_{2}|(t-s)e^{-\int_{0}^{s}p(s^{\prime},N_{1}(t-s+s^{\prime}))\mathrm{d}s^{\prime}}\mathrm{d}s$
$\displaystyle A_{3}$ $\displaystyle=$
$\displaystyle\int_{0}^{t}N_{2}(t-s)\left|e^{-\int_{0}^{s}p(s^{\prime},N_{1}(t-s+s^{\prime}))\mathrm{d}s^{\prime}}-e^{-\int_{0}^{s}p(s^{\prime},N_{2}(t-s+s^{\prime}))\mathrm{d}s^{\prime}}\right|\mathrm{d}s.$

We proceed by estimating each term. For simplicity we write
$n_{1}=n[N_{1}],n_{2}=n[N_{2}]$, so that for $A_{1}$ we have

$A_{1}\leq\int_{0}^{\infty}n^{0}(s)\int_{0}^{t}|p(s+t^{\prime},N_{1}(t^{\prime}))-p(s+t^{\prime},N_{2}(t^{\prime}))|\,\mathrm{d}t^{\prime}\,\mathrm{d}s\leq
T\|n^{0}\|_{1}\|\partial_{N}p\|_{\infty}\|N_{1}-N_{2}\|_{\infty},$

while for $A_{2}$ we have

$A_{2}\leq\int_{0}^{t}|N_{1}-N_{2}|(t-s)\mathrm{d}s\leq
T\|N_{1}-N_{2}\|_{\infty},$

and for $A_{3}$ we get

$\begin{split}A_{3}&\leq\|p\|_{\infty}\|n^{0}\|_{1}\int_{0}^{t}\int_{0}^{s}|p(s^{\prime},N_{1}(t-s+s^{\prime}))-p(s^{\prime},N_{2}(t-s+s^{\prime}))|\mathrm{d}s^{\prime}\,\mathrm{d}s\\\
&\leq\frac{T^{2}}{2}\|p\|_{\infty}\|n^{0}\|_{1}\|\partial_{N}p\|_{\infty}\|N_{1}-N_{2}\|_{\infty}.\end{split}$

By combining these inequalities we have finally

$\sup_{t\in[0,T]}\|n_{1}(t,\cdot)-n_{2}(t,\cdot)\|_{1}\leq\left(\|n^{0}\|_{1}\|\partial_{N}p\|_{\infty}+\frac{T}{2}\|p\|_{\infty}\|n^{0}\|_{1}\|\partial_{N}p\|_{\infty}+1\right)T\|N_{1}-N_{2}\|_{\infty}.$
(17)

Since the mass is conserved, we apply Lemma 2.2 to get

$\sup_{t\in[0,T]}|\psi(n^{1}(t,\cdot))-\psi(n^{2}(t,\cdot))|\leq\frac{\|p\|_{\infty}}{1-\gamma\|n^{0}\|_{1}}\sup_{t\in[0,T]}\|n^{1}(t,\cdot)-n^{2}(t,\cdot)\|_{1},$
(18)

and therefore, we deduce from (17) and (18) that the following estimate holds
for $H$

$\|H(N_{1})-H(N_{2})\|_{\infty}\leq\frac{\|p\|_{\infty}T}{1-\gamma\|n^{0}\|_{1}}\left(\frac{T}{2}\|p\|_{\infty}\|n^{0}\|_{1}\|\partial_{N}p\|_{\infty}+\|n^{0}\|_{1}\|\partial_{N}p\|_{\infty}+1\right)\|N_{1}-N_{2}\|_{\infty},$

so that for $T>0$ small enough we get a unique fixed point of $H$ by
contraction principle, implying the existence of a unique solution
$n\in\mathcal{C}\left([0,T],L^{1}(\mathbb{R}^{+})\right)$ with
$N\in\mathcal{C}[0,T]$ of Equation (1). Since the mass is conserved, we can
iterate this argument to conclude the solution is indeed defined for all
$T>0$. ∎

###  2.2 Numerical scheme for ITM

In this subsection, a first-order explicit upwind scheme to approach ITM
equation (1) is introduced based on the finite-volume framework [24, 25, 26].
We consider an uniform discretization of $\Omega=[0,T]\times[0,+\infty)$ with
cells $I_{j}=[s_{j-\frac{1}{2}},s_{j+\frac{1}{2}})$, interface points
$s_{j+\frac{1}{2}}=j\Delta s>0$ and centers $s_{j}=(j-\frac{1}{2})\Delta s$,
$j\in\mathbb{N}$ such that $[0,+\infty)=\cup_{j\in\mathbb{N}}I_{j}$ and
$I^{m}=[t^{m-1},t^{m}]$, $t^{m}=m\Delta t>0$, $m\in\mathbb{N}$, $T=M\Delta t$
such that $[0,T]=\cup_{m=0}^{M}I^{m}$. Let $n_{j}^{m}=\frac{1}{\Delta
s}\int_{I_{j}}n(t^{m},s)\mathrm{d}s$ be the cell average of $n(t,s)$ at time
$t^{m}$ in the cell $I_{j}$, then applying an explicit finite-volume
approximation with an upwind discretization for the convective term in
Equation (1), we obtain

$\frac{n_{j}^{m+1}-n_{j}^{m}}{\Delta t}+\frac{n_{j}^{m}-n_{j-1}^{m}}{\Delta
s}+p(s_{j},N(t^{m}))n_{j}^{m}=0,\quad j\in\mathbb{N},\,m\in\mathbb{N},$ (19)

and

$N(t^{m})=n(t^{m},0)\approx\Delta
s\sum_{j\in\mathbb{N}}p(s_{j},N(t^{m}))n_{j}^{m},\quad m\in\mathbb{N}.$

Now, if we define $N^{m}:=N(t^{m})$, the solution of the partial differential
equation in (1) can be solved by the explicit upwind scheme

$n_{j}^{m+1}=n_{j}^{m}-\frac{\Delta t}{\Delta s}(n_{j}^{m}-n_{j-1}^{m})-\Delta
tp(s_{j},N^{m})n_{j}^{m},\quad j\in\mathbb{N},\,m\in\mathbb{N}.$ (20)

In particular for $j=1$ we have

$n_{1}^{m+1}=n_{1}^{m}-\frac{\Delta t}{\Delta s}(n_{1}^{m}-N^{m})-\Delta
tp(s_{1},N^{m})n_{1}^{m}.$

The explicit upwind scheme (20) is stable if the CFL condition hold

$1-\Delta t\left(\frac{1}{\Delta s}+p(s_{j},N^{m})\right)\geq 0,\quad
j\in\mathbb{N},\,m\in\mathbb{N}.$

In that regard, the numerical scheme for $n_{j}^{m}$ and $N^{m}$ can be
summarized in the next algorithm.

######  Algorithm 2.1.

ITM numerical scheme

  * Input: Approximate initial data $\\{n^{0}_{j}\\}_{j\in\mathbb{N}}$

  * Solve for $N^{0}$

$N^{0}=\sum_{j\in\mathbb{N}}\Delta sp(s_{j},N^{0})n_{j}^{0}.$ (21)

  * Choose $\Delta t$ such that

$\Delta t\leq\left(\frac{1}{\Delta s}+\|p\|_{\infty}\right)^{-1}.$ (22)

  * For $m\in\mathbb{N}_{0}$ do

  * For $j\in\mathbb{N}$ do

$n_{j}^{m+1}\leftarrow\begin{cases}n_{1}^{m}-\frac{\Delta t}{\Delta
s}(n_{1}^{m}-N^{m})-\Delta tp(s_{1},N^{m})n_{1}^{m}&j=1,\\\
n_{j}^{m}-\frac{\Delta t}{\Delta s}(n_{j}^{m}-n_{j-1}^{m})-\Delta
tp(s_{j},N^{m})n_{j}^{m}&j>1.\end{cases}$

  * end

  * Solve for $N^{m+1}$

$N^{m+1}=\sum_{j\in\mathbb{N}}\Delta sp(s_{j},N^{m+1})n_{j}^{m+1}.$ (23)

  * end

  * Output: Approximate solution $\\{n_{j}^{m+1}\\}_{j\in\mathbb{N}}$ and $N^{m+1}$ at time $t^{m+1}=(m+1)\Delta t$

The solution of Equation (23) for $N^{m+1}$ can be solved with different
numerical methods such as Newton-Raphson, bisection or inverse quadratic
interpolation. In particular for the inhibitory case and the weak
interconnections regime, the solution of Equation (23) $N^{m+1}$ can be
approximated in terms of $N^{m}$ through the following formula if $\Delta t$
is small enough:

$N^{m+1}=\sum_{j\in\mathbb{N}}\Delta sp(s_{j},N^{m})n_{j}^{m+1}.$

For simplicity in the estimates we assume that we can compute the solution of
Equation (23) exactly, but the results remain valid if we take into account an
specific method to get an approximation.

######  Remark 2.2.

Analogously to Remark 1.1, we will prove in this section that Equations (21)
and (23) have a unique solution in the inhibitory case and the weak
interconnections regime. In the excitatory case, we may have multiple
solutions for $N^{0}$ that lead to different branches of numerical solutions
for the ITM equation (1) defined in some interval of time. Depending on how we
calculate the solution of the fixed point problem (23), the numerical method
will approximate one of the multiple possible solutions.

In order to prove the convergence of the upwind scheme we follow the ideas of
the previous subsection on well-posedness and we prove a BV-estimate that will
be the crucial in the analysis. For simplicity we assume that initial data
$n^{0}$ is compactly supported, but the theoretical results still hold when
the initial data $n^{0}$ vanishes at infinity.

We start with some lemmas that will be useful in the sequel.

######  Lemma 2.3.

($L^{1}$-norm) Numerical approximation obtained with the Algorithm 2.1
satisfies

$\|n^{m}\|_{1}:=\sum_{j\in\mathbb{N}}\Delta sn_{j}^{m}=\|n^{0}\|_{1},\quad
m\in\mathbb{N}.$

###### Proof.

Multiply equation (20) by $\Delta s$ and sum over $j\in\mathbb{N}$ and using
the boundary condition we obtain

$\sum_{j\in\mathbb{N}}\Delta sn_{j}^{m+1}=\sum_{j\in\mathbb{N}}\Delta
sn_{j}^{m}-\Delta t\left(N^{m}-\sum_{j\in\mathbb{N}}\Delta
sp(s_{j},N^{m})n_{j}^{m}\right).$

∎

######  Lemma 2.4.

($L^{\infty}$-norm) Assume the initial data $n^{0}\in
L^{\infty}(\mathbb{R}^{+})$ is non-negative. Then under CFL restriction (22)
the numerical solution obtained with the Algorithm 2.1 satisfies

$\displaystyle 0\leq n^{m}_{j}\leq\|n^{0}\|_{\infty},$ $\displaystyle 0\leq
N^{m}\leq\|p\|_{\infty}\|n^{0}\|_{1}$ $\displaystyle\textrm{for
all}\>j,m\in\mathbb{N}.$

###### Proof.

By the CFL condition, observe that

$0\leq n_{j}^{1}=\left(1-\Delta t\left(\frac{1}{\Delta
s}+p(s_{j},N^{0})\right)\right)n_{j}^{0}+\frac{\Delta t}{\Delta
s}n_{j-1}^{0}\leq\|n^{0}\|_{\infty}\quad\textrm{for all}\>j\in\mathbb{N},$

and by induction in $m$ we conclude that $0\leq
n^{m}_{j}\leq\|n^{0}\|_{\infty}$ for all $j,m\in\mathbb{N}$. Now for the
estimates involving $N^{m}$, we apply the previous lemma to conclude the
following inequality

$0\leq N^{m}\leq\|p\|_{\infty}\sum_{j\in\mathbb{N}}\Delta
sn^{m}_{j}=\|p\|_{\infty}\sum_{j\in\mathbb{N}}\Delta
sn^{0}_{j}=\|p\|_{\infty}\|n^{0}\|_{1},$

and we get the desired result. ∎

We now prove that for each iteration in $m$ of the numerical method. Equation
(23) has a indeed a unique solution $N^{m}$ in the inhibitory case and the
weak interconnection regime. The following result corresponds to the discrete
version of Lemma 2.2.

######  Lemma 2.5.

Consider a discretization $n^{0}=(n^{0}_{j})_{j\in\mathbb{N}}$ of non-negative
terms. Let $\gamma\coloneqq\sup_{s,N}\partial_{N}p(s,N)$ and assume that
$\gamma\|n^{0}\|_{1}<1$. Then there exists a unique solution for $N$ of the
equation

$N=F(N)\coloneqq\sum_{j\in\mathbb{N}}\Delta sp(s_{j},N)n^{0}_{j},$

that we call $n^{0}\coloneqq\psi(n^{0})$. Moreover the map $\psi$ satisfies
the following estimate

$|\psi(n^{1})-\psi(n^{2})|\leq\frac{\|p\|_{\infty}}{1-\gamma\|n^{0}\|_{1}}\sum_{j\in\mathbb{N}}\Delta
s|n^{1}_{j}-n^{2}_{j}|$ (24)

for sequences
$n^{1}=(n^{1}_{j})_{j\in\mathbb{N}},\,n^{2}=(n^{2}_{j})_{j\in\mathbb{N}}$ of
non-negative terms with $\|n^{1}\|_{1}=\|n^{2}\|_{1}=\|n^{0}\|_{1}$.

###### Proof.

The proof is similar to that of Lemma 2.2 by the replacing the integral terms
with discrete summations. ∎

######  Remark 2.3.

In the excitatory case, the ideas of the previous lemma can be applied as long
as the following invertibility condition holds

$\Psi(N,n):=1-\Delta
s\sum_{j\in\mathbb{N}}\partial_{N}p(s_{j},N)n_{j}\mathrm{d}s\neq 0,$ (25)

so that we can extend the numerical solution through the implicit function
theorem in order to approximate a continuous solution of the ITM equation (10)
in some interval $[0,T]$. This is the analog to the invertibility condition
(13).

We now establish some BV-lemmas on the discretization given by the upwind
scheme to prove that the numerical method approximates the solution of
Equation (1) when $\Delta t$ and $\Delta s$ converge to zero. For the
discretization $n=(n_{j})_{j\in\mathbb{N}}$ we define the total variation as

$TV(n)\coloneqq\sum_{j=0}^{\infty}|n_{j+1}-n_{j}|.$ (26)

In this context, we prove the following key lemma.

######  Lemma 2.6.

(BV-estimate) Assume that $TV(n^{0})<\infty$ and the CFL condition (22). Then
there exist constants $C_{1},C_{2}>0$ (depending only $p$ and the norms of
$n^{0}$) such that for $m\in\mathbb{N}$ we have

$TV(n^{m})\leq e^{C_{1}T}TV(n^{0})+C_{2}(e^{C_{1}T}-1),$ (27)

with $T=m\Delta t$ and
$TV(n^{m})=\displaystyle\sum_{j=0}^{\infty}|n_{j+1}^{m}-n_{j}^{m}|$.

###### Proof.

Using the notation $\Delta^{+}n_{j}^{m}=n_{j+1}^{m}-n_{j}^{m}$, we have

$\displaystyle\Delta^{+}n_{j}^{m+1}$ $\displaystyle=$
$\displaystyle\Delta^{+}n_{j}^{m}-\frac{\Delta t}{\Delta
s}\left(\Delta^{+}n_{j}^{m}-\Delta^{+}n_{j-1}^{m}\right)-\Delta
t\left(p(s_{j+1},N^{m})n_{j+1}^{m}-p(s_{j},N^{m})n_{j}^{m}\right)$
$\displaystyle=$ $\displaystyle\Delta^{+}n_{j}^{m}-\frac{\Delta t}{\Delta
s}\left(\Delta^{+}n_{j}^{m}-\Delta^{+}n_{j-1}^{m}\right)-\Delta
tp(s_{j+1},N^{m})\Delta^{+}n_{j}^{m}$ $\displaystyle-\Delta
tn_{j}^{m}(p(s_{j+1},N^{m})-p(s_{j},N^{m})).$

now applying (8) and taking absolute value and using the CFL condition (22) we
obtain

$|\Delta^{+}n_{j}^{m+1}|\leq\left(1-\Delta t\left(\frac{1}{\Delta
s}+p(s_{j+1},N^{m})\right)\right)|\Delta^{+}n_{j}^{m}|+\frac{\Delta t}{\Delta
s}|\Delta^{+}n_{j-1}^{m}|+\Delta t\Delta
s\|\partial_{s}p\|_{\infty}n_{j}^{m}.$

Now by summing over all $j\geq 1$, we have

$\begin{split}\sum_{j=1}^{\infty}|\Delta^{+}n_{j}^{m+1}|&\leq\sum_{j=1}^{\infty}|\Delta^{+}n_{j}^{m}|+\frac{\Delta
t}{\Delta s}|\Delta^{+}n_{0}^{m}|-\Delta
t\sum_{j\in\mathbb{N}}^{\infty}p(s_{j+1},N^{m})|\Delta^{+}n_{j}^{m}|\\\
&\qquad+\Delta t\|\partial_{s}p\|_{\infty}\sum_{j=1}^{\infty}\Delta
sn_{j}^{m}\end{split}$

and from the mass conservation we deduce

$\sum_{j=1}^{\infty}|\Delta^{+}n_{j}^{m+1}|\leq\sum_{j=1}^{\infty}|\Delta^{+}n_{j}^{m}|+\frac{\Delta
t}{\Delta s}|\Delta^{+}n_{0}^{m}|+\Delta
t\|\partial_{s}p\|_{\infty}\|n^{0}\|_{1}.$ (28)

On the other hand, by assuming $\Delta t$ small enough we have

$\begin{split}|\Delta^{+}n_{0}^{m+1}|&=|n_{1}^{m+1}-N^{m+1}|\\\
&\leq\left(1-\frac{\Delta t}{\Delta
s}\right)|n_{1}^{m}-N^{m}|+|N^{m+1}-N^{m}|+\Delta tp(s_{1},N^{m})n^{m}_{1}\\\
&\leq\left(1-\frac{\Delta t}{\Delta
s}\right)|\Delta^{+}n_{0}^{m}|+|N^{m+1}-N^{m}|+\Delta
t\|p\|_{\infty}\|n^{0}\|_{\infty}.\end{split}$ (29)

Furthermore, from Lemma 2.5 there exists $C_{1}>0$ depending only on $p$ such
that

$\displaystyle|N^{m+1}-N^{m}|$ $\displaystyle=$
$\displaystyle\left|\psi(n^{m+1})-\psi(n^{m})\right|$ (30) $\displaystyle\leq$
$\displaystyle C_{1}\sum_{j\in\mathbb{N}}\Delta s|n_{j}^{m+1}-n_{j}^{m}|$
$\displaystyle\leq$ $\displaystyle C_{1}\Delta
t\left(\sum_{j\in\mathbb{N}}|n_{j}^{m}-n_{j-1}^{m}|+\|p\|_{\infty}\sum_{j\in\mathbb{N}}\Delta
sn_{j}^{m}\right),$

and replacing (30) in (29), we obtain from the mass conservation

$\begin{split}|\Delta^{+}n_{0}^{m+1}|&\leq\left(1-\frac{\Delta t}{\Delta
s}\right)|\Delta^{+}n_{0}^{m}|+C_{1}\Delta
t\left(\sum_{j\in\mathbb{N}}|\Delta^{+}n^{m}_{j-1}|+\|p\|_{\infty}\|n^{0}\|_{1}\right)\\\
&\qquad+\Delta t\|p\|_{\infty}\|n^{0}\|_{\infty}.\end{split}$ (31)

Now, summing (28) and (31), we deduce

$\sum_{j=0}^{\infty}|\Delta^{+}n_{j}^{m+1}|\leq\left(1+C_{1}\Delta
t\right)\sum_{j=0}^{\infty}|\Delta^{+}n_{j}^{m}|+C_{2}\Delta t.$ (32)

with
$C_{2}\coloneqq\|\partial_{s}p\|_{\infty}\|n^{0}\|_{1}+\|p\|_{\infty}C_{1}\|n^{0}\|_{1}+\|p\|_{\infty}\|n^{0}\|_{\infty}$.

Finally, proceeding recursively on $m$, we obtain

$\sum_{j=0}^{\infty}|\Delta^{+}n_{j}^{m}|\leq(1+C_{1}\Delta
t)^{m}\sum_{j=0}^{\infty}|\Delta^{+}n_{j}^{0}|+\frac{C_{2}}{C_{1}}\left((1+C_{1}\Delta
t)^{m}-1\right),$

and the estimate (27) readily follows. ∎

######  Remark 2.4.

The previous lemma is also valid for the case when the hazard rate is of the
form

$p(s,N)=\varphi(N)\chi_{\\{s>\sigma(N)\\}}$

with $\varphi$ and $\sigma$ Lipschitz bounded functions and satisfying the
analogous hypothesis on $\partial_{N}p$. Indeed, the inequality in (28) is
replaced by

$\begin{split}\sum_{j=1}^{\infty}|\Delta^{+}n_{j}^{m+1}|&\leq\sum_{j=1}^{\infty}|\Delta^{+}n_{j}^{m}|+\frac{\Delta
t}{\Delta s}|\Delta^{+}n_{0}^{m}|+\Delta
t\varphi(N^{m})\sum_{j\in\mathbb{N}}n^{m}_{j}\left|\chi_{\\{s_{j+1}>\sigma(N^{m})\\}}-\chi_{\\{s_{j}>\sigma(N^{m})\\}}\right|\\\
&\leq\sum_{j=1}^{\infty}|\Delta^{+}n_{j}^{m}|+\frac{\Delta t}{\Delta
s}|\Delta^{+}n_{0}^{m}|+\Delta t\|p\|_{\infty}n^{m}_{j_{m}}\\\
&\leq\sum_{j=1}^{\infty}|\Delta^{+}n_{j}^{m}|+\frac{\Delta t}{\Delta
s}|\Delta^{+}n_{0}^{m}|+\Delta t\|p\|_{\infty}\|n^{0}\|_{\infty},\end{split}$

where $j_{m}\coloneqq\min\\{j\in N\colon
s_{j}\geq\sigma(N^{m})\\}=\left\lceil\frac{\sigma(N^{m})}{\Delta
s}+\frac{1}{2}\right\rceil$ and the rest of proof is analogous.

Now we prove that the numerical approximation of the solution of Equation (1)
$n(t,s)$, which is constructed by a simple piece-wise linear interpolation,
has a limit when the time step $\Delta t$ and age step $\Delta s$ converge to
$0$. For simplicity we assume that the initial data $n^{0}\in
BV(\mathbb{R}^{+})$ and with compact support.

######  Lemma 2.7.

Assume that $n^{0}\in BV(\mathbb{R}^{+})$ is compactly supported and the rate
$p$ satisfies the hypothesis of Theorem 2.1. Consider the function $n_{\Delta
t,\Delta s}\in\mathcal{C}\left([0,T],L^{1}(\mathbb{R}^{+})\right)$ defined by

$n_{\Delta t,\Delta s}(t,s)\coloneqq\frac{t^{m}-t}{\Delta
t}\sum_{j\in\mathbb{N}}n_{j}^{m-1}\chi_{[s_{j-\frac{1}{2}},s_{j+\frac{1}{2}}]}(s)+\frac{t-t^{m-1}}{\Delta
t}\sum_{j\in\mathbb{N}}n_{j}^{m}\chi_{[s_{j-\frac{1}{2}},s_{j+\frac{1}{2}}]}(s)\;\textrm{if}\;t\in[t^{m-1},t^{m}].$
(33)

Then there exists a sub-sequence $(\Delta t_{k},\Delta s_{k})\to(0,0)$ when
$k\to\infty$ and a function $\overline{n}(t,s)$ such that $n_{\Delta
t_{k},\Delta s_{k}}\to\overline{n}$ in
$\mathcal{C}\left([0,T],L^{1}(\mathbb{R}^{+})\right)$. Moreover, if we define
the function $N_{\Delta t,\Delta s}(t)$ as the unique solution of the equation

$N(t)=\int_{0}^{\infty}p(s,N(t))n_{\Delta t,\Delta s}(t,s)\,\mathrm{d}s,$

then there exists $\overline{N}$ such that $N_{\Delta t_{k},\Delta
s_{k}}\to\overline{N}$ in $\mathcal{C}[0,T]$ and $\overline{N}$ is a solution
of the equation

$\overline{N}(t)=\int_{0}^{\infty}p(s,\overline{N}(t))\overline{n}(t,s)\,\mathrm{d}s.$
(34)

###### Proof.

The idea is to apply the compactness criterion in
$\mathcal{C}([0,T],L^{1}(\mathbb{R}^{+}))$ in order to extract a convergent
sub-sequence of $n_{\Delta t,\Delta s}$ when $\Delta t$ and $\Delta s$
converge to $0$. From Lemma 2.3, we deduce that

$\|n_{\Delta t,\Delta s}(t,\cdot)\|_{1}=\|n^{0}\|_{1},\qquad\forall
t\in[0,T].$

We prove that sequence $n_{\Delta t,\Delta s}$ has a modulus of continuity in
the $L^{1}$ in both variables. In the variable $s$ we have that for
$t\in[t^{m-1},t^{m}]$ and $|h|<\varepsilon$ the following estimate holds

$\begin{split}\int_{0}^{\infty}|n_{\Delta t,\Delta s}(t,s+h)-n_{\Delta
t,\Delta s}(t,s)|\,\mathrm{d}s&\leq\frac{t^{m}-t}{\Delta
t}|h|\sum_{j\in\mathbb{N}}\Delta sn_{j}^{m-1}+\frac{t-t^{m-1}}{\Delta
t}|h|\sum_{j\in\mathbb{N}}\Delta sn_{j}^{m}\\\
&\leq\varepsilon\|n^{0}\|_{1}\end{split}$

Now we prove we have modulus of continuity in the variable $t$. Consider
$t_{1},t_{2}\in[0,T]$ and without loss of generality assume that
$t_{1},t_{2}\in[t^{m-1},t^{m}]$ for some $m\in\mathbb{N}$. Then from Lemma 2.6
we have following estimate

$\begin{split}\int_{0}^{\infty}|n_{\Delta t,\Delta s}(t_{1},s)-n_{\Delta
t,\Delta
s}(t_{2},s)|\,\mathrm{d}s&\leq|t_{1}-t_{2}|\sum_{j\in\mathbb{N}}\Delta
s\frac{|n^{m}_{j}-n^{m-1}_{j}|}{\Delta t}\\\
&\leq|t_{1}-t_{2}|\left(\sum_{j\in\mathbb{N}}|n^{m-1}_{j}-n^{m-1}_{j-1}|+\Delta
sp(s_{j},N^{m-1})n^{m-1}_{j}\right)\\\
&\leq|t_{1}-t_{2}|\left(C_{T}\,TV(n^{0})+\|p\|_{\infty}\|n^{0}\|_{1}\right),\end{split}$
(35)

thus we have the modulus of continuity in time.

Since $n^{0}$ has its support contained in some interval $[0,R]$, there exists
$K\in\mathbb{N}$ such that $n^{0}_{j}=0$ for $j\geq K$ and
$s_{K}=(K-\frac{1}{2})\Delta s\geq R$. From the numerical scheme we deduce
that $n^{m}_{j}=0$ for $j\geq K+M$, which implies that $n_{\Delta t,\Delta s}$
vanishes for

$s\geq s_{K+M}=(K+M-\frac{1}{2})\Delta s\geq R+M\Delta{t}\,\frac{\Delta
s}{\Delta t}\geq R+T,$

so $n_{\Delta t,\Delta s}$ has also compact support. From the estimates on the
modulus of continuity, we can apply the compactness criterion in
$\mathcal{C}([0,T],L^{1}(\mathbb{R}^{+}))$ in order to extract a convergent
sub-sequence of $n_{\Delta t,\Delta s}$ when $\Delta t$ and $\Delta s$
converge to $0$.

Observe now that $N_{\Delta t,\Delta s}\in\mathcal{C}[0,T]$. Indeed from Lemma
2.2, for each $t\in[0,T]$ we get that

$N_{\Delta t,\Delta s}(t)=\psi(n_{\Delta t,\Delta s})$

and from the continuity of $\psi$ and $n_{\Delta t,\Delta
s}\in\mathcal{C}([0,T],L^{1}(\mathbb{R}^{+})$ we obtain that $N_{\Delta
t,\Delta s}\in\mathcal{C}[0,T]$. Since $\|n_{\Delta t,\Delta
s}(t,\cdot)\|_{1}=\|n^{0}\|_{1},$ for all $t\in[0,T]$ we have the following
estimate

$\|N_{\Delta t,\Delta s}\|_{\infty}\leq\|p\|_{\infty}\|n^{0}\|_{1},$

so that $N_{\Delta t,\Delta s}$ is uniformly bounded.

We now prove that the family $N_{\Delta t,\Delta s}$ is equicontinuous. For
$t_{1},t_{2}\in[0,T]$ we deduce from Lemma (2.2) and estimate (35) the
following inequality

$\begin{split}\left|N_{\Delta t,\Delta s}(t_{1})-N_{\Delta t,\Delta
s}(t_{2})\right|&=|\psi(n_{\Delta t,\Delta s})(t_{1})-\psi(n_{\Delta t,\Delta
s})(t_{2})|\\\ &\leq C\int_{0}^{\infty}|n_{\Delta t,\Delta
s}(t_{1},s)-n_{\Delta t,\Delta s}(t_{2},s)|\,\mathrm{d}s\\\
&\leq|t_{1}-t_{2}|\left(C_{T}\,TV(n^{0})+\|p\|_{\infty}\|n^{0}\|_{1}\right),\end{split}$

where $C$ is a constant independent of $\Delta t$ and $\Delta s$. Therefore
the family $N_{\Delta t,\Delta s}$ is equicontinuous in $[0,T]$ and we can
extract a convergent sub-sequence by applying Arzelà-Ascoli Theorem. Finally,
by passing to the limit in $\Delta t$ and $\Delta s$ we obtain Equation (34).
∎

With the previous lemmas, we are now ready to the prove the following theorem
on convergence of the numerical scheme for the ITM equation (1).

######  Theorem 2.2 (Convergence of the numerical scheme).

Assume that $n^{0}\in BV(\mathbb{R}^{+})$ is compactly supported and the rate
$p$ satisfies the hypothesis of Theorem 2.1. Then for all $T>0$, the numerical
approximation (33) given by the upwind scheme converges to the unique weak
solution $n\in\mathcal{C}([0,T],L^{1}(\mathbb{R}^{+}))$ of the ITM equation
(1).

###### Proof.

Consider the functions $n_{\Delta t,\Delta s}$ and $N_{\Delta t,\Delta s}$
defined in Lemma 2.7. From this result we get Equation (2.7). Now we take
$\varphi\in\mathcal{C}^{1}_{c}([0,T)\times[0,\infty))$ a test function. If we
multiply Equation (19) by $\varphi_{j}^{m}\coloneqq\varphi(t^{m},s_{j})$ and
compute the discrete integral we get

$\sum_{m=0}^{M}\sum_{j\in\mathbb{N}}\Delta
s(n_{j}^{m+1}-n_{j}^{m})\varphi_{j}^{m}+\sum_{m=0}^{M}\sum_{j\in\mathbb{N}}\Delta
t(n_{j}^{m}-n_{j-1}^{m})\varphi_{j}^{m}+\sum_{m=0}^{M}\sum_{j\in\mathbb{N}}\Delta
t\Delta sp(s_{j},N^{m})n_{j}^{m}\varphi_{j}^{m}=0,$ (36)

We study each term of Equation (36). From summation by parts we have the
following inequality for the first term

$\sum_{m=0}^{M}\sum_{j\in\mathbb{N}}\Delta
s(n_{j}^{m+1}-n_{j}^{m})\varphi_{j}^{m}=-\sum_{j\in\mathbb{N}}\Delta
s\varphi^{0}_{j}n^{0}_{j}-\sum_{m=1}^{M}\sum_{j\in\mathbb{N}}\Delta
sn^{m}_{j}(\varphi^{m}_{j}-\varphi^{m-1}_{j}),$

thus applying Lemma 2.7, we get the following limit when $(\Delta t,\Delta
s)\to 0$

$\sum_{m=0}^{M}\sum_{j\in\mathbb{N}}\Delta
s(n_{j}^{m+1}-n_{j}^{m})\varphi_{j}^{m}\to-\int_{0}^{\infty}\varphi(0,s)n^{0}(s)\,\mathrm{d}s-\int_{0}^{T}\int_{0}^{\infty}\overline{n}(t,s)\partial_{t}\varphi(t,s)\,\mathrm{d}s\,\mathrm{d}t.$

Similarly, for the second term of Equation (36) the following equality holds

$\sum_{m=0}^{M}\sum_{j\in\mathbb{N}}\Delta
t(n_{j}^{m}-n_{j-1}^{m})\varphi_{j}^{m}=-\Delta
t\varphi_{0}^{0}n^{0}_{0}-\sum_{m=1}^{M}\Delta
t\varphi_{0}^{m}N^{m}-\sum_{m=0}^{M}\sum_{j\in\mathbb{N}}\Delta
tn^{m}_{j}(\varphi^{m}_{j}-\varphi^{m}_{j-1}),$

and by passing to the limit in $(\Delta t,\Delta s)$ we get

$\sum_{m=0}^{M}\sum_{j\in\mathbb{N}}\Delta
t(n_{j}^{m}-n_{j-1}^{m})\varphi_{j}^{m}\to-\int_{0}^{T}\varphi(t,0)\overline{N}(t)\,\mathrm{d}t-\int_{0}^{T}\int_{0}^{\infty}\overline{n}(t,s)\partial_{s}\varphi(t,s)\,\mathrm{d}s\,\mathrm{d}t,$

and in the same way

$\sum_{m=0}^{M}\sum_{j\in\mathbb{N}}\Delta t\Delta
s\,p(s_{j},N^{m})n_{j}^{m}\varphi_{j}^{m}\to\int_{0}^{T}\int_{0}^{\infty}p(s,\overline{N}(t))\overline{n}(t,s)\varphi(t,s)\,\mathrm{d}s\,\mathrm{d}t.$

Therefore $\overline{n}(t,s)$ is the weak solution of the ITM equation (1). ∎

##  3 Distributed Delay Model (DDM)

###  3.1 Well-posedness of DDM

In this subsection we prove the well-posedness of the DDM equation (10). As in
Section 2, we improve the proof of [8] by extending existence and uniqueness
for a more general calls of hazard rates $p$. We essentially follow the same
ideas from of the previous section with some slight modifications.

######  Theorem 3.1.

Assume that $p\in W^{1,\infty}(\mathbb{R}^{+}\times\mathbb{R}^{+})$ and
$\alpha\in L^{1}(\mathbb{R}^{+})$ is bounded. Then for a non-negative
$n^{0}\in L^{1}(\mathbb{R}^{+})$, Equation (10) has a unique solution
$n\in\mathcal{C}\left([0,\infty),L^{1}(\mathbb{R}^{+})\right)$ and
$N,X\in\mathcal{C}[0,\infty)$.

For the proof we need the following lemma, which is the analogous of Lemma
2.2.

######  Lemma 3.1.

Consider a non-negative function
$n\in\mathcal{C}\left([0,T],L^{1}(\mathbb{R}^{+})\right)$. Then for $T$ small
enough, the exists a unique solution $X\in\mathcal{C}([0,T])$ of the integral
equation

$X(t)=F(X(t))\coloneqq\int_{0}^{t}\int_{0}^{\infty}\alpha(t-\tau)p(s,X(\tau))n(\tau,s)\,\mathrm{d}s\,\mathrm{d}\tau,$
(37)

that we call $X\coloneqq\psi(n)$, where the map
$\psi\colon\mathcal{C}\left([0,T],L^{1}(\mathbb{R}^{+})\right)\mapsto\mathcal{C}([0,T])$
satisfies the following estimate

$\|\psi(n^{1})-\psi(n^{2})\|_{\infty}\leq\tfrac{T\|\alpha\|_{\infty}\|p\|_{\infty}}{1-T\|\alpha\|_{\infty}\|\partial_{X}p\|_{\infty}\|n^{0}\|_{1}}\sup_{t\in[0,T]}\|n^{1}(t,\cdot)-n^{2}(t,\cdot)\|_{1}$
(38)

for non-negative integrable functions $n^{1},n^{2}$ with
$\|n^{1}(t,\cdot)\|_{1}=\|n^{2}(t,\cdot)\|_{1}=\|n^{0}\|_{1}$ for all
$t\in[0,T]$.

###### Proof.

Observe that the map $F\colon\mathcal{C}[0,T]\mapsto\mathcal{C}[0,T]$
satisfies the following estimate

$\|F(X_{1})-F(X_{2})\|_{\infty}\leq
T\|\alpha\|_{\infty}\|\partial_{X}p\|_{\infty}\|n\|\|X_{1}-X_{2}\|_{\infty},$

with $\|n\|=\sup_{t\in[0,T]}\|n(t,\cdot)\|_{1}$. Hence for $T>0$ such that
$T\|\alpha\|_{\infty}\|\partial_{X}p\|_{\infty}\|n\|<1$ the map $F$ is a
contraction and then the map $\psi$ is well-defined.

Let $n^{1},n^{2}\in\mathcal{C}([0,T],L^{1}(\mathbb{R}^{+}))$ with
$\|n^{1}(t,\cdot)\|_{1}=\|n^{2}(t,\cdot)\|_{1}=\|n^{0}\|_{1}$ for all
$t\in[0,T]$. Then we have the following inequality

$\begin{split}|X_{1}-X_{2}|(t)&\leq\int_{0}^{t}\int_{0}^{\infty}\alpha(t-\tau)\left(|p(s,X_{1}(\tau))-p(s,X_{2}(\tau))|n^{1}(\tau,s)+p(s,X_{2}(\tau))|n^{1}-n^{2}|(\tau,s)\right)\mathrm{d}s\mathrm{d}\tau\\\
&\leq
T\|\alpha\|_{\infty}\|\partial_{X}p\|_{\infty}\|n^{0}\|_{1}\|X_{1}-X_{2}\|_{\infty}+T\|\alpha\|_{\infty}\|p\|_{\infty}\sup_{t\in[0,T]}\|n^{1}(t,\cdot)-n^{2}(t,\cdot)\|_{1},\end{split}$

and therefore for $T>0$ such that
$T\|\alpha\|_{\infty}\|\partial_{X}p\|_{\infty}\|n^{0}\|_{1}<1$ we get

$\|X_{1}-X_{2}\|_{\infty}\leq\frac{T\|\alpha\|_{\infty}\|p\|_{\infty}}{1-T\|\alpha\|_{\infty}\|\partial_{X}p\|_{\infty}\|n^{0}\|_{1}}\sup_{t\in[0,T]}\|n^{1}(t,\cdot)-n^{2}(t,\cdot)\|_{1},$

and estimate (38) holds. ∎

With this lemma, we continue the proof of the Theorem 3.1.

###### Proof.

Consider $T>0$ and a given non-negative $X\in\mathcal{C}[0,T]$. Like in the
proof of Theorem 2.1, we define
$n[X]\in\mathcal{C}([0,T],L^{1}(\mathbb{R}^{+}))$ as the weak solution of the
linear equation

$\begin{cases}\partial_{t}n+\partial_{s}n+p(s,X(t))n=0&t>0,\>s>0,\\\
n(t,s=0)=N(t)=\int_{0}^{+\infty}p(s,X(t))n(t,s)\mathrm{d}s&t>0,\\\
n(0,s)=n^{0}(s)\geq 0&s\geq 0.\end{cases}$

which can be expressed trough the method of characteristics

$n(t,s)=n^{0}(s-t)e^{-\int_{0}^{t}p(s-t+t^{\prime},X(t^{\prime}))\mathrm{d}t^{\prime}}\chi_{\\{s>t\\}}+N(t-s)e^{-\int_{0}^{s}p(s^{\prime},X(t-s+s^{\prime}))\mathrm{d}s^{\prime}}\chi_{\\{t>s\\}}.$

From the mass conservation property we know that
$\|n[X](t)\|_{1}=\|n^{0}\|_{1}$.

On the other hand, for a given $n\in\mathcal{C}([0,T],L^{1}(\mathbb{R}^{+}))$
and for $T>0$ small enough, let
$\psi\colon\mathcal{C}\left([0,T],L^{1}(\mathbb{R}^{+})\right)\mapsto\mathcal{C}[0,T]$
the unique solution $X\in\mathcal{C}[0,T]$ of Equation (37). Under this
setting, we get a solution of the non-linear equation (10) if only if we find
$X\in\mathcal{C}[0,T]$ such that is a solution of the equation

$X=H(X)\coloneqq\psi\left(n[X]\right),$

where $H\colon\mathcal{C}[0,T]\mapsto\mathcal{C}[0,T]$.

We assert that for $T$ small enough the map $H$ is a contraction following the
proof Theorem 2.1. For non-negative functions
$X_{1},X_{2}\in\mathcal{C}[0,T]$, we get from the method of characteristics we
get the following inequality

$\int_{0}^{\infty}|n[X_{1}]-n[X_{2}]|(t,s)\mathrm{d}s\leq A_{1}+A_{2}+A_{3},$

where $A_{1},A_{2},A_{3}$ are given by

$\displaystyle A_{1}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}n^{0}(s)\left|e^{-\int_{0}^{t}p(s+t^{\prime},X_{1}(t^{\prime}))\mathrm{d}t^{\prime}}-e^{-\int_{0}^{t}p(s+t^{\prime},X_{2}(t^{\prime}))\mathrm{d}t^{\prime}}\right|\,\mathrm{d}s$
$\displaystyle A_{2}$ $\displaystyle=$
$\displaystyle\int_{0}^{t}|N_{1}-N_{2}|(t-s)e^{-\int_{0}^{s}p(s^{\prime},X_{1}(t-s+s^{\prime}))\mathrm{d}s^{\prime}}\mathrm{d}s$
$\displaystyle A_{3}$ $\displaystyle=$
$\displaystyle\int_{0}^{t}N_{2}(t-s)\left|e^{-\int_{0}^{s}p(s^{\prime},X_{1}(t-s+s^{\prime}))\mathrm{d}s^{\prime}}-e^{-\int_{0}^{s}p(s^{\prime},X_{2}(t-s+s^{\prime}))\mathrm{d}s^{\prime}}\right|\mathrm{d}s.$

We proceed by estimating each term. For simplicity we write
$n_{1}=n[X_{1}],n_{2}=n[X_{2}]$, so that for $A_{1}$ we have

$A_{1}\leq\int_{0}^{\infty}n^{0}(s)\int_{0}^{t}|p(s+t^{\prime},X_{1}(t^{\prime}))-p(s+t^{\prime},X_{2}(t^{\prime}))|\,\mathrm{d}t^{\prime}\,\mathrm{d}s\leq
T\|n^{0}\|_{1}\|\partial_{N}p\|_{\infty}\|X_{1}-X_{2}\|_{\infty},$

while for $A_{2}$ we have

$\begin{split}A_{2}&\leq\int_{0}^{t}\int_{0}^{\infty}|p(s^{\prime},X_{1}(s))-p(s^{\prime},X_{2}(s))|n_{1}(s,s^{\prime})\mathrm{d}s^{\prime}\,\mathrm{d}s+\int_{0}^{t}\int_{0}^{\infty}p(s^{\prime},X_{2}(s))|n_{1}-n_{2}|(s,s^{\prime})\mathrm{d}s^{\prime}\,\mathrm{d}s,\\\
&\leq
T\|\partial_{X}p\|\|n^{0}\|_{1}+T\|p\|_{\infty}\sup_{t\in[0,T]}\|n^{1}(t,\cdot)-n^{2}(t,\cdot)\|_{1},\end{split}$

and for $A_{3}$ we get

$\begin{split}A_{3}&\leq\|p\|_{\infty}\|n^{0}\|_{1}\int_{0}^{t}\int_{0}^{s}|p(s^{\prime},X_{1}(t-s+s^{\prime}))-p(s^{\prime},X_{2}(t-s+s^{\prime}))|\mathrm{d}s^{\prime}\,\mathrm{d}s\\\
&\leq\frac{T^{2}}{2}\|p\|_{\infty}\|n^{0}\|_{1}\|\partial_{N}p\|_{\infty}\|X_{1}-X_{2}\|_{\infty}.\end{split}$

By combining these estimates, we get for $T<\frac{1}{\|p\|_{\infty}}$ the
following inequality

$\sup_{t\in[0,T]}\|n_{1}(t,\cdot)-n_{2}(t,\cdot)\|_{1}\leq\tfrac{T}{1-T\|p\|_{\infty}}\left(2\|n^{0}\|_{1}\|\partial_{N}p\|_{\infty}+\frac{T}{2}\|p\|_{\infty}\|n^{0}\|_{1}\|\partial_{N}p\|_{\infty}\right)\|X_{1}-X_{2}\|_{\infty}.$
(39)

From the mass conservation property, we get from (38) and (39) that the
following estimate holds for $H$

$\begin{split}\|H(X_{1})-H(X_{2})\|_{\infty}&\leq\tfrac{T^{2}\|p\|_{\infty}^{2}\|\alpha\|_{\infty}}{(1-T\|\alpha\|_{\infty}\|n^{0}\|_{1}\|\partial_{X}p\|_{\infty})(1-T\|p\|_{\infty})}\\\
&\qquad\left(2\|n^{0}\|_{1}\|\partial_{N}p\|_{\infty}+\frac{T}{2}\|p\|_{\infty}\|n^{0}\|_{1}\|\partial_{N}p\|_{\infty}\right)\|X_{1}-X_{2}\|_{\infty}.\end{split}$

For $T>0$ small enough we get a unique fixed point of $H$ by contraction
principle, implying the existence of a unique solution
$n\in\mathcal{C}\left([0,T],L^{1}(\mathbb{R}^{+})\right)$ with
$N,X\in\mathcal{C}[0,T]$ of Equation (10). In order to extend the solution for
all times, we split the integral involving the distributed delay

$X(t)=\int_{0}^{T}\alpha(t-\tau)N(\tau)\mathrm{d}\tau+\int_{T}^{t}\alpha(t-\tau)N(\tau)\mathrm{d}\tau.$

Since the first term is already known and $T$ is independent of the initial
data, we can reapply the argument on existence to have the solution of
Equation (10) defined for all $t\in[T,2T]$. By iterating the splitting
argument involving $X(t)$, we conclude that the solution of Equation (10) is
defined for all $t>0$. ∎

######  Remark 3.1.

As in Theorem 2.1, the regularity of the rate $p$ is not fundamental for the
proof and Theorem 2.1 is still valid for rates $p$ of the form

$p(s,X)=\varphi(X)\chi_{\\{s>\sigma(X)\\}},$

with $\varphi$ and $\sigma$ Lipschitz bounded functions. Moreover, the proof
can be adapted when $p$ is not necessarily bounded, but the continuous
solution might not be defined for all $t>0$. We can extend the solution as
long as $X(t)<\infty$.

###  3.2 Numerical scheme for DDM

In this section we make the respective numerical analysis for DDM equation
(10). Using the same discretization and notation from Section 2, Equation (10)
can be solved numerically by the explicit scheme given by

$n_{j}^{m+1}=n_{j}^{m}-\frac{\Delta t}{\Delta s}(n_{j}^{m}-n_{j-1}^{m})-\Delta
tp(s_{j},X(t^{m}))n_{j}^{m},\quad j\in\mathbb{N}.$ (40)

In particular for $j=1$

$n_{1}^{m+1}=n_{1}^{m}-\frac{\Delta t}{\Delta s}(n_{1}^{m}-N(t^{m}))-\Delta
tp(s_{1},X(t^{m}))n_{1}^{m},$

where

$N(t^{m})=n(t^{m},0)\approx\Delta
s\sum_{j\in\mathbb{N}}p(s_{j},X(t^{m}))n_{j}^{m}$ (41)

and using a trapezoidal quadrature rule we have

$X(t^{m})=\int_{0}^{t^{m}}N(t^{m}-s)\alpha(s)\,\mathrm{d}s\approx\frac{\Delta
t}{2}\sum_{k=0}^{m}N(t^{k})\alpha_{m-k},$ (42)

where $\alpha_{m}=\alpha(t^{m})$ and we denote
$\|\alpha\|_{1}\coloneqq\sum_{k\in\mathbb{N}}\Delta t|\alpha_{k}|$. If we set
$X^{m}\coloneqq X(t^{m})$ and $N^{m}\coloneqq N(t^{m})$, we can combine the
equations (41) and (41) and to solve for $X^{m}$

$X^{m}=\frac{\Delta t}{2}\left(\Delta
s\sum_{j\in\mathbb{N}}p(s_{j},X^{m})n_{j}^{m}\alpha_{0}+\sum_{k=0}^{m-1}N^{k}\alpha_{m-k},\right)$
(43)

and then we obtain $N^{m}$ from (41). In that regard, the numerical method is
given as follows.

######  Algorithm 3.1.

DDM numerical scheme

  * Input: Approximate initial data $\\{n_{j}^{0}\\}_{j\in\mathbb{N}}$

  * $X^{0}\leftarrow 0,\quad N^{0}\leftarrow\Delta s\sum_{j\in\mathbb{N}}p(s_{j},X^{0})n_{j}^{0}.$

  * Choose $\Delta t$ such that

$\Delta t<\min\left\\{\left(\frac{1}{\Delta
s}+\|p\|_{\infty}\right)^{-1},\frac{2}{\alpha_{0}\|\partial_{X}p\|_{\infty}\|n^{0}\|_{1}}\right\\}$
(44)

  * For $m\in\mathbb{N}_{0}$

  * Do $j\in\mathbb{N}$,

$n_{j}^{m+1}\leftarrow\begin{cases}n_{1}^{m}-\frac{\Delta t}{\Delta
s}(n_{1}^{m}-N^{m})-\Delta tp(s_{1},N^{m})n_{1}^{m}&j=1\\\
n_{j}^{m}-\frac{\Delta t}{\Delta s}(n_{j}^{m}-n_{j-1}^{m})-\Delta
tp(s_{j},N^{m})n_{j}^{m}&j>1\end{cases}$

  * end

  * Solve for $X^{m+1}$

$X^{m+1}=\frac{\Delta t}{2}\left(\Delta
s\sum_{j\in\mathbb{N}}p(s_{j},X^{m+1})n_{j}^{m+1}\alpha_{0}+\sum_{k=0}^{m}N^{k}\alpha_{m-k}.\right)$
(45)

  * $N^{m+1}\leftarrow\Delta s\sum_{j\in\mathbb{N}}p(s_{j},X^{m+1})n_{j}^{m+1}$

  * end

  * Output: approximate solution $\\{n_{j}^{m+1}\\}_{j\in\mathbb{N}}$ and $X^{m+1},\,N^{m+1}$ at time $t^{m+1}=(m+1)\Delta t$.

Analogously to the numerical scheme for the ITM equation. The solution of
Equation (45) for $X^{m+1}$ can be solved with different numerical methods.
Unlike the ITM equation, there no restriction on the rate $p$ to have unique
solution of Equation (10). Hence, by following the idea of Lemma 3.1 and the
contraction principle, the solution of Equation (23) $N^{m+1}$ can be
approximated in terms of $N^{m}$ through the following formula if $\Delta t$
is small enough:

$X^{m+1}=\frac{\Delta t}{2}\left(\Delta
s\sum_{j\in\mathbb{N}}p(s_{j},X^{m})n_{j}^{m+1}\alpha_{0}+\sum_{k=0}^{m}N^{k}\alpha_{m-k}.\right)$

For simplicity in the estimates we assume that we can compute the solution of
Equation (45) exactly, but the results remain valid if we take into account an
specific method to get an approximation.

In order to prove the convergence of the upwind scheme we follow the ideas of
the previous subsection on well-posedness and we prove a BV-estimate that will
be the crucial in the analysis. For simplicity we assume that initial data
$n^{0}$ is compactly supported, but the theoretical results still hold when
the initial data $n^{0}$ vanishes at infinity.

As in the case of the ITM equation, we get the corresponding lemmas on $L^{1}$
and $L^{\infty}$ norms.

######  Lemma 3.2.

($L^{1}$-norm) Numerical approximation obtained with the Algorithm 3.1
satisfies

$\|n^{m}\|_{1}:=\sum_{j\in\mathbb{N}}\Delta sn_{j}^{m}=\|n^{0}\|_{1},\quad
m\in\mathbb{N}.$

###### Proof.

The proof is the same as Lemma 2.3. ∎

######  Lemma 3.3.

($L^{\infty}$-norm) Assume that $n^{0}\in L^{\infty}(\mathbb{R}^{+})$ is non-
negative, then under the condition (44) Equation (45) has a unique solution
and the numerical solution obtained by Algorithm 3.1, satisfies the following
estimates

$\displaystyle 0\leq n_{j}^{m}\leq\|n^{0}\|_{\infty},$ $\displaystyle 0\leq
X^{m}\leq\|p\|_{\infty}\|n^{0}\|_{1}\|\alpha\|_{1},$ $\displaystyle 0\leq
N^{m}\leq\|p\|_{\infty}\|n^{0}\|_{1}\quad\textrm{for all}\>j,m\in\mathbb{N}.$

###### Proof.

First, observe that $X^{0}=0$ and $N^{0}=\Delta
s\sum_{j\in\mathbb{N}}p(s_{j},0)n_{j}^{0}\leq\|p\|_{\infty}$. Now, using the
CFL condition we have

$0\leq n_{j}^{1}=n_{j}^{0}\left(1-\Delta t\left(\frac{1}{\Delta
s}+p(s_{j},0)\right)\right)+\frac{\Delta t}{\Delta
s}n_{j-1}^{0}\leq\|n^{0}\|_{\infty},$

for $j\geq 1$. In order to proof existence of $X^{1}$ which satisfy equation
(43), we consider the function

$F(X)=\frac{\Delta t}{2}\left(\Delta
s\sum_{j\in\mathbb{N}}p(s_{j},X)n_{j}^{1}\alpha_{0}+N^{0}\alpha_{1}.\right),$

and observe that

$\displaystyle 0\leq F(X)\leq\|p\|_{\infty}\|n^{0}\|_{1}\|\alpha\|_{1},$
$\displaystyle\text{for all }X\geq 0,$
$\displaystyle|F(X_{1})-F(X_{2})|\leq\frac{\Delta
t}{2}\|n^{0}\|_{1}\alpha_{0}\|\partial_{X}p\|_{\infty}|X_{1}-X_{2}|,$
$\displaystyle\text{for all }X_{1},X_{2}\geq 0.$

From Condition (44) we have

$\frac{\Delta t}{2}\alpha_{0}\|\partial_{X}p\|_{\infty}\|n^{0}\|_{1}<1,$

so that from contraction principle that there exists a unique $X^{1}$ such
that

$X^{1}=F(X^{1})=\frac{\Delta t}{2}\left(\Delta
s\sum_{j\in\mathbb{N}}p(s_{j},X^{1})n_{j}^{1}\alpha_{0}+N^{0}\alpha_{1}.\right)$

and we have that $0\leq X^{1}\leq\|p\|_{\infty}\|n^{0}\|_{1}\|\alpha\|_{1}$.
Moreover we get the following estimate

$N^{1}:=\Delta
s\sum_{j\in\mathbb{N}}p(s_{j},X^{1})n_{j}^{1}\leq\|p\|_{\infty}\|n^{0}\|_{1},$

and we conclude the desired result by iterating this argument for all
$m\in\mathbb{N}$. ∎

We note that in Condition (44), the first term in right-hand side corresponds
to the CFL condition of the explicit scheme and the second term ensures that
the Equation (45) has a unique solution so that $X^{m}$ is well-defined.

Next, we proceed with the corresponding BV-estimate that gives the necessary
compactness to prove the convergence of the numerical scheme.

######  Lemma 3.4.

(BV-estimate) Assume that $TV(n^{0})<\infty$ and the CFL condition (22). Then
there exist constants $C_{1},C_{2}>0$ (depending only $p,\alpha$ and the norms
of $n^{0}$) such that for $m\in\mathbb{N}$ we have

$TV(n^{m})\leq e^{C_{1}T}TV(n^{0})+C_{2}(e^{C_{1}T}-1),$ (46)

with $T=m\Delta t$ and
$TV(n^{m})=\displaystyle\sum_{j=0}^{\infty}|n_{j+1}^{m}-n_{j}^{m}|$.

###### Proof.

using notation $\Delta^{+}n_{j}^{m}=n_{j+1}^{m}-n_{j}^{m}$, we have

$\begin{split}\Delta^{+}n_{j}^{m+1}&=\Delta^{+}n_{j}^{m}-\frac{\Delta
t}{\Delta s}\left(\Delta^{+}n_{j}^{m}-\Delta^{+}n_{j-1}^{m}\right)-\Delta
tp(s_{j+1},X^{m})\Delta^{+}n_{j}^{m}\\\ &\qquad-\Delta
tn_{j}^{m}(p(s_{j+1},X^{m})-p(s_{j},X^{m})).\end{split}$

Next, from the the CFL condition (22) we obtain

$|\Delta^{+}n_{j}^{m+1}|\leq\left(1-\Delta t\left(\frac{1}{\Delta
s}+p(s_{j+1},X^{m})\right)\right)|\Delta^{+}n_{j}^{m}|+\frac{\Delta t}{\Delta
s}|\Delta^{+}n_{j-1}^{m}|+\Delta t\Delta
s\|\partial_{s}p\|_{\infty}n_{j}^{m}.$

By summing over all $j\geq 1$ we deduce that

$\sum_{j=1}^{\infty}|\Delta^{+}n_{j}^{m+1}|\leq\sum_{j=1}^{\infty}|\Delta^{+}n_{j}^{m}|+\frac{\Delta
t}{\Delta s}|\Delta^{+}n_{0}^{m}|+\Delta
t\|\partial_{s}p\|_{\infty}\|n^{0}\|_{1}.$ (47)

We now take into account the boundary term. From the numerical scheme we have

$\begin{split}|\Delta^{+}n_{0}^{m+1}|&=|n_{1}^{m+1}-N^{m+1}|\\\
&\leq\left(1-\frac{\Delta t}{\Delta
s}\right)|n_{1}^{m}-N^{m}|+|N^{m+1}-N^{m}|+\Delta tp(s_{1},N^{m})n^{m}_{1}\\\
&\leq\left(1-\frac{\Delta t}{\Delta
s}\right)|\Delta^{+}n_{0}^{m}|+|N^{m+1}-N^{m}|+\Delta
t\|p\|_{\infty}\|n^{0}\|_{\infty}.\end{split}$

Next we estimate the second term of the last inequality. First note that

$\begin{split}N^{m+1}-N^{m}&=\Delta
s\sum_{j\in\mathbb{N}}p(s_{j},X^{m+1})n_{j}^{m+1}-\Delta
s\sum_{j\in\mathbb{N}}p(s_{j},X^{m})n_{j}^{m}\\\ &=\Delta
s\sum_{j\in\mathbb{N}}(p(s_{j},X^{m+1})-p(s_{j},X^{m}))n_{j}^{m+1}+\Delta
s\sum_{j\in\mathbb{N}}p(s_{j},X^{m})(n_{j}^{m+1}-n_{j}^{m}),\end{split}$

thus we have

$\begin{split}|N^{m+1}-N^{m}|&\leq\|\partial_{X}p\|_{\infty}\|n^{0}\|_{1}|X^{m+1}-X^{m}|+\|p\|_{\infty}\sum_{j\in\mathbb{N}}\Delta
s|n^{m+1}_{j}-n^{m}_{j}|\\\
&\leq\|\partial_{X}p\|_{\infty}\|n^{0}\|_{1}|X^{m+1}-X^{m}|+\|p\|_{\infty}\Delta
t\left(\sum_{j\in\mathbb{N}}|n_{j}^{m}-n_{j-1}^{m}|+\|p\|_{\infty}\sum_{j\in\mathbb{N}}\Delta
sn_{j}^{m}\right)\\\
&\leq\|\partial_{X}p\|_{\infty}\|n^{0}\|_{1}|X^{m+1}-X^{m}|+\|p\|_{\infty}\Delta
t\sum_{j\in\mathbb{N}}|\Delta^{+}n^{m}_{j-1}|+\Delta
t\|p\|_{\infty}^{2}\|n^{0}\|_{1}\end{split}$ (48)

And similarly for $|X^{m+1}-X^{m}|$ we get by using Equation (43)

$\begin{split}|X^{m+1}-X^{m}|&\leq\frac{\Delta
t}{2}\alpha_{0}\|\partial_{X}p\|_{\infty}\|n^{0}\|_{1}|X^{m+1}-X^{m}|+\frac{\Delta
t}{2}\alpha_{0}\|p\|_{\infty}\sum_{j\in\mathbb{N}}\Delta
s|n^{m+1}_{j}-n^{m}_{j}|\\\ &\qquad+\frac{\Delta
t}{2}\left|\sum_{k=0}^{m}N^{k}\alpha_{m+1-k}-\sum_{k=0}^{m-1}N^{k}\alpha_{m-k}\right|,\end{split}$

so we get

$\begin{split}|X^{m+1}-X^{m}|&\leq\frac{\Delta
t}{2}\alpha_{0}\left(\|\partial_{X}p\|_{\infty}\|n^{0}\|_{1}|X^{m+1}-X^{m}|+\Delta
t\|p\|_{\infty}\sum_{j\in\mathbb{N}}|\Delta^{+}n^{m}_{j-1}|+\Delta
t\|p\|_{\infty}^{2}\|n^{0}\|_{1}\right)\\\ &\qquad+\frac{\Delta
t}{2}N^{m}\alpha_{1}+\frac{\Delta
t}{2}\|p\|_{\infty}\|n^{0}\|_{1}\sum_{k=0}^{m}|\alpha_{m+1-k}-\alpha_{m-k}|,\end{split}$

and therefore we have the following estimate

$\begin{split}|X^{m+1}-X^{m}|&\leq\frac{\Delta t}{2(1-\frac{\Delta
t}{2}\alpha_{0}\|\partial_{X}p\|_{\infty}\|n^{0}\|_{1})}\left(\Delta
t\alpha_{0}\|p\|_{\infty}\sum_{j\in\mathbb{N}}|\Delta^{+}n^{m}_{j-1}|+\Delta
t\alpha_{0}\|p\|_{\infty}^{2}\|n^{0}\|_{1}\right.\\\
&\left.\phantom{\sum_{j\in\mathbb{N}}}+\|p\|_{\infty}\|n^{0}\|_{1}\|\alpha\|_{\infty}+\|p\|_{\infty}\|n^{0}\|_{1}TV(\alpha)\right).\end{split}$
(49)

By plugging (49) into (48) we obtain for $\Delta t$ small

$|N^{m+1}-N^{m}|\leq A_{1}\Delta
t\sum_{j\in\mathbb{N}}|\Delta^{+}n^{m}_{j-1}|+B_{1}\Delta t,$

where $A_{1},B_{1}>0$ are constants depending on $p,\alpha$ and the norms of
$n^{0}$. So that the boundary term is estimated as follows

$\begin{split}|\Delta^{+}n_{0}^{m+1}|&\leq\left(1-\frac{\Delta t}{\Delta
s}\right)|\Delta^{+}n_{0}^{m}|+A_{1}\Delta
t\sum_{j\in\mathbb{N}}|\Delta^{+}n^{m}_{j-1}|+B_{2}\Delta t,\end{split}$ (50)

and by adding (50) with (47), we obtain

$\sum_{j=0}^{\infty}|\Delta^{+}n_{j}^{m+1}|\leq(1+C_{1}\Delta
t)\sum_{j=0}^{\infty}|\Delta^{+}n_{j}^{m}|+C_{2}\Delta t,$ (51)

where $C_{1},C_{2}>0$ constants depending only on $p,\alpha$ and the norms of
$n^{0}$.

Finally, proceeding recursively on $m$, we obtain

$\sum_{j=0}^{\infty}|\Delta^{+}n_{j}^{m}|\leq(1+C_{1}\Delta
t)^{m}\sum_{j=0}^{\infty}|\Delta^{+}n_{j}^{0}|+\frac{C_{2}}{C_{1}}\left((1+C_{1}\Delta
t)^{m}-1\right),$

and the estimate (46) readily follows. ∎

As we did in Section 2 for the ITM equation. We now prove that the numerical
approximation of the solution of Equation (10) $n(t,s)$, which is constructed
by a simple piece-wise linear interpolation, has a limit when the time step
$\Delta t$ and age step $\Delta s$ converge to $0$.

######  Lemma 3.5.

Assume that $n^{0}\in BV(\mathbb{R}^{+})$ is compactly supported and the rate
$p$ satisfies the hypothesis of Theorem 3.1. Consider the function $n_{\Delta
t,\Delta s}(t,s)\in\mathcal{C}\left([0,T],L^{1}(\mathbb{R}^{+})\right)$
defined by

$n_{\Delta t,\Delta s}(t,s)\coloneqq\frac{t^{m}-t}{\Delta
t}\sum_{j\in\mathbb{N}}n_{j}^{m-1}\chi_{[s_{j-\frac{1}{2}},s_{j+\frac{1}{2}}]}(s)+\frac{t-t^{m-1}}{\Delta
t}\sum_{j\in\mathbb{N}}n_{j}^{m}\chi_{[s_{j-\frac{1}{2}},s_{j+\frac{1}{2}}]}(s)\;\textrm{if}\;t\in[t^{m-1},t^{m}].$

Then there exists a sub-sequence $(\Delta t_{k},\Delta s_{k})\to(0,0)$ when
$k\to\infty$ and a function $\overline{n}(t,s)$ such that $n_{\Delta
t_{k},\Delta s_{k}}\to\overline{n}$ in
$\mathcal{C}\left([0,T],L^{1}(\mathbb{R}^{+})\right)$. Moreover, if we define
the function $X_{\Delta t,\Delta s}(t)\in\mathcal{C}[0,T]$ as the unique
solution of the integral equation

$X(t)=\int_{0}^{t}\int_{0}^{\infty}\alpha(t-\tau)p(s,X(\tau))n_{\Delta
t,\Delta s}(\tau,s)\,\mathrm{d}s\,\mathrm{d}\tau,$

then there exists $\overline{X}\in\mathcal{C}[0,T]$ such that $X_{\Delta
t_{k},\Delta s_{k}}\to\overline{X}$ in $\mathcal{C}[0,T]$ and $\overline{X}$
is a solution of the equation

$\overline{X}(t)=\int_{0}^{t}\int_{0}^{\infty}\alpha(t-\tau)p(s,\overline{X}(\tau))\overline{n}(\tau,s)\,\mathrm{d}s\,\mathrm{d}\tau,$
(52)

and similarly $N_{\Delta t,\Delta s}\in\mathcal{C}[0,T]$ defined as

$N_{\Delta t,\Delta s}(t)=\int_{0}^{\infty}p\left(s,X_{\Delta t,\Delta
s}(t)\right)n_{\Delta t,\Delta s}\,\mathrm{d}s$

converges respectively to $\overline{N}\in\mathcal{C}[0,T]$, where
$\overline{N}$ satisfies the equality

$\overline{N}(t)=\int_{0}^{\infty}p(s,\overline{X}(t))\overline{n}(t,s)\,\mathrm{d}s.$
(53)

###### Proof.

The proof on the compactness of $n_{\Delta t,\Delta s}$ is the same as Lemma
2.7 by using Lemma 3.4. Hence there exists a sub-sequence $(\Delta
t_{k},\Delta s_{k})\to(0,0)$ when $k\to\infty$ and a function
$\overline{n}(t,s)$ such that $n_{\Delta t_{k},\Delta s_{k}}\to\overline{n}$
in $\mathcal{C}\left([0,T],L^{1}(\mathbb{R}^{+})\right)$. For the sequence
$X_{\Delta t,\Delta s}$ observe that

$\|X_{\Delta t,\Delta
s}\|_{\infty}\leq\|\alpha\|_{1}\|p\|_{\infty}\|n^{0}\|_{1}$

and for the derivative we get

$\begin{split}\frac{d}{dt}X_{\Delta t,\Delta
s}(t)&=\int_{0}^{\infty}\alpha(0)p(s,X_{\Delta t,\Delta s}(t))n_{\Delta
t,\Delta s}(\tau,s)\,\mathrm{d}s\,\mathrm{d}\tau\\\
&\qquad+\int_{0}^{t}\int_{0}^{\infty}\alpha^{\prime}(t-\tau)p(s,X_{\Delta
t,\Delta s}(\tau))n_{\Delta t,\Delta
s}(\tau,s)\,\mathrm{d}s\,\mathrm{d}\tau,\end{split}$

thus we have the following estimate

$\left\|\frac{d}{dt}X_{\Delta t,\Delta
s}\right\|_{\infty}\leq\alpha(0)\|p\|_{\infty}\|n^{0}\|_{1}+T\|\alpha^{\prime}\|_{\infty}\|p\|_{\infty}\|n^{0}\|_{1},$

and from Arzelà-Ascoli theorem we conclude that $X_{\Delta t,\Delta s}$
converges to some $\overline{X}$ in $\mathcal{C}[0,T]$ for a sub-sequence. By
passing to the limit, Equations (52) and (53) readily follow. ∎

With this result, we finally we get the result on convergence of the numerical
scheme for the DDM equation (10).

######  Theorem 3.2 (Convergence of the numerical scheme).

Assume that $n^{0}\in BV(\mathbb{R}^{+})$ is compactly supported and the rate
$p$ satisfies the hypothesis of Theorem 3.1. Then for all $T>0$, the numerical
scheme converges to the unique weak solution
$n\in\mathcal{C}([0,T],L^{1}(\mathbb{R}^{+}))$ of the DDM equation (10).

###### Proof.

The proof is the same as Theorem 2.2. ∎

######  Remark 3.2.

The previous results are also valid for the case when the rate $p$ is of the
form

$p(s,X)=\varphi(X)\chi_{\\{s>\sigma(X)\\}},$

with $\varphi$ and $\sigma$ Lipschitz bounded functions.

##  4 Numerical Results

In order to illustrate the theoretical results of the previous sections, we
present in different scenarios for the dynamics of the ITM equation (1) and
DDM equation (10) which are solved by the finite volume method described in
Algorithms (2.1) and (3.1) respectively, where the non-linear problems (23)
and (45) where solved for $N$ or $X$ using the Newton-Raphson iterative method
with a relative error less than $10^{-12}$. For all numerical tests, we
consider the prototypical rate $p$ with absolute refractory period $\sigma>0$
given by

$p(s,N)=\varphi(N)\chi_{\\{s>\sigma\\}}(s),$

where $\varphi(N)$ and $\sigma$ are specified in each example. For the DDM
equation we will consider for the delay kernel $\alpha(t)$ the following
examples

$\alpha_{1}(t)=\frac{e^{-t/\lambda}}{\lambda}\quad\text{or}\quad\alpha_{2}(t)=\frac{1}{\sqrt{2\pi}\lambda}e^{-\frac{1}{2}(\frac{t-d}{\lambda})^{2}},\quad\text{with
}\lambda=10^{-3}.$

For this choice of $\lambda$ we essentially consider the approximation
$\alpha_{1}(t)\approx\delta(t)$, where we are interested in comparing both ITM
and DDM equations when we are close to this limit case. Similarly for the
second kernel we get $\alpha_{2}(t)\approx\delta(t-d)$ and we study the the
behavior of DDM equation (10) with different values of the parameter $d$.

### Example 1: A strongly inhibitory case

We start with an inhibitory hazard rate, i.e. $\varphi^{\prime}(N)<0$, given
by

$\varphi(N)=e^{-9N},\quad\sigma=\frac{1}{2},\quad
n^{0}(s)=\frac{1}{2}e^{-(s-1)^{+}}.$ (54)

With these choice of parameters Equation (1) has a unique steady state with
$N^{*}\approx 0.1800$, by solving Equation (15). For the ITM equation we
display in Figures 1(a-b) the numerical solution for $n(t,s)$ with
$(t,s)\in[0,30]\times[0,40]$ and $N(t)$ with $t\in[0,10]$, where we observe
that the total activity $N(t)$ converge to the unique steady state $N^{*}$,
while for $n(t,s)$ the initial condition moves to the right (with respect to
age $s$) as it decays exponentially and approaches to the equilibrium density
given by Equation 6. Moreover, this example is clearly consistent with the
theory on the convergence to the equilibrium for the strongly excitatory case
studied in [14].

(a) | (b)
---|---
|
(c) | (d)
|
Figure 1: An inhibitory case. (a-b) Density $n(t,s)$ and discharging flux
$N(t)$ for ITM. (c-d) Density $n(t,s)$, discharging flux $N(t)$ and total
activity $X(t)$ for DDM with $\alpha_{2}(t)\approx\delta(t-\frac{1}{2})$.

On the other hand, we solve the DDM (10) with the parameters in (54) and we
choose $d=\frac{1}{2}$ so that $\alpha_{2}(t)\approx\delta(t-\frac{1}{2})$. In
Figs 1(c,d) we display numerical solution for $n(t,s)$ for
$(t,s)\in[0,15]\times[0,20]$, and $N(t)$ and $X(t)$ for $t\in[0,20]$. Unlike
the ITM, for the DDM equation the solutions for $N$ and $X$ converge to a
periodic profile that we conjecture to be of $2d$-periodic and they tend to
differ by a period of time equal to $d$, which means that
$|N(t-d)-X(t)|\approx 0$ for $t$ large. The density $n(t,s)$ is asymptotic to
the its respective periodic profile in time. In this case, we observe a
periodic solution induced by a negative feedback delay, which is a classical
behavior in the context of delay differential equations. This negative
feedback corresponds to the inhibition determined by the rate $p$, so that
when combined with the delay it induces cycles of increase and decrease in the
discharging flux $N$ and total activity $X$, which are favorable to form
periodic solutions (see [27] for a reference).

### Example 2: An excitatory case with a unique steady state

Now we consider an excitatory case, i.e. $\varphi^{\prime}(N)>0$, given by

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

This example was previously studied in [14] for the ITM equation. We know that
under this choice of parameters Equation (1) has a unique steady state with
$N^{*}\approx 0.8186$.

For the ITM equation, in Figure 2(a) we display the numerical solution for
$n(t,s)$ for $(t,s)\in[0,20]\times[0,4]$ and in the blue curve of Figure 2(b)
we show $N(t)$ for $t\in[0,20]$ where the solution is asymptotic periodic
pattern with jump discontinuities. This is due to the invertibility condition
$\Psi(N,n)$ in (13) is close to zero as we show in Figure 2(c), where we plot
$\Psi(N(t),n(t,\cdot))$ for $t\in[0,10]$. We observe that when a discontinuity
arises for $N(t)$ in Figure 2(b), then $\Psi(N,n)$ is close to zero in Figure
2(c). This means that the invertibility condition (13) is a key criterion that
determine the existence and continuity of solutions.

For the DDM equation with $\alpha_{1}(t)\approx\delta(t)$ we observe in the
red curve of Figure 2(b) that the respective discharging flux is a smooth
approximation of the activity of $N(t)$ for the ITM equation. This is due to
the regularizing effect of the delay kernel $\alpha$ through the convolution.

(a) | (b)
---|---
|
(c)
Figure 2: An excitatory case with periodic patterns. (a) Density $n(t,s)$ for
ITM, (b) Comparison of $N(t)$ for both ITM and DDM equations with
$\alpha_{1}(t)\approx\delta(t)$, (c) Invertibility condition $\Psi(N,n)$ for
ITM.

Next we take $d=1$ so that $\alpha_{2}(t)\approx\delta(t-1)$ In Figure 3(a) we
display $n(t,s)$ for $(t,s)\in[0,6]\times[0,15]$ and in Figure 3(b) we display
the graphics of $N(t)$ and $X(t)$ for $t\in[0,15]$, where we observe an
asymptotic periodic pattern for both discharging flux $N$ and total activity
$X$, which we conjecture to be $d$-periodic. In this case we observe a a
synchronization phenomena which means that $|N(t)-X(t)|\approx 0$ for large
$t$, unlike the inhibitory case shown in Fig 1(d) where they tend to differ in
time by $d$. In this excitatory case we conjecture that periodic solutions are
due to effect of the refractory period $\sigma$ as it was studied in [14] and
the solution are continuous due to the regularizing effect of the kernel
$\alpha$. Therefore we observe that periodic solution that may arise in the
inhibitory and excitatory are of different nature.

(a) | (b)
---|---
|
Figure 3: A periodic solution for the DDM equation. (a-b) Density $n(t,s)$,
discharging flux $N(t)$ and total activity $X(t)$ with $d=1$ and
$\alpha_{2}(t)\approx\delta(t-1),\,X(t)\approx N(t-1)$.

### Example 3: An excitatory case with multiples steady states

Next, we consider a case where $\varphi^{\prime}(N)>0$ with parameters given
by

$\varphi(N)=\frac{1}{1+e^{-9N+3.5}},\quad\sigma=\frac{1}{2},\quad
n^{0}(s)=e^{-(s-\frac{1}{2})^{+}}\chi_{\\{s>\frac{1}{2}\\}}(s).$ (56)

This example was previously studied in [14] for the ITM equation. Under this
choice of parameters, we have three different solutions for $N(0)$ according
to Equation (4), that are given by $N^{0}_{1}\approx 0.0410$,
$N^{0}_{2}\approx 0.3650$ and $N^{0}_{3}\approx 0.6118$. These values
determine three different branches of local continuous solutions. For the ITM
equation we display the numerical solution for $n(t,s)$ and $N(t)$ in Figures
4. The dynamics for the discharging flux $N$ is determined by the initial
condition $N(0)$ and thus it also determines the dynamics for $n$. In this
case, we observe three different types of numerical approach for $N(t)$, which
converge to two different equilibrium points $N^{*}_{1}$ and $N^{*}_{2}$.

(a) | (b)
---|---
|
(c) | (d)
|
(e) | (f)
|
Figure 4: An excitatory problems with multiples solutions for the ITM equation
with different initial approximations for $N^{0}$ in (23). (a-b) Density
$n(t,s)$ and discharging flux $N(t)$ for $N^{0}_{1}\approx 0.0281$, (c-d) for
$N^{0}_{2}\approx 0.4089$ and (e-f) for $N^{0}_{3}\approx 0.7114$

In the DDM equation for both $\alpha_{1}(t)\approx\delta(t)$ and
$\alpha_{2}(t)\approx\delta(t-d)$, we observe in Figure 5 that $N(t)$ converge
to the first equilibrium of the ITM equation. Moreover, when the ITM equation
has multiple branches of solutions for the same initial condition, i.e.
multiples solutions for $N(0)$, we conjecture that when $\alpha(t)$ converges
$\delta(t)$ in the sense of distributions, the total activity of $X(t)$ in the
DDM equation converges a.e. to the solution of $N(t)$ in the ITM equation,
whose value of $N(0)$ is the closest one to zero and we expect
$L^{1}$-convergence for the corresponding probability densities. This is due
to the condition $X(0)=0$ imposed for the DDM equation. We observe in Figure
5(b) that $X(t)$ and $N(t)$ follow the same behavior of the Figure 4(b) and in
particular the total activity $X(t)$ grows fast from $X(0)=0$ until it
approaches to the solution of $N(t)$.

Similarly when $\alpha(t)$ converges to $\delta(t-d)$, we conjecture that
total activity $X(t)$ in the DDM equation converges to the solution of $N(t)$
of Equation 12 that satisfies $N(t)\equiv 0$ for $t\in[-d,0]$, as it is
suggested by the numerical solution in Figure 5(d) since we would formally get
$X(t)=N(t-d)$ and $X(t)=0$ for $t\in[0,d]$.

(a) | (b)
---|---
|
(c) | (d)
|
Figure 5: An excitatory case for the DDM equation. (a-b) Density $n(t,s)$,
discharging flux $N(t)$ and total activity $X(t)$ for the DDM equation with
$\alpha(t)\approx\delta(t)$. (c-d) Same variables of the system with
$\alpha(t)\approx\delta(t-1)$.

### Example 4: A variable refractory period [8]

Based on Example 2 in [8], we consider a hazard rate with variable refractory
period as in Equation (9) for the the DDM equation with parameters given by

$p(s,X)=\chi_{\\{s>\sigma(X)\\}}(s),\quad\sigma(X)=2-\frac{X^{4}}{X^{4}+1},\quad\alpha(t)=J\alpha_{1}(t),\quad
n^{0}(s)=e^{-(s-1)}\chi_{\\{s>1\\}}(s),$ (57)

where $J>0$ is the connectivity parameter of the network. As it was studied in
[8], the system has different behaviors depending value of $J$. When $J$ is
small the network is weakly connected and the dynamics are close to the linear
case, while if $J$ is large the network is strongly connected and different
asymptotic behaviors are possible. We recall that when $\alpha(t)=\delta_{0}$,
we formally obtain the ITM equation where

$X(t)=JN(t),\quad\text{and},\quad p(s,N)=\chi_{\\{s>\sigma(JN)\\}}(s).$ (58)

Taking $J=2.5$ as in [8], in Figure 6 we compare the numerical approximation
of $N(t)$ with $t\in[0,14]$, for both ITM and DDM equations. In Figure 6(a) we
observe that the solution of the ITM equation is asymptotic to a periodic
pattern with jump discontinuities, where this type of solutions were also
observed in [14]. We also observe that when a jump discontinuity arises for
$N(t)$ in Figure 6(a), the function $\Psi(N,n)$ is close to zero as we see in
Figure 6(b), verifying numerically the invertibility condition (13) that
ensures the continuity of solutions.

In Figure 6(c) we observe that the discharging flux $N(t)$ in the DDM equation
is a smooth approximation of the solution observed in Figure 6(a) and we see
the same phenomenon for $\tilde{X}(t)\coloneqq X(t)/J$, corresponding to a
normalization of the total activity in order to the compare these quantities.
Finally in Figure 6(d) we display the corresponding numerical approximation of
the DDM equation with $n(t,s)$ for $(t,s)\in[0,14]\times[0,8]$, which also
follows a periodic pattern.

(a) | (b)
---|---
|
(c) | (d)
|
Figure 6: Hazard rate with variable refractory period. Comparing numerical
approximation of $N(t)$ for both ITM and DDM equations. (a) Discharging flux
$N(t)$ for ITM equation, (b) Invertibility condition $\Psi(N,n)$ for ITM. (c)
$N(t)$ and $X(t)/J$ for DDM equation, (d) Density $n(t,s)$ for DDM.

## Conclusion and perspectives

In this article we managed to improve proof of [8, 11] on well-posedness for
both ITM and DDM equations, allowing to extend the theory for wider types of
hazard rates $p$. The key idea is to apply the implicit function theorem to
the correct fixed point problem and the arguments can be extended when this
rate is not necessarily bounded. This motivates the study of elapsed time
model when the activity of neurons may increase to infinity and blow-up or
other special phenomena might arise.

Another interesting question is convergence of the delay kernel $\alpha(t)$ to
the Dirac’s mass $\delta(t)$ in order to compare both ITM and DDM equations.
We conjecture that the total activity $X(t)$ of the DDM equation converges
almost everywhere (or in some norm) to the discharging flux $N(t)$ of the ITM
equation. In particular we believe that the convergence holds for every $t>0$
except when $N(t)$ has a jump discontinuity. This motivates to determine if
the assumption of instantaneous transmission is actually a good approximation
of the neural dynamics, which have indeed a certain delay. Similarly, when the
delay kernel $\alpha(t)$ converges to the Dirac’s mass $\delta(t-d)$ in the
sense of distributions, we conjecture that solutions of the DDM equations
converge to the solutions of Equation (12).

From a numerical point of view, we proved the convergence of the explicit
upwind scheme for the elapse time model relying on the mass-preserving
property, the analysis of the fixed point equations (23), (45) and the key BV-
estimate, from where we obtain the compactness to conclude the result. We can
extend the analysis of the elapsed time model by considering implicit or semi-
discrete schemes, but a more detailed analysis on the mass conservation the
estimates must be considered. Other possible discretizations to solve the
equations include for the example high-order Runge-Kutta-WENO methods (see for
example [24, 25, 28, 26]). These alternatives might be useful to analyse
numerically the time elapsed equation when the rate $p$ is not bounded, which
implies that the total activity may be also unbounded. Furthermore, this
numerical analysis can be considered for other extensions of the elapsed time
equation such as the model with fragmentation [10], spatial dependence [15],
the multiple-renewal equation [16] and the model with leaky memory variable in
[17] or other type of structured equations.

#### Acknowledgements

This work has been supported by ANID project ECOS200018. MS and LMV was
partially supported by ANID-Chile through Centro de Modelamiento Matemático
(FB210005) and INRIA Associated team ANACONDA. MS was supported by Fondecyt-
ANID project 1220869, and Jean d’Alembert fellowship program, Université de
Paris-Saclay. NT was supported by the grant Juan de la Cierva FJC2021-046894-I
funded by MCIN/AEI and the European Union NextGenerationEU/PRTR.

#### ORCID iDs

Mauricio Sepulveda: https://orcid.org/0000-0001-8463-3830.
Nicolas Torres: https://orcid.org/0000-0001-6059-9754.
Luis Miguel Villada: https://orcid.org/0000-0002-4860-4431.

## References

  * [1] Gerstner W, Kistler WM.  Spiking neuron models: Single neurons, populations, plasticity.  Cambridge university press; 2002.
  * [2] Pham J, Pakdaman K, Champagnat J, Vibert JF.  Activity in sparsely connected excitatory neural networks: effect of connectivity.  Neural Networks. 1998;11(3):415-34.
  * [3] Ly C, Tranchina D.  Spike train statistics and dynamics with synaptic input from any renewal process: a population density approach.  Neural Computation. 2009;21(2):360-96.
  * [4] Chevallier J, Cáceres MJ, Doumic M, Reynaud-Bouret P.  Microscopic approach of a time elapsed neural model.  Mathematical Models and Methods in Applied Sciences. 2015;25(14):2669-719.
  * [5] Chevallier J.  Mean-field limit of generalized Hawkes processes.  Stochastic Processes and their Applications. 2017;127(12):3870-912.
  * [6] Quiñinao C.  A microscopic spiking neuronal network for the age-structured model.  Acta Applicandae Mathematicae. 2016;146:29-55.
  * [7] Schwalger T, Chizhov AV.  Mind the last spike—firing rate models for mesoscopic populations of spiking neurons.  Current opinion in neurobiology. 2019;58:155-66.
  * [8] Pakdaman K, Perthame B, Salort D.  Dynamics of a structured neuron population.  Nonlinearity. 2009;23(1):55.
  * [9] Pakdaman K, Perthame B, Salort D.  Relaxation and self-sustained oscillations in the time elapsed neuron network model.  SIAM Journal on Applied Mathematics. 2013;73(3):1260-79.
  * [10] Pakdaman K, Perthame B, Salort D.  Adaptation and fatigue model for neuron networks and large time asymptotics in a nonlinear fragmentation equation.  The Journal of Mathematical Neuroscience. 2014;4:1-26.
  * [11] Cañizo JA, Yoldaş H.  Asymptotic behaviour of neuron population models structured by elapsed-time.  Nonlinearity. 2019;32(2):464.
  * [12] Mischler S, Weng Q.  Relaxation in time elapsed neuron network models in the weak connectivity regime.  Acta Applicandae Mathematicae. 2018;157(1):45-74.
  * [13] Mischler S, Quiñinao C, Weng Q.  Weak and strong connectivity regimes for a general time elapsed neuron network model.  Journal of Statistical Physics. 2018;173(1):77-98.
  * [14] Torres N, Cáceres MJ, Perthame B, Salort D.  An elapsed time model for strongly coupled inhibitory and excitatory neural networks.  Physica D: Nonlinear Phenomena. 2021:132977.
  * [15] Torres N, Salort D.  Dynamics of neural networks with elapsed time model and learning processes.  Acta Applicandae Mathematicae. 2020;170(1):1065-99.
  * [16] Torres N, Perthame B, Salort D.  A multiple time renewal equation for neural assemblies with elapsed time model.  Nonlinearity. 2022;35(10):5051.
  * [17] Fonte C, Schmutz V.  Long Time Behavior of an Age-and Leaky Memory-Structured Neuronal Population Equation.  SIAM Journal on Mathematical Analysis. 2022;54(4):4721-56.
  * [18] Abia L, Angulo O, López-Marcos J.  Age-structured population models and their numerical solution.  Ecological modelling. 2005;188(1):112-36.
  * [19] Sulsky D.  Numerical solution of structured population models.  Journal of Mathematical Biology. 1994;32(5):491-514.
  * [20] Carrillo JA, Gwiazda P, Ulikowska A.  Splitting-particle methods for structured population models: convergence and applications.  Mathematical Models and Methods in Applied Sciences. 2014;24(11):2171-97.
  * [21] Gwiazda P, Lorenz T, Marciniak-Czochra A.  A nonlinear structured population model: Lipschitz continuity of measure-valued solutions with respect to model ingredients.  Journal of Differential Equations. 2010;248(11):2703-35.
  * [22] Ackleh AS, Lyons R, Saintier N.  Finite difference schemes for a structured population model in the space of measures.  Mathematical biosciences and engineering: MBE. 2019;17(1):747-75.
  * [23] Perthame B.  Transport equations in biology.  Frontiers in Mathematics. Basel: Springer Science & Business Media; 2006\.
  * [24] LeVeque RJ, Leveque RJ.  Numerical methods for conservation laws. vol. 132.  Springer; 1992.
  * [25] Godlewski E, Raviart PA.  Numerical approximation of hyperbolic systems of conservation laws. vol. 118.  Springer; 1996.
  * [26] Bouchut F.  Nonlinear stability of finite Volume Methods for hyperbolic conservation laws: And Well-Balanced schemes for sources.  Springer Science & Business Media; 2004.
  * [27] Diekmann O, Van Gils SA, Lunel SMV, Walther HO.  Delay equations: functional-, complex-, and nonlinear analysis. vol. 110\.  Springer Science & Business Media; 2012.
  * [28] Cockburn B, Shu CW, Johnson C, Tadmor E, Shu CW.  Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws.  Springer; 1998.