# A multigroup approach to delayed prion production

Mostafa Adimy, Abdennasser Chekroun, Laurent Pujo-Menjouet, Mattia Sensi

► To cite this version:

Mostafa Adimy, Abdennasser Chekroun, Laurent Pujo-Menjouet, Mattia Sensi. A multigroup approach to delayed prion production. Discrete and Continuous Dynamical Systems - Series B, 2023, 10.3934/dcdsb.2023209. hal-04162164

HAL Id: hal-04162164

https://inria.hal.science/hal-04162164

Submitted on 8 Jan 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# A multigroup approach to delayed prion production

Mostafa Adimy, Abdennasser Chekroun, Laurent Pujo-Menjouet, Mattia Sensi

► To cite this version:

Mostafa Adimy, Abdennasser Chekroun, Laurent Pujo-Menjouet, Mattia Sensi. A multigroup approach to delayed prion production. Discrete and Continuous Dynamical Systems - Series B, 2024, 29 (7), pp.2972-2998. 10.3934/dcdsb.2023209. hal-04162164

HAL Id: hal-04162164

https://inria.hal.science/hal-04162164v1

Submitted on 8 Jan 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# A multigroup approach to delayed prion production

Mostafa Adimy<fnref n="1" />, Abdennasser Chekroun<fnref n="2" />, Laurent Pujo-Menjouet<fnref n="3" />, Mattia Sensi<fnref n="4" />,<fnref n="5" />,*

January 8, 2024

## Abstract

We generalize the model proposed in [Adimy, Babin, Pujo-Menjouet, *SIAM Journal on Applied Dynamical Systems* (2022)] for prion infection to a network of neurons. We do so by applying a so-called *multigroup approach* to the system of Delay Differential Equations (DDEs) proposed in the aforementioned paper. We derive the classical threshold quantity $R_0$, i.e. the basic reproduction number, exploiting the fact that the DDEs of our model qualitatively behave like Ordinary Differential Equations (ODEs) when evaluated at the Disease Free Equilibrium. We prove analytically that the disease naturally goes extinct when $R_0 < 1$, whereas it persists when $R_0 > 1$. We conclude with some selected numerical simulations of the system, to illustrate our analytical results.

# 1 Introduction

Prion is a protein involved in neurodegenerative diseases and more particularly the transmissible spongiform encephalopathies such as scrapie for sheep, bovine spongiform encephalopathy, also known as mad cow disease in cattle, and the Creutzfeldt-Jakob disease in humans [1, 2]. Produced by the cells, this protein in its normal form is called $PrP^C$ (for Prion Protein Cellular) and appears to be protective [2]. However, it becomes harmful and fatal when its shape changes. This misfolded pathological conformation also known as $PrP^{Sc}$ (for Prion Protein Scrapie) can be acquired either through transmission (this was the case for instance under the mad cow disease spread in the 1990s), or spontaneously, mostly above 75 years old for humans [3].

Even if extensively studied in the past decades, the action of this protein on the neurons leading to a fatal issue remains unclear. However, some recent discoveries may bring possible explanations and open new therapeutic strategies. This mechanism also known as Unfolded Protein Response (or UPR) [4, 5, 6, 7, 8] can be described as follows.

First, when produced by the cell, the $PrP^C$ proteins remain anchored to its membrane, unless misfolded $PrP^{Sc}$ in the extracellular matrix forces it to set it free and to join the pathological cohort. It is important to remind here that, by contact, a $PrP^{Sc}$ protein allows the normal form $PrP^C$ to change its conformation and to become misconformed. Once in this state, the proteins have the ability to polymerize, that is to tie together. They can easily reach very large sizes, stay in the neighbourhood of the cell or diffuse in the extracellular matrix to seed other neurons (see Fig. 1).

<footnote n="1">¹Inria, Univ Lyon, Université Claude Bernard Lyon 1, CNRS UMR 5208, Institut Camille Jordan, F-69603 Villeurbanne, France</footnote>

<footnote n="2">²Laboratoire d'Analyse Nonlinéaire et Mathématiques Appliquées, University of Tlemcen, Tlemcen 13000, Algeria</footnote>

<footnote n="3">³Univ Lyon, Inria, Université Claude Bernard Lyon 1, CNRS UMR 5208, Institut Camille Jordan, F-69603 Villeurbanne, France</footnote>

<footnote n="4">⁴MathNeuro Team, Inria at Université Côte d'Azur, 2004 Rte des Lucioles, 06410 Biot, France,</footnote>

<footnote n="5">⁵Department of Mathematical Sciences “G. L. Lagrange”, Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino Italy</footnote>

<footnote>*Corresponding author: mattia.sensi@polito.it</footnote>
# A multigroup approach to delayed prion production

Mostafa Adimy¹, Abdennasser Chekroun², Laurent Pujo-Menjouet³, Mattia Sensi⁴,⁵,*

¹Inria, Univ Lyon, Université Claude Bernard Lyon 1, CNRS UMR 5208, Institut Camille Jordan, F-69603 Villeurbanne, France

²Laboratoire d'Analyse Nonlinéaire et Mathématiques Appliquées, University of Tlemcen, Tlemcen 13000, Algeria

³Univ Lyon, Inria, Université Claude Bernard Lyon 1, CNRS UMR 5208, Institut Camille Jordan, F-69603 Villeurbanne, France

⁴MathNeuro Team, Inria at Université Côte d'Azur, 2004 Rte des Lucioles, 06410 Biot, France,

⁵Department of Mathematical Sciences “G. L. Lagrange”, Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino Italy

*Corresponding author: mattia.sensi@polito.it

January 8, 2024

## Abstract

We generalize the model proposed in [Adimy, Babin, Pujo-Menjouet, *SIAM Journal on Applied Dynamical Systems* (2022)] for prion infection to a network of neurons. We do so by applying a so-called *multigroup approach* to the system of Delay Differential Equations (DDEs) proposed in the aforementioned paper. We derive the classical threshold quantity $R_0$, i.e. the basic reproduction number, exploiting the fact that the DDEs of our model qualitatively behave like Ordinary Differential Equations (ODEs) when evaluated at the Disease Free Equilibrium. We prove analytically that the disease naturally goes extinct when $R_0 < 1$, whereas it persists when $R_0 > 1$. We conclude with some selected numerical simulations of the system, to illustrate our analytical results.

## 1 Introduction

Prion is a protein involved in neurodegenerative diseases and more particularly the transmissible spongiform encephalopathies such as scrapie for sheep, bovine spongiform encephalopathy, also known as mad cow disease in cattle, and the Creutzfeldt-Jakob disease in humans [1, 2]. Produced by the cells, this protein in its normal form is called $PrP^C$ (for Prion Protein Cellular) and appears to be protective [2]. However, it becomes harmful and fatal when its shape changes. This misfolded pathological conformation also known as $PrP^{Sc}$ (for Prion Protein Scrapie) can be acquired either through transmission (this was the case for instance under the mad cow disease spread in the 1990s), or spontaneously, mostly above 75 years old for humans [3].

Even if extensively studied in the past decades, the action of this protein on the neurons leading to a fatal issue remains unclear. However, some recent discoveries may bring possible explanations and open new therapeutic strategies. This mechanism also known as Unfolded Protein Response (or UPR) [4, 5, 6, 7, 8] can be described as follows.

First, when produced by the cell, the $PrP^C$ proteins remain anchored to its membrane, unless misfolded $PrP^{Sc}$ in the extracellular matrix forces it to set it free and to join the pathological cohort. It is important to remind here that, by contact, a $PrP^{Sc}$ protein allows the normal form $PrP^C$ to change its conformation and to become misconformed. Once in this state, the proteins have the ability to polymerize, that is to tie together. They can easily reach very large sizes, stay in the neighbourhood of the cell or diffuse in the extracellular matrix to seed other neurons (see Fig. 1).
<figure><img src="image_2.png" /><figcaption>Figure 1: schematic view of the PrP<sup>C</sup> protein production (in blue) by two neurons (green). The PrP<sup>C</sup> protein can aggregate and form pathological PrP<sup>Sc</sup> (pink and orange). The PrP<sup>Sc</sup> proteins diffuse and a certain amount can reach the neighbourhood of another neuron (the orange ones can reach the neighbourhood of neuron 1, while the pink ones can reach neuron 2. We refer to Section 2 for a complete description of this case and of the parameters and variables involved.</figcaption></figure>

If for some reason, such as an over-expression of PrP<sup>C</sup> or a slow diffusion, they accumulate in the neuron proximity, this latter feels it and under this induced stress shuts down almost all its activities except the vital ones.

This global shutdown, created by a high concentration of PrP<sup>Sc</sup> in the neuron surrounding, causes the neuron to stop producing PrP<sup>C</sup>, not vital for the cell (see Fig. 2).

This break ends only if these proteins move away by diffusion or degradation. When the zone is clear, the cell starts again its protein production and the process continues until the next stress period.

Still under investigation, the detailed UPR mechanism remains to be fully understood, even if several papers may be referred to the reader [5, 6, 7, 8]. Besides, the link between UPR, PrP<sup>Sc</sup> has been put in evidence [9, 10, 11, 12, 13, 14].

Because of its complexity, the UPR *modus operandi* has already been the object of mathematical models, from a gene regulatory point of view [15, 16, 17, 18, 19] or through regulation of UPR intra- or extra-cellular pathways [17, 18]. Our goal here is to generalise the pioneering mathematical model [20] dealing explicitly with prion, and investigating the parameters causing the oscillating neural activity. In [20], the authors investigated the case of one neuron only, and for two neurons they gave analytical results specifically when both cells would exhibit the exact same behaviour. In this paper, we briefly remind the model with two neurons and give new theoretical results to complete the ones of [20], then we extend the construction to any neuron number $n \in \mathbb{N}$, $n \ge 2$.

We exploit the formulation of the Delay Differential Equation (DDE) system set up in detail in [20], where the delay is only present in the infectious/infected variables, in order to apply a classical tool of Ordinary Differential Equations (ODE) epidemic models, namely the Next Generation Matrix. This technique was first introduced in [21], then generalized in [22] (see also [23]). Through an appropriate decomposition of the Jacobian matrix evaluated at the Disease Free Equilibrium, we are able to provide a formulation for the Basic Reproduction Number $R_0$ of the $n$-dimensional system, under biologically acceptable conditions.

Then, we apply the definition of the threshold quantity $R_0$ to prove either global stability of the Disease Free
<figure><img src="image_2.png" /><figcaption>Figure 1: schematic view of the $PrP^C$ protein production (in blue) by two neurons (green). The $PrP^C$ protein can aggregate and form pathological $PrP^{Sc}$ (pink and orange). The $PrP^{Sc}$ proteins diffuse and a certain amount can reach the neighbourhood of another neuron (the orange ones can reach the neighbourhood of neuron 1, while the pink ones can reach neuron 2. We refer to Section 2 for a complete description of this case and of the parameters and variables involved.</figcaption></figure>

If for some reason, such as an over-expression of $PrP^C$ or a slow diffusion, they accumulate in the neuron proximity, this latter feels it and under this induced stress shuts down almost all its activities except the vital ones.

This global shutdown, created by a high concentration of $PrP^{Sc}$ in the neuron surrounding, causes the neuron to stop producing $PrP^C$, not vital for the cell (see Fig. 2).

This break ends only if these proteins move away by diffusion or degradation. When the zone is clear, the cell starts again its protein production and the process continues until the next stress period.

Still under investigation, the detailed UPR mechanism remains to be fully understood, even if several papers may be referred to the reader [5, 6, 7, 8]. Besides, the link between UPR, $PrP^{Sc}$ has been put in evidence [9, 10, 11, 12, 13, 14].

Because of its complexity, the UPR *modus operandi* has already been the object of mathematical models, from a gene regulatory point of view [15, 16, 17, 18, 19] or through regulation of UPR intra- or extra-cellular pathways [17, 18]. Our goal here is to generalise the pioneering mathematical model [20] dealing explicitly with prion, and investigating the parameters causing the oscillating neural activity. In [20], the authors investigated the case of one neuron only, and for two neurons they gave analytical results specifically when both cells would exhibit the exact same behaviour. In this paper, we briefly remind the model with two neurons and give new theoretical results to complete the ones of [20], then we extend the construction to any neuron number $n \in \mathbb{N}$, $n \ge 2$.

We exploit the formulation of the Delay Differential Equation (DDE) system set up in detail in [20], where the delay is only present in the infectious/infected variables, in order to apply a classical tool of Ordinary Differential Equations (ODE) epidemic models, namely the Next Generation Matrix. This technique was first introduced in [21], then generalized in [22] (see also [23]). Through an appropriate decomposition of the Jacobian matrix evaluated at the Disease Free Equilibrium, we are able to provide a formulation for the Basic Reproduction Number $R_0$ of the $n$-dimensional system, under biologically acceptable conditions.

Then, we apply the definition of the threshold quantity $R_0$ to prove either global stability of the Disease Free
<figure><img src="image_3.png" /><figcaption>Figure 2: schematic representation of a neuron (neuron 1) under Unfolded Protein Response (UPR). Stressed by the overcrowded amount of $PrP^{Sc}$ in its neighbourhood, neuron 1 shuts down its activities (except the vital ones). No $PrP^C$ protein is then produced, and the population of $PrP^{Sc}$ pathological proteins diffuse out the neuron surroundings.</figcaption></figure>

Equilibrium (when $R_0 < 1$) or permanence of the system (when $R_0 > 1$). Moreover, under slightly stricter conditions, we are able to prove the existence of at least one Endemic Equilibrium.

The paper is structured as follows. In Section 2, we recall the 2 neurons model introduced in [20]. In Section 3, we generalize this construction to a network of $n$ neurons, with $n \ge 2$; moreover, we show two other results: first, how the fully connected and fully homogeneous case can be qualitatively reduced to a single neuron model (but with a different $R_0$) proposed in [20], and second how the case of one-way direction connection of several neurons behaves like a single one. In Section 4, we prove global stability of the Disease Free Equilibrium when $R_0 < 1$. In Section 5, we show a condition for the existence (but not uniqueness) of the Endemic Equilibrium. In Section 6, we show the persistence of the system when $R_0 > 1$. In Section 7, we provide extensive numerical simulations of the model proposed in Section 3. Lastly, in Section 8 we conclude.

## 2 System with 2 neurons

We begin by recalling the system of 2 neurons from [20]. It describes the dynamics of the $PrP^C$ protein associated with neuron 1 and neuron 2, respectively $x_1$ and $x_2$, as well as the $PrP^{Sc}$ concentrations in the environment of neuron 1 and neuron 2, $y_1$ and $y_2$. Due to their biological interpretation, we only consider $x_i, y_i \ge 0$. This model is represented,
<figure><img src="image_3.png" /><figcaption>Figure 2: schematic representation of a neuron (neuron 1) under Unfolded Protein Response (UPR). Stressed by the overcrowded amount of $PrP^{Sc}$ in its neighbourhood, neuron 1 shuts down its activities (except the vital ones). No $PrP^C$ protein is then produced, and the population of $PrP^{Sc}$ pathological proteins diffuse out the neuron surroundings.</figcaption></figure>

Equilibrium (when $R_0 < 1$) or permanence of the system (when $R_0 > 1$). Moreover, under slightly stricter conditions, we are able to prove the existence of at least one Endemic Equilibrium.

The paper is structured as follows. In Section 2, we recall the 2 neurons model introduced in [20]. In Section 3, we generalize this construction to a network of $n$ neurons, with $n \ge 2$; moreover, we show two other results: first, how the fully connected and fully homogeneous case can be qualitatively reduced to a single neuron model (but with a different $R_0$) proposed in [20], and second how the case of one-way direction connection of several neurons behaves like a single one. In Section 4, we prove global stability of the Disease Free Equilibrium when $R_0 < 1$. In Section 5, we show a condition for the existence (but not uniqueness) of the Endemic Equilibrium. In Section 6, we show the persistence of the system when $R_0 > 1$. In Section 7, we provide extensive numerical simulations of the model proposed in Section 3. Lastly, in Section 8 we conclude.

## 2 System with 2 neurons

We begin by recalling the system of 2 neurons from [20]. It describes the dynamics of the $PrP^C$ protein associated with neuron 1 and neuron 2, respectively $x_1$ and $x_2$, as well as the $PrP^{Sc}$ concentrations in the environment of neuron 1 and neuron 2, $y_1$ and $y_2$. Due to their biological interpretation, we only consider $x_i, y_i \ge 0$. This model is represented,
for $t > 0$, by the following system

$$
\begin{align}
\frac{dx_1}{dt} &= K_1\beta(y_1(t-T_1)) - \mu_1 x_1(t) - dx_1(t)(y_1(t) + \kappa\alpha_2 y_2(t)), \nonumber \\
\frac{dx_2}{dt} &= K_2\beta(y_2(t-T_2)) - \mu_2 x_2(t) - dx_2(t)(y_2(t) + \kappa\alpha_1 y_1(t)), \nonumber \\
\frac{dy_1}{dt} &= dx_1(t)(y_1(t) + \kappa\alpha_2 y_2(t)) - \alpha_1 y_1(t), \nonumber \\
\frac{dy_2}{dt} &= dx_2(t)(y_2(t) + \kappa\alpha_1 y_1(t)) - \alpha_2 y_2(t), \tag{1}
\end{align}
$$

where $K_i > 0$ ($i = 1$ or 2) represents the PrP$^C$ production rate of the neuron $i$ and $d > 0$ characterizes the force of the interaction between PrP$^C$ and PrP$^{Sc}$. The terms $dx_1(t)$ ($y_1(t) + \kappa\alpha_2 y_2(t)$) and $dx_2(t)$ ($y_2(t) + \kappa\alpha_1 y_1(t)$) stand for the new PrP$^{Sc}$ produced. The parameter $\mu_i$ represents the degradation rate of PrP$^C$ produced by the neuron $i$ and $\alpha_i$ is the rate at which PrP$^{Sc}$ proteins are lost through degradation or diffusion. The factor $\kappa$ indicates the interaction between proteins from different neurons. The parameter $T_i$ is the time required for a neuron $i$ to process the PrP$^C$ protein synthesis. Due to the UPR effect, increasing the amount of PrP$^{Sc}$ around a neuron decreases its activity and consequently the PrP$^C$ production. The contribution of PrP$^{Sc}$ concentration to PrP$^C$ production is therefore given through the decreasing Hill function (negative feedback [20])

$$ \beta(y) = \frac{1}{1 + (y/y_c)^p}, \quad (2) $$

where $p > 0$ is the sensitivity of PrP$^{Sc}$ production to PrP$^{Sc}$ overload. The parameter $y_c > 0$ is the PrP$^{Sc}$ threshold beyond which the neuron stops PrP$^C$ production.

Compared to the notation in [20], to avoid confusion we write $\beta$ instead of $\beta_n$ since $n$ will represent the number of neurons in the system from Section 3 onward.

The Disease Free Equilibrium corresponding to the system (1) is

$$ (x_1, x_2, y_1, y_2) = \left( \frac{K_1}{\mu_1}, \frac{K_2}{\mu_2}, 0, 0 \right). \quad (3) $$

The linearized version of system (1) around the Disease Free Equilibrium (3) is given by

$$
\begin{align}
\frac{dx_1}{dt} &= K_1\beta'(0)y_1(t-T_1) - \mu_1 x_1(t) - d\frac{K_1}{\mu_1}(y_1(t) + \kappa\alpha_2 y_2(t)), \nonumber \\
\frac{dx_2}{dt} &= K_2\beta'(0)y_2(t-T_2) - \mu_2 x_2(t) - d\frac{K_2}{\mu_2}(y_2(t) + \kappa\alpha_1 y_1(t)), \nonumber \\
\frac{dy_1}{dt} &= d\frac{K_1}{\mu_1}(y_1(t) + \kappa\alpha_2 y_2(t)) - \alpha_1 y_1(t), \nonumber \\
\frac{dy_2}{dt} &= d\frac{K_2}{\mu_2}(y_2(t) + \kappa\alpha_1 y_1(t)) - \alpha_2 y_2(t). \tag{4}
\end{align}
$$

Note that $\beta'(0) = 0$. Then, the linearized system (4) becomes an ordinary differential system and its Jacobian matrix is

$$ J_{\text{DFE}} = \begin{pmatrix} -\mu_1 & 0 & -d\frac{K_1}{\mu_1} & -d\kappa\alpha_2\frac{K_1}{\mu_1} \\ 0 & -\mu_2 & -d\kappa\alpha_1\frac{K_2}{\mu_2} & -d\frac{K_2}{\mu_2} \\ 0 & 0 & d\frac{K_1}{\mu_1} - \alpha_1 & d\kappa\alpha_2\frac{K_1}{\mu_1} \\ 0 & 0 & d\kappa\alpha_1\frac{K_2}{\mu_2} & d\frac{K_2}{\mu_2} - \alpha_2 \end{pmatrix}. $$

Now, we use the Next Generation Matrix method, firstly introduced in [21], then generalized in [22] (see also [23]) to obtain the basic reproduction number $R_0$ of the system (1). In order to do so, we need to write $J_{\text{DFE}}$ as $J_{\text{DFE}} = M - V$, with $M$ having non-negative entries and $V$ invertible. One possible choice is the following

$$ M = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & d\frac{K_1}{\mu_1} & d\kappa\alpha_2\frac{K_1}{\mu_1} \\ 0 & 0 & d\kappa\alpha_1\frac{K_2}{\mu_2} & d\frac{K_2}{\mu_2} \end{pmatrix} \quad \text{and} \quad V = \begin{pmatrix} \mu_1 & 0 & d\frac{K_1}{\mu_1} & d\kappa\alpha_2\frac{K_1}{\mu_1} \\ 0 & \mu_2 & d\kappa\alpha_1\frac{K_2}{\mu_2} & d\frac{K_2}{\mu_2} \\ 0 & 0 & \alpha_1 & 0 \\ 0 & 0 & 0 & \alpha_2 \end{pmatrix}. $$
for $t > 0$, by the following system

$$ 
\begin{aligned}
\frac{dx_1}{dt} &= K_1\beta(y_1(t-T_1)) - \mu_1 x_1(t) - dx_1(t)(y_1(t) + \kappa\alpha_2 y_2(t)), \\
\frac{dx_2}{dt} &= K_2\beta(y_2(t-T_2)) - \mu_2 x_2(t) - dx_2(t)(y_2(t) + \kappa\alpha_1 y_1(t)), \\
\frac{dy_1}{dt} &= dx_1(t)(y_1(t) + \kappa\alpha_2 y_2(t)) - \alpha_1 y_1(t), \\
\frac{dy_2}{dt} &= dx_2(t)(y_2(t) + \kappa\alpha_1 y_1(t)) - \alpha_2 y_2(t),
\end{aligned}
\quad (1) $$

where $K_i > 0$ ($i = 1$ or 2) represents the PrP<sup>C</sup> production rate of the neuron i and $d > 0$ characterizes the force of the interaction between PrP<sup>C</sup> and PrP<sup>Sc</sup>. The terms $dx_1(t)$ ($y_1(t) + \kappa\alpha_2 y_2(t)$) and $dx_2(t)$ ($y_2(t) + \kappa\alpha_1 y_1(t)$) stand for the new PrP<sup>Sc</sup> produced. The parameter $\mu_i$ represents the degradation rate of PrP<sup>C</sup> produced by the neuron i and $\alpha_i$ is the rate at which PrP<sup>Sc</sup> proteins are lost through degradation or diffusion. The factor $\kappa$ indicates the interaction between proteins from different neurons. The parameter $T_i$ is the time required for a neuron i to process the PrP<sup>C</sup> protein synthesis. Due to the UPR effect, increasing the amount of PrP<sup>Sc</sup> around a neuron decreases its activity and consequently the PrP<sup>C</sup> production. The contribution of PrP<sup>Sc</sup> concentration to PrP<sup>C</sup> production is therefore given through the decreasing Hill function (negative feedback [20])

$$ \beta(y) = \frac{1}{1 + (y/y_c)^p}, \quad (2) $$

where $p > 0$ is the sensitivity of PrP<sup>Sc</sup> production to PrP<sup>Sc</sup> overload. The parameter $y_c > 0$ is the PrP<sup>Sc</sup> threshold beyond which the neuron stops PrP<sup>C</sup> production.

Compared to the notation in [20], to avoid confusion we write $\beta$ instead of $\beta_n$ since n will represent the number of neurons in the system from Section 3 onward.

The Disease Free Equilibrium corresponding to the system (1) is

$$ (x_1, x_2, y_1, y_2) = \left( \frac{K_1}{\mu_1}, \frac{K_2}{\mu_2}, 0, 0 \right). \quad (3) $$

The linearized version of system (1) around the Disease Free Equilibrium (3) is given by

$$ 
\begin{aligned}
\frac{dx_1}{dt} &= K_1\beta'(0)y_1(t-T_1) - \mu_1 x_1(t) - d\frac{K_1}{\mu_1}(y_1(t) + \kappa\alpha_2 y_2(t)), \\
\frac{dx_2}{dt} &= K_2\beta'(0)y_2(t-T_2) - \mu_2 x_2(t) - d\frac{K_2}{\mu_2}(y_2(t) + \kappa\alpha_1 y_1(t)), \\
\frac{dy_1}{dt} &= d\frac{K_1}{\mu_1}(y_1(t) + \kappa\alpha_2 y_2(t)) - \alpha_1 y_1(t), \\
\frac{dy_2}{dt} &= d\frac{K_2}{\mu_2}(y_2(t) + \kappa\alpha_1 y_1(t)) - \alpha_2 y_2(t).
\end{aligned}
\quad (4) $$

Note that $\beta'(0) = 0$. Then, the linearized system (4) becomes an ordinary differential system and its Jacobian matrix is

$$ J_{\text{DFE}} = \begin{pmatrix} -\mu_1 & 0 & -d\frac{K_1}{\mu_1} & -d\kappa\alpha_2\frac{K_1}{\mu_1} \\ 0 & -\mu_2 & -d\kappa\alpha_1\frac{K_2}{\mu_2} & -d\frac{K_2}{\mu_2} \\ 0 & 0 & d\frac{K_1}{\mu_1} - \alpha_1 & d\kappa\alpha_2\frac{K_1}{\mu_1} \\ 0 & 0 & d\kappa\alpha_1\frac{K_2}{\mu_2} & d\frac{K_2}{\mu_2} - \alpha_2 \end{pmatrix}. $$

Now, we use the Next Generation Matrix method, firstly introduced in [21], then generalized in [22] (see also [23]) to obtain the basic reproduction number $R_0$ of the system (1). In order to do so, we need to write $J_{\text{DFE}}$ as $J_{\text{DFE}} = M - V$, with M having non-negative entries and V invertible. One possible choice is the following

$$ M = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & d\frac{K_1}{\mu_1} & d\kappa\alpha_2\frac{K_1}{\mu_1} \\ 0 & 0 & d\kappa\alpha_1\frac{K_2}{\mu_2} & d\frac{K_2}{\mu_2} \end{pmatrix} \quad \text{and} \quad V = \begin{pmatrix} \mu_1 & 0 & d\frac{K_1}{\mu_1} & d\kappa\alpha_2\frac{K_1}{\mu_1} \\ 0 & \mu_2 & d\kappa\alpha_1\frac{K_2}{\mu_2} & d\frac{K_2}{\mu_2} \\ 0 & 0 & \alpha_1 & 0 \\ 0 & 0 & 0 & \alpha_2 \end{pmatrix}. $$
The basic reproduction number $\mathcal{R}_0$ is exactly $\rho(MV^{-1})$. Remark that, in order to compute this spectral radius, we implicitly assumed that $\alpha_1, \alpha_2 \neq 0$. This means that each neuron receives a strictly positive amount of infection from the other. Recall that in [20] the basic reproduction number of neuron $i$ was computed as

$$R_{0i} = d \frac{K_i}{\mu_i \alpha_i}.$$

Since the first two rows of $M$ are 0, it suffices to observe the matrix

$$F = \begin{pmatrix} d \frac{K_1}{\mu_1 \alpha_1} & d \frac{K_1}{\mu_1 \alpha_1} \\ d \frac{K_2}{\mu_2 \alpha_2} & d \frac{K_2}{\mu_2 \alpha_2} \end{pmatrix} = \begin{pmatrix} R_{01} & \kappa \alpha_1 R_{01} \\ \kappa \alpha_2 R_{02} & R_{02} \end{pmatrix},$$

which has eigenvalues

$$\lambda_{\pm} = \frac{R_{01} + R_{02} \pm \sqrt{(R_{01} - R_{02})^2 + 4\kappa^2 \alpha_1 \alpha_2 R_{01} R_{02}}}{2}, \quad (5)$$

with $\lambda_+ = \rho(F)$ being the new $\mathcal{R}_0$ of the 2-neurons system. We remark that the connectivity between the two neurons $\kappa$ plays a fundamental role in the dynamics: even if both $R_{0i} < 1$, with $\kappa$ large enough the disease could remain endemic. However, we recall that due to its biological interpretation, the relevant region we should consider is $\kappa \in [0, 1]$.

In the next section, we generalize this construction to a network of $n \in \mathbb{N}_{\ge 2}$ neurons.

## 3 System with $n$ neurons

The construction from the previous section can be generalized to a $n$ neurons case by similarly constructing the matrices $M$ and $V$. In this case, we would generally obtain $\mathcal{R}_0$ implicitly, as the spectral radius of a $2n \times 2n$ matrix. However, such a matrix can be reduced to $n \times n$ as in the previous section since $M$ will only have the lower-right quarter of non-zero entries.

We consider the following system of Delay Differential Equations (DDEs)

$$\begin{aligned} \frac{dx_i}{dt} &= K_i \beta(y_i(t-T_i)) - \mu_i x_i(t) - dx_i(t) \left( y_i(t) + \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i} y_j(t) \right), \\ \frac{dy_i}{dt} &= dx_i(t) \left( y_i(t) + \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i} y_j(t) \right) - \left( \sum_{j \neq i} \alpha_{i \to j} \right) y_i(t). \end{aligned} \quad (6)$$

