Well-posedness and numerical analysis of an elapsed time model with
  strongly coupled neural networks
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
nonlinearity, 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.
Mauricio Sepulveda, Nicolas Torres, Luis Miguel Villada
Well-posedness and numerical analysis of an elapsed time model with strongly coupled neural networks

###### 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

\[\text{(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} \tag{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, \tag{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, \tag{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. \tag{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,\\ N=n(s=0)\coloneqq\int_{0}^{\infty}p(s,N)n(s)\,\mathrm{d}s,\\ \int_{0}^{\infty}n(s)\,\mathrm{d}s=1,\quad n(s)\geq 0.\end{cases} \tag{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}. \tag{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}. \tag{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\}}, \tag{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)\}}. \tag{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

\[\text{(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} \tag{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} \tag{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
\[\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} \tag{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 Canizo 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).
Instantaneous Transmission Model (ITM)

### 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], \tag{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} \tag{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, \tag{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 \tag{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

\[A_{1} = \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\] \[A_{2} = \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\] \[A_{3} = \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

\[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}.\]
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}. \tag{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}, \tag{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\). 

### 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}, \tag{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}. \tag{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_{j}^{0}\}_{j\in\mathbb{N}}\)_

_Solve for \(N^{0}\)_

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

_Choose \(\Delta t\) such that_

\[\Delta t\leq\left(\frac{1}{\Delta s}+\|p\|_{\infty}\right)^{-1}. \tag{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}. \tag{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_

\[0\leq n_{j}^{m}\leq\|n^{0}\|_{\infty},\ \ 0\leq N^{m}\leq\|p\|_{\infty}\|n^{0} \|_{1}\ \ \text{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}\ \ \ \text{for all}\,j\in\mathbb{N},\]

and by induction in \(m\) we conclude that \(0\leq n_{j}^{m}\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_{j}^{m}=\|p\|_{ \infty}\sum_{j\in\mathbb{N}}\Delta sn_{j}^{0}=\|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_{j}^{0})_{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_{j}^{0},\]

_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_{j}^{1}-n_{j}^{2}| \tag{24}\]

_for sequences \(n^{1}=(n_{j}^{1})_{j\in\mathbb{N}}\), \(n^{2}=(n_{j}^{2})_{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, \tag{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):=\sum_{j=0}^{\infty}|n_{j+1}-n_{j}|. \tag{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), \tag{27}\]

_with \(T=m\Delta t\) and \(TV(n^{m})=\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

\[\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 t\left(p(s_{j+1},N^{m})n_{j+1}^{m}-p(s_ {j},N^{m})n_{j}^{m}\right)\] \[= \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}\] \[-\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

\[\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}|\]

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}. \tag{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_{1}^{m}\\ &\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} \tag{29}\]

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

\[|N^{m+1}-N^{m}| = \left|\psi(n^{m+1})-\psi(n^{m})\right| \tag{30}\] \[\leq C_{1}\sum_{j\in\mathbb{N}}\Delta s|n_{j}^{m+1}-n_{j}^{m}|\] \[\leq 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_{j-1}^{m}|+\|p\|_{\infty}\|n^{0}\|_{1}\right)\\ &\qquad+\Delta t\|p\|_{\infty}\|n^{0}\|_{\infty}.\end{split} \tag{31}\]

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

\[\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. \tag{32}\]

with \(C_{2}:=\|\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_
\[\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_{j}^{m}\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_{j_{m}}^{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},\]

_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)\;\text{if}\;t\in[t^{m-1},t^{m}]. \tag{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. \tag{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

\[\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}\]
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

\[\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_{j}^{m}-n _{j}^{m-1}|}{\Delta t}\] \[\leq|t_{1}-t_{2}|\left(\sum_{j\in\mathbb{N}}|n_{j}^{m-1}-n_{j-1} ^{m-1}|+\Delta sp(s_{j},N^{m-1})n_{j}^{m-1}\right)\] \[\leq|t_{1}-t_{2}|\left(C_{T}\,TV(n^{0})+\|p\|_{\infty}\|n^{0}\|_{ 1}\right), \tag{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_{j}^{0}=0\) for \(j\geq K\) and \(s_{K}=(K-\frac{1}{2})\Delta s\geq R\). From the numerical scheme we deduce that \(n_{j}^{m}=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

\[|N_{\Delta t,\Delta s}(t_{1})-N_{\Delta t,\Delta s}(t_{2})| =|\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),\]

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 Arzela-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, \tag{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_{j}^{0}n_{j}^{0}-\sum_{m=1}^{M}\sum_ {j\in\mathbb{N}}\Delta sn_{j}^{m}(\varphi_{j}^{m}-\varphi_{j}^{m-1}),\]

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_{j}^{m}(\varphi_{j}^{m}- \varphi_{j-1}^{m}),\]

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)

### 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, \tag{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} \tag{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

\[|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},\]

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

\[A_{1} = \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\] \[A_{2} = \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\] \[A_{3} = \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

\[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},\]

and for \(A_{3}\) we get

\[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}.\]

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\frac{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}. \tag{39}\]
From the mass conservation property, we get from (38) and (39) that the following estimate holds for \(H\)

\[\|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\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}.\]

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\)._