Due to their biological interpretation, we only consider $x_i, y_i \ge 0$. The parameter $\alpha_{i \to j}$ represents the fraction of prions produced by neuron $i$ and moving towards neuron $j$. It also includes prion degradation. In other words, $\alpha_{i \to j}$ describes the diffusive property (including degradation) of $\text{PrP}^{Sc}$ to the neuron $j \neq i$. This includes both prions which die while moving away and prions which actually reach neuron $j$. The interactions between $\text{PrP}^C$ from neuron $i$ with $\text{PrP}^{Sc}$ of another neuron $j \neq i$ is given by the factor $\kappa_{ji}$ (it characterizes the difference between prion species); hence, $\sum_{j \neq i} \kappa_{ji} \le 1$ for all $j$, since this sum represents the fraction of prions “orbiting” neuron $i$ (neuron has a number of prions it can spread to others) which does not die and manages to spread to other neurons.

For ease of notation, let

$$\alpha_i := \sum_{j \neq i} \alpha_{i \to j}$$

denote the total rate of migration of prions from neuron $i$, which can result in either the death of the prion or contact with any other neuron $j \neq i$.

Let $C := C([-T, 0], \mathbb{R})$, $T := \max_{i=1,...,n} T_i$, be the space of continuous functions on $[-T, 0]$ and $C^+ := C([-T, 0], \mathbb{R}^+)$ be the space of nonnegative continuous functions on $[-T, 0]$. We assume throughout this paper that the initial conditions for the system (6), i.e. $(x_{i0}, \varphi_i) \in \mathbb{R}^+ \times C^+$, for $i = 1, ..., n$. The existence and uniqueness of nonnegative solutions of (6) can be obtained by using the theory of functional differential equations.
The basic reproduction number $\mathcal{R}_0$ is exactly $\rho(MV^{-1})$. Remark that, in order to compute this spectral radius, we implicitly assumed that $\alpha_1, \alpha_2 \neq 0$. This means that each neuron receives a strictly positive amount of infection from the other. Recall that in [20] the basic reproduction number of neuron $i$ was computed as

$$R_{0i} = d \frac{K_i}{\mu_i \alpha_i}.$$

Since the first two rows of $M$ are 0, it suffices to observe the matrix

$$F = \begin{pmatrix} d \frac{K_1}{\mu_1 \alpha_1} & d \frac{K_1}{\mu_1 \alpha_1} \\ d \frac{K_2}{\mu_2 \alpha_2} & d \frac{K_2}{\mu_2 \alpha_2} \end{pmatrix} = \begin{pmatrix} R_{01} & \kappa \alpha_1 R_{01} \\ \kappa \alpha_2 R_{02} & R_{02} \end{pmatrix},$$

which has eigenvalues

$$\lambda_{\pm} = \frac{R_{01} + R_{02} \pm \sqrt{(R_{01} - R_{02})^2 + 4\kappa^2 \alpha_1 \alpha_2 R_{01} R_{02}}}{2}, \quad (5)$$

with $\lambda_+ = \rho(F)$ being the new $\mathcal{R}_0$ of the 2-neurons system. We remark that the connectivity between the two neurons $\kappa$ plays a fundamental role in the dynamics: even if both $R_{0i} < 1$, with $\kappa$ large enough the disease could remain endemic. However, we recall that due to its biological interpretation, the relevant region we should consider is $\kappa \in [0, 1]$.

In the next section, we generalize this construction to a network of $n \in \mathbb{N}_{\ge 2}$ neurons.

## 3 System with $n$ neurons

The construction from the previous section can be generalized to a $n$ neurons case by similarly constructing the matrices $M$ and $V$. In this case, we would generally obtain $\mathcal{R}_0$ implicitly, as the spectral radius of a $2n \times 2n$ matrix. However, such a matrix can be reduced to $n \times n$ as in the previous section since $M$ will only have the lower-right quarter of non-zero entries.

We consider the following system of Delay Differential Equations (DDEs)

$$\frac{dx_i}{dt} = K_i \beta(y_i(t - T_i)) - \mu_i x_i(t) - dx_i(t) \left( y_i(t) + \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i} y_j(t) \right), \quad (6)$$

$$\frac{dy_i}{dt} = dx_i(t) \left( y_i(t) + \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i} y_j(t) \right) - \left( \sum_{j \neq i} \alpha_{i \to j} \right) y_i(t).$$

Due to their biological interpretation, we only consider $x_i, y_i \ge 0$. The parameter $\alpha_{i \to j}$ represents the fraction of prions produced by neuron $i$ and moving towards neuron $j$. It also includes prion degradation. In other words, $\alpha_{i \to j}$ describes the diffusive property (including degradation) of $\text{PrP}^{Sc}$ to the neuron $j \neq i$. This includes both prions which die while moving away and prions which actually reach neuron $j$. The interactions between $\text{PrP}^C$ from neuron $i$ with $\text{PrP}^{Sc}$ of another neuron $j \neq i$ is given by the factor $\kappa_{ji}$ (it characterizes the difference between prion species); hence, $\sum_{j \neq i} \kappa_{ji} \le 1$ for all $j$, since this sum represents the fraction of prions “orbiting” neuron $i$ (neuron has a number of prions it can spread to others) which does not die and manages to spread to other neurons.

For ease of notation, let

$$\alpha_i := \sum_{j \neq i} \alpha_{i \to j}$$

denote the total rate of migration of prions from neuron $i$, which can result in either the death of the prion or contact with any other neuron $j \neq i$.

Let $C := C([-T, 0], \mathbb{R})$, $T := \max_{i=1,...,n} T_i$, be the space of continuous functions on $[-T, 0]$ and $C^+ := C([-T, 0], \mathbb{R}^+)$ be the space of nonnegative continuous functions on $[-T, 0]$. We assume throughout this paper that the initial conditions for the system (6), i.e. $(x_{i0}, \varphi_i) \in \mathbb{R}^+ \times C^+$, for $i = 1, ..., n$. The existence and uniqueness of nonnegative solutions of (6) can be obtained by using the theory of functional differential equations.
Since the delay is discrete, the continuity of $\beta$ is sufficient to ensure the existence and uniqueness of the solution (see, [24, 25]). We call the history function each function $u_t \in C$, for $t \ge 0$ and $u \in C([-T, +\infty), \mathbb{R})$ satisfying $u_t(\theta) = u(t + \theta)$ for $\theta \in [-T, 0]$. Now, we show the nonnegativity and boundedness of solutions of the system (6).

**Proposition 1.** All solutions of the system (6) with nonnegative initial conditions remain nonnegative and bounded.

*Proof.* We prove nonnegativity by applying the Theorem 3.4 of [26]. In fact, for $i = 1, ..., n$, if $x_i(t) = 0$ then

$$ \frac{dx_i}{dt} = K_i\beta(y_i(t - T_i)) \ge 0, \quad \text{for } y_i \in C^+, $$

and if $y_i(t) = 0$ then

$$ \frac{dy_i}{dt} = dx_i(t) \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i} y_j(t) \ge 0, \quad \text{for } x_i \in \mathbb{R}^+, y_j \in C^+. $$

Then, by Theorem 3.4 of [26], we get $x_i(t) \ge 0$ and $y_i(t) \ge 0$ for $t \ge 0$.

Now, by adding both equation of $x_i$ and $y_i$, we get, for $t \ge 0$,

$$ \frac{dx_i}{dt} + \frac{dy_i}{dt} = K_i\beta(y_i(t - T_i)) - \mu_i x_i(t) - \left( \sum_{j \neq i} \alpha_{i \to j} \right) y_i(t). $$

This implies that, for $t \ge 0$,

$$ \frac{d(x_i + y_i)}{dt} \le K_i\beta(0) - \min\{\mu_i, \alpha_i\}(x_i + y_i). $$

This means that

$$ \limsup_{t \to +\infty} x_i(t) + y_i(t) \le \frac{K_i\beta(0)}{\min\{\mu_i, \alpha_i\}}. $$

Therefore, the solution should be necessarily bounded.

We now introduce a formula for the Basic Reproduction Number (BRN) $\mathcal{R}_0$ of the system (6), given as the spectral radius of an $n \times n$ matrix. We do so by applying the Next Generation Matrix method [21, 22, 23]. We remark that this method was developed specifically for systems of ODEs. However, when evaluated in its Disease Free Equilibrium, namely

$$ (x_1, \dots, x_n, y_1, \dots, y_n) = \left( \frac{K_1}{\mu_1}, \dots, \frac{K_n}{\mu_n}, 0, \dots, 0 \right), \qquad (7) $$

the system (6) does not exhibit any form of delay, and qualitatively reduces to a system of ODEs. We focus on the Jacobian on the system evaluated in this equilibrium, obtaining a reliable threshold quantity.

**Proposition 2.** Recall from [20] that the “Basic Reproduction Number of neuron i” is $R_{0i} = d \frac{K_i}{\mu_i \alpha_i}$.

The Basic Reproduction Number $\mathcal{R}_0$ of the system (6) is given by the spectral radius of the matrix $F \in \mathbb{R}^{n \times n}$ defined as

$$ (F)_{ij} = \begin{cases} R_{0i} & \text{if } j = i, \\ \kappa_{ji} \alpha_{j \to i} R_{0i} & \text{if } j \neq i. \end{cases} \qquad (8) $$

*Proof.* For ease of notation, we use diag(·) to indicate diag(·)$_{1 \le i \le n}$, since all the diagonal matrices we consider are of dimension $n \times n$.

We compute the Jacobian $J$ of (6), dropping the explicit dependence on $(t)$ everywhere for ease of notation. By “splitting” the system into $x$ and $y$, we can write

$$ J = \begin{pmatrix} J_{11} & J_{12} \\ J_{21} & J_{22} \end{pmatrix}, $$

where

$$ J_{11} = \operatorname{diag} \left( -\mu_i - d \left( y_i + \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i} y_j \right) \right), \quad (J_{12})_{ij} = \begin{cases} K_i \beta'(y_i) - dx_i & \text{if } j = i, \\ -dx_i \kappa_{ji} \alpha_{j \to i} & \text{if } j \neq i, \end{cases} $$
Since the delay is discrete, the continuity of $\beta$ is sufficient to ensure the existence and uniqueness of the solution (see, [24, 25]). We call the history function each function $u_t \in C$, for $t \ge 0$ and $u \in C([-T, +\infty), \mathbb{R})$ satisfying $u_t(\theta) = u(t + \theta)$ for $\theta \in [-T, 0]$. Now, we show the nonnegativity and boundedness of solutions of the system (6).

**Proposition 1.** All solutions of the system (6) with nonnegative initial conditions remain nonnegative and bounded.

*Proof.* We prove nonnegativity by applying the Theorem 3.4 of [26]. In fact, for $i = 1, ..., n$, if $x_i(t) = 0$ then

$$ \frac{dx_i}{dt} = K_i\beta(y_i(t - T_i)) \ge 0, \quad \text{for } y_i \in C^+, $$

and if $y_i(t) = 0$ then

$$ \frac{dy_i}{dt} = dx_i(t) \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i} y_j(t) \ge 0, \quad \text{for } x_i \in \mathbb{R}^+, y_j \in C^+. $$

Then, by Theorem 3.4 of [26], we get $x_i(t) \ge 0$ and $y_i(t) \ge 0$ for $t \ge 0$.

Now, by adding both equation of $x_i$ and $y_i$, we get, for $t \ge 0$,

$$ \frac{dx_i}{dt} + \frac{dy_i}{dt} = K_i\beta(y_i(t - T_i)) - \mu_i x_i(t) - \left( \sum_{j \neq i} \alpha_{i \to j} \right) y_i(t). $$

This implies that, for $t \ge 0$,

$$ \frac{d(x_i + y_i)}{dt} \le K_i\beta(0) - \min\{\mu_i, \alpha_i\}(x_i + y_i). $$

This means that

$$ \limsup_{t \to +\infty} x_i(t) + y_i(t) \le \frac{K_i\beta(0)}{\min\{\mu_i, \alpha_i\}}. $$

Therefore, the solution should be necessarily bounded.

We now introduce a formula for the Basic Reproduction Number (BRN) $\mathcal{R}_0$ of the system (6), given as the spectral radius of an $n \times n$ matrix. We do so by applying the Next Generation Matrix method [21, 22, 23]. We remark that this method was developed specifically for systems of ODEs. However, when evaluated in its Disease Free Equilibrium, namely

$$ (x_1, \dots, x_n, y_1, \dots, y_n) = \left( \frac{K_1}{\mu_1}, \dots, \frac{K_n}{\mu_n}, 0, \dots, 0 \right), \quad (7) $$

the system (6) does not exhibit any form of delay, and qualitatively reduces to a system of ODEs. We focus on the Jacobian on the system evaluated in this equilibrium, obtaining a reliable threshold quantity.

**Proposition 2.** Recall from [20] that the “Basic Reproduction Number of neuron i” is $R_{0i} = d \frac{K_i}{\mu_i \alpha_i}$.

The Basic Reproduction Number $\mathcal{R}_0$ of the system (6) is given by the spectral radius of the matrix $F \in \mathbb{R}^{n \times n}$ defined as

$$ (F)_{ij} = \begin{cases} R_{0i} & \text{if } j = i, \\ \kappa_{ji} \alpha_{j \to i} R_{0i} & \text{if } j \neq i. \end{cases} \quad (8) $$

*Proof.* For ease of notation, we use $\text{diag}(\cdot)$ to indicate $\text{diag}(\cdot)_{1 \le i \le n}$, since all the diagonal matrices we consider are of dimension $n \times n$.

We compute the Jacobian $J$ of (6), dropping the explicit dependence on $(t)$ everywhere for ease of notation. By “splitting” the system into $x$ and $y$, we can write

$$ J = \begin{pmatrix} J_{11} & J_{12} \\ J_{21} & J_{22} \end{pmatrix}, $$

where

$$ J_{11} = \text{diag} \left( -\mu_i - d \left( y_i + \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i} y_j \right) \right), \quad (J_{12})_{ij} = \begin{cases} K_i \beta'(y_i) - dx_i & \text{if } j = i, \\ -dx_i \kappa_{ji} \alpha_{j \to i} & \text{if } j \neq i, \end{cases} $$
$$J_{21} = \operatorname{diag} \left( d \left( y_i + \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i} y_j \right) \right) \quad \text{and} \quad (J_{22})_{ij} = \begin{cases} dx_i - \alpha_i & \text{if } j = i, \\ dx_i \kappa_{ji} \alpha_{j \to i} & \text{if } j \neq i. \end{cases}$$

We now evaluate the Jacobian in the Disease Free Equilibrium of System (6), given in (7).

Recall that $\beta'(0) = 0$. We obtain

$$J_{11,DFE} = \operatorname{diag}(-\mu_i), \quad J_{21,DFE} = 0,$$

$$(J_{12,DFE})_{ij} = \begin{cases} -d \frac{K_i}{\mu_i} & \text{if } j = i, \\ -d \frac{K_i}{\mu_i} \kappa_{ji} \alpha_{j \to i} & \text{if } j \neq i, \end{cases} \quad \text{and} \quad (J_{22,DFE})_{ij} = \begin{cases} d \frac{K_i}{\mu_i} - \alpha_i & \text{if } j = i, \\ d \frac{K_i}{\mu_i} \kappa_{ji} \alpha_{j \to i} & \text{if } j \neq i. \end{cases}$$

Finally, we decompose $J_{DFE} = M - V$, with

$$M = \begin{pmatrix} 0 & 0 \\ 0 & M_{22} \end{pmatrix} \quad \text{and} \quad V = \begin{pmatrix} V_{11} & V_{12} \\ 0 & V_{22} \end{pmatrix},$$

where

$$(M_{22})_{ij} = \begin{cases} d \frac{K_i}{\mu_i} & \text{if } j = i, \\ d \frac{K_i}{\mu_i} \kappa_{ji} \alpha_{j \to i} & \text{if } j \neq i, \end{cases}$$

and

$$V_{22} = \operatorname{diag}(\alpha_i).$$

We do not write $V_{11}$ and $V_{12}$ explicitly, since they are not needed for our computations. Then, the basic reproduction number of the whole system is the spectral radius $\mathcal{R}_0 = \rho(M_{22}V_{22}^{-1}) = \rho(F)$, with $F \in \mathbb{R}^{n \times n}$ defined as

$$(F)_{ij} = \begin{cases} R_{0i} & \text{if } j = i, \\ \kappa_{ji} \alpha_{j \to i} R_{0i} & \text{if } j \neq i. \end{cases}$$

□

Notice that, in order to compute this spectral radius, we implicitly assumed that $\alpha_i \neq 0$ for all $i = 1, 2, \dots, n$. This means that we assume that each neuron receives some infection from at least one of its neighbours. We comment more on this in Section 7. We derived the Basic Reproduction Number of System (6) similarly to how we proceeded on page 2 for the 2 neurons case. Here, the influence of the various $\kappa_{ji}$ is less obvious, and we shall investigate it numerically, except for the case $\mathcal{R}_0 < 1$, for which we analytically prove global convergence towards the Disease Free Equilibrium (7) in Section 4.

## 3.1 Fully homogeneous case

Recall that we are interested in the number of neurons $n \ge 2$, so the divisions we make in this section by $n-1$ are not problematic. Assume now that the system is fully homogeneous, and that all the neurons are connected to each other. This is clearly an unrealistic setting, however, it is instructive to obtain an intuition of what the role of $n$, the number of neurons, is in the spread of the prion.

Full homogeneity in this setting means that in the system (6) the parameters are $K_i = K$, $\mu_i = \mu$ for all $i = 1, \dots, n$, and $\kappa_{ij} = \kappa$, $\alpha_{i \to j} = \alpha/(n-1)$ for all $i, j = 1, \dots, n$. Then, $\alpha_i = \alpha$; moreover, for each neuron the local Basic Reproduction Number is

$$R_{0i} = \frac{dK}{\mu\alpha} =: R_0,$$

and the matrix $F$ defining the global Basic Reproduction Number $\mathcal{R}_0$ is given by

$$(F)_{ij} = \begin{cases} \frac{R_0}{\kappa\alpha} & \text{if } j = i \\ \frac{R_0}{n-1} & \text{if } j \neq i \end{cases} = R_0 \begin{cases} 1 & \text{if } j = i \\ \frac{\kappa\alpha}{n-1} & \text{if } j \neq i \end{cases} = R_0 \left( \left(1 - \frac{\kappa\alpha}{n-1}\right) I_n + \frac{\kappa\alpha}{n-1} I_n \right), \quad (9)$$
$$J_{21} = \operatorname{diag} \left( d \left( y_i + \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i} y_j \right) \right) \quad \text{and} \quad (J_{22})_{ij} = \begin{cases} dx_i - \alpha_i & \text{if } j = i, \\ dx_i \kappa_{ji} \alpha_{j \to i} & \text{if } j \neq i. \end{cases}$$

We now evaluate the Jacobian in the Disease Free Equilibrium of System (6), given in (7).

Recall that $\beta'(0) = 0$. We obtain

$$J_{11,DFE} = \operatorname{diag}(-\mu_i), \quad J_{21,DFE} = 0,$$

$$(J_{12,DFE})_{ij} = \begin{cases} -d \frac{K_i}{\mu_i} & \text{if } j = i, \\ -d \frac{K_i}{\mu_i} \kappa_{ji} \alpha_{j \to i} & \text{if } j \neq i, \end{cases} \quad \text{and} \quad (J_{22,DFE})_{ij} = \begin{cases} d \frac{K_i}{\mu_i} - \alpha_i & \text{if } j = i, \\ d \frac{K_i}{\mu_i} \kappa_{ji} \alpha_{j \to i} & \text{if } j \neq i. \end{cases}$$

Finally, we decompose $J_{DFE} = M - V$, with

$$M = \begin{pmatrix} 0 & 0 \\ 0 & M_{22} \end{pmatrix} \quad \text{and} \quad V = \begin{pmatrix} V_{11} & V_{12} \\ 0 & V_{22} \end{pmatrix},$$

where

$$(M_{22})_{ij} = \begin{cases} d \frac{K_i}{\mu_i} & \text{if } j = i, \\ d \frac{K_i}{\mu_i} \kappa_{ji} \alpha_{j \to i} & \text{if } j \neq i, \end{cases}$$

and

$$V_{22} = \operatorname{diag}(\alpha_i).$$

We do not write $V_{11}$ and $V_{12}$ explicitly, since they are not needed for our computations. Then, the basic reproduction number of the whole system is the spectral radius $\mathcal{R}_0 = \rho(M_{22}V_{22}^{-1}) = \rho(F)$, with $F \in \mathbb{R}^{n \times n}$ defined as

$$(F)_{ij} = \begin{cases} R_{0i} & \text{if } j = i, \\ \kappa_{ji} \alpha_{j \to i} R_{0i} & \text{if } j \neq i. \end{cases}$$

□

Notice that, in order to compute this spectral radius, we implicitly assumed that $\alpha_i \neq 0$ for all $i = 1, 2, \dots, n$. This means that we assume that each neuron receives some infection from at least one of its neighbours. We comment more on this in Section 7. We derived the Basic Reproduction Number of System (6) similarly to how we proceeded on page 2 for the 2 neurons case. Here, the influence of the various $\kappa_{ji}$ is less obvious, and we shall investigate it numerically, except for the case $\mathcal{R}_0 < 1$, for which we analytically prove global convergence towards the Disease Free Equilibrium (7) in Section 4.

## 3.1 Fully homogeneous case

Recall that we are interested in the number of neurons $n \ge 2$, so the divisions we make in this section by $n-1$ are not problematic. Assume now that the system is fully homogeneous, and that all the neurons are connected to each other. This is clearly an unrealistic setting, however, it is instructive to obtain an intuition of what the role of $n$, the number of neurons, is in the spread of the prion.

Full homogeneity in this setting means that in the system (6) the parameters are $K_i = K$, $\mu_i = \mu$ for all $i = 1, \dots, n$, and $\kappa_{ij} = \kappa$, $\alpha_{i \to j} = \alpha/(n-1)$ for all $i, j = 1, \dots, n$. Then, $\alpha_i = \alpha$; moreover, for each neuron the local Basic Reproduction Number is

$$R_{0i} = \frac{dK}{\mu\alpha} =: R_0,$$

and the matrix $F$ defining the global Basic Reproduction Number $\mathcal{R}_0$ is given by

$$(F)_{ij} = \begin{cases} \frac{R_0}{\kappa\alpha} & \text{if } j = i \\ \frac{R_0}{n-1} & \text{if } j \neq i \end{cases} = R_0 \begin{cases} 1 & \text{if } j = i \\ \frac{\kappa\alpha}{n-1} & \text{if } j \neq i \end{cases} = R_0 \left( \left(1 - \frac{\kappa\alpha}{n-1}\right) I_n + \frac{\kappa\alpha}{n-1} I_n \right), \quad (9)$$
where $I_n$ is the $n \times n$ identity matrix, and $1_n$ is the $n \times n$ matrix with 1 in all its entries.

Then, $1_n$ has one eigenvalue $n$ (its trace) and $n-1$ zero eigenvalues (since it has rank 1), whereas the matrix

$$ \left( 1 - \frac{\kappa\alpha}{n-1} \right) I_n, $$

clearly has $n$ eigenvalues equal to $1 - \kappa\alpha/(n-1)$.

Recall that, if a matrix $A$ has eigenvalues $\lambda_1, \dots, \lambda_n$, then the matrix $cI_n + bA$ has eigenvalues $c+b\lambda_1, \dots, c+b\lambda_n$, for any $b, c \in \mathbb{R}$, since any eigenvector $v$ of $A$ will also satisfy $cIv = cv$.

Hence, the sum (9) (ignoring for a moment the scalar coefficient $R_0$ in front of the brackets) has one eigenvalue equal to $\kappa\alpha + 1$ (its spectral radius) and $n-1$ eigenvalues equal to $1 - \kappa\alpha/(n-1)$. Consequently,

$$ \mathcal{R}_0 = \rho(F) = R_0 (\kappa\alpha + 1). \qquad (10) $$

This value is clearly strictly greater than $R_0$, and independent on $n$. This means that, as long as the network is fully connected and fully homogeneous, the number of neurons has no direct impact on the dynamics of the system, according to our model. However, the connectivity between the neurons, expressed by $\kappa$, plays a fundamental role in the dynamics in the sense that even if all $R_0 < 1$, with $\kappa$ big enough the global system might have a Basic Reproduction Number $\mathcal{R}_0 > 1$, and the disease could hence remain endemic. Figure 4a in Section 7 is an example of a fully connected network with $n=3$.

<figure><img src="image_4.png" /><figcaption>Figure 3: The three networks we consider in the sections 3.2, 3.3 and 3.4. We only show the forward case from Section 3.4, as the backward case would look exactly the same but with each arrow reversed.</figcaption></figure>

## 3.2 Case of ring network

Next, we consider a case of a ring network with homogeneous coefficients. In this scenario, as illustrated in Figure 3a, each neuron $i$ is only connected to the previous one ($i-1$) and the following one ($i+1$), as well as the $n$-th neuron being connected to the first. This leads then to the following expression

$$ (F)_{ij} = R_0 \begin{cases} 1 & \text{if } j=i, \\ \frac{\kappa\alpha}{2} & \text{if } j=i-1, 2 \le i \le n, \\ \frac{\kappa\alpha}{2} & \text{if } j=i+1, 1 \le i \le n-1, \\ \frac{\kappa\alpha}{2} & \text{if } j=1, i=n, \\ \frac{\kappa\alpha}{2} & \text{if } j=n, i=1, \\ 0 & \text{otherwise.} \end{cases} $$
where $I_n$ is the $n \times n$ identity matrix, and $1_n$ is the $n \times n$ matrix with 1 in all its entries.

Then, $1_n$ has one eigenvalue $n$ (its trace) and $n-1$ zero eigenvalues (since it has rank 1), whereas the matrix

$$ \left( 1 - \frac{\kappa\alpha}{n-1} \right) I_n, $$

clearly has $n$ eigenvalues equal to $1 - \kappa\alpha/(n-1)$.