### 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}. \tag{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} \tag{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}, \tag{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) \tag{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}\gets 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\} \tag{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) \tag{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_

\[0\leq n_{j}^{m}\leq\|n^{0}\|_{\infty},\ \ 0\leq X^{m}\leq\|p\|_{\infty}\|n^{0} \|_{1}\|\alpha\|_{1},\ \ 0\leq N^{m}\leq\|p\|_{\infty}\|n^{0}\|_{1}\quad\text{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

\[0\leq F(X)\leq\|p\|_{\infty}\|n^{0}\|_{1}\|\alpha\|_{1}, \text{for all}\;X\geq 0,\] \[|F(X_{1})-F(X_{2})|\leq\tfrac{\Delta t}{2}\|n^{0}\|_{1}\alpha_{0 }\|\partial_{X}p\|_{\infty}|X_{1}-X_{2}|, \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), \tag{46}\]

_with \(T=m\Delta t\) and \(TV(n^{m})=\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

\[\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})).\]

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}. \tag{47}\]

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

\[|\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_{1}^{m}\] \[\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}.\]

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

\[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}),\]
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_{j}^{m+1}-n_{j}^{m}|\\ &\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_{j-1}^{m}|+\Delta t\|p\|_ {\infty}^{2}\|n^{0}\|_{1}\end{split} \tag{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\chi 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_{j}^{m+1}-n_{j}^{m} |\\ &\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_{j-1}^{m}|+\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_{j-1}^{m}|+ \Delta t\alpha_{0}\|p\|_{\infty}^{2}\|n^{0}\|_{1}\right.\\ &\qquad+\|p\|_{\infty}\|n^{0}\|_{1}\|\alpha\|_{\infty}+\|p\|_{ \infty}\|n^{0}\|_{1}TV(\alpha)\right).\end{split} \tag{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_{j-1}^{m}|+ 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

\[|\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_{j-1}^{m}|+B_{2} \Delta t, \tag{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, \tag{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)\;\text{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, \tag{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. \tag{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

\[\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,\]

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 Arzela-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)^ {+}}. \tag{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].

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

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})\).
\(|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). \tag{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.

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.

### 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). \tag{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_{1}^{0}\approx 0.0410\), \(N_{2}^{0}\approx 0.3650\) and \(N_{3}^{0}\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
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)\).

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.
observe three different types of numerical approach for \(N(t)\), which converge to two different equilibrium points \(N_{1}^{*}\) and \(N_{2}^{*}\).

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]\).

### 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), \tag{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). \tag{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.
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_{1}^{0}\approx 0.0281\), (c-d) for \(N_{2}^{0}\approx 0.4089\) and (e-f) for \(N_{3}^{0}\approx 0.7114\)
## 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).

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)\).
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.

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.
### Acknowledgements

This work has been supported by ANID project ECOS200018. MS and LMV was partially supported by ANID-Chile through Centro de Modelamiento Matematico (FB210005) and INRIA Associated team ANACONDA. MS was supported by Fondecyt-ANID project 1220869, and Jean d'Alembert fellowship program, Universite 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](https://orcid.org/0000-0001-8463-3830).

Nicolas Torres: [https://orcid.org/0000-0001-6059-9754](https://orcid.org/0000-0001-6059-9754).

Luis Miguel Villada: [https://orcid.org/0000-0002-4860-4431](https://orcid.org/0000-0002-4860-4431).