Recall that, if a matrix $A$ has eigenvalues $\lambda_1, \dots, \lambda_n$, then the matrix $cI_n + bA$ has eigenvalues $c+b\lambda_1, \dots, c+b\lambda_n$, for any $b, c \in \mathbb{R}$, since any eigenvector $v$ of $A$ will also satisfy $cIv = cv$.

Hence, the sum (9) (ignoring for a moment the scalar coefficient $R_0$ in front of the brackets) has one eigenvalue equal to $\kappa\alpha + 1$ (its spectral radius) and $n-1$ eigenvalues equal to $1 - \kappa\alpha/(n-1)$. Consequently,

$$ \mathcal{R}_0 = \rho(F) = R_0 (\kappa\alpha + 1). \qquad (10) $$

This value is clearly strictly greater than $R_0$, and independent on $n$. This means that, as long as the network is fully connected and fully homogeneous, the number of neurons has no direct impact on the dynamics of the system, according to our model. However, the connectivity between the neurons, expressed by $\kappa$, plays a fundamental role in the dynamics in the sense that even if all $R_0 < 1$, with $\kappa$ big enough the global system might have a Basic Reproduction Number $\mathcal{R}_0 > 1$, and the disease could hence remain endemic. Figure 4a in Section 7 is an example of a fully connected network with $n=3$.

<figure><img src="image_4.png" /><figcaption>Figure 3: The three networks we consider in the sections 3.2, 3.3 and 3.4. We only show the forward case from Section 3.4, as the backward case would look exactly the same but with each arrow reversed.</figcaption></figure>

## 3.2 Case of ring network

Next, we consider a case of a ring network with homogeneous coefficients. In this scenario, as illustrated in Figure 3a, each neuron $i$ is only connected to the previous one ($i-1$) and the following one ($i+1$), as well as the $n$-th neuron being connected to the first. This leads then to the following expression

$$ (F)_{ij} = R_0 \begin{cases} 1 & \text{if } j = i, \\ \frac{\kappa\alpha}{2} & \text{if } j = i-1, 2 \le i \le n, \\ \frac{\kappa\alpha}{2} & \text{if } j = i+1, 1 \le i \le n-1, \\ \frac{\kappa\alpha}{2} & \text{if } j = 1, i = n, \\ \frac{\kappa\alpha}{2} & \text{if } j = n, i = 1, \\ 0 & \text{otherwise.} \end{cases} $$
$F$ is a symmetric circulant matrix (see [27]); its eigenvalues are given by

$$ \lambda_k = R_0 \left( 1 + \kappa \alpha \cos \left( \frac{2k\pi}{n} \right) \right), \quad k = 0, \dots, n-1. $$

Then, the spectral radius of the matrix $F$ is given by

$$ \mathcal{R}_0 = \rho(F) = R_0 (1 + \kappa \alpha). $$

As with a fully connected network, connectivity between neurons plays an important role in the dynamics of a ring network. Figure 4c in Section 7 is a particular case of ring network with $n = 5$.

### 3.3 Case of line network

Then, we consider a case of a line network with homogeneous coefficients. This corresponds to the ring network considered in Section 3.2, from which we remove the connection between neuron $n$ and neuron 1. This scenario is illustrated in Figure 3b. This leads then to the following expression

$$ (F)_{ij} = R_0 \begin{cases} 1 & \text{if } j = i, \\ \frac{\kappa\alpha}{2} & \text{if } j = i - 1, 2 \le i \le n, \\ \frac{\kappa\alpha}{2} & \text{if } j = i + 1, 1 \le i \le n - 1, \\ 0 & \text{otherwise.} \end{cases} $$

We obtain a tridiagonal symmetric Toeplitz matrix (see [28]). Then, the eigenvalues of $F$ are given by

$$ \lambda_k = R_0 \left( 1 + \kappa \alpha \cos \left( \frac{k\pi}{n+1} \right) \right), \quad k = 1, \dots, n, $$

and the spectral radius is

$$ \mathcal{R}_0 = \rho(F) = R_0 \left( 1 + \kappa \alpha \cos \left( \frac{\pi}{n+1} \right) \right). $$

We see here that not only the connectivity $\kappa$ between neurons plays a fundamental role, but also the number $n$ of neurons involved. This network acts like a ring network when $n$ is very large ($n \to +\infty$). Figures 4b and 4d in Section 7 are particular cases of line networks with $n = 5$ and $n = 9$, respectively.

### 3.4 Case of one-way direction

Lastly, we consider a case in which each neuron is only connected to the ones preceding (or following) it, meaning neuron $i$ only spreads the infection to neurons $j < i$ (or $j > i$), with homogeneous coefficients. This leads then to the following expressions,

$$ (F)_{ij} = R_0 \begin{cases} 1 & \text{if } i = j, \\ \frac{\kappa\alpha}{n-j} & \text{if } i \ge j + 1, j \le n - 1, \\ 0 & \text{otherwise,} \end{cases} \quad \text{or} \quad (F)_{ij} = R_0 \begin{cases} 1 & \text{if } i = j, \\ \frac{\kappa\alpha}{j-1} & \text{if } i \le j - 1, j \ge 2, \\ 0 & \text{otherwise.} \end{cases} $$

In this case, since the matrix $F$ is triangular, we obtain

$$ \mathcal{R}_0 = \rho(F) = R_0. $$

The former case is illustrated in Figure 3c. Unlike fully connected networks, ring networks and line networks, the global Basic Reproduction Number $\mathcal{R}_0$ of the one-way direction network is independent of the connectivity between neurons.
$F$ is a symmetric circulant matrix (see [27]); its eigenvalues are given by

$$ \lambda_k = R_0 \left( 1 + \kappa \alpha \cos \left( \frac{2k\pi}{n} \right) \right), \quad k = 0, \dots, n-1. $$

Then, the spectral radius of the matrix $F$ is given by

$$ \mathcal{R}_0 = \rho(F) = R_0 (1 + \kappa \alpha). $$

As with a fully connected network, connectivity between neurons plays an important role in the dynamics of a ring network. Figure 4c in Section 7 is a particular case of ring network with $n = 5$.

### 3.3 Case of line network

Then, we consider a case of a line network with homogeneous coefficients. This corresponds to the ring network considered in Section 3.2, from which we remove the connection between neuron $n$ and neuron 1. This scenario is illustrated in Figure 3b. This leads then to the following expression

$$ (F)_{ij} = R_0 \begin{cases} 1 & \text{if } j = i, \\ \frac{\kappa \alpha}{2} & \text{if } j = i - 1, 2 \le i \le n, \\ \frac{\kappa \alpha}{2} & \text{if } j = i + 1, 1 \le i \le n - 1, \\ 0 & \text{otherwise.} \end{cases} $$

We obtain a tridiagonal symmetric Toeplitz matrix (see [28]). Then, the eigenvalues of $F$ are given by

$$ \lambda_k = R_0 \left( 1 + \kappa \alpha \cos \left( \frac{k\pi}{n+1} \right) \right), \quad k = 1, \dots, n, $$

and the spectral radius is

$$ \mathcal{R}_0 = \rho(F) = R_0 \left( 1 + \kappa \alpha \cos \left( \frac{\pi}{n+1} \right) \right). $$

We see here that not only the connectivity $\kappa$ between neurons plays a fundamental role, but also the number $n$ of neurons involved. This network acts like a ring network when $n$ is very large ($n \to +\infty$). Figures 4b and 4d in Section 7 are particular cases of line networks with $n = 5$ and $n = 9$, respectively.

### 3.4 Case of one-way direction

Lastly, we consider a case in which each neuron is only connected to the ones preceding (or following) it, meaning neuron $i$ only spreads the infection to neurons $j < i$ (or $j > i$), with homogeneous coefficients. This leads then to the following expressions,

$$ (F)_{ij} = R_0 \begin{cases} 1 & \text{if } i = j, \\ \frac{\kappa \alpha}{n-j} & \text{if } i \ge j + 1, j \le n - 1, \\ 0 & \text{otherwise,} \end{cases} \quad \text{or} \quad (F)_{ij} = R_0 \begin{cases} 1 & \text{if } i = j, \\ \frac{\kappa \alpha}{j-1} & \text{if } i \le j - 1, j \ge 2, \\ 0 & \text{otherwise.} \end{cases} $$

In this case, since the matrix $F$ is triangular, we obtain

$$ \mathcal{R}_0 = \rho(F) = R_0. $$

The former case is illustrated in Figure 3c. Unlike fully connected networks, ring networks and line networks, the global Basic Reproduction Number $\mathcal{R}_0$ of the one-way direction network is independent of the connectivity between neurons.
# 4 Global stability of the Disease Free Equilibrium

In this section, we prove the global stability of the Disease Free Equilibrium (7) for System (6) when $\mathcal{R}_0 = \rho(F) < 1$. In order to do so, we proceed similarly to [29, Thm. 5].

**Theorem 1.** The Disease Free Equilibrium (7) of System (6) is globally asymptotically stable when $\mathcal{R}_0 = \rho(F) < 1$.

*Proof.* Recall from (2) that $\beta(x) \le 1$ for all $x \ge 0$. Then, we can bound the first $n$ DDEs of System (6) from above by

$$ \frac{dx_i}{dt} = K_i\beta(y_i(t-T_i)) - \mu_i x_i(t) - dx_i \left( y_i(t) + \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i} y_j(t) \right) \le K_i - \mu_i x_i(t). $$

Consider the auxiliary system

$$ \frac{dz_i}{dt} = K_i - \mu_i z_i(t), \quad i = 1, 2, \dots, n. \qquad (11) $$

Clearly, the first $n$ entries of Disease Free Equilibrium (7) form a point which is globally asymptotically stable for (11). Then, for any $\varepsilon > 0$, there exists a $\bar{t}_i > 0$ such that, for $t \ge \bar{t}_i$,

$$ x_i(t) \le \frac{K_i}{\mu_i} + \varepsilon. $$

Take $\bar{t} = \max_i \bar{t}_i$. Then, for $t \ge \bar{t}$, the second $n$ ODEs of the system (6) can be bound from above by

$$ \frac{dy_i}{dt} \le d \left( \frac{K_i}{\mu_i} + \varepsilon \right) \left( y_i(t) + \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i} y_j(t) \right) - \alpha_i y_i(t). $$

Consider the second auxiliary system

$$ \frac{dw_i}{dt} = d \left( \frac{K_i}{\mu_i} + \varepsilon \right) \left( w_i(t) + \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i} w_j(t) \right) - \alpha_i w_i(t). $$

This system is linear in $w = (w_1, w_2, \dots, w_n)$, and can be rewritten as

$$ \frac{dw}{dt} = (M_{22}(\varepsilon) - V_{22})w, $$

where

$$ (M_{22}(\varepsilon))_{ij} = \begin{cases} d\left(\frac{K_i}{\mu_i} + \varepsilon\right) & \text{if } j = i, \\ d\left(\frac{K_i}{\mu_i} + \varepsilon\right) \kappa_{ji} \alpha_{j \to i} & \text{if } j \neq i, \end{cases} $$

meaning the matrix $M_{22}$ used in the definition of $\mathcal{R}_0$ is actually $M_{22}(0)$, and

$$ V_{22} = \mathrm{diag}(\alpha_i), $$

as above. For $\varepsilon > 0$ small enough, as a consequence of our assumption $\mathcal{R}_0 < 1$, we can have $\rho(M_{22}(\varepsilon)V_{22}^{-1}) < 1$. We then use the following lemma:

**Lemma 1** ([30], Lemma 2). If $M$ is non-negative and $V$ is a non-singular $M$-matrix, then $\mathcal{R}_0 = \rho(MV^{-1}) < 1$ if and only if all eigenvalues of $(M - V)$ have negative real parts.
# 4 Global stability of the Disease Free Equilibrium

In this section, we prove the global stability of the Disease Free Equilibrium (7) for System (6) when $\mathcal{R}_0 = \rho(F) < 1$. In order to do so, we proceed similarly to [29, Thm. 5].

**Theorem 1.** The Disease Free Equilibrium (7) of System (6) is globally asymptotically stable when $\mathcal{R}_0 = \rho(F) < 1$.

*Proof.* Recall from (2) that $\beta(x) \le 1$ for all $x \ge 0$. Then, we can bound the first $n$ DDEs of System (6) from above by

$$ \frac{dx_i}{dt} = K_i\beta(y_i(t-T_i)) - \mu_i x_i(t) - dx_i \left( y_i(t) + \sum_{j \neq i} \kappa_{ji}\alpha_{j \to i}y_j(t) \right) \le K_i - \mu_i x_i(t). $$

Consider the auxiliary system

$$ \frac{dz_i}{dt} = K_i - \mu_i z_i(t), \quad i = 1, 2, \dots, n. \qquad (11) $$

Clearly, the first $n$ entries of Disease Free Equilibrium (7) form a point which is globally asymptotically stable for (11). Then, for any $\varepsilon > 0$, there exists a $\bar{t}_i > 0$ such that, for $t \ge \bar{t}_i$,

$$ x_i(t) \le \frac{K_i}{\mu_i} + \varepsilon. $$

Take $\bar{t} = \max_i \bar{t}_i$. Then, for $t \ge \bar{t}$, the second $n$ ODEs of the system (6) can be bound from above by

$$ \frac{dy_i}{dt} \le d \left( \frac{K_i}{\mu_i} + \varepsilon \right) \left( y_i(t) + \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i} y_j(t) \right) - \alpha_i y_i(t). $$

Consider the second auxiliary system

$$ \frac{dw_i}{dt} = d \left( \frac{K_i}{\mu_i} + \varepsilon \right) \left( w_i(t) + \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i} w_j(t) \right) - \alpha_i w_i(t). $$

This system is linear in $w = (w_1, w_2, \dots, w_n)$, and can be rewritten as

$$ \frac{dw}{dt} = (M_{22}(\varepsilon) - V_{22})w, $$

where

$$ (M_{22}(\varepsilon))_{ij} = \begin{cases} d\left(\frac{K_i}{\mu_i} + \varepsilon\right) & \text{if } j = i, \\ d\left(\frac{K_i}{\mu_i} + \varepsilon\right) \kappa_{ji} \alpha_{j \to i} & \text{if } j \neq i, \end{cases} $$

meaning the matrix $M_{22}$ used in the definition of $\mathcal{R}_0$ is actually $M_{22}(0)$, and

$$ V_{22} = \mathrm{diag}(\alpha_i), $$

as above. For $\varepsilon > 0$ small enough, as a consequence of our assumption $\mathcal{R}_0 < 1$, we can have $\rho(M_{22}(\varepsilon)V_{22}^{-1}) < 1$. We then use the following lemma:

**Lemma 1** ([30], Lemma 2). If $M$ is non-negative and $V$ is a non-singular $M$-matrix, then $\mathcal{R}_0 = \rho(MV^{-1}) < 1$ if and only if all eigenvalues of $(M - V)$ have negative real parts.
This means that, if $\rho(M_{22}(\varepsilon)V_{22}^{-1}) < 1$, then

$$ \lim_{t \to +\infty} w_i(t) = 0 $$

for all $i = 1, 2, \dots, n$, which implies

$$ \lim_{t \to +\infty} y_i(t) = 0. $$

Thus, for any $\delta > 0$, there exists $t^* > 0$ such that, for all $t \ge t^*$ and for all $i = 1, 2, \dots, n$, we have $y_i(t) \le \delta$. Hence, introducing for ease of notation $T = \max_i T_i$, for $t \ge t^* + T$, we have

$$ \beta(y_i(t - T_i)) = \frac{1}{1 + (y_i(t - T_i)/y_c)^p} \ge \frac{1}{1 + (\delta/y_c)^p} = \beta(\delta). $$

Notice that $\beta(\delta) \to 1$ as $\delta \to 0$. We can then bound the first $n$ DDEs of System (6) from below by

$$ \frac{dx_i}{dt} \ge K_i\beta(\delta) - \mu_i x_i(t) - dx_i(t) \left( \delta + \sum_{j \ne i} \kappa_{ji} \alpha_{j \to i} \delta \right). $$

Consider the final auxiliary system

$$ \frac{dv_i}{dt} = K_i\beta(\delta) - \mu_i v_i(t) - dv_i(t) \left( \delta + \sum_{j \ne i} \kappa_{ji} \alpha_{j \to i} \delta \right). $$

Clearly, each for each $i$ we have

$$ \lim_{t \to +\infty} v_i(t) = \frac{K_i \beta(\delta)}{\mu_i + \delta d \left(1 + \sum_{j \ne i} \kappa_{ji} \alpha_{j \to i}\right)}. $$

Hence, for each $i = 1, 2, \dots, n$ and for all $\varepsilon, \delta > 0$, we have the following lower and upper bounds:

$$ \frac{K_i \beta(\delta)}{\mu_i + \delta d \left(1 + \sum_{j \ne i} \kappa_{ji} \alpha_{j \to i}\right)} \le \liminf_{t \to +\infty} x_i(t) \le \limsup_{t \to +\infty} x_i(t) \le \frac{K_i}{\mu_i} + \varepsilon. $$

Letting $\varepsilon, \delta \to 0$ concludes the proof.

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

**Corollary 1.** *The Disease Free Equilibrium is locally unstable when $R_0 > 1$.*

*Proof.* Direct consequence of Theorem 1 and [30, Thm. 1].

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

We remark that, for all our results thus far, the only assumptions on the function $\beta(\cdot)$ are: $\beta(0) = 1$, $\beta'(0) = 0$ and $\beta(x)$ decreasing in $x$. Our specific choice (2) was made for consistency with [20] and because it appears biologically relevant. However, other choices might lead to interesting results. We comment more on this in Section 8.

# 5 Existence of an endemic equilibrium

We now prove, under stronger assumptions than $R_0 > 1$ (but weaker than $R_{0i} > 1$), that System (6) admits at least one Endemic Equilibrium (EE), i.e. an equilibrium such that $y_i > 0$ for all $i$.

**Theorem 2.** *Assume that the matrix $F(8)$ is such that the minimum row sum is strictly bigger than 1. Then, System (6) admits at least one Endemic Equilibrium.*
This means that, if $\rho(M_{22}(\varepsilon)V_{22}^{-1}) < 1$, then

$$ \lim_{t \to +\infty} w_i(t) = 0 $$

for all $i = 1, 2, \dots, n$, which implies

$$ \lim_{t \to +\infty} y_i(t) = 0. $$

Thus, for any $\delta > 0$, there exists $t^* > 0$ such that, for all $t \ge t^*$ and for all $i = 1, 2, \dots, n$, we have $y_i(t) \le \delta$. Hence, introducing for ease of notation $T = \max_i T_i$, for $t \ge t^* + T$, we have

$$ \beta(y_i(t - T_i)) = \frac{1}{1 + (y_i(t - T_i)/y_c)^p} \ge \frac{1}{1 + (\delta/y_c)^p} = \beta(\delta). $$

Notice that $\beta(\delta) \to 1$ as $\delta \to 0$. We can then bound the first $n$ DDEs of System (6) from below by

$$ \frac{dx_i}{dt} \ge K_i\beta(\delta) - \mu_i x_i(t) - dx_i(t) \left( \delta + \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i} \delta \right). $$

Consider the final auxiliary system

$$ \frac{dv_i}{dt} = K_i\beta(\delta) - \mu_i v_i(t) - dv_i(t) \left( \delta + \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i} \delta \right). $$

Clearly, each for each $i$ we have

$$ \lim_{t \to +\infty} v_i(t) = \frac{K_i \beta(\delta)}{\mu_i + \delta d \left(1 + \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i}\right)}. $$

Hence, for each $i = 1, 2, \dots, n$ and for all $\varepsilon, \delta > 0$, we have the following lower and upper bounds:

$$ \frac{K_i \beta(\delta)}{\mu_i + \delta d \left(1 + \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i}\right)} \le \liminf_{t \to +\infty} x_i(t) \le \limsup_{t \to +\infty} x_i(t) \le \frac{K_i}{\mu_i} + \varepsilon. $$

Letting $\varepsilon, \delta \to 0$ concludes the proof.

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

**Corollary 1.** *The Disease Free Equilibrium is locally unstable when $\mathcal{R}_0 > 1$.*

*Proof.* Direct consequence of Theorem 1 and [30, Thm. 1].

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

We remark that, for all our results thus far, the only assumptions on the function $\beta(\cdot)$ are: $\beta(0) = 1$, $\beta'(0) = 0$ and $\beta(x)$ decreasing in $x$. Our specific choice (2) was made for consistency with [20] and because it appears biologically relevant. However, other choices might lead to interesting results. We comment more on this in Section 8.

# 5 Existence of an endemic equilibrium

We now prove, under stronger assumptions than $\mathcal{R}_0 > 1$ (but weaker than $R_{0i} > 1$), that System (6) admits at least one Endemic Equilibrium (EE), i.e. an equilibrium such that $y_i > 0$ for all $i$.

**Theorem 2.** *Assume that the matrix $F(8)$ is such that the minimum row sum is strictly bigger than 1. Then, System (6) admits at least one Endemic Equilibrium.*
Proof. We know that

min row/column sum of $F \le \rho(F) \le \max$ row/column sum of $F$,

hence under our assumption, $\rho(F) = \mathcal{R}_0 > 1$.

We begin by noticing that an equilibrium of the system (6) necessarily satisfies

$$x_i = \frac{K_i \beta(y_i)}{\mu_i + d \left(y_i + \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i} y_j\right)} \quad (12)$$

Substituting (12) in the ODEs for $y_i$ and equating them to 0, we obtain

$$0 = \frac{dK_i\beta(y_i)\left(y_i + \sum_{j \neq i} \kappa_{ji}\alpha_{j \to i}y_j\right)}{\mu_i + d\left(y_i + \sum_{j \neq i} \kappa_{ji}\alpha_{j \to i}y_j\right)} - \alpha_i y_i \quad (13)$$

Notice that, for $y_i$ large enough, the right hand side (RHS) of (13) is clearly negative. Let us denote with $M_i$ a large number such that

$$\left( \frac{dK_i\beta(y_i)\left(y_i + \sum_{j \neq i} \kappa_{ji}\alpha_{j \to i}y_j\right)}{\mu_i + d\left(y_i + \sum_{j \neq i} \kappa_{ji}\alpha_{j \to i}y_j\right)} - \alpha_i y_i \right) \Big|_{y_i=M_i} < 0,$$

for all non-negative values of $y_j, j \neq i$. Moreover, let us denote with $M = \max M_i$.

If we find a value $\varepsilon > 0$ such that the RHS of (13) is positive for all $i = 1, 2, \dots, n$, we can apply the Poincaré-Miranda theorem [31, 32] (qualitatively, a higher dimensional version of the intermediate value theorem) to conclude the existence of at least one Endemic Equilibrium of System (6).

Let us evaluate the RHS of (13) at $y_i = \varepsilon$ for all $i = 1, 2, \dots, n$, and study its sign. We have

$$\frac{dK_i\beta(\varepsilon)\left(\varepsilon + \varepsilon \sum_{j \neq i} \kappa_{ji}\alpha_{j \to i}\right)}{\mu_i + d\left(\varepsilon + \varepsilon \sum_{j \neq i} \kappa_{ji}\alpha_{j \to i}\right)} - \alpha_i\varepsilon > 0.$$

In fact, we can divide by $\varepsilon > 0$ on both sides, obtaining

$$\frac{dK_i\beta(\varepsilon)\left(1 + \sum_{j \neq i} \kappa_{ji}\alpha_{j \to i}\right)}{\mu_i + d\left(\varepsilon + \varepsilon \sum_{j \neq i} \kappa_{ji}\alpha_{j \to i}\right)} - \alpha_i > 0. \quad (14)$$

Recall that $\beta(0) = 1$. Then, for $\varepsilon = 0$, (14) coincides with the $i$-th row sum of $F$ being strictly greater than 1. Since by assumption the minimum of the row sums (hence, all the row sums) is greater than 1, by continuity there exists a small $\varepsilon_i > 0$ such that the RHS of (13) is strictly positive. Let us denote with $\varepsilon = \min \varepsilon_i$.

Applying the Poincaré-Miranda theorem on the set $[\varepsilon, M]^n$ allows us to conclude the existence of at least one Endemic Equilibrium, i.e. with $0 < \varepsilon < y_i < M$ for $i = 1, 2, \dots, n$. $\square$

We conjecture the following, based on our extensive numerical simulations:

**Conjecture 1.** System (6) admits at least one Endemic Equilibrium when $\mathcal{R}_0 > 1$.

Our proof of Theorem 2 relies heavily on the assumption on the minimum row sum being strictly bigger than 1, hence it fails so for a generic matrix $F$, if we only assume $\rho(F) > 1$. However, the application of the Poincaré-Miranda theorem might not be necessary to prove this result.
Proof. We know that

min row/column sum of $F \le \rho(F) \le \max$ row/column sum of $F$,

hence under our assumption, $\rho(F) = \mathcal{R}_0 > 1$.

We begin by noticing that an equilibrium of the system (6) necessarily satisfies

$$x_i = \frac{K_i \beta(y_i)}{\mu_i + d \left(y_i + \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i} y_j\right)} \quad (12)$$

Substituting (12) in the ODEs for $y_i$ and equating them to 0, we obtain

$$0 = \frac{dK_i\beta(y_i)\left(y_i + \sum_{j \neq i} \kappa_{ji}\alpha_{j \to i}y_j\right)}{\mu_i + d\left(y_i + \sum_{j \neq i} \kappa_{ji}\alpha_{j \to i}y_j\right)} - \alpha_i y_i \quad (13)$$

Notice that, for $y_i$ large enough, the right hand side (RHS) of (13) is clearly negative. Let us denote with $M_i$ a large number such that

$$\left( \frac{dK_i\beta(y_i)\left(y_i + \sum_{j \neq i} \kappa_{ji}\alpha_{j \to i}y_j\right)}{\mu_i + d\left(y_i + \sum_{j \neq i} \kappa_{ji}\alpha_{j \to i}y_j\right)} - \alpha_i y_i \right) \Big|_{y_i=M_i} < 0,$$

for all non-negative values of $y_j, j \neq i$. Moreover, let us denote with $M = \max M_i$.

If we find a value $\varepsilon > 0$ such that the RHS of (13) is positive for all $i = 1, 2, \dots, n$, we can apply the Poincaré-Miranda theorem [31, 32] (qualitatively, a higher dimensional version of the intermediate value theorem) to conclude the existence of at least one Endemic Equilibrium of System (6).

Let us evaluate the RHS of (13) at $y_i = \varepsilon$ for all $i = 1, 2, \dots, n$, and study its sign. We have

$$\frac{dK_i\beta(\varepsilon)\left(\varepsilon + \varepsilon \sum_{j \neq i} \kappa_{ji}\alpha_{j \to i}\right)}{\mu_i + d\left(\varepsilon + \varepsilon \sum_{j \neq i} \kappa_{ji}\alpha_{j \to i}\right)} - \alpha_i\varepsilon > 0.$$

In fact, we can divide by $\varepsilon > 0$ on both sides, obtaining

$$\frac{dK_i\beta(\varepsilon)\left(1 + \sum_{j \neq i} \kappa_{ji}\alpha_{j \to i}\right)}{\mu_i + d\left(\varepsilon + \varepsilon \sum_{j \neq i} \kappa_{ji}\alpha_{j \to i}\right)} - \alpha_i > 0. \quad (14)$$

Recall that $\beta(0) = 1$. Then, for $\varepsilon = 0$, (14) coincides with the $i$-th row sum of $F$ being strictly greater than 1. Since by assumption the minimum of the row sums (hence, all the row sums) is greater than 1, by continuity there exists a small $\varepsilon_i > 0$ such that the RHS of (13) is strictly positive. Let us denote with $\varepsilon = \min \varepsilon_i$.

Applying the Poincaré-Miranda theorem on the set $[\varepsilon, M]^n$ allows us to conclude the existence of at least one Endemic Equilibrium, i.e. with $0 < \varepsilon < y_i < M$ for $i = 1, 2, \dots, n$. $\square$

We conjecture the following, based on our extensive numerical simulations:

**Conjecture 1.** System (6) admits at least one Endemic Equilibrium when $\mathcal{R}_0 > 1$.

Our proof of Theorem 2 relies heavily on the assumption on the minimum row sum being strictly bigger than 1, hence it fails so for a generic matrix $F$, if we only assume $\rho(F) > 1$. However, the application of the Poincaré-Miranda theorem might not be necessary to prove this result.
# 6 Persistence of solutions

In this section, we treat the long-term behavior of the system (6) when $\mathcal{R}_0 > 1$. We start with the following proposition:

**Proposition 3.** Consider a fixed $i = \tilde{i} \in \{1, \dots, n\}$. If $R_{0\tilde{i}} = dK_{\tilde{i}}/(\mu_{\tilde{i}}\alpha_{\tilde{i}}) > 1$, then $\mathcal{R}_0 > 1$ and there exists a constant $\varepsilon_{\tilde{i}} > 0$ such that

$$ \limsup_{t \to +\infty} y_{\tilde{i}}(t) > \varepsilon_{\tilde{i}}, \quad \text{with } \varphi_{\tilde{i}} \in C([-T, 0], \mathbb{R}_{+}), \quad \varphi_{\tilde{i}}(0) \neq 0. $$

*Proof.* By considering the second equation of System (6) for $i = \tilde{i}$, we have, for $t > 0$,

$$ \frac{dy_{\tilde{i}}}{dt} \geq dx_{\tilde{i}}(t)y_{\tilde{i}}(t) - \alpha_i y_{\tilde{i}}(t). $$

We have also, for $t > 0$,

$$ \frac{d(x_{\tilde{i}} + y_{\tilde{i}})}{dt} = K_{\tilde{i}}\beta(y_{\tilde{i}}(t - T_1)) - \mu_{\tilde{i}}x_{\tilde{i}}(t) - \alpha_{\tilde{i}}y_{\tilde{i}}(t). $$

We suppose by contradiction that $\limsup_{t \to +\infty} y_{\tilde{i}}(t) \le \varepsilon_{\tilde{i}}$, for any small $\varepsilon_{\tilde{i}} > 0$. By the boundedness of solutions, we consider $\liminf_{t \to +\infty} x_{\tilde{i}}(t) = x_{\tilde{i}\infty}$ and $\liminf_{t \to +\infty} y_{\tilde{i}}(t) = y_{\tilde{i}\infty} = 0$. Then, there exists a sequence $t_k \to +\infty$ as $k \to +\infty$, such that $x_{\tilde{i}}(t_k) \to x_{\tilde{i}\infty}$, $y_{\tilde{i}}(t_k) \to 0$, $y'_{\tilde{i}}(t_k) \to 0$ and $x'_{\tilde{i}}(t_k) \to 0$. This yields to

$$ 0 \geq K_{\tilde{i}}\beta(0) - \mu_{\tilde{i}}x_{\tilde{i}\infty} \Rightarrow x_{\tilde{i}\infty} \geq \frac{K_{\tilde{i}}\beta(0)}{\mu_{\tilde{i}}} = \frac{K_{\tilde{i}}}{\mu_{\tilde{i}}}. $$

For a very large time $t$, the ODE of $y_{\tilde{i}}$ then satisfies

$$ \frac{dy_{\tilde{i}}}{dt} \geq \frac{dK_{\tilde{i}}}{\mu_{\tilde{i}}}y_{\tilde{i}}(t) - \alpha_i y_{\tilde{i}}(t). $$

By using $R_{0\tilde{i}} = dK_{\tilde{i}}/(\mu_{\tilde{i}}\alpha_{\tilde{i}}) > 1$, then $\lim_{t \to +\infty} y_{\tilde{i}}(t) = y_{\tilde{i}\infty} = +\infty$, which contradicts the hypothesis and clashes with the results derived earlier on the boundedness of solutions.

Clearly, if $R_{0\tilde{i}} > 1$, then $\mathcal{R}_0 > 1$. This is a consequence of Theorem 1: if $\mathcal{R}_0 < 1$, then the solution approach zero in every $y_i$ component, which is not the case for $i = \tilde{i}$. $\square$

Next, we show the weak persistence of each $y_i$, $i = 1, \dots, n$, in the following proposition.

**Proposition 4.** Suppose that $\mathcal{R}_0 > 1$. Then, there exists a constant $\varepsilon > 0$ such that, for any initial condition $(x_{i0}, \varphi_i) \in \mathbb{R}_+ \times C([-T, 0], \mathbb{R}_+)$, for $i = 1, \dots, n$, we have

$$ \limsup_{t \to +\infty} y_i(t) > \varepsilon, \quad \varphi_i(0) \neq 0. $$

*Proof.* We suppose by contradiction that $\limsup_{t \to +\infty} y_i(t) \le \varepsilon$, for $i = 1, \dots, n$ and for any small $\varepsilon > 0$. Then, there exists a sufficiently large $t_{1\varepsilon} > 0$ such that $y_i(t) \le \varepsilon$, for all $t \ge t_{1\varepsilon}$. Hence, we get for all $t \ge t_{1\varepsilon}$,

$$ \frac{dx_i}{dt} \geq K_i \beta(\varepsilon) - \mu_i x_i(t) - dx_i \left( \varepsilon + \varepsilon \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i} \right). $$

We denote $\liminf_{t \to +\infty} x_i(t) = x_{i\infty}$, for $i = 1, \dots, n$. Then, there exists a sequence $t_m \to +\infty$ as $m \to +\infty$, such that $x_i(t_m) \to x_{i\infty}$ and $x'_i(t_m) \to 0$ (see Lemma A.14 of [33]). This yields

$$ 0 \geq K_i \beta(\varepsilon) - \mu_i x_{i\infty} - dx_{i\infty} \left( \varepsilon + \varepsilon \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i} \right). $$

Then, we have

$$ x_{i\infty} \geq \frac{K_i \beta(\varepsilon)}{\mu_i x_{i\infty} + d \left( \varepsilon + \varepsilon \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i} \right)} =: x_{i\varepsilon}. $$
# 6 Persistence of solutions

In this section, we treat the long-term behavior of the system (6) when $\mathcal{R}_0 > 1$. We start with the following proposition:

**Proposition 3.** Consider a fixed $i = \tilde{i} \in \{1, \dots, n\}$. If $R_{0\tilde{i}} = dK_{\tilde{i}}/(\mu_{\tilde{i}}\alpha_{\tilde{i}}) > 1$, then $\mathcal{R}_0 > 1$ and there exists a constant $\varepsilon_{\tilde{i}} > 0$ such that

$$ \limsup_{t \to +\infty} y_{\tilde{i}}(t) > \varepsilon_{\tilde{i}}, \quad \text{with } \varphi_{\tilde{i}} \in C([-T, 0], \mathbb{R}_{+}), \quad \varphi_{\tilde{i}}(0) \neq 0. $$

*Proof.* By considering the second equation of System (6) for $i = \tilde{i}$, we have, for $t > 0$,

$$ \frac{dy_{\tilde{i}}}{dt} \geq dx_{\tilde{i}}(t)y_{\tilde{i}}(t) - \alpha_i y_{\tilde{i}}(t). $$

We have also, for $t > 0$,

$$ \frac{d(x_{\tilde{i}} + y_{\tilde{i}})}{dt} = K_{\tilde{i}}\beta(y_{\tilde{i}}(t - T_1)) - \mu_{\tilde{i}}x_{\tilde{i}}(t) - \alpha_{\tilde{i}}y_{\tilde{i}}(t). $$

We suppose by contradiction that $\limsup_{t \to +\infty} y_{\tilde{i}}(t) \le \varepsilon_{\tilde{i}}$, for any small $\varepsilon_{\tilde{i}} > 0$. By the boundedness of solutions, we consider $\liminf_{t \to +\infty} x_{\tilde{i}}(t) = x_{\tilde{i}\infty}$ and $\liminf_{t \to +\infty} y_{\tilde{i}}(t) = y_{\tilde{i}\infty} = 0$. Then, there exists a sequence $t_k \to +\infty$ as $k \to +\infty$, such that $x_{\tilde{i}}(t_k) \to x_{\tilde{i}\infty}$, $y_{\tilde{i}}(t_k) \to 0$, $y'_{\tilde{i}}(t_k) \to 0$ and $x'_{\tilde{i}}(t_k) \to 0$. This yields to

$$ 0 \geq K_{\tilde{i}}\beta(0) - \mu_{\tilde{i}}x_{\tilde{i}\infty} \Rightarrow x_{\tilde{i}\infty} \geq \frac{K_{\tilde{i}}\beta(0)}{\mu_{\tilde{i}}} = \frac{K_{\tilde{i}}}{\mu_{\tilde{i}}}. $$

For a very large time $t$, the ODE of $y_{\tilde{i}}$ then satisfies

$$ \frac{dy_{\tilde{i}}}{dt} \geq \frac{dK_{\tilde{i}}}{\mu_{\tilde{i}}}y_{\tilde{i}}(t) - \alpha_i y_{\tilde{i}}(t). $$

By using $R_{0\tilde{i}} = dK_{\tilde{i}}/(\mu_{\tilde{i}}\alpha_{\tilde{i}}) > 1$, then $\lim_{t \to +\infty} y_{\tilde{i}}(t) = y_{\tilde{i}\infty} = +\infty$, which contradicts the hypothesis and clashes with the results derived earlier on the boundedness of solutions.

Clearly, if $R_{0\tilde{i}} > 1$, then $\mathcal{R}_0 > 1$. This is a consequence of Theorem 1: if $\mathcal{R}_0 < 1$, then the solution approach zero in every $y_i$ component, which is not the case for $i = \tilde{i}$. $\square$

Next, we show the weak persistence of each $y_i$, $i = 1, \dots, n$, in the following proposition.

**Proposition 4.** Suppose that $\mathcal{R}_0 > 1$. Then, there exists a constant $\varepsilon > 0$ such that, for any initial condition $(x_{i0}, \varphi_i) \in \mathbb{R}_{+} \times C([-T, 0], \mathbb{R}_{+})$, for $i = 1, \dots, n$, we have

$$ \limsup_{t \to +\infty} y_i(t) > \varepsilon, \quad \varphi_i(0) \neq 0. $$

*Proof.* We suppose by contradiction that $\limsup_{t \to +\infty} y_i(t) \le \varepsilon$, for $i = 1, \dots, n$ and for any small $\varepsilon > 0$. Then, there exists a sufficiently large $t_{1\varepsilon} > 0$ such that $y_i(t) \le \varepsilon$, for all $t \ge t_{1\varepsilon}$. Hence, we get for all $t \ge t_{1\varepsilon}$,

$$ \frac{dx_i}{dt} \geq K_i \beta(\varepsilon) - \mu_i x_i(t) - dx_i \left( \varepsilon + \varepsilon \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i} \right). $$

We denote $\liminf_{t \to +\infty} x_i(t) = x_{i\infty}$, for $i = 1, \dots, n$. Then, there exists a sequence $t_m \to +\infty$ as $m \to +\infty$, such that $x_i(t_m) \to x_{i\infty}$ and $x'_i(t_m) \to 0$ (see Lemma A.14 of [33]). This yields

$$ 0 \geq K_i \beta(\varepsilon) - \mu_i x_{i\infty} - dx_{i\infty} \left( \varepsilon + \varepsilon \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i} \right). $$

Then, we have

$$ x_{i\infty} \geq \frac{K_i \beta(\varepsilon)}{\mu_i x_{i\infty} + d \left( \varepsilon + \varepsilon \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i} \right)} =: x_{i\varepsilon}. $$
Hence, for every small $\nu > 0$, there exists a sufficiently large $t_{2\nu} > 0$ such that, for $t \ge t_{2\nu}$,

$$x_i(t) \ge x_{i\epsilon} - \nu := x_{i\epsilon}^{\nu}.$$

Then, for a significant large time, we get

$$\frac{dy_i}{dt} \ge dx_{i\epsilon}^{\nu} \left( y_i(t) + \sum_{j \ne i} \kappa_{ji} \alpha_{j \to i} y_j(t) \right) - \alpha_i y_i(t).$$

As in the proof of Theorem 1, we consider the following system of ODEs:

$$\frac{dw_i}{dt} = dx_{i\epsilon}^{\nu} \left( w_i(t) + \sum_{j \ne i} \kappa_{ji} \alpha_{j \to i} w_j(t) \right) - \alpha_i w_i(t).$$

This system is linear in $w = (w_1, w_2, \dots, w_n)$, and can be rewritten as

$$\frac{dw}{dt} = (M_{22}(\epsilon, \nu) - V_{22})w,$$

where

$$(M_{22}(\epsilon, \nu))_{ij} = \begin{cases} dx_{i\epsilon}^{\nu} & \text{if } j = i, \\ dx_{i\epsilon}^{\nu} \kappa_{ji} \alpha_{j \to i} & \text{if } j \ne i, \end{cases} \quad V_{22} = \text{diag}(\alpha_i).$$

Using the hypothesis that $\mathcal{R}_0 = \rho(M_{22}(0,0)V_{22}^{-1}) > 1$, we can consider $\epsilon$ and $\nu$ sufficiently small such that

$$\mathcal{R}_0^{\epsilon, \nu} := \rho(M_{22}(\epsilon, \nu)V_{22}^{-1}) > 1.$$

We then use Lemma 1 to conclude that at least one eigenvalue of $(M_{22}(\epsilon, \nu) - V_{22})$ has positive real part. This leads to a contradiction with the assumption $\limsup_{t \to +\infty} y_i(t) \le \epsilon$, for all $i = 1, \dots, n$. As a consequence, there exists at least one $i = \tilde{i}$ such that

$$\limsup_{t \to +\infty} y_i(t) = y_{\tilde{i}} > \epsilon.$$

This is sufficient to conclude the weak persistence for each $i \in \{1, \dots, n\}$.

Suppose this is not true, and for some $i \in \{1, \dots, n\}$ and $i \ne \tilde{i}$ we have $\limsup_{t \to +\infty} y_i(t) = 0$. This means since solutions remain non-negative, that

$$\lim_{t \to +\infty} y_i(t) = 0.$$

By the boundedness of solutions and as a consequence of Barbalat's Lemma [34, 35], we obtain

$$\lim_{t \to +\infty} y_i'(t) = 0.$$

We can then choose a sequence $t_m \to +\infty$ as $m \to +\infty$, such that $y_{\tilde{i}}(t_m) \to y_{\tilde{i}\infty}$. Recall the equation of $y_i$, for $t > 0$,

$$\frac{dy_i}{dt} = dx_i(t) \left( y_i(t) + \sum_{j \ne i} \kappa_{ji} \alpha_{j \to i} y_j(t) \right) - \alpha_i y_i(t).$$

Therefore, by letting $t_m \to +\infty$ we obtain

$$0 \ge dx_{i\infty} \kappa_{i\tilde{i}} \alpha_{\tilde{i} \to i} y_{\tilde{i}\infty} > 0.$$

This contradiction completes the proof. $\square$

Using the boundedness of the solution (see Proposition 1) and the fact that $\beta$ is nonincreasing, we can show easily the following result.
Hence, for every small $\nu > 0$, there exists a sufficiently large $t_{2\nu} > 0$ such that, for $t \ge t_{2\nu}$,

$$x_i(t) \ge x_{i\epsilon} - \nu := x_{i\epsilon}^{\nu}.$$

Then, for a significant large time, we get

$$\frac{dy_i}{dt} \ge dx_{i\epsilon}^{\nu} \left( y_i(t) + \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i} y_j(t) \right) - \alpha_i y_i(t).$$

As in the proof of Theorem 1, we consider the following system of ODEs:

$$\frac{dw_i}{dt} = dx_{i\epsilon}^{\nu} \left( w_i(t) + \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i} w_j(t) \right) - \alpha_i w_i(t).$$

This system is linear in $w = (w_1, w_2, \dots, w_n)$, and can be rewritten as

$$\frac{dw}{dt} = (M_{22}(\epsilon, \nu) - V_{22})w,$$

where

$$(M_{22}(\epsilon, \nu))_{ij} = \begin{cases} dx_{i\epsilon}^{\nu} & \text{if } j = i, \\ dx_{i\epsilon}^{\nu} \kappa_{ji} \alpha_{j \to i} & \text{if } j \neq i, \end{cases} \quad V_{22} = \text{diag}(\alpha_i).$$

Using the hypothesis that $\mathcal{R}_0 = \rho(M_{22}(0, 0)V_{22}^{-1}) > 1$, we can consider $\epsilon$ and $\nu$ sufficiently small such that

$$\mathcal{R}_0^{\epsilon, \nu} := \rho(M_{22}(\epsilon, \nu)V_{22}^{-1}) > 1.$$

We then use Lemma 1 to conclude that at least one eigenvalue of $(M_{22}(\epsilon, \nu) - V_{22})$ has positive real part. This leads to a contradiction with the assumption $\limsup_{t \to +\infty} y_i(t) \le \epsilon$, for all $i = 1, \dots, n$. As a consequence, there exists at least one $i = \tilde{i}$ such that

$$\limsup_{t \to +\infty} y_i(t) = y_{\tilde{i}} > \epsilon.$$

This is sufficient to conclude the weak persistence for each $i \in \{1, \dots, n\}$.

Suppose this is not true, and for some $i \in \{1, \dots, n\}$ and $i \neq \tilde{i}$ we have $\limsup_{t \to +\infty} y_i(t) = 0$. This means since solutions remain non-negative, that

$$\lim_{t \to +\infty} y_i(t) = 0.$$

By the boundedness of solutions and as a consequence of Barbalat's Lemma [34, 35], we obtain

$$\lim_{t \to +\infty} y_i'(t) = 0.$$

We can then choose a sequence $t_m \to +\infty$ as $m \to +\infty$, such that $y_{\tilde{i}}(t_m) \to y_{\tilde{i}\infty}$. Recall the equation of $y_i$, for $t > 0$,

$$\frac{dy_i}{dt} = dx_i(t) \left( y_i(t) + \sum_{j \neq i} \kappa_{ji} \alpha_{j \to i} y_j(t) \right) - \alpha_i y_i(t).$$

Therefore, by letting $t_m \to +\infty$ we obtain

$$0 \ge dx_{i\infty} \kappa_{i\tilde{i}} \alpha_{\tilde{i} \to i} y_{\tilde{i}\infty} > 0.$$

This contradiction completes the proof. $\square$

Using the boundedness of the solution (see Proposition 1) and the fact that $\beta$ is nonincreasing, we can show easily the following result.
**Proposition 5.** There exists a constant $\tilde{\varepsilon} > 0$ such that, for any initial condition $(x_{i0}, \varphi_i) \in \mathbb{R}_+ \times C([-T, 0], \mathbb{R}_+)$, for $i = 1, ..., n$, we have

$$ \liminf_{t \to +\infty} x_i(t) > \tilde{\varepsilon}. $$

Now, we can establish the following result stating the strong uniform persistence of System (6) when $\mathcal{R}_0 > 1$.

**Theorem 3.** Suppose that $\mathcal{R}_0 > 1$. Then, there exists a constant $\varepsilon > 0$ such that, for any initial condition $(x_{i0}, \varphi_i) \in \mathbb{R}_+ \times C([-T, 0], \mathbb{R}_+)$, for $i = 1, ..., n$, we have

$$ \liminf_{t \to +\infty} y_i(t) > \varepsilon, \quad \varphi_i(0) \neq 0. $$

The proof can be adapted from the demonstration of Theorem 1 of [36] and it follows that the uniform weak persistence implies the uniform (strong) persistence (see also Theorem 7.3 of [37]).

# 7 Numerical simulations

In this section, we provide an extensive, but not exhaustive numerical exploration of the system (6), for various values of the parameters involved. Specifically, we simulate the model (6) for the four choices of networks depicted in Fig. 4: fully connected network with $n = 3$ neurons; line networks with $n = 5$ and $n = 9$ neurons; and ring network with $n = 5$ neurons. Our selection of initial conditions is influenced by a potential relevance from a biological standpoint. In the context of initiating an experiment, there is an inherent interest in commencing with all neurons exhibiting similar behavior. This entails the production of an identical amount of monomers and, through a protein misconfiguration process, the generation of oligomers in approximately the same quantities. This initial production may be set at zero in certain simulations, as illustrated in Figure 7, Figure 8, Figure 9, or at a non-zero value, as observed in Figure 5. Other choices of initial conditions can be explored, but as the analytical study has shown, the asymptotic behaviour of the solution will remain unchanged, so no further simulation is required.

We start by considering two fully connected, and fully homogeneous scenarios, which could correspond to the figure (a) in Fig. 4.

An illustration for the case $\mathcal{R}_0 < 1$ with $n = 3$ is given in Fig. 5. In this configuration, we observe an extinction of the disease, as analytically expected and proven in Section 4.

In Fig. 6, instead, we illustrate the case $\mathcal{R}_0 > 1$ with $n = 3$. These figure showcase the role of $\mu_i$, i.e. degradation rate of PrP<sup>C</sup> produced by the neuron i, particularly its effects on the basic reproduction number and the spread of the infection, in a non-trivial case.

We are also interested in the fully homogeneous case ($\kappa_{ij} := \kappa$ and $\alpha_{i \to j} = \alpha/(n-1)$) with a line connection between neurons as depicted in Fig. 4b, 4c and 4d. In this case, the system has the following form, for $t \ge 0$,

$$ \frac{dx_1}{dt} = K\beta(y_1(t-T)) - \mu x_1(t) - dx_1(t)y_1(t), $$

$$ \frac{dy_1}{dt} = dx_1(t)y_1(t) - \alpha y_1(t), $$

$$ \frac{dx_2}{dt} = K\beta(y_2(t-T)) - \mu x_2(t) - dx_2(t) \left(y_2(t) + \frac{\kappa\alpha}{n-1}y_1(t)\right), $$

$$ \frac{dy_2}{dt} = dx_2(t) \left(y_2(t) + \frac{\kappa\alpha}{n-1}y_1(t)\right) - \alpha y_2(t), $$

$$ \frac{dx_i}{dt} = K\beta(y_i(t-T)) - \mu x_i(t) - dx_i(t) \left(y_i(t) + \frac{\kappa\alpha}{n-1}y_{i-1}(t)\right), $$

$$ \frac{dy_i}{dt} = dx_i(t) \left(y_i(t) + \frac{\kappa\alpha}{n-1}y_{i-1}(t)\right) - \alpha y_i(t), $$
**Proposition 5.** There exists a constant $\tilde{\epsilon} > 0$ such that, for any initial condition $(x_{i0}, \varphi_i) \in \mathbb{R}_+ \times C([-T, 0], \mathbb{R}_+)$, for $i = 1, ..., n$, we have

$$ \liminf_{t \to +\infty} x_i(t) > \tilde{\epsilon}. $$

Now, we can establish the following result stating the strong uniform persistence of System (6) when $\mathcal{R}_0 > 1$.

**Theorem 3.** Suppose that $\mathcal{R}_0 > 1$. Then, there exists a constant $\epsilon > 0$ such that, for any initial condition $(x_{i0}, \varphi_i) \in \mathbb{R}_+ \times C([-T, 0], \mathbb{R}_+)$, for $i = 1, ..., n$, we have

$$ \liminf_{t \to +\infty} y_i(t) > \epsilon, \quad \varphi_i(0) \neq 0. $$

The proof can be adapted from the demonstration of Theorem 1 of [36] and it follows that the uniform weak persistence implies the uniform (strong) persistence (see also Theorem 7.3 of [37]).

# 7 Numerical simulations

In this section, we provide an extensive, but not exhaustive numerical exploration of the system (6), for various values of the parameters involved. Specifically, we simulate the model (6) for the four choices of networks depicted in Fig. 4: fully connected network with $n = 3$ neurons; line networks with $n = 5$ and $n = 9$ neurons; and ring network with $n = 5$ neurons. Our selection of initial conditions is influenced by a potential relevance from a biological standpoint. In the context of initiating an experiment, there is an inherent interest in commencing with all neurons exhibiting similar behavior. This entails the production of an identical amount of monomers and, through a protein misconfiguration process, the generation of oligomers in approximately the same quantities. This initial production may be set at zero in certain simulations, as illustrated in Figure 7, Figure 8, Figure 9, or at a non-zero value, as observed in Figure 5. Other choices of initial conditions can be explored, but as the analytical study has shown, the asymptotic behaviour of the solution will remain unchanged, so no further simulation is required.

We start by considering two fully connected, and fully homogeneous scenarios, which could correspond to the figure (a) in Fig. 4.

An illustration for the case $\mathcal{R}_0 < 1$ with $n = 3$ is given in Fig. 5. In this configuration, we observe an extinction of the disease, as analytically expected and proven in Section 4.

In Fig. 6, instead, we illustrate the case $\mathcal{R}_0 > 1$ with $n = 3$. These figure showcase the role of $\mu_i$, i.e. degradation rate of PrP<sup>C</sup> produced by the neuron i, particularly its effects on the basic reproduction number and the spread of the infection, in a non-trivial case.

We are also interested in the fully homogeneous case ($\kappa_{ij} := \kappa$ and $\alpha_{i \to j} = \alpha/(n-1)$) with a line connection between neurons as depicted in Fig. 4b, 4c and 4d. In this case, the system has the following form, for $t \ge 0$,

$$ \frac{dx_1}{dt} = K\beta(y_1(t-T)) - \mu x_1(t) - dx_1(t)y_1(t), $$

$$ \frac{dy_1}{dt} = dx_1(t)y_1(t) - \alpha y_1(t), $$

$$ \frac{dx_2}{dt} = K\beta(y_2(t-T)) - \mu x_2(t) - dx_2(t) \left(y_2(t) + \frac{\kappa\alpha}{n-1}y_1(t)\right), $$

$$ \frac{dy_2}{dt} = dx_2(t) \left(y_2(t) + \frac{\kappa\alpha}{n-1}y_1(t)\right) - \alpha y_2(t), $$

$\vdots$

$$ \frac{dx_i}{dt} = K\beta(y_i(t-T)) - \mu x_i(t) - dx_i(t) \left(y_i(t) + \frac{\kappa\alpha}{n-1}y_{i-1}(t)\right), $$

$$ \frac{dy_i}{dt} = dx_i(t) \left(y_i(t) + \frac{\kappa\alpha}{n-1}y_{i-1}(t)\right) - \alpha y_i(t), $$

$\vdots$
<figure><img src="image_7.png" /></figure>

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

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

<figure><img src="image_10.png" /><figcaption>Figure 4: The four networks we consider in our numerical simulations. Notice that only the network with $n = 3$ has double arrows on each edge, representing the fully connected network. The three remaining networks all have unidirectional edges.</figcaption></figure>

and the matrix (8) becomes

$$ (F)_{ij} = R_0 \begin{cases} 1 & \text{if } j = i, \\ \frac{\kappa\alpha}{n-1} & \text{if } j = i-1, i \ge 2, \\ 0 & \text{otherwise.} \end{cases} $$

Hence, we obtain

$$ \mathcal{R}_0 = \rho(F) = R_0 = \frac{dK}{\mu\alpha}. $$

Note that in the non-fully homogeneous case, and non-fully connected like in Fig. 4 (b), (c) or (d) (where parameters can be different from one neuron to another) we cannot apply this definition of $\mathcal{R}_0$, since some of the $\alpha_i$'s may end up to be equal to 0. We need then to go back to the more general theory.

In Fig. 7, a fully homogeneous line network, after some initial “wobbling”, the system approaches an Endemic Equilibrium. In particular, we notice that neurons further down the line (i.e., couples $(x_i, y_i)$ with larger $i$'s ($i = 4$, or 5) approach equilibrium with a higher value for the infected compartment $y_i$.

In Fig. 8-11, we examine the behavior of the neural network when the parameter $\kappa$ is varied. In fact, we obtain similar results, see Fig. 12, by varying the time delay $T$.

In Fig. 8 we “close” the line, forming a close ring with the five neurons. Without changing the other parameters, this additional link completely destabilizes the system. Each neuron approaches a limit cycle. Proving the existence of such a limit cycle analytically remains to prove and will be the object of future works. This result was proved for two neurons [20], by a Hopf bifurcation with respect to the parameters $\kappa$ and $T$. It can naturally be extended to a network of $n$ neurons. It is also interesting to notice that in the case of the close ring, all neurons oscillate at the same time, with the same period.
<figure><img src="image_7.png" /></figure>

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

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

<figure><img src="image_10.png" /><figcaption>Figure 4: The four networks we consider in our numerical simulations. Notice that only the network with $n = 3$ has double arrows on each edge, representing the fully connected network. The three remaining networks all have unidirectional edges.</figcaption></figure>

and the matrix (8) becomes

$$ (F)_{ij} = R_0 \begin{cases} 1 & \text{if } j = i, \\ \frac{\kappa\alpha}{n-1} & \text{if } j = i-1, i \ge 2, \\ 0 & \text{otherwise.} \end{cases} $$

Hence, we obtain

$$ \mathcal{R}_0 = \rho(F) = R_0 = \frac{dK}{\mu\alpha}. $$

Note that in the non-fully homogeneous case, and non-fully connected like in Fig. 4 (b), (c) or (d) (where parameters can be different from one neuron to another) we cannot apply this definition of $\mathcal{R}_0$, since some of the $\alpha_i$'s may end up to be equal to 0. We need then to go back to the more general theory.

In Fig. 7, a fully homogeneous line network, after some initial “wobbling”, the system approaches an Endemic Equilibrium. In particular, we notice that neurons further down the line (i.e., couples $(x_i, y_i)$ with larger $i$'s ($i = 4$, or 5) approach equilibrium with a higher value for the infected compartment $y_i$.

In Fig. 8-11, we examine the behavior of the neural network when the parameter $\kappa$ is varied. In fact, we obtain similar results, see Fig. 12, by varying the time delay $T$.

In Fig. 8 we “close” the line, forming a close ring with the five neurons. Without changing the other parameters, this additional link completely destabilizes the system. Each neuron approaches a limit cycle. Proving the existence of such a limit cycle analytically remains to prove and will be the object of future works. This result was proved for two neurons [20], by a Hopf bifurcation with respect to the parameters $\kappa$ and $T$. It can naturally be extended to a network of $n$ neurons. It is also interesting to notice that in the case of the close ring, all neurons oscillate at the same time, with the same period.
<figure><img src="image_11.png" /><figcaption>Figure 5: the fully connected network of $n = 3$ neurons (recall Fig. 4a), showing the stability of disease free equilibrium when $R_0 = \rho(F) = 0.8194 < 1$. In this case $R_{0i} = 0.6944 < 1$ for $i = 1, 2, 3$. The parameters are: $\alpha_{i \to j} = 0.9$, $\kappa_{ij} = 0.1$, $p = 5$, $d = 0.015$, $y_c = 60$, $T = 0.17$, $K_i = 1500$ and $\mu_i = 18$.</figcaption></figure>

In Fig. 9, we revisit the $n = 5$ line structure, illustrating a case in which neurons 4 and 5 exhibit instability in their asymptotic behaviour, whereas neuron 1 converges to an equilibrium. This is possible because neuron 1 only spreads the infection, and does not receive feedback from the remaining neurons in the network. Moreover, we provide more detail on the parameter $\kappa$, which incorporates crucial information, namely the interconnectivity between neurons. We provide the bifurcation values of $\kappa$, assuming all the other parameters to be fixed, for which each neuron of the system destabilizes. From a mathematical point of view, the oscillating neuron acts as a force on the next neuron, which in turn oscillates, and the oscillations thus propagate to all the following neurons.

In Fig. 10, we explored the case of Fig. 9 a little further in the following sense: beginning our simulations from the stable case, we increased $\kappa$, one of the key parameters for the Hopf bifurcations and managed to compute the exact value of $\kappa$ at which the first neuron of the line would oscillate. For instance, at $\kappa = 0.07$ only the last (the fifth) one is destabilized, while for $\kappa = 0.077$ only the first one is stable, while all the others oscillate. We observe that this process is non-linear. Predicting the number of neurons destabilized with respect to $\kappa$ analytically is also an open problem that we keep for future work.

In Fig. 11, we showcase how a continuous variation of the parameter $\kappa$ impacts the asymptotic value of each $x_i$ and $y_i$, again for the $n = 5$ line network. Increasing $\kappa$ causes more neurons to destabilize; for each of them, we plot the maximum and minimum values assumed by each variable as they asymptotically approach a limit cycle. Similarly, in Fig. 12 we showcase what the influence of the delay $T$ is on these oscillations.

Finally, in Fig. 13, we explore the effect of increasing the number of neurons in the line network in the case of our choice of $\kappa_{ij} = \kappa = 0.125$ set up in Section 3.1. Using similar parameter values as in Fig. 9 (except for $\kappa_{ij}$ and $\alpha_{i \to j}$, to keep them biologically feasible). We manage to show that adding more neurons (going from 5 to 9 neurons here) can damp the oscillations, and lead the system of neurons to stabilize again. This process has to be explored biologically to be confirmed experimentally.

We remark that our results in Section 3.1 only concern $R_0$, from which however we can only predict extinction or permanence of the disease, and nothing on the asymptotic stability of orbits. Hence, further analytical results in this direction are a promising research outlook.
<figure><img src="image_11.png" /><figcaption>Figure 5: the fully connected network of $n = 3$ neurons (recall Fig. 4a), showing the stability of disease free equilibrium when $\mathcal{R}_0 = \rho(F) = 0.8194 < 1$. In this case $R_{0i} = 0.6944 < 1$ for $i = 1, 2, 3$. The parameters are: $\alpha_{i \to j} = 0.9$, $\kappa_{ij} = 0.1$, $p = 5$, $d = 0.015$, $y_c = 60$, $T = 0.17$, $K_i = 1500$ and $\mu_i = 18$.</figcaption></figure>

In Fig. 9, we revisit the $n = 5$ line structure, illustrating a case in which neurons 4 and 5 exhibit instability in their asymptotic behaviour, whereas neuron 1 converges to an equilibrium. This is possible because neuron 1 only spreads the infection, and does not receive feedback from the remaining neurons in the network. Moreover, we provide more detail on the parameter $\kappa$, which incorporates crucial information, namely the interconnectivity between neurons. We provide the bifurcation values of $\kappa$, assuming all the other parameters to be fixed, for which each neuron of the system destabilizes. From a mathematical point of view, the oscillating neuron acts as a force on the next neuron, which in turn oscillates, and the oscillations thus propagate to all the following neurons.

In Fig. 10, we explored the case of Fig. 9 a little further in the following sense: beginning our simulations from the stable case, we increased $\kappa$, one of the key parameters for the Hopf bifurcations and managed to compute the exact value of $\kappa$ at which the first neuron of the line would oscillate. For instance, at $\kappa = 0.07$ only the last (the fifth) one is destabilized, while for $\kappa = 0.077$ only the first one is stable, while all the others oscillate. We observe that this process is non-linear. Predicting the number of neurons destabilized with respect to $\kappa$ analytically is also an open problem that we keep for future work.

In Fig. 11, we showcase how a continuous variation of the parameter $\kappa$ impacts the asymptotic value of each $x_i$ and $y_i$, again for the $n = 5$ line network. Increasing $\kappa$ causes more neurons to destabilize; for each of them, we plot the maximum and minimum values assumed by each variable as they asymptotically approach a limit cycle. Similarly, in Fig. 12 we showcase what the influence of the delay $T$ is on these oscillations.

Finally, in Fig. 13, we explore the effect of increasing the number of neurons in the line network in the case of our choice of $\kappa_{ij} = \kappa = 0.125$ set up in Section 3.1. Using similar parameter values as in Fig. 9 (except for $\kappa_{ij}$ and $\alpha_{i \to j}$, to keep them biologically feasible). We manage to show that adding more neurons (going from 5 to 9 neurons here) can damp the oscillations, and lead the system of neurons to stabilize again. This process has to be explored biologically to be confirmed experimentally.

We remark that our results in Section 3.1 only concern $\mathcal{R}_0$, from which however we can only predict extinction or permanence of the disease, and nothing on the asymptotic stability of orbits. Hence, further analytical results in this direction are a promising research outlook.
<figure><img src="image_12.png" /><figcaption>Figure 6: the fully connected network of $n = 3$ neurons (depicted by Fig. 4a), showing the stability of endemic equilibrium when $R_0 = \rho(F) = 1.1346 > 1$. In this case $R_{0i} = 0.9615 < 1$ for $i = 1, 2, 3$. The parameters are: $\alpha_{i \to j} = 0.9$, $\kappa_{ij} = 0.1$, $p = 5$, $d = 0.015$, $y_c = 60$, $T = 0.17$, $K_i = 1500$ and $\mu_i = 13$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 7: line network of $n = 5$ neurons (depicted in Fig. 4b). This case shows that cutting the connection showed stabilization. The parameters are: $\alpha_{i \to j} = 2.5$, $\kappa_{ij} = 0.17$ (when considered), $p = 10$, $d = 0.15$, $y_c = 50$, $K_i = 1500$, $\mu_i = 20$ and $T = 0.15$.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 6: the fully connected network of $n = 3$ neurons (depicted by Fig. 4a), showing the stability of endemic equilibrium when $R_0 = \rho(F) = 1.1346 > 1$. In this case $R_{0i} = 0.9615 < 1$ for $i = 1, 2, 3$. The parameters are: $\alpha_{i \to j} = 0.9$, $\kappa_{ij} = 0.1$, $p = 5$, $d = 0.015$, $y_c = 60$, $T = 0.17$, $K_i = 1500$ and $\mu_i = 13$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 7: line network of $n = 5$ neurons (depicted in Fig. 4b). This case shows that cutting the connection showed stabilization. The parameters are: $\alpha_{i \to j} = 2.5$, $\kappa_{ij} = 0.17$ (when considered), $p = 10$, $d = 0.15$, $y_c = 50$, $K_i = 1500$, $\mu_i = 20$ and $T = 0.15$.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 8: circle (ring) network of $n = 5$ neurons (depicted in Fig. 4c). This case shows that linking the connection in Fig. 7 showed destabilization of the system. The parameters are: $\alpha_{i \to j} = 2.5$, $\kappa_{ij} = 0.17$ (when considered), $p = 10$, $d = 0.15$, $y_c = 50$, $K_i = 1500$, $\mu_i = 20$ and $T = 0.15$. Similar results can be obtained by varying the time delay $T$ and fixing all other parameters.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 9: line network of $n = 5$ neurons (depicted in Fig. 4b). For this figure showing the oscillation of only some neurons ($n = 4, 5$, the last ones), we took $\kappa = 0.071$. Parameters are: $p = 10$, $d = 0.15$, $y_c = 60$, $K_i = 1800$, $\mu_i = 50$, $\alpha_{i \to j} = 0.9$ ($\alpha_i = 3.6$) and $T = 0.15$. Similar results can be produced by changing the delay $T$.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 8: circle (ring) network of $n = 5$ neurons (depicted in Fig. 4c). This case shows that linking the connection in Fig. 7 showed destabilization of the system. The parameters are: $\alpha_{i \to j} = 2.5$, $\kappa_{ij} = 0.17$ (when considered), $p = 10$, $d = 0.15$, $y_c = 50$, $K_i = 1500$, $\mu_i = 20$ and $T = 0.15$. Similar results can be obtained by varying the time delay $T$ and fixing all other parameters.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 9: line network of $n = 5$ neurons (depicted in Fig. 4b). For this figure showing the oscillation of only some neurons ($n = 4, 5$, the last ones), we took $\kappa = 0.071$. Parameters are: $p = 10$, $d = 0.15$, $y_c = 60$, $K_i = 1800$, $\mu_i = 50$, $\alpha_{i \to j} = 0.9$ ($\alpha_i = 3.6$) and $T = 0.15$. Similar results can be produced by changing the delay $T$.</figcaption></figure>
<figure><img src="image_16.png" /></figure>

Figure 10: line network of $n = 5$ neurons (depicted in Fig. 4b). Our starting point was the case where the system is stable and we have increased the bifurcation parameter $\kappa$ and noted each value for which a neuron is destabilized. The value of bifurcations are 0.07, 0.0709, 0.728 and 0.768. The neurons are destabilized one by one from the last and going up to the second. The last neurons always look destabilized almost at the same time (approximately the same bifurcation values). Parameters are: $p = 10$, $d = 0.15$, $y_c = 60$, $K_i = 1800$, $\mu_i = 50$, $\alpha_{i \to j} = 0.9$ ($\alpha_i = 3.6$) and $T = 0.15$. A curve similar to this one can be obtained as a function of the time delay $T$.

# 8 Conclusions and outlook

In this paper, we presented a model for the delayed spread of prion in a network of $n$ neurons, building on the 1 neuron model proposed in [20]. We studied its analytical properties and provided extensive numerical simulations to illustrate various scenarios.

Due to the high dimension of the system, and of the analytical complexity of systems of DDEs, many questions remain unanswered. How can we overcome the requirement that $\alpha_i > 0$ in the definition of $R_0$? The system has a clear biological interpretation even when this condition is not satisfied; hence, we would like to find a threshold quantity in such a scenario. Moreover, does the Endemic Equilibrium exist for all systems with $R_0 > 1$? Is it unique? If yes, when is it stable? In our numerical exploration, we found both convergences to Endemic Equilibrium and sustained oscillations, indicating a possible stable limit cycle arising in the system. It would be of interest to understand which relations between the parameters of the system lead to the former or the latter.

Finally, we should point out here also that spatial structure has not been taken into account. Indeed, considering diffusion may appear challenging for the following reason: in the case of Alzheimer's disease, oligomers diffuse randomly in the brain tissue since the $A\beta$ monomers are no longer anchored to the cell membrane. On the contrary, for the prion disease, pathological $PrP^{Sc}$ proteins spread following the axon (and thus the cell membrane) where the source of non-pathological $PrP^C$ proteins are attached (thanks to a GPI anchor). Thus diffusion cannot be represented in the same way depending on the neurodegenerative disease studied. Furthermore, $PrP^{Sc}$ might not be produced immediately after the contact with previous neurons. A time lag may be needed, and thus this could involve some de-synchronisation and perhaps some chaotic behaviour. We leave these fundamental and other interesting questions as the possible outlook for future work.

**Acknowledgements.** Mattia Sensi was supported by the Italian Ministry for University and Research (MUR) through the PRIN 2020 project "Integrated Mathematical Approaches to Socio-Epidemiological Dynamics" (No. 2020JLWP23, CUP: E15F21005420006).

Abdennasser Chekroun was supported by the grant PRFU: C00L03UN29012022002, from DGRSDT of Algeria.
Laurent Pujo-Menjouet was supported by ANR grant PrionDiff ANR-21-CE15-0011-02.
<figure><img src="image_16.png" /></figure>

Figure 10: line network of $n = 5$ neurons (depicted in Fig. 4b). Our starting point was the case where the system is stable and we have increased the bifurcation parameter $\kappa$ and noted each value for which a neuron is destabilized. The value of bifurcations are 0.07, 0.0709, 0.728 and 0.768. The neurons are destabilized one by one from the last and going up to the second. The last neurons always look destabilized almost at the same time (approximately the same bifurcation values). Parameters are: $p = 10$, $d = 0.15$, $y_c = 60$, $K_i = 1800$, $\mu_i = 50$, $\alpha_{i \to j} = 0.9$ ($\alpha_i = 3.6$) and $T = 0.15$. A curve similar to this one can be obtained as a function of the time delay $T$.

# 8 Conclusions and outlook

In this paper, we presented a model for the delayed spread of prion in a network of $n$ neurons, building on the 1 neuron model proposed in [20]. We studied its analytical properties and provided extensive numerical simulations to illustrate various scenarios.

Due to the high dimension of the system, and of the analytical complexity of systems of DDEs, many questions remain unanswered. How can we overcome the requirement that $\alpha_i > 0$ in the definition of $R_0$? The system has a clear biological interpretation even when this condition is not satisfied; hence, we would like to find a threshold quantity in such a scenario. Moreover, does the Endemic Equilibrium exist for all systems with $R_0 > 1$? Is it unique? If yes, when is it stable? In our numerical exploration, we found both convergences to Endemic Equilibrium and sustained oscillations, indicating a possible stable limit cycle arising in the system. It would be of interest to understand which relations between the parameters of the system lead to the former or the latter.

Finally, we should point out here also that spatial structure has not been taken into account. Indeed, considering diffusion may appear challenging for the following reason: in the case of Alzheimer's disease, oligomers diffuse randomly in the brain tissue since the $A\beta$ monomers are no longer anchored to the cell membrane. On the contrary, for the prion disease, pathological $PrP^{Sc}$ proteins spread following the axon (and thus the cell membrane) where the source of non-pathological $PrP^C$ proteins are attached (thanks to a GPI anchor). Thus diffusion cannot be represented in the same way depending on the neurodegenerative disease studied. Furthermore, $PrP^{Sc}$ might not be produced immediately after the contact with previous neurons. A time lag may be needed, and thus this could involve some de-synchronisation and perhaps some chaotic behaviour. We leave these fundamental and other interesting questions as the possible outlook for future work.

**Acknowledgements.** Mattia Sensi was supported by the Italian Ministry for University and Research (MUR) through the PRIN 2020 project "Integrated Mathematical Approaches to Socio-Epidemiological Dynamics" (No. 2020JLWP23, CUP: E15F21005420006).

Abdennasser Chekroun was supported by the grant PRFU: C00L03UN29012022002, from DGRSDT of Algeria.
Laurent Pujo-Menjouet was supported by ANR grant PrionDiff ANR-21-CE15-0011-02.
<figure><img src="image_17.png" /><figcaption>Figure 11: line network of $n = 5$ neurons (recall Fig. 4b). The amplitude of the oscillations as a function of $\kappa \in [0.055, 0.1]$. The amplitude of the oscillations becomes almost the same for higher values of $\kappa$. Parameters are: $p = 10$, $d = 0.15$, $y_c = 60$, $K_i = 1800$, $\mu_i = 50$, $\alpha_{i \to j} = 0.9$ ($\alpha_i = 3.6$) and $T = 0.15$. The red, blue, green and yellow curves are associated respectively with $(x_2, y_2)$, $(x_3, y_3)$, $(x_4, y_4)$ and $(x_5, y_5)$. The same behavior can be obtained as a function of time delay $T$, see Figure 12.</figcaption></figure>

References

[1] S. B. Prusiner, M. R. Scott, S. J. DeArmond, and F. E. Cohen. Prion Protein Biology. *Cell*, 93(3):337-348, May 1998.

[2] X. Roucou and A. C. LeBlanc. Cellular prion protein neuroprotective function: implications in prion diseases. *Journal of Molecular Medicine*, 83(1):3-11, January 2005.

[3] Stanley B. Prusiner. Prions. *Proceedings of the National Academy of Sciences*, 95(23):13363-13383, 1998.

[4] J. C. Genereux and R. L. Wiseman. Regulating extracellular proteostasis capacity through the unfolded protein response. *Prion*, 9(1):10-21, January 2015.

[5] C. Hetz and B. Mollereau. Disturbance of endoplasmic reticulum proteostasis in neurodegenerative diseases. *Nature Reviews Neuro-science*, 15(4):233-249, April 2014.

[6] C. Hetz and S. Saxena. ER stress and the unfolded protein response in neurodegeneration. *Nature Reviews Neurology*, 13(8):477-491, August 2017.

[7] C. Hetz, K. Zhang, and R. J. Kaufman. Mechanisms, regulation and functions of the unfolded protein response. *Nature Reviews Molecular Cell Biology*, 21(8):421-438, August 2020.

[8] H. L. Smith and G. R. Mallucci. The unfolded protein response: mechanisms and therapy of neurodegeneration. *Brain*, 139(8):2113-2121, August 2016.

[9] C. Hetz, M. Russelakis-Carneiro, K. Maundrell, J. Castilla, and C. Soto. Caspase-12 and endoplasmic reticulum stress mediate neurotoxicity of pathological prion protein. *The EMBO Journal*, 22(20):5435-5445, October 2003.

[10] M. Tanaka, T. Yamasaki, R. Hasebe, A. Suzuki, and M. Horiuchi. Enhanced phosphorylation of PERK in primary cultured neurons as an autonomous neuronal response to prion infection. *PLOS ONE*, 15(6):e0234147, June 2020.

[11] B. Schneider, A. Baudry, M. Pietri, A. Alleaume-Butaux, C. Bizingre, P. Nioche, O. Kellermann, and J.-M. Launay. The Cellular Prion Protein—ROCK Connection: Contribution to Neuronal Homeostasis and Neurodegenerative Diseases. *Frontiers in Cellular Neuroscience*, 15:101, 2021.
<figure><img src="image_17.png" /><figcaption>Figure 11: line network of $n = 5$ neurons (recall Fig. 4b). The amplitude of the oscillations as a function of $\kappa \in [0.055, 0.1]$. The amplitude of the oscillations becomes almost the same for higher values of $\kappa$. Parameters are: $p = 10$, $d = 0.15$, $y_c = 60$, $K_i = 1800$, $\mu_i = 50$, $\alpha_{i \to j} = 0.9$ ($\alpha_i = 3.6$) and $T = 0.15$. The red, blue, green and yellow curves are associated respectively with $(x_2, y_2)$, $(x_3, y_3)$, $(x_4, y_4)$ and $(x_5, y_5)$. The same behavior can be obtained as a function of time delay $T$, see Figure 12.</figcaption></figure>

References

[1] S. B. Prusiner, M. R. Scott, S. J. DeArmond, and F. E. Cohen. Prion Protein Biology. *Cell*, 93(3):337-348, May 1998.

[2] X. Roucou and A. C. LeBlanc. Cellular prion protein neuroprotective function: implications in prion diseases. *Journal of Molecular Medicine*, 83(1):3-11, January 2005.

[3] Stanley B. Prusiner. Prions. *Proceedings of the National Academy of Sciences*, 95(23):13363-13383, 1998.

[4] J. C. Genereux and R. L. Wiseman. Regulating extracellular proteostasis capacity through the unfolded protein response. *Prion*, 9(1):10-21, January 2015.

[5] C. Hetz and B. Mollereau. Disturbance of endoplasmic reticulum proteostasis in neurodegenerative diseases. *Nature Reviews Neuro-science*, 15(4):233-249, April 2014.

[6] C. Hetz and S. Saxena. ER stress and the unfolded protein response in neurodegeneration. *Nature Reviews Neurology*, 13(8):477-491, August 2017.

[7] C. Hetz, K. Zhang, and R. J. Kaufman. Mechanisms, regulation and functions of the unfolded protein response. *Nature Reviews Molecular Cell Biology*, 21(8):421-438, August 2020.

[8] H. L. Smith and G. R. Mallucci. The unfolded protein response: mechanisms and therapy of neurodegeneration. *Brain*, 139(8):2113-2121, August 2016.

[9] C. Hetz, M. Russelakis-Carneiro, K. Maundrell, J. Castilla, and C. Soto. Caspase-12 and endoplasmic reticulum stress mediate neurotoxicity of pathological prion protein. *The EMBO Journal*, 22(20):5435-5445, October 2003.

[10] M. Tanaka, T. Yamasaki, R. Hasebe, A. Suzuki, and M. Horiuchi. Enhanced phosphorylation of PERK in primary cultured neurons as an autonomous neuronal response to prion infection. *PLOS ONE*, 15(6):e0234147, June 2020.

[11] B. Schneider, A. Baudry, M. Pietri, A. Alleaume-Butaux, C. Bizingre, P. Nioche, O. Kellermann, and J.-M. Launay. The Cellular Prion Protein—ROCK Connection: Contribution to Neuronal Homeostasis and Neurodegenerative Diseases. *Frontiers in Cellular Neuroscience*, 15:101, 2021.
<figure><img src="image_18.png" /><figcaption>Figure 12: line network of $n = 5$ neurons (recall Fig. 4b). The amplitude of the oscillations as a function of $T \in [0.15, 0.4]$. All the cycles except for the blue one are almost overlapping for every value of $T$. Parameters are: $p = 10$, $d = 0.15$, $y_c = 60$, $K_i = 1800$, $\mu_i = 50$, $\alpha_{i \to j} = 0.9$ ($\alpha_i = 3.6$) and $\kappa = 0.07$. The blue, red, green, purple and yellow curves are associated respectively with $(x_1, y_1)$, $(x_2, y_2)$, $(x_3, y_3)$, $(x_4, y_4)$ and $(x_5, y_5)$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 13: line network of $n = 9$ neurons (recall Fig. 4d). We only plot the time series of $(x_i, y_i)$ with $i = 1, 2, 5, 8, 9$, respectively the first two, the central one and the last two neurons in the line network. By increasing the number $n$ of neurons from 5 to 9 with the same parameters as in Fig. 10 and 11 and with $\kappa = 0.125$, the system becomes stable. Recall parameters: $p = 10$, $d = 0.15$, $y_c = 60$, $K_i = 1800$, $\mu_i = 50$, $\alpha_{i \to j} = 0.45$ ($\alpha_i = 3.6$) and $T = 0.15$.</figcaption></figure>

<footnote>[12] H. L. Smith, O. J. Freeman, A. J. Butcher, S. Holmqvist, I. Humoud, T. Schätzl, D. T. Hughes, N. C. Verity, D. P. Swinden, J. Hayes, L. de Weerd, D. H. Rowitch, R. J. M. Franklin, and G. R. Mallucci. Astrocyte Unfolded Protein Response Induces a Specific Reactivity State that Causes Non-Cell-Autonomous Neuronal Degeneration. *Neuron*, 105(5):855-866.e5, March 2020.</footnote>

<footnote>[13] J. A. Moreno, H. Radford, D. Peretti, J. R. Steinert, N. Verity, M. G. Martin, M. Halliday, J. Morgan, D. Dinsdale, C. A. Ortori, D. A. Barrett, P. Tsaytler, A. Bertolotti, A. E. Willis, M. Bushell, and G. R. Mallucci. Sustained translational repression by eIF2α-P mediates prion neurodegeneration. *Nature*, 485(7399):507-511, May 2012.</footnote>
<figure><img src="image_18.png" /><figcaption>Figure 12: line network of $n = 5$ neurons (recall Fig. 4b). The amplitude of the oscillations as a function of $T \in [0.15, 0.4]$. All the cycles except for the blue one are almost overlapping for every value of $T$. Parameters are: $p = 10$, $d = 0.15$, $y_c = 60$, $K_i = 1800$, $\mu_i = 50$, $\alpha_{i \to j} = 0.9$ ($\alpha_i = 3.6$) and $\kappa = 0.07$. The blue, red, green, purple and yellow curves are associated respectively with $(x_1, y_1)$, $(x_2, y_2)$, $(x_3, y_3)$, $(x_4, y_4)$ and $(x_5, y_5)$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 13: line network of $n = 9$ neurons (recall Fig. 4d). We only plot the time series of $(x_i, y_i)$ with $i = 1, 2, 5, 8, 9$, respectively the first two, the central one and the last two neurons in the line network. By increasing the number $n$ of neurons from 5 to 9 with the same parameters as in Fig. 10 and 11 and with $\kappa = 0.125$, the system becomes stable. Recall parameters: $p = 10$, $d = 0.15$, $y_c = 60$, $K_i = 1800$, $\mu_i = 50$, $\alpha_{i \to j} = 0.45$ ($\alpha_i = 3.6$) and $T = 0.15$.</figcaption></figure>

<footnote>[12] H. L. Smith, O. J. Freeman, A. J. Butcher, S. Holmqvist, I. Humoud, T. Schätzl, D. T. Hughes, N. C. Verity, D. P. Swinden, J. Hayes, L. de Weerd, D. H. Rowitch, R. J. M. Franklin, and G. R. Mallucci. Astrocyte Unfolded Protein Response Induces a Specific Reactivity State that Causes Non-Cell-Autonomous Neuronal Degeneration. *Neuron*, 105(5):855-866.e5, March 2020.</footnote>

<footnote>[13] J. A. Moreno, H. Radford, D. Peretti, J. R. Steinert, N. Verity, M. G. Martin, M. Halliday, J. Morgan, D. Dinsdale, C. A. Ortori, D. A. Barrett, P. Tsaytler, A. Bertolotti, A. E. Willis, M. Bushell, and G. R. Mallucci. Sustained translational repression by eIF2α-P mediates prion neurodegeneration. *Nature*, 485(7399):507-511, May 2012.</footnote>
[14] M. Torres, K. Castillo, R. Armisén, A. Stutzin, C. Soto, and C. Hetz. Prion Protein Misfolding Affects Calcium Homeostasis and Sensitizes Cells to Endoplasmic Reticulum Stress. *PLOS ONE*, 5(12):e15658, December 2011.

[15] K. L. Cook, P. A. G. Clarke, J. Parmar, R. Hu, J. L. Schwartz-Roberts, M. Abu-Asab, A. Wärri, W. T. Baumann, and R. Clarke. Knockdown of estrogen receptor-α induces autophagy and inhibits antiestrogen-mediated unfolded protein response activation, promoting ROS-induced breast cancer cell death. *The FASEB Journal*, 28(9):3891-3905, September 2014.

[16] S. Schnell. A Model of the Unfolded Protein Response: Pancreatic β-Cell as a Case Study. *Cellular Physiology and Biochemistry*, 23(4-6):233-244, 2009.

[17] A. Trusina, F. R. Papa, and C. Tang. Rationalizing translation attenuation in the network architecture of the unfolded protein response. *Proc Natl Acad Sci USA*, 105(51):20280, December 2008.

[18] A. Trusina and C. Tang. The unfolded protein response and translation attenuation: a modelling approach. *Diabetes, Obesity and Metabolism*, 12(s2):27-31, October 2010.

[19] R. Luke Wiseman, Evan T. Powers, Joel N. Buxbaum, Jeffery W. Kelly, and William E. Balch. An Adaptable Standard for Protein Export from the Endoplasmic Reticulum. *Cell*, 131(4):809-821, November 2007.

[20] M. Adimy, L. Babin, and L. Pujo-Menjouet. Neuron Scale Modeling of Prion Production with the Unfolded Protein Response. *SIAM Journal on Applied Dynamical Systems*, 21(4):2487-2517, 2022.

[21] O. Diekmann, J. A. P. Heesterbeek, and J. A. J. Metz. On the definition and the computation of the basic reproduction ratio R 0 in models for infectious diseases in heterogeneous populations. *Journal of mathematical biology*, 28(4):365-382, 1990.

[22] P. Van den Driessche and J. Watmough. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. *Mathematical Biosciences*, 180(1):29-48, 2002.

[23] O. Diekmann, J. A. P. Heesterbeek, and M. G. Roberts. The construction of next-generation matrices for compartmental epidemic models. *Journal of the royal society interface*, 7(47):873-885, 2010.

[24] J. K. Hale and S. M. Verduyn Lunel. *Introduction to Functional Differential Equations*. Springer, 1993.

[25] Y. Kuang. *Delay Differential Equations: With Applications in Population Dynamics*. Academic Press, 1993.

[26] H. Smith. *An Introduction to Delay Differential Equations with Applications to the Life Sciences*. Texts in Applied Mathematics. Springer, 2011.

[27] Oscar Rojo and Héctor Rojo. Some results on symmetric circulant matrices and on symmetric centrosymmetric matrices. *Linear Algebra and its Applications*, 392:211-233, 2004.

[28] Kaveh Mousavand Dan Kucerovsky and Aydin Sarraf. On some properties of toeplitz matrices. *Cogent Mathematics*, 3(1):1154705, 2016.

[29] S. Ottaviano, M. Sensi, and S. Sottile. Global stability of multi-group SAIRS epidemic models. *Mathematical Methods in the Applied Sciences*, 2023.

[30] P. Van den Driessche and J. Watmough. Further notes on the basic reproduction number. *Mathematical epidemiology*, pages 159-178, 2008.

[31] W. Kulpa. The Poincaré-Miranda theorem. *The American Mathematical Monthly*, 104(6):545-550, 1997.

[32] J. Mawhin. Simple proofs of the Hadamard and Poincaré-Miranda theorems using the Brouwer fixed point theorem. *The American Mathematical Monthly*, 126(3):260-263, 2019.

[33] H. L Smith and H. R Thieme. *Dynamical Systems and Population Persistence*. Americal Mathematical Society, 2011.

[34] H. K. Khalil. Nonlinear systems third edition. *Patience Hall*, 115, 2002.

[35] Z. Sun. A gathering of Barbalat's lemmas and their (unsung) cousins. *arXiv preprint arXiv:2301.00466*, 2023.

[36] H. I. Freedman and P. Moson. Persistence definitions and their connections. *Proceedings of the american mathematical society*, 109(4):1025-1033, 1990.

[37] M. Adimy, A. Chekroun, and C. P. Ferreira. Global dynamics of a differential-difference system: a case of Kermack-McKendrick SIR model with age-structured protection phase. *Mathematical Biosciences and Engineering*, 17(2):1329-1354, 2020.
[14] M. Torres, K. Castillo, R. Armisén, A. Stutzin, C. Soto, and C. Hetz. Prion Protein Misfolding Affects Calcium Homeostasis and Sensitizes Cells to Endoplasmic Reticulum Stress. *PLOS ONE*, 5(12):e15658, December 2011.

[15] K. L. Cook, P. A. G. Clarke, J. Parmar, R. Hu, J. L. Schwartz-Roberts, M. Abu-Asab, A. Wärri, W. T. Baumann, and R. Clarke. Knockdown of estrogen receptor-α induces autophagy and inhibits antiestrogen-mediated unfolded protein response activation, promoting ROS-induced breast cancer cell death. *The FASEB Journal*, 28(9):3891-3905, September 2014.

[16] S. Schnell. A Model of the Unfolded Protein Response: Pancreatic β-Cell as a Case Study. *Cellular Physiology and Biochemistry*, 23(4-6):233-244, 2009.

[17] A. Trusina, F. R. Papa, and C. Tang. Rationalizing translation attenuation in the network architecture of the unfolded protein response. *Proc Natl Acad Sci USA*, 105(51):20280, December 2008.

[18] A. Trusina and C. Tang. The unfolded protein response and translation attenuation: a modelling approach. *Diabetes, Obesity and Metabolism*, 12(s2):27-31, October 2010.

[19] R. Luke Wiseman, Evan T. Powers, Joel N. Buxbaum, Jeffery W. Kelly, and William E. Balch. An Adaptable Standard for Protein Export from the Endoplasmic Reticulum. *Cell*, 131(4):809-821, November 2007.

[20] M. Adimy, L. Babin, and L. Pujo-Menjouet. Neuron Scale Modeling of Prion Production with the Unfolded Protein Response. *SIAM Journal on Applied Dynamical Systems*, 21(4):2487-2517, 2022.

[21] O. Diekmann, J. A. P. Heesterbeek, and J. A. J. Metz. On the definition and the computation of the basic reproduction ratio R 0 in models for infectious diseases in heterogeneous populations. *Journal of mathematical biology*, 28(4):365-382, 1990.

[22] P. Van den Driessche and J. Watmough. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. *Mathematical Biosciences*, 180(1):29-48, 2002.

[23] O. Diekmann, J. A. P. Heesterbeek, and M. G. Roberts. The construction of next-generation matrices for compartmental epidemic models. *Journal of the royal society interface*, 7(47):873-885, 2010.

[24] J. K. Hale and S. M. Verduyn Lunel. *Introduction to Functional Differential Equations*. Springer, 1993.

[25] Y. Kuang. *Delay Differential Equations: With Applications in Population Dynamics*. Academic Press, 1993.

[26] H. Smith. *An Introduction to Delay Differential Equations with Applications to the Life Sciences*. Texts in Applied Mathematics. Springer, 2011.

[27] Oscar Rojo and Héctor Rojo. Some results on symmetric circulant matrices and on symmetric centrosymmetric matrices. *Linear Algebra and its Applications*, 392:211-233, 2004.

[28] Kaveh Mousavand Dan Kucerovsky and Aydin Sarraf. On some properties of toeplitz matrices. *Cogent Mathematics*, 3(1):1154705, 2016.

[29] S. Ottaviano, M. Sensi, and S. Sottile. Global stability of multi-group SAIRS epidemic models. *Mathematical Methods in the Applied Sciences*, 2023.

[30] P. Van den Driessche and J. Watmough. Further notes on the basic reproduction number. *Mathematical epidemiology*, pages 159-178, 2008.

[31] W. Kulpa. The Poincaré-Miranda theorem. *The American Mathematical Monthly*, 104(6):545-550, 1997.

[32] J. Mawhin. Simple proofs of the Hadamard and Poincaré-Miranda theorems using the Brouwer fixed point theorem. *The American Mathematical Monthly*, 126(3):260-263, 2019.

[33] H. L Smith and H. R Thieme. *Dynamical Systems and Population Persistence*. Americal Mathematical Society, 2011.

[34] H. K. Khalil. Nonlinear systems third edition. *Patience Hall*, 115, 2002.

[35] Z. Sun. A gathering of Barbalat's lemmas and their (unsung) cousins. *arXiv preprint arXiv:2301.00466*, 2023.

[36] H. I. Freedman and P. Moson. Persistence definitions and their connections. *Proceedings of the american mathematical society*, 109(4):1025-1033, 1990.

[37] M. Adimy, A. Chekroun, and C. P. Ferreira. Global dynamics of a differential-difference system: a case of Kermack-McKendrick SIR model with age-structured protection phase. *Mathematical Biosciences and Engineering*, 17(2):1329-1354, 2020.