Mohammad Ghousein
email: mohammad.ghousein@univ-poitiers.fr

Emmanuel Moulay
email: emmanuel.moulay@univ-poitiers.fr

| Patrick Coirault

Patrick Coirault
email: patrick.coirault@univ-poitiers.fr

Semi-Markov models of epidemics over networks with time delays

Keywords: Epidemic models, semi-Markov theory, time-delay systems, multi-agent systems

In this work, we extend the Markov models describing the susceptible-infectedsusceptible (SIS) epidemics over undirected networks to take into account the virus minimum incubation period and the minimum recovery period of an infected individual. We represent both periods as time delays in the states of the extended model. However, due to the addition of time delays, the process loses its Markovain property. We use the generalized semi-Markov theory to introduce both the incubation and recovery delays to the probabilistic dynamical models. Hence, in this paper, we propose a time-delay version of the two principal models of the SIS epidemics over undirected networks: the exact 2 𝑁 -state model and the approximated 𝑁-intertwined model. Finally, using Lyapunov analysis, we give sufficient conditions that guarantee the global exponential stability of the time-delay 𝑁-intertwined model.

INTRODUCTION

Since the seminal work of Kermack and McKendrick 1 introducing the SIS model, where S stands for susceptible and I for infectious, a lot of research has been done for developing mathematical epidemic models that studies the spread of an infectious disease in a certain population. The SIS class can model the propagation of specific viruses in which immunity is not acquired upon recovery from infection, and individuals become susceptible again, for example the cold flu and influenza. The most famous SIS models are the compartmental models which assume homogeneous mixing of the population [START_REF] Kermack | Contributions to the mathematical theory of epidemics. II.-The problem of endemicity[END_REF] . These models do not allow for heterogeneous contacts between individuals, a limitation that is overcomed by the so called degree-based mean field (DBMF) models [START_REF] Pastor-Satorras | Evolution and structure of the Internet: A statistical physics approach[END_REF][START_REF] Barrat | Dynamical processes on complex networks[END_REF] . The DBMF models divide the population into a set of groups. Individuals belonging to the same group have the same number of connections and equivalently have the same probability to become infected. The drawback of the compartmental and DBMF models is that they do not precisely compute the probability of infection of each individual in the population. Knowing every individual infection probability helps in using mobile applications that give each one specific isolation instructions, instead of imposing centralized lockdowns which are economically heavy [START_REF] Ferretti | Quantifying SARS-CoV-2 transmission suggests epidemic control with digital contact tracing[END_REF] . The epidemic models that serve this goal are called multi-agent epidemic models. They can describe the spread of infections between countries, provinces, cities, towns, districts, computers, individuals and many more [START_REF] Paré | Modeling, estimation, and analysis of epidemics over networks: An overview[END_REF] . Their advantages is to consider the network heterogeneity on a graph of nodes and to assign to each class of the model a probability of infection. Ganesh et al. [START_REF] Ganesh | The effect of network topology on the spread of epidemics[END_REF] have proposed a multi-agent probablistic model for the SIS class over graphs using the theory of Markov chains. The model is in continuous-time and it is of 2 𝑁 states, 𝑁 is the population size. Their analysis is mainly focused on estimating how fast or slow an epidemic will die out. Afterwards, Van Mieghem et al. 7 have extended the work of Ganesh et receiving an infection min. recovery period min. incubation period FIGURE 1 Sample path of the stochastic process 𝑋 𝑖 (𝑡) which is the state of agent 𝑖 at time 𝑡. In an SIS model, 𝑋 𝑖 (𝑡) is either susceptible 𝑋 𝑖 (𝑡) = 𝑆 or infected 𝑋 𝑖 (𝑡) = 𝐼. 𝜏 0 and 𝜏 1 are the constant time delays for the minimum incubation and recovery periods respectively. al. [START_REF] Ganesh | The effect of network topology on the spread of epidemics[END_REF] by approximating the exact model with an 𝑁-intertwined model of 𝑁 states using a mean field approximation. Since each virus has a incubation period, which is the time to be infectious after being infected, there is a minimum intrinsic time delay 𝜏 0 before transitioning from the susceptible state S to the infected state I. After becoming infected, there is also a minimum time for recovery 𝜏 1 to return back susceptible, see Figure 1. The models of Ganesh et al. [START_REF] Ganesh | The effect of network topology on the spread of epidemics[END_REF] and Van Mieghem et al. 7 do not consider this time-delay aspect in their modeling. Because of the time delays, the waiting time in the states is not any more exponentially distributed and the SIS process loses its Markovian property as shown in Figure 1. The main goal of this paper is to introduce the minimum incubation time delay 𝜏 0 and the minimum recovery time delay 𝜏 1 into the SIS multi-agent model of Ganesh et al. [START_REF] Ganesh | The effect of network topology on the spread of epidemics[END_REF] and Van Mieghem et al. 7 . Note that we refer to 𝜏 0 as the minimum incubation period and to 𝜏 1 as the minimum recovery period. Actually, the incubation period and the recovery period vary between each individual, and hence, 𝜏 0 and 𝜏 1 are random variables and are not constants. However, in this paper we define them as the minimum periods for incubation and recovery and they are considered as constants.

Introducing time delays to the continuous-time SIS-SIR models, where R stands for recovery, is well studied by the research community. Hethcote et al. [START_REF] Hethcote | An SIS epidemic model with variable population size and a delay[END_REF] , J.Arino et al. [START_REF] Arino | Time delays in epidemic models[END_REF]Chapter 13 , Tien et al. [START_REF] Tien | Multiple transmission pathways and disease dynamics in a waterborne pathogen model[END_REF] , Bretta et al. [START_REF] Beretta | Global stability of an SIR epidemic model with time delays[END_REF] and Cooke et al. [START_REF] Cooke | Analysis of an SEIRS epidemic model with two delays[END_REF] , not to mention many others, have considered introducing time delays to the compartmental models SIS, SIR, SEIRS, etc. They have assumed general probabilistic distributions for the sojourn times in the states and they have encapsulated different time delays between the different compartments. The authors have also studied the existence of solutions of these models, calculated their equilibrium points and analyses the epidemic stability. The models are based on homogeneous contacts between individuals, i.e. they are based on the assumption that all individuals are equally likely to be infected or to be recovered. Another type of models for SIS-SIR classes that could incorporate time delays is modeling using the piecewise deterministic Markov process (PDMP). Clancy et al. [START_REF] Clancy | SIR epidemic models with general infectious period distribution[END_REF][START_REF] Gómez-Corral | On SIR epidemic models with generally distributed infectious periods: Number of secondary cases and probability of infection[END_REF] use the theory of PDMP developed by Davis et al. [START_REF] Davis | Piecewise-deterministic Markov processes: a general class of non-diffusion stochastic models[END_REF] to model the SIR epidemic. The idea is to consider infectious periods that are not necessarily exponentially distributed. Finally, in the work of Huang et al. [START_REF] Huang | Stability analysis of SIR model with distributed delay on complex networks[END_REF][START_REF] Wang | Global stability analysis of an SIR epidemic model with demographics and time delay on networks[END_REF] the time delays are encapsulated in the DBMF models developed by Pastor et al. [START_REF] Pastor-Satorras | Evolution and structure of the Internet: A statistical physics approach[END_REF][START_REF] Barrat | Dynamical processes on complex networks[END_REF] , and their results is focused in formulating the stability conditions in terms of the basic reproduction number 𝑅 0 . These models predict the infection probability of groups that have the same number of connections. We should also mention that the models of Huang et al. [START_REF] Huang | Stability analysis of SIR model with distributed delay on complex networks[END_REF][START_REF] Wang | Global stability analysis of an SIR epidemic model with demographics and time delay on networks[END_REF] consider only time delays for the incubation period. The main difference between all these works and our work is in modeling the contacts between the individuals. We model contacts between individuals over general graphs where we consider that each individual in the population is a node in the graph and therefore it has its specific number of neighbors and has specific number of contacts. Hence, our modeling introduces the contacts heterogeneity of the network through graphs and we give a probability of infection for each individual of the population at each instant of time. Now, we look on multi-agent models that could incorporate time delays. Bovenkamp et al. [START_REF] Cator | Susceptible-infected-susceptible epidemics on networks with general infection and cure times[END_REF] have proposed a steady state model for the SIS infection process with general curing and infection times over undirected graphs with 𝑁 nodes. Ogura et al. [START_REF] Ogura | Stability of SIS spreading processes in networks with non-Markovian transmission and recovery[END_REF] have also explored the exponential decay rate of the SIS processes toward the disease free equilibrium over undirected graphs. They consider non-markovian processes and allows for a "phase-type" distributions for the waiting times. Kiss et al. [START_REF] Kiss | Generalization of pairwise models to non-Markovian epidemics on networks[END_REF] have proposed a delayed SIS model for pairwise correlated networks. The time delay is considered from the class S to the class I. Hang-Hyun et al. [START_REF] Jo | Analytically solvable model of spreading dynamics with non-Poissonian processes[END_REF] have modeled the (SI) infection process over a temporal network of infinite size. The goal is to study the bursty nature of some dynamical processes that appear in nature such as earthquake statistics, financial interactions, etc. In these processes, the distribution of the inter-event times follows a power-law and not the homogeneous Poisson process. Furthermore, Doerr et al. [START_REF] Doerr | Lognormal infection times of online information spread[END_REF] have modeled the distribution of the spreading time 𝑇 of information in social networks as a convolution of two lognormally distributed random variables.

Our contribution in this paper is to propose a multi-agent model for the SIS epidemics that includes time delays for the incubation and the recovery periods. The modeling is done over a finite-size graph composed of 𝑁 agents and we consider no correlations between the nodes. It is a dynamical model that predicts at each instant of time a probability of infection for each individual in the population. We approximate the exact model which has 2 𝑁 states with an 𝑁-intertwined model of 𝑁 states using a mean field approximation. While the multi-agent dynamical model for SIS epidemics is presented by Ganesh et al. [START_REF] Ganesh | The effect of network topology on the spread of epidemics[END_REF] and Van Mieghem et al. 7 , their model is delay free and assumes the Markovian property. We upgrade this model to incorporate time delays by using the semi-Markovian analysis and the resulting model coincides with the model of Van Mieghem et al. 7 when the delays are zero. Using Lyapunov analysis, we then obtain sufficient conditions for the epidemic to die out. The stability conditions of the disease free equilibrium point (DFE) coincides with the results of Ganesh et al. [START_REF] Ganesh | The effect of network topology on the spread of epidemics[END_REF] , Wang et al. [START_REF] Wang | Epidemic spreading in real networks: An eigenvalue viewpoint[END_REF] and Paré et al. [START_REF] Paré | Epidemic processes over time-varying networks[END_REF] for the delay-free case, more precisely Condition 1 in the article of Ganesh et al. [START_REF] Ganesh | The effect of network topology on the spread of epidemics[END_REF] , Theorem 1 in the article of Paré et al. [START_REF] Paré | Epidemic processes over time-varying networks[END_REF] and Corollary 1 in the article of Wang et al. [START_REF] Wang | Epidemic spreading in real networks: An eigenvalue viewpoint[END_REF] ).

The paper is organized as follows. After some preliminaries on semi-Markov processes in Section 2, we present the exact 2 𝑁 states semi-Markov model of the SIS class in Section 3 where 𝑁 is the size of the population. In Section 4, the exact model is approximated by the 𝑁-intertwined model with 𝑁 states. Stability results are provided in Section 5. Section 6 is devoted to numerical simulations and the paper is concluded in Section 7.

Notations

{𝑍(𝑡), 𝑡 ≥ 0} Ω is a stochastic process defined over the discrete state space Ω. 𝑃 𝑟(𝑋 = 𝑥) is the probability that the random variable 𝑋 takes the value 𝑥. 𝐑 is the set of real numbers. Let 𝑓 (𝑡) ∶ [0, +∞) → 𝐑 be a piece-wise continuous function, then 𝑓 𝑙 (𝑠) = 𝐿{𝑓 (𝑡)} is the Laplace transform of 𝑓 (𝑡). 𝐻(𝑢) is the Heavise function, i.e. 𝐻(𝑢) = 1 if 𝑢 ≥ 0 and 𝐻(𝑢) = 0 if 𝑢 < 0. 𝛿 𝑖𝑗 is the Kronecker delta function, 𝛿 𝑖𝑗 = 1 if 𝑖 = 𝑗 and 𝛿 𝑖𝑗 = 0 otherwise. (𝑣) 𝑛×1 denotes a vector in 𝐑 𝑛 . 𝑑𝑖𝑎𝑔(𝑣) is the 𝑛 × 𝑛 diagonal matrix formed with elements of 𝑣. 𝟏 𝑛 is the vector of all ones of length 𝑛. 𝐈 𝑛×𝑛 is the identity matrix. The indicator function is defined by 1 𝑥 = 1 if 𝑥 is true and 1 𝑥 = 0 if 𝑥 is false.

PRELIMINARIES ON SEMI-MARKOV PROCESSES

Let us begin by introducing some necessary concepts on semi-Markov processes that we shall need later in our analysis. The reader may refer to the chapter 5 of the book of Ross [START_REF] Ross | Applied probability models with optimization applications[END_REF] for more details on semi-Markov processes. A stochastic process {𝑆(𝑡), 𝑡 ≥ 0} Ω that makes a transition from one state to the other while satisfying the Markov property at the instants of the jumps, but in which the time spent in any state is random, is called a semi-Markov process. More precisely, consider the following definition.

Definition 1.

Let Ω = {1, … , 𝑁} be a discrete state space. A stochastic process {𝑆(𝑡), 𝑡 ≥ 0} Ω is called a time-homogeneous semi-Markov process on the probability space with finite state space Ω if the following conditions hold:

1. the sample paths {𝑆(𝑡), 𝑡 ≥ 0} Ω are right-continuous step functions and the probability of the left-handed limits is equal to one;

2. denote by 0 = 𝑇 0 < 𝑇 1 < … < 𝑇 𝑘 < … the jump times of (𝑆 𝑡 ) 𝑡≥0 and by 𝜏 𝑘 = 𝑇 𝑘 -𝑇 𝑘-1 the sojourn times in the states, the process {𝑆(𝑡), 𝑡 ≥ 0} Ω satisfies the Markov property at each jump instant 𝑇 𝑘 , 𝑘 = 0, 1, 2, …

𝑃 𝑟 ( 𝑆(𝑇 𝑘+1 ) = 𝑗|𝑆(𝑇 𝑘 ) = 𝑖, 𝑆(𝑇 𝑘-1 ) = 𝑚, … , 𝑆(𝑇 0 ) = 𝑙 ) = 𝑃 𝑟 ( 𝑆(𝑇 𝑘+1 ) = 𝑗|𝑆(𝑇 𝑘 ) = 𝑖 ) , 𝑖, 𝑗, 𝑚, … , 𝑙 ∈ Ω 3. 𝐹 𝑖𝑗 (𝑡) = 𝑃 𝑟 ( 𝜏 𝑘+1 ≤ 𝑡 | 𝑆(𝑇 𝑘 ) = 𝑖, 𝑆(𝑇 𝑘+1 ) = 𝑗 )
does not depend on the jump instant 𝑇 𝑘 .

The time-homogeneous semi-Markov process {𝑆(𝑡), 𝑡 ≥ 0} Ω can be described by the following two notions:

-the transition probability matrix (𝑃 ) 𝑁×𝑁 such that:

𝑝 𝑖𝑗 = 𝑃 𝑟 ( 𝑆(𝑇 𝑘+1 ) = 𝑗 | 𝑆(𝑇 𝑘 ) = 𝑖 ) , 𝑖, 𝑗 ∈ Ω
where 𝑝 𝑖𝑗 is the conditional probability that the process will make a transition from the state 𝑖 to the state 𝑗 at the instant 𝑇 𝑘+1 , 𝑘 ≥ 0, and ∑ 𝑗∈Ω 𝑝 𝑖𝑗 = 1, 𝑖 ∈ Ω; -the holding time density 𝑓 𝑖𝑗 (𝑡) in a state 𝑖 given that the next transition is to state

𝑗 𝑓 𝑖𝑗 (𝑡)𝑑𝑡 = 𝑃 𝑟 ( 𝑡 ≤ 𝜏 𝑘+1 ≤ 𝑡 + 𝑑𝑡 | 𝑆(𝑇 𝑘 ) = 𝑖, 𝑆(𝑇 𝑘+1 ) = 𝑗 )
where 𝑡 ≥ 0 and 𝑖, 𝑗 ∈ Ω. Hence, one can think of generating a semi-Markov process as follows: after having entered state 𝑖, make a random draw to determine the successor state 𝑗 based on the set of transition probabilities 𝑝 𝑖𝑗 , 𝑖, 𝑗 ∈ Ω, then make another draw to determine the waiting time in 𝑖 before transiting to 𝑗 based on the holding densities 𝑓 𝑖𝑗 (𝑡), 𝑖, 𝑗 ∈ Ω.

Another useful notion in semi-Markov processes is the transition density functions defined by

ℎ 𝑖𝑗 (𝑡) = 𝑝 𝑖𝑗 𝑓 𝑖𝑗 (𝑡), 𝑖, 𝑗 ∈ Ω (1)
where

ℎ 𝑖𝑗 (𝑡)𝑑𝑡 = 𝑃 𝑟 ( 𝑡 ≤ 𝜏 𝑘+1 ≤ 𝑡 + 𝑑𝑡, 𝑆(𝑇 𝑘+1 ) = 𝑗 | 𝑆(𝑇 𝑘 ) = 𝑖 )
is the probability of a transition to state 𝑗 between 𝑡 and 𝑡 + 𝑑𝑡 given that the last transition was to state 𝑖. Let 𝑤 𝑖 (𝑡) be the unconditional waiting time density in a state 𝑖 ∈ Ω,

𝑤 𝑖 (𝑡) = ∑ 𝑗∈Ω ℎ 𝑖𝑗 (𝑡) (2)
then 𝑤 𝑖 (𝑡)𝑑𝑡 is the probability of leaving state 𝑖 between 𝑡 and 𝑡 + 𝑑𝑡. After recalling the basic features of a semi-Markov process, we can now directly give the dynamics of the state transition probabilities 𝜙 𝑖𝑗 (𝑡) = 𝑃 𝑟(𝑆(𝑡) = 𝑗 | 𝑆(0) = 𝑖), for 𝑡 ≥ 0 and 𝑖, 𝑗 ∈ Ω,

𝜙 𝑖𝑗 (𝑡) = 𝛿 𝑖𝑗 𝑊 𝑖 (𝑡) + 𝑁 ∑ 𝑘=1 𝑡 ∫ 0 ℎ 𝑖𝑘 (𝑣)𝜙 𝑘𝑗 (𝑡 -𝑣)𝑑𝑣 ( 3
)
where 𝛿 𝑖𝑗 = 1 if 𝑖 = 𝑗 and 𝛿 𝑖𝑗 = 0 if 𝑖 ≠ 𝑗, and 𝑊 𝑖 (𝑡) = 1 -∫ 𝑡 0 𝑤 𝑖 (𝑣)𝑑𝑣. By (3), 𝜙 𝑖𝑗 (𝑡), which is the probability that the process is in state 𝑗 at time 𝑡, given that it entered state 𝑖 at time zero, is a sum of the following two terms:

1. the non-zero term 𝑊 𝑖 (𝑡) when 𝑖 = 𝑗, which represents the probability that the process may have not left state 𝑖 by time 𝑡; 2. the integral over all 𝑣 between 0 and 𝑡 of ℎ 𝑖𝑘 (𝑣)𝑑𝑣 (probability of a transition to state 𝑘 between 𝑣 and 𝑣+𝑑𝑣) multiplied by 𝜙 𝑘𝑗 (𝑡 -𝑗) (probability that the process will reach state 𝑗 in the remaining time 𝑡 -𝑣 after having entered state 𝑘) summed all over the states 𝑘.

In what follows, we use the semi-Markov theory to model the SIS epidemic over undirected networks while taking into account the time delays.

2 𝑁 -STATE SEMI-MARKOV MODEL

We model the spread of an infection in a certain population by using an undirected graph 𝐺(𝑁, 𝐸) with a fixed adjacency matrix 𝐴. 𝐴 is a binary matrix, i.e. 𝑎 𝑖𝑗 = 1 if individual 𝑖 is connected to individual 𝑗 otherwise 𝑎 𝑖𝑗 = 0, and 𝑎 𝑖𝑖 = 0. For details on graph theory, the reader may refer to Chapter 3 of the book of Bullo et al. [START_REF] Bullo | Lectures on network systems[END_REF] . The arrival of an infection to a link and the curing process of an infected node is modeled using two independent Poisson processes of rates 𝛽 and 𝛿 respectively. Each node in the graph represents an agent in the population whom can be either healthy or infected at time 𝑡. This type of modeling belongs to the class of SIS models in continuous time. Define Ω = {0, 1, 2, … , 2 𝑁 -1} to represent all the possible states of the population at time 𝑡 which are detailled in Table 1. State 0 means that all the agents in the population are healthy, state 2 𝑁 -1 means that all the agents are infected and so on, see Page 3 of Van Mieghem et al. 7 .

⋯ 𝑥 2 𝑥 1 0 0 0 ⋯ 0 0 1 0 0 ⋯ 0 1 ⋮ ⋮ 2 𝑁 -1 1 1 ⋯ 1 1
Without time delays, the susceptible-infection process defines a continous-time Markov chain defined on the discrete state space Ω. The infinitesimal rate matrix of this Markov chain is defined, similarly to equation (4) of Van Mieghem et al. 7 , by:

𝑞 𝑖𝑗 = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝑞 𝛿 𝑖𝑗 = 𝛿 if 𝑖 = 𝑗 + 2 𝑚-1 , 𝑚 = 1, … , 𝑁, 𝑥 𝑚 = 1 𝑞 𝐴 𝑖𝑗 = 𝛽 ∑ 𝑁 𝑘=1 𝑎 𝑚𝑘 𝑥 𝑘 if 𝑖 = 𝑗 -2 𝑚-1 , 𝑚 = 1, … , 𝑁, 𝑥 𝑚 = 0 - ∑ 𝑁 𝑘=1;𝑘≠𝑗 𝑞 𝑘𝑗 if 𝑖 = 𝑗 0 else (4)
where 7 . The novelty of our approach with respect to the work of the authors Van Mieghem et al. 7 is to consider that when an agent 𝑖 receives an infection, it cannot be infected or infectious instantly, there is a minimum time 𝜏 0 before transiting from 𝑆 to 𝐼 because of the disease incubation period as shown on Figure 1. Moreover, when an agent 𝑖 is infected, it cannot recover instantly and return back to be susceptible, there is a minimum time for recovery to transition from 𝐼 to 𝑆 represented by the minimum time 𝜏 1 . Adding a minimum time to jump from 𝑆 to 𝐼 or from 𝐼 to 𝑆 does not affect the Markov property at the instants of the jumps, but affects the waiting time in a given state that becomes not exponentially distributed, see Figure 1. Hence, in our case, the process defined over the discrete state space Ω is not a continuous-time Markov process but a semi-Markov process. Let {𝑍(𝑡), 𝑡 ≥ 0} Ω be the semi-Markov process defined over the discrete state space Ω. 𝑍(𝑡) models the evolution of the SIS epidemic with time delays over the graph 𝐺. Recall that the objective is to find the probability of infection for each agent as a function of time in the presence of the time delays. We have mentioned in Section 2 that in order to generate the semi-Markov process {𝑍(𝑡), 𝑡 ≥ 0} Ω , the first step is to specify the transition probability matrix (𝑃 ) 2 𝑁 ×2 𝑁 and the holding time densities 𝑓 𝑖𝑗 (𝑡), 𝑖, 𝑗 ∈ Ω. Let us start by computing (𝑃 ) 2 𝑁 ×2 𝑁 . Since {𝑍(𝑡), 𝑡 ≥ 0} Ω is Markovian at the instants of the jumps, the transition probability matrix (𝑃 ) 2 𝑁 ×2 𝑁 is the same as that in the continuous-time Markov model given in Van Mieghem et al. 7 . Hence, we can define the transition probability matrix (𝑃 ) 2 𝑁 ×2 𝑁 for the semi-Markov process {𝑍(𝑡), 𝑡 ≥ 0} Ω as:

𝑖 = ∑ 𝑁 𝑘=1 2 𝑘-1 𝑥 𝑘 ,
𝑝 𝑖𝑗 = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ - 𝑞 𝛿 𝑖𝑗 𝑞 𝑖𝑖 if 𝑖 = 𝑗 + 2 𝑚-1 , 𝑚 = 1, … , 𝑁, 𝑥 𝑚 = 1 - 𝑞 𝐴 𝑖𝑗 𝑞 𝑖𝑖 if 𝑖 = 𝑗 -2 𝑚-1 , 𝑚 = 1, … , 𝑁, 𝑥 𝑚 = 0 0 else
Note that 𝑝 𝑖𝑗 is just the ratio of the rate to be cured or to be infected over the total rate of leaving a given state. One can notice also that 𝑝 𝑖𝑖 = 0, and this means that we only allow real transitions from 𝑖 to 𝑗 where 𝑖 ≠ 𝑗. The second step to generate the semi-Markov process {𝑍(𝑡), 𝑡 ≥ 0} Ω is to define the densities 𝑓 𝑖𝑗 (𝑡). We will directly define ℎ 𝑖𝑗 (𝑡) as it will be used later in the probability dynamics (𝑓 𝑖𝑗 (𝑡) can be directly calculated from (1) as well as 𝑤 𝑖 (𝑡) can be calculated from (2)):

ℎ 𝑖𝑗 (𝑡) = ⎧ ⎪ ⎨ ⎪ ⎩ 𝑞 𝛿 𝑖𝑗 𝑒 𝑞 𝑖𝑖 (𝑡-𝜏 1 ) 𝐻(𝑡 -𝜏 1 ) if 𝑖 = 𝑗 + 2 𝑚-1 , , 𝑚 = 1, … , 𝑁, 𝑥 𝑚 = 1 𝑞 𝐴 𝑖𝑗 𝑒 𝑞 𝑖𝑖 (𝑡-𝜏 0 ) 𝐻(𝑡 -𝜏 0 ) if 𝑖 = 𝑗 -2 𝑚-1 , 𝑚 = 1, … , 𝑁, 𝑥 𝑚 = 0 0 else (5)
The equations of ℎ 𝑖𝑗 (𝑡) in (5) mean that the holding time in a state 𝑖 before transitioning to a state 𝑗 is a delayed exponential random variable. The process cannot transition from a state 𝑖 to a state 𝑗 before the minimum time delays of the incubation period or the recovery period are ended, as explained in Figure 1. We have set all the distribution properties of the semi-Markov process {𝑍(𝑡), 𝑡 ≥ 0} Ω . We start now computing the dynamics of the state transition probabilities 𝜙 𝑖𝑗 (𝑡). Take the Laplace transform of (3) to get:

𝜙 𝑙 𝑖𝑗 (𝑠) = 𝛿 𝑖𝑗 ( 1 -𝑤 𝑙 𝑖 (𝑠) 𝑠 ) + 2 𝑁 -1 ∑ 𝑘=0 ℎ 𝑙 𝑖𝑘 (𝑠)𝜙 𝑙 𝑘𝑗 (𝑠) (6)
Now, we compute the Laplace transform ℎ 𝑖𝑗 (𝑡) in ( 5) and 𝑤 𝑖 (𝑡) in (2) as:

ℎ 𝑙 𝑖𝑗 (𝑠) = ⎧ ⎪ ⎨ ⎪ ⎩ 𝑒 -𝑠𝜏 1 𝑠-𝑞 𝑖𝑖 𝑞 𝛿 𝑖𝑗 if 𝑖 = 𝑗 + 2 𝑚-1 , , 𝑚 = 1, … , 𝑁, 𝑥 𝑚 = 1 𝑒 -𝑠𝜏 0 𝑠-𝑞 𝑖𝑖 𝑞 𝐴 𝑖𝑗 if 𝑖 = 𝑗 -2 𝑚-1 , 𝑚 = 1, … , 𝑁, 𝑥 𝑚 = 0 0 else (7) 𝑤 𝑙 𝑖 (𝑠) = 𝑒 -𝑠𝜏 0 𝑠 -𝑞 𝑖𝑖 ∑ 𝑗 𝑞 𝐴 𝑖𝑗 + 𝑒 -𝑠𝜏 1 𝑠 -𝑞 𝑖𝑖 ∑ 𝑗 𝑞 𝛿 𝑖𝑗 (8)
For 𝑖 = 𝑗, equation ( 6) becomes:

𝜙 𝑙 𝑖𝑖 (𝑠) = 1 𝑠 (1 -𝑤 𝑙 𝑖 (𝑠)) + 2 𝑁 -1 ∑ 𝑘=0 ℎ 𝑙 𝑖𝑘 (𝑠)𝜙 𝑙 𝑘𝑖 (𝑠) (9)
By substituting ( 7) and ( 8) in ( 9) one gets:

𝜙 𝑙 𝑖𝑖 (𝑠) = 1 𝑠 ( 1 - 𝑒 -𝑠𝜏 0 𝑠 -𝑞 𝑖𝑖 ∑ 𝑗 𝑞 𝐴 𝑖𝑗 - 𝑒 -𝑠𝜏 1 𝑠 -𝑞 𝑖𝑖 ∑ 𝑗 𝑞 𝛿 𝑖𝑗 ) + 𝑒 -𝑠𝜏 0 𝑠 -𝑞 𝑖𝑖 ∑ 𝑘 𝑞 𝐴 𝑖𝑘 𝜙 𝑙 𝑘𝑖 (𝑠) + 𝑒 -𝑠𝜏 1 𝑠 -𝑞 𝑖𝑖 ∑ 𝑘 𝑞 𝛿 𝑖𝑘 𝜙 𝑙 𝑘𝑖 (𝑠) (10)
Multiplying both sides of ( 10) by (𝑠 -𝑞 𝑖𝑖 ) leads to:

𝑠𝜙 𝑙 𝑖𝑖 (𝑠) -𝑞 𝑖𝑖 𝜙 𝑙 𝑖𝑖 (𝑠) = 1 - 𝑞 𝑖𝑖 𝑠 - 𝑒 -𝑠𝜏 0 𝑠 ∑ 𝑗 𝑞 𝐴 𝑖𝑗 - 𝑒 -𝑠𝜏 1 𝑠 ∑ 𝑗 𝑞 𝛿 𝑖𝑗 + ∑ 𝑗 𝑞 𝐴 𝑖𝑗 𝜙 𝑙 𝑗𝑖 (𝑠)𝑒 -𝑠𝜏 0 + ∑ 𝑗 𝑞 𝛿 𝑖𝑗 𝜙 𝑙 𝑗𝑖 (𝑠)𝑒 -𝑠𝜏 1 (11)
If 𝑖 ≠ 𝑗, then (6) becomes:

𝜙 𝑙 𝑖𝑗 (𝑠) = 2 𝑁 -1 ∑ 𝑘=0 ℎ 𝑙 𝑖𝑘 (𝑠)𝜙 𝑙 𝑘𝑗 (𝑠) = 𝑒 -𝑠𝜏 0 𝑠 -𝑞 𝑖𝑖 ∑ 𝑘 𝑞 𝐴 𝑖𝑘 𝜙 𝑙 𝑘𝑗 (𝑠) + 𝑒 -𝑠𝜏 1 𝑠 -𝑞 𝑖𝑖 ∑ 𝑘 𝑞 𝛿 𝑖𝑘 𝜙 𝑙 𝑘𝑗 (𝑠) (12)
Multiplying both sides of ( 12) by (𝑠 -𝑞 𝑖𝑖 ) leads directly to:

𝑠𝜙 𝑙 𝑖𝑗 (𝑠) -𝑞 𝑖𝑖 𝜙 𝑙 𝑖𝑗 (𝑠) = ∑ 𝑘 𝑞 𝐴 𝑖𝑘 𝜙 𝑙 𝑘𝑗 (𝑠)𝑒 -𝑠𝜏 0 + ∑ 𝑘 𝑞 𝛿 𝑖𝑘 𝜙 𝑙 𝑘𝑗 (𝑠)𝑒 -𝑠𝜏 1 (13)
Define the matrix (Φ(𝑡)) 2 𝑁 ×2 𝑁 with elements 𝜙 𝑖𝑗 (𝑡), 𝑖, 𝑗 ∈ Ω. By taking the inverse Laplace transform of (11) one gets

φ𝑖𝑖 (𝑡) + 𝛿(𝑡)𝜙 𝑖𝑖 (0) -𝑞 𝑖𝑖 𝜙 𝑖𝑖 (𝑡) =𝛿(𝑡) -𝑞 𝑖𝑖 -𝐻(𝑡 -𝜏 0 ) ∑ 𝑗 𝑞 𝐴 𝑖𝑗 -𝐻(𝑡 -𝜏 1 ) ∑ 𝑗 𝑞 𝛿 𝑖𝑗 + 𝐻(𝑡 -𝜏 0 ) ∑ 𝑘 𝑞 𝐴 𝑖𝑘 𝜙 𝑘𝑖 (𝑡 -𝜏 0 ) + 𝐻(𝑡 -𝜏 1 ) ∑ 𝑘 𝑞 𝛿 𝑖𝑘 𝜙 𝑘𝑖 (𝑡 -𝜏 1 ) (14)
where 𝛿(𝑡) is the Dirac function. Since 𝜙 𝑖𝑖 (0) = 1 which is the initial condition, ( 14) becomes:

φ𝑖𝑖 (𝑡) = 𝑞 𝑖𝑖 𝜙 𝑖𝑖 (𝑡) - ( 𝑞 𝑖𝑖 + 𝐻(𝑡 -𝜏 0 ) ∑ 𝑗 𝑞 𝐴 𝑖𝑗 + 𝐻(𝑡 -𝜏 1 ) ∑ 𝑗 𝑞 𝛿 𝑖𝑗 ) + 𝐻(𝑡 -𝜏 0 ) ∑ 𝑘 𝑞 𝐴 𝑖𝑘 𝜙 𝑘𝑖 (𝑡 -𝜏 0 ) + 𝐻(𝑡 -𝜏 1 ) ∑ 𝑘 𝑞 𝛿 𝑖𝑘 𝜙 𝑘𝑖 (𝑡 -𝜏 1 ) (15)
Now, take the inverse Laplace transform of (13) to get:

φ𝑖𝑗 (𝑡) + 𝛿(𝑡)𝜙 𝑖𝑗 (0) -𝑞 𝑖𝑖 𝜙 𝑖𝑗 (𝑡) = 𝐻(𝑡 -𝜏 0 ) ∑ 𝑘 𝑞 𝐴 𝑖𝑘 𝜙 𝑘𝑗 (𝑡 -𝜏 0 ) + 𝐻(𝑡 -𝜏 1 ) ∑ 𝑘 𝑞 𝛿 𝑖𝑘 𝜙 𝑘𝑗 (𝑡 -𝜏 1 ) (16)
Since 𝜙 𝑖𝑗 = 0 which is the initial condition, ( 16) becomes:

φ𝑖𝑗 (𝑡) = 𝑞 𝑖𝑖 𝜙 𝑖𝑗 (𝑡) + 𝐻(𝑡 -𝜏 0 ) ∑ 𝑘 𝑞 𝐴 𝑖𝑘 𝜙 𝑘𝑗 (𝑡 -𝜏 0 ) + 𝐻(𝑡 -𝜏 1 ) ∑ 𝑘 𝑞 𝛿 𝑖𝑘 𝜙 𝑘𝑗 (𝑡 -𝜏 1 ) (17)
Combining ( 15) and ( 17) in the matrix form leads directly to the dynamics of the state transition probability matrix Φ(𝑡):

Φ(𝑡) = 𝑄 𝑑 Φ(𝑡) + 𝐻(𝑡 -𝜏 0 )𝑄 𝐴 Φ(𝑡 -𝜏 0 ) + 𝐻(𝑡 -𝜏 1 )𝑄 𝛿 Φ(𝑡 -𝜏 1 ) - ( 𝑄 𝑑 + 𝐻(𝑡 -𝜏 0 )𝑑𝑖𝑎𝑔(𝑄 𝐴 𝟏 2 𝑁 ) + 𝐻(𝑡 -𝜏 1 )𝑑𝑖𝑎𝑔(𝑄 𝛿 𝟏 2 𝑁 ) ) (18)
where Φ(0) = 𝐈 2 𝑁 ×2 𝑁 is the matrix of the initial conditions.

Remark 1. In the delay-free case, i.e. 𝜏 0 = 𝜏 1 = 0, equation ( 18) becomes:

Φ(𝑡) = 𝑄Φ(𝑡) (19)
which is the Kolmogorov backward differential equations for continuous-time Markov chains over discrete state space. When 𝜏 0 = 𝜏 1 = 0, the waiting time in any state is exponentially distributed and the semi-Markov process {𝑍(𝑡), 𝑡 ≥ 0} becomes a continuous-time Markov chain. In this case, the dynamics of the state probability transitions in (19) are the same as the ones obtained by the authors in Van Mieghem et al. 7 for the exact 2 𝑁 -state SIS Markov models without time delays.

After obtaining the dynamics of Φ(𝑡) in ( 18), we can proceed to calculate the infection probabilities 𝑝 𝑗 (𝑡) = 𝑃 𝑟(agent 𝑗 is infected at time 𝑡). Let 𝑠 𝑖 (𝑡) = 𝑃 𝑟(𝑍(𝑡) = 𝑖) be an element of the state probability vector (𝑠(𝑡)) 2 𝑁 ×1 for all 𝑖 ∈ Ω, and 𝑝 𝑗 (𝑡) is an element of the infection probability vector (𝑌 (𝑡)) 𝑁×1 for all 𝑗 = 1, … , 𝑁. 𝑌 (𝑡) is computed as follows:

𝑌 (𝑡) = 𝑀 𝑇 𝑠(𝑡) ( 20
)
where 𝑠(𝑡) = Φ(𝑡)𝑠(0) and (𝑀) 2 𝑁 ×𝑁 is defined in Van Mieghem et al. 7 by

𝑀 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 0 ⋯ 0 1 0 0 ⋯ 0 0 1 0 ⋯ 0 1 1 0 ⋯ 0 ⋮ ⋮ ⋮ ⋮ 1 1 1 ⋯ 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ .
The time-delay system (18) with the output (20) models the probability of infection of each agent taking into account the incubation and recovery time delays. The main problem of this exact 2 𝑁 -state semi-Markov model is that the number of states 2 𝑁 increases rapidly with 𝑁 which makes it computationally heavy and difficult to implement especially for large populations. Hence, an approximation for this model is necessary and this is the subject of the next section.

𝑁-INTERTWINED SEMI-MARKOV MODEL

The exact model in Section 3 computes the probabilities of all the possible combinations of infected-susceptible agents over the network. We aim in this section at approximating this model by looking at each node directly. But before modeling the infection process at the node level, we show how to obtain the infection probability with time delays for only one isolated agent (no connections with the network), then we upgrade this model by introducing the connections with the other agents.

Semi-Markov process with time delays for one agent

Let {𝑆(𝑡), 𝑡 ≥ 0} Ω 𝑆 be a semi-Markov process with two states Ω 𝑆 = {0, 1}, and with minimum state waiting times 𝜏 0 and 𝜏 1 . This configuration is the same as the model in Section 3 for the case 𝑁 = 1 with the following 𝑄 matrix:

𝑄 = [ -𝑞 𝐴 01 𝑞 𝐴 01 𝑞 𝛿 10 -𝑞 𝛿 10 ]
where 𝑄 𝐴 = 𝑞 𝐴 01 = 𝜆 0 and 𝑄 𝛿 = 𝑞 𝛿 10 = 𝜆 1 , 𝜆 0 and 𝜆 1 are the post-delay transition rates from state 0 to state 1 and vice-versa. By inserting 𝑄 in (18), one gets the following dynamics for the state transition probabilities:

φ00 (𝑡) + 𝜆 0 𝜙 00 (𝑡) = 𝜆 0 𝐻(𝑡 -𝜏 0 )𝜙 10 (𝑡 -𝜏 0 ) + 𝜆 0 (1 -𝐻(𝑡 -𝜏 0 )) (21) φ10 (𝑡) + 𝜆 1 𝜙 10 (𝑡) = 𝜆 1 𝐻(𝑡 -𝜏 1 )𝜙 00 (𝑡 -𝜏 1 ) (22) φ01 (𝑡) + 𝜆 0 𝜙 01 (𝑡) = 𝜆 0 𝐻(𝑡 -𝜏 0 )𝜙 11 (𝑡 -𝜏 0 ) (23) φ11 (𝑡) + 𝜆 1 𝜙 11 (𝑡) = 𝜆 1 𝐻(𝑡 -𝜏 1 )𝜙 01 (𝑡 -𝜏 1 ) + 𝜆 1 (1 -𝐻(𝑡 -𝜏 1 )) (24)
with 𝜙 00 (0) = 1, 𝜙 10 (0) = 0, 𝜙 01 (0) = 0 and 𝜙 11 (0) = 1 the initial conditions. Our goal is to compute directly the dynamics of the probability of infection 𝑝(𝑡) = 𝑃 𝑟(𝑆(𝑡) = 1) for 𝑡 ≥ 0. We recall that by definition, 𝑝(𝑡) is given by:

𝑝(𝑡) = 𝜙 01 (𝑡)(1 -𝑝 0 ) + 𝜙 11 (𝑡)𝑝 0 (25)
where 𝑝(0) = 𝑝 0 is the probability of being in state 1 at time 𝑡 = 0. By taking the Laplace transform of ( 23)-( 24) one gets:

𝜙 𝑙 01 (𝑠) = 𝜆 0 𝑒 -𝑠𝜏 0 + 𝜆 0 𝜆 1 𝑠 𝑒 -𝑠𝜏 0 - 𝜆 0 𝜆 1 𝑠 𝑒 -𝑠(𝜏 0 +𝜏 1 ) 𝑠 2 + (𝜆 0 + 𝜆 1 )𝑠 + 𝜆 0 𝜆 1 (1 -𝑒 -𝑠(𝜏 0 +𝜏 1 ) ) ( 26
)
𝜙 𝑙 11 (𝑠) = 𝑠 + 𝜆 0 + 𝜆 1 (𝑠+𝜆 0 ) 𝑠 - 𝜆 1 (𝑠+𝜆 0 ) 𝑠 𝑒 -𝑠𝜏 1 𝑠 2 + (𝜆 0 + 𝜆 1 )𝑠 + 𝜆 0 𝜆 1 (1 -𝑒 -𝑠(𝜏 0 +𝜏 1 ) ) (27)
Substitute ( 26) and ( 27) in (25) to get:

𝑠𝑝 𝑙 (𝑠) + (𝜆 0 + 𝜆 1 )𝑝 𝑙 (𝑠) + 𝜆 0 𝜆 1 𝑠 (1 -𝑒 -𝑠(𝜏 0 +𝜏 1 ) )𝑝 𝑙 (𝑠) = 𝜆 0 𝑠 𝑒 -𝑠𝜏 0 + 𝜆 0 𝜆 1 𝑠 2 𝑒 -𝑠𝜏 0 - 𝜆 0 𝜆 1 𝑠 2 𝑒 -𝑠(𝜏 0 +𝜏 1 ) + ( 1 + 𝜆 0 𝑠 (1 -𝑒 -𝑠𝜏 0 ) + 𝜆 1 𝑠 (1 -𝑒 -𝑠𝜏 1 ) + 𝜆 0 𝜆 1 𝑠 2 ( 1 -𝑒 -𝑠𝜏 0 -𝑒 -𝑠𝜏 1 + 𝑒 -𝑠(𝜏 0 +𝜏 1 ) ) ) 𝑝 0 (28)
where 𝑝 𝑙 (𝑠) is the Laplace transform of 𝑝(𝑡). Now by taking the inverse Laplace transform of each term in (28), for example,

𝐿 -1 {𝑠𝑝 𝑙 (𝑠)} = ṗ(𝑡) + 𝛿(𝑡)𝑝 0 , 𝐿 -1 { 𝜆 0 𝜆 1 𝑠 (1 -𝑒 -𝑠(𝜏 0 +𝜏 1 ) )𝑝 𝑙 (𝑠)} = 𝜆 0 𝜆 1 (𝑥 𝐼 (𝑡) -𝐻(𝑡 -𝜏 0 -𝜏 1 )𝑥 𝐼 (𝑡 -𝜏 0 -𝜏 1 )
), where ẋ𝐼 (𝑡) = 𝑝(𝑡), and so on, we obtain the explicit dynamics of 𝑝(𝑡):

ẋ𝐼 (𝑡) = 𝑝(𝑡) ṗ(𝑡) = 𝜆 0 (𝐻(𝑡 -𝜏 0 ) -𝑝(𝑡)) -𝜆 1 𝑝(𝑡) -𝜆 0 𝜆 1 ( 𝑥 𝐼 (𝑡) -𝐻(𝑡 -𝐷)𝑥 𝐼 (𝑡 -𝐷) ) + 𝜆 0 𝜆 1 ( 𝐻(𝑡 -𝜏 0 )(𝑡 -𝜏 0 ) -𝐻(𝑡 -𝐷)(𝑡 -𝐷) ) + ( 𝜆 0 ( 1 -𝐻(𝑡 -𝜏 0 ) ) + 𝜆 1 ( 1 -𝐻(𝑡 -𝜏 1 ) ) + 𝜆 0 𝜆 1 ( 𝑡 -𝐻(𝑡 -𝜏 0 )(𝑡 -𝜏 0 ) -𝐻(𝑡 -𝜏 1 )(𝑡 -𝜏 1 ) + 𝐻(𝑡 -𝐷)(𝑡 -𝐷) )) 𝑝 0 (29)
where 𝐷 = 𝜏 0 + 𝜏 1 , 𝑝(0) = 𝑝 0 and 𝑥 𝐼 (0) = 0 are the system initial conditions. 𝑥 𝐼 (𝑡) represents the history of the infection. 𝑥 𝐼 (𝑡) contains information on the number of times agent 𝑖 becomes infected in the time interval [0, 𝑡]. Notice that 𝑥 𝐼 (𝑡) disappears in the delay-free case because the process is Markovian.It is known that for Markovian process, the history of the process has no effect on the future predictions. System (29) is the time-delay model of the infection probability for one agent not connected to the rest of the network. Now, we will introduce the network aspect into (29).

Remark 2. In the delay-free case, i.e. 𝜏 0 = 𝜏 1 = 0, the dynamics of 𝑝(𝑡) following ( 29) become:

ṗ(𝑡) = 𝜆 0 (1 -𝑝(𝑡)) -𝜆 1 𝑝(𝑡) ( 30
)
which are the dynamics of infection in a 2-state continuous-time Markov chain, because as mentioned in Remark 1, in the delay-free case, the semi-Markov process becomes a continuous-time Markov chain.

Multi-agent 𝑁-intertwined semi-Markov model

Let {𝑋 𝑖 (𝑡), 𝑡 ≥ 0} {0,1} be the stochastic process attached to an agent 𝑖. At any time 𝑡, 𝑋 𝑖 (𝑡) can take only two values, 𝑋 𝑖 (𝑡) = 1 if agent 𝑖 is infected and 𝑋 𝑖 (𝑡) = 0 if agent 𝑖 is healthy. If agent 𝑖 is susceptible i.e. 𝑋 𝑖 (𝑡) = 0, it can receive an infection from the neighbors with a rate 𝜆 0,𝑖 = 𝛽 ∑ 𝑁 𝑗=1 𝑎 𝑖𝑗 1 {𝑋 𝑗 (𝑡)=1} , where 𝛽 is the infection rate. After receiving the infection, the process 𝑋 𝑖 (𝑡) waits for the minimum incubation period 𝜏 0 , then it transition towards the state 1. When 𝑋 𝑖 (𝑡) enters the state 1, it stays for a minimum recovery time 𝜏 1 . After 𝜏 1 , it can transition to state 0 with a curing rate 𝜆 1,𝑖 = 𝛿. Agent 𝑖 is coupled with the rest of the network through the rate 𝜆 0,𝑖 which is a random variable that makes the process doubly-stochastic. Due to the randomness in the value of 𝜆 0,𝑖 , the process {𝑋 𝑖 (𝑡), 𝑡 ≥ 0} is no more semi-Markovian. The random nature of 𝜆 0,𝑖 can be removed by conditioning to all the possible combinations of the states 𝑋 𝑗 (𝑡) = 1 (and 𝑋 𝑗 (𝑡) = 0) which are neighbors of node 𝑖. But we will end up with the 2 𝑁 -state semi-Markov process defined earlier in Section 3. This is a model of a large number of states (2 𝑁 states). Therefore, we approximate 𝜆 0,𝑖 by its expectation λ0,𝑖 = 𝐸{𝜆 0,𝑖 } = 𝛽 ∑ 𝑁 𝑗=1 𝑎 𝑖𝑗 𝑃 𝑟(𝑋 𝑗 (𝑡) = 1) = 𝛽 ∑ 𝑁 𝑗=1 𝑎 𝑖𝑗 𝑝 𝑗 (𝑡), similar to the mean field approximation done in Van Mieghem et al. 7 . As a result, the process {𝑋 𝑖 (𝑡), 𝑡 ≥ 0} is now semi-Markovian with post-delay rates λ0,𝑖 and 𝜆 1,𝑖 . Let 𝑝 𝑖 (𝑡), 𝑖 = 1, … , 𝑁, be the probability of infection of agent 𝑖 at time 𝑡. Hence, using (29), the time-delayed dynamics of 𝑝 𝑖 (𝑡) become:

ẋ𝐼,𝑖 (𝑡) =𝑝 𝑖 (𝑡) ṗ𝑖 (𝑡) =𝛽 𝑁 ∑ 𝑗=1 𝑎 𝑖𝑗 𝑝 𝑗 (𝑡)(𝐻(𝑡 -𝜏 0 ) -𝑝 𝑖 (𝑡)) -𝛿𝑝 𝑖 (𝑡) -𝛽𝛿 𝑁 ∑ 𝑗=1 𝑎 𝑖𝑗 𝑝 𝑗 (𝑡) ( 𝑥 𝐼,𝑖 (𝑡) -𝐻(𝑡 -𝐷)𝑥 𝐼,𝑖 (𝑡 -𝐷) -𝐻(𝑡 -𝜏 0 )(𝑡 -𝜏 0 ) + 𝐻(𝑡 -𝐷)(𝑡 -𝐷) ) + ( 𝛽 𝑁 ∑ 𝑗=1 𝑎 𝑖𝑗 𝑝 𝑗 (𝑡) ( 1 -𝐻(𝑡 -𝜏 0 ) ) + 𝛿 ( 1 -𝐻(𝑡 -𝜏 1 )
)

+ 𝛽𝛿 𝑁 ∑ 𝑗=1 𝑎 𝑖𝑗 𝑝 𝑗 (𝑡) ( 𝑡 -𝐻(𝑡 -𝜏 0 )(𝑡 -𝜏 0 ) -𝐻(𝑡 -𝜏 1 )(𝑡 -𝜏 1 ) + 𝐻(𝑡 -𝐷)(𝑡 -𝐷) ) ) 𝑝 0,𝑖 (31)
where 𝑝 𝑖 (0) = 𝑝 0,𝑖 and 𝑥 𝐼,𝑖 (0) = 0 are the initial conditions. Unlike model ( 18)-( 20), model ( 31) is composed of 𝑁 states.

Remark 3. If 𝜏 0 = 𝜏 1 = 0, the model (31) reduces to

ṗ𝑖 (𝑡) = 𝛽 𝑁 ∑ 𝑗=1 𝑎 𝑖𝑗 𝑝 𝑗 (𝑡)(1 -𝑝 𝑖 (𝑡)) -𝛿𝑝 𝑖 (𝑡), 𝑖 = 1, … , 𝑁 (32)
where 𝑝 𝑖 (0) = 𝑝 0,𝑖 are the initial conditions. The model (32) is the same as the 𝑁-interwined Markov model obtained by the authors in Van Mieghem et al. 7 describing the evolution of SIS epidemics over networks.

STABILITY RESULTS

In practice, it is vital to know under which conditions the epidemic will terminate. For example, it is important to know how many individuals a random individual can meet, or how fast medications and also vaccinations can stop an epidemic. This is crucial to decrease the number of hospitalisations caused by the disease to protect the health system. All these questions are answered by doing a stability analysis of the model, see for instance the work of Giordano et al. [START_REF] Giordano | Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy[END_REF] . In this section, we study the stability of the disease free equilibrium (DFE) 𝑝 𝑒𝑞 𝑖 = 0 for 𝑖 = 1, … , 𝑁, of the 𝑁-intertwined model with time delays given in (31).

The first step is to prove that the probability of infection predicted by the approximated time-delay system (31) remains bounded between 0 and 1 for all times and is really a probability, as in Paré et al. [START_REF] Paré | Epidemic processes over time-varying networks[END_REF] for the delay free case. This is given in the following theorem.

Theorem 1. Consider system (31). If 0 ≤ 𝑝 0,𝑖 ≤ 1 then we have 0 ≤ 𝑝 𝑖 (𝑡) ≤ 1 for all 𝑡 ≥ 0 and 𝑖 = 1, … , 𝑁.

Proof. We assume without any loss of generality that 𝜏 0 ≤ 𝜏 1 . Due to the time delays, the time domain is divided into four intervals.

Interval 1: 0 ≤ 𝑡 < 𝜏 0 . One can notice that 𝑝 𝑖 (𝑡) = 𝑝 0,𝑖 , 𝑖 = 1, … , 𝑁, is a solution of (31) in the interval 0 ≤ 𝑡 < 𝜏 0 . This result is obvious, since wherever the process starts it should wait at least 𝜏 0 amount of time before jumping to the other state, and hence the probability in the interval 0 ≤ 𝑡 < 𝜏 0 remains constant and equal to the initial condition 𝑝 0,𝑖 . Since 0 ≤ 𝑝 0,𝑖 ≤ 1, then we have 0 ≤ 𝑝 𝑖 (𝑡) ≤ 1 for 0 ≤ 𝑡 < 𝜏 0 and 𝑖 = 1, … , 𝑁.

Interval 2: 𝜏 0 ≤ 𝑡 < 𝜏 1 . The dynamics in (31) become:

ṗ𝑖 (𝑡) =𝛽 𝑁 ∑ 𝑗=1 𝑎 𝑖𝑗 𝑝 𝑗 (𝑡)(1 -𝑝 𝑖 (𝑡)) -𝛿𝑝 𝑖 (𝑡) + 𝛿𝑝 0,𝑖 + 𝛽𝛿 𝑁 ∑ 𝑗=1 𝑎 𝑖𝑗 𝑝 𝑗 (𝑡) 𝑡 ∫ 𝜏 0 (1 -𝑝 𝑖 (𝜏))𝑑𝜏 (33)
If 𝑝 0,𝑖 = 1, then by (33) we obtain ṗ𝑖 (𝜏 0 ) = 0, and hence the probability of infection in the interval 𝜏 0 ≤ 𝑡 < 𝜏 1 remains constant and equals to 1. This is also clear from the fact that if the initial start is in the state 1, the process is obliged to stay there for a time 𝜏 1 before it is able to transition. If 𝑝 0,𝑖 = 0 for some 𝑖, and 𝑝 0,𝑗 > 0, for 𝑗 ≠ 𝑖, 𝑖, 𝑗 = 1, … , 𝑁, then by (33) we obtain ṗ𝑖 (𝜏 0 ) = 𝛽 ∑ 𝑁 𝑗=1 𝑎 𝑖𝑗 𝑝 𝑗 (𝜏 0 ) > 0; therefore, 𝑝 𝑖 (𝑡) increases when it touches zero at 𝜏 0 . If 0 ≤ 𝑝 0,𝑖 < 1, 𝑖 = 1, … , 𝑁, let 𝜖 𝑖,0 and 𝜖 𝑖,1 be the two time instants between 𝜏 0 and 𝜏 1 such that 𝑝 𝑖 (𝜖 𝑖,0 ) = 0 and 𝑝 𝑖 (𝜖 𝑖,1 ) = 1 for the first time in the interval 𝜏 0 ≤ 𝑡 < 𝜏 1 . By (33), one obtains:

ṗ𝑖 (𝜖 𝑖,0 ) =𝛽 𝑁 ∑ 𝑗=1 𝑎 𝑖𝑗 𝑝 𝑗 (𝜖 𝑖,0 ) ⎛ ⎜ ⎜ ⎝ 1 + 𝛿 𝜖 𝑖,0 ∫ 𝜏 0 (1 -𝑝 𝑖 (𝜏))𝑑𝜏 ⎞ ⎟ ⎟ ⎠ + 𝛿𝑝 0,𝑖 (34)
ṗ𝑖 (𝜖 𝑖,1 ) = -𝛿 ⎛ ⎜ ⎜ ⎝ (1 -𝑝 𝑖,0 ) -𝛽 𝑁 ∑ 𝑗=1 𝑎 𝑖𝑗 𝑝 𝑗 (𝜖 𝑖,1 ) 𝜖 𝑖,1 ∫ 𝜏 0 (1 -𝑝 𝑖 (𝜏))𝑑𝜏 ⎞ ⎟ ⎟ ⎠ (35)
We should emphasis that the chronological order at which the probability states reach 0 or reach 1 does not affect the proof, as we can always consider that when 𝑝 𝑖 reaches 0 or 1 at a given time, all the other states 𝑝 𝑗 , 𝑖 ≠ 𝑗, do not reach 0 or 1 by that time, then afterwards, the proof can then be repeated similarly. The term between the brackets in (35) is the probability of being infected at the time 𝜖 𝑖,0 , which is equal to the probability to be infected initially 𝑝 𝑖,0 summed with the probability to not be infected in the time interval [ 𝜏 0 , 𝜖 𝑖,0 ] and then to be infected from the neighbors at 𝜖 𝑖,0 . Hence, we have ṗ𝑖 (𝜖 𝑖,0 ) = 0 and therefore, when 𝑝 𝑖 touches 1 in [𝜏 0 , 𝜖 𝑖,0 ], it remains constant since the time delay 𝜏 1 has not yet been reached for the process to be able to transition. Now consider (34), we have ṗ𝑖 (𝜖 𝑖,0 ) ≥ 0 and therefore by the continuity of the solution in the interval [ 𝜏 0 , 𝜖 𝑖,0 ] , 𝑝 𝑖 (𝑡) stays bounded between 0 and 1 in the interval

[ 𝜏 0 , 𝜖 𝑖,0 ]
. The same analysis is repeated on the complete interval [𝜏 0 , 𝜏 1 ]. Interval 3: 𝜏 1 ≤ 𝑡 < 𝐷. The dynamics in (31) read as:

ṗ𝑖 (𝑡) = 𝛽 𝑁 ∑ 𝑗=1 𝑎 𝑖𝑗 𝑝 𝑗 (𝑡)(1 -𝑝 𝑖 (𝑡)) -𝛿𝑝 𝑖 (𝑡) + 𝛽𝛿 𝑁 ∑ 𝑗=1 𝑎 𝑖𝑗 𝑝 𝑗 (𝑡) ( -𝑥 𝐼,𝑖 (𝑡) + 𝑡 -𝜏 0 -𝑝 0,𝑖 (𝑡 -𝐷) ) (36)
If 𝑝 𝑖 (𝜏 1 ) = 0, then by (36) and the boundedness of 𝑝 𝑖 (𝑡) in the interval [0, 𝜏 1 ], we have ṗ𝑖 (𝜏

∫ 𝜏 0 (1 -𝑝 𝑖 (𝜏))𝑑𝜏 ⎞ ⎟ ⎟ ⎠ (38)
The last term in (38) is the probability of not being infected in the interval [𝜏 0 , 𝜖 𝑖,1 ] and then to be infected from the neighbors at time 𝜖 𝑖,1 . Because 𝑝 𝑖 (𝑡) is a probability in [0,𝜖 𝑖,1 ], the last term in (38) is less than 1 and therefore ṗ𝑖 (𝜖 𝑖,1 ) ≤ 0. Consider now equation (37). Since 𝑝 𝑖 (𝑡) stays less than 1 before the states touch zero, then we have ṗ𝑖 (𝜖 𝑖,0 ) ≥ 0. By repeating the same procedure over the time instants at which the solution touch 0 or 1 and by using the continuity of the dynamics, we obtain that 𝑝 𝑖 (𝑡) remains bounded between 0 and 1 in the interval [𝜏 1 , 𝐷]. Interval 4: 𝑡 ≥ 𝐷. The system has the following dynamics:

ṗ𝑖 (𝑡) = 𝛽 𝑁 ∑ 𝑗=1 𝑎 𝑖𝑗 𝑝 𝑗 (𝑡)(1 -𝑝 𝑖 (𝑡)) -𝛿𝑝 𝑖 (𝑡) + 𝛽𝛿 𝑁 ∑ 𝑗=1 𝑎 𝑖𝑗 𝑝 𝑗 (𝑡) ( 𝜏 1 -𝑥 𝐼,𝑖 (𝑡) + 𝑥 𝐼,𝑖 (𝑡 -𝐷) ) (39)
We argue as before. Let 𝜖 𝑖,0 ≥ 𝐷 and 𝜖 𝑖,1 ≥ 𝐷 be the two time instants such that 𝑝 𝑖 (𝜖 𝑖,0 ) = 0 and 𝑝 𝑖 (𝜖 𝑖,1 ) = 1 for the first time in the interval 𝑡 ≥ 𝐷. By (39) we have:

ṗ𝑖 (𝜖 𝑖,0 ) = 𝛽 𝑁 ∑ 𝑗=1 𝑎 𝑖𝑗 𝑝 𝑗 (𝜖 𝑖,0 ) ⎛ ⎜ ⎜ ⎜ ⎝ 1 + 𝛿𝜏 1 -𝛿 𝜖 𝑖,0 ∫ 𝜖 𝑖,0 -𝐷 𝑝 𝑖 (𝜏)𝑑𝜏 ⎞ ⎟ ⎟ ⎟ ⎠ (40) ṗ𝑖 (𝜖 𝑖,1 ) = -𝛿 ⎛ ⎜ ⎜ ⎜ ⎝ 1 + 𝛽𝜏 0 𝑁 ∑ 𝑗=1 𝑎 𝑖𝑗 𝑝 𝑗 (𝜖 𝑖,1 ) -𝛽 𝑁 ∑ 𝑗=1 𝑎 𝑖𝑗 𝑝 𝑗 (𝜖 𝑖,1 ) 𝜖 𝑖,1 ∫ 𝜖 𝑖,1 -𝐷 (1 -𝑝 𝑖 (𝜏))𝑑𝜏 ⎞ ⎟ ⎟ ⎟ ⎠ (41)
One can notice that the last term in (40) is the probability to be infected in the interval [𝜖 𝑖,0 -𝐷,𝜖 𝑖,0 ] and then to be cured at time 𝜖 𝑖,0 . Since 𝑝 𝑖 (𝑡) is a probability in the time interval [0, 𝜖 𝑖,0 ], then the last term in (40) is less than 1 and hence ṗ𝑖 (𝜖 𝑖,0 ) ≥ 0. Doing the same analysis for the last term in (41), one can directly deduce that ṗ𝑖 (𝜖 𝑖,1 ) ≤ 0. Repeating the same argument over all the time instants at which the solution intersects 0 or 1 and by using the continuity of the dynamics, we deduce that 0 ≤ 𝑝 𝑖 (𝑡) ≤ 1 for 𝑡 ≥ 𝐷 and the proof is complete.

The second step is to compute the steady state points of the 𝑁-intertwined model (31). The interval of interest is Interval 4 for 𝑡 ≥ 𝐷. We begin by writing (39) in the matrix form:

Ẋ𝐼 (𝑡) = 𝑝(𝑡) ṗ(𝑡) = ( 𝛽(1 + 𝛿𝜏 1 )𝐴 -𝛽𝑃 𝑑 (𝑡)𝐴 -𝛿𝐈 𝐍×𝐍 -𝛽𝛿 ( 𝑋 𝑑 (𝑡) -𝑋 𝑑 (𝑡 -𝐷) ) 𝐴 ) 𝑝(𝑡) ( 42
)
where

𝑝(𝑡) = [ 𝑝 1 (𝑡), … , 𝑝 𝑁 (𝑡) ] 𝑇 , 𝑋 𝐼 (𝑡) = [ 𝑥 𝐼,1 (𝑡), … , 𝑥 𝐼,𝑁 (𝑡)
] 𝑇 , 𝑃 𝑑 (𝑡) = 𝑑𝑖𝑎𝑔(𝑝(𝑡)) and 𝑋 𝑑 (𝑡) = 𝑑𝑖𝑎𝑔(𝑋 𝐼 (𝑡)). In the steady state when 𝑡 goes to +∞, one can obtain the following algebraic steady state equation from (42):

0 = 𝛽(1 + 𝛿𝜏 1 ) 𝑁 ∑ 𝑗=1 𝑎 𝑖𝑗 𝑝 ∞ 𝑗 -𝛽𝑝 ∞ 𝑖 𝑁 ∑ 𝑗=1 𝑎 𝑖𝑗 𝑝 ∞ 𝑗 -𝛿𝑝 ∞ 𝑖 ( 43
)
where

𝑝 ∞ = [ 𝑝 ∞ 1 , … , 𝑝 ∞ 𝑁
] 𝑇 is the steady state vector of the infection probabilities. We can rewrite (43) in a more compact way as:

𝑝 ∞ 𝑖 = (1 + 𝛿𝜏 1 ) ( 1 - 1 1 + 𝛽 𝛿 ∑ 𝑁 𝑗=1 𝑎 𝑖𝑗 𝑝 ∞ 𝑗 ) (44)
Obviously from (44), 𝑝 ∞ 𝑖 = 0 for all 𝑖 is a trivial solution. If 𝛿 = 0, then 𝑝 ∞ 𝑖 = 1 for all 𝑖, which is also obvious from the fact that eventually all the nodes will become infected if the curing rate is zero. Now for 𝛽 ≠ 0 and 𝛿 ≠ 0 and if the graph is connected, we have that either all the 𝑝 ∞ 𝑖 = 0 for all 𝑖 or non of the 𝑝 ∞ 𝑖 is zero. Since (44) is a nonlinear equation, this gives rise to a second solution called the meta-stable solution, i.e. to have a long live epidemic. An idea on the magnitude of this solution is given in the following lemma: Lemma 1. If 𝛽∕𝛿 ≥ 0, then the non-zero solution 𝑝 ∞ 𝑖 can be expressed as a continued fraction

𝑝 ∞ 𝑖 = (1 + 𝛿𝜏 1 ) ⎛ ⎜ ⎜ ⎜ ⎝ 1 - 1 1 + 𝛽 𝛿 (1 + 𝛿𝜏 1 )𝑑 𝑖 -𝛽 𝛿 ∑ 𝑁 𝑗=1 𝑎 𝑖𝑗 1+ 𝛽 𝛿 (1+𝛿𝜏 1 )𝑑 𝑗 -𝛽 𝛿 ∑ 𝑁 𝑘=1 𝑎 𝑗𝑘 1+ 𝛽 𝛿 (1+𝛿𝜏 1 )𝑑 𝑘 - 𝛽 𝛿 ∑ 𝑁 𝑙=1 𝑎 𝑘𝑙 ... ⎞ ⎟ ⎟ ⎟ ⎠ ( 45
)
where 𝑑 𝑖 = ∑ 𝑗 𝑎 𝑖𝑗 is the degree of node 𝑖. Furthermore, 𝑝 ∞ 𝑖 is bounded as follows:

0 ≤ 𝑝 ∞ 𝑖 ≤ min { (1 + 𝛿𝜏 1 ) ( 1 - 1 1 + 𝛽 𝛿 𝑑 𝑖 ) , 1 } (46)
Proof. Consider equation ( 44). We have that:

𝑝 ∞ 𝑖 = (1 + 𝛿𝜏 1 ) ( 1 - 1 1 + 𝛽 𝛿 ∑ 𝑁 𝑗=1 𝑎 𝑖𝑗 𝑝 ∞ 𝑗 ) = (1 + 𝛿𝜏 1 ) ( 1 - 1 1 + 𝛽 𝛿 𝑑 𝑖 -𝛽 𝛿 ∑ 𝑁 𝑗=1 𝑎 𝑖𝑗 (1 -𝑝 ∞ 𝑗 ) ) ≤ (1 + 𝛿𝜏 1 ) ( 1 - 1 1 + 𝛽 𝛿 𝑑 𝑖 ) Since 𝛽 𝛿 ∑ 𝑁 𝑗=1 𝑎 𝑖𝑗 (1 -𝑝 ∞ 𝑗 ) ≥ 0 as a consequence of 0 ≤ 𝑝 ∞ 𝑗 ≤ 1 for all 𝑗 proved in Theorem 1.
The proof of (46) is complete. Furthermore, by doing consecutive iterations of (44), one can find directly the infinite continued fraction expansion in (45).

We should remark that the calculation of the steady state solution is similar to Section 4-B of Van Mieghem et al. 7 . Hence, we will not proceed to the proof of existence of the steady state as it is exactly similar to Van Mieghem et al. 7 . However, it is useful to comment on the upper bound of the steady state probability given in (46). The upper bound depends on the degree of connection of each node 𝑑 𝑖 . If the node 𝑖 is isolated, i.e. 𝑑 𝑖 = 0, then individual 𝑖 has no neighbors and cannot be infected (𝑝 ∞ 𝑖 = 0). Otherwise, if a node has many connections then it is more likely to be infected. Similarly, if the rate of infection 𝛽 is small then the steady state probability is small, otherwise the probability of infection increases with the increase of 𝛽. The bound is linear in 𝜏 1 , so if 𝜏 1 is very large (i.e. 𝛿 is very small) then the upper bound is no more significant, this is clear from the fact that 𝜏 1 is the time spent in the infected state for any individual. If this time is large, the probability of infection will be large.

The third step is to analyse the stability of the DFE 𝑝 𝑒𝑞 𝑖 = 0, 𝑖 = 1, … , 𝑁 for the 𝑁-intertwined model (31). This is given in the following theorem: Theorem 2. Consider system (42) for 𝑡 ≥ 𝐷. If there exists 𝛾 > 0 such that

𝜆 𝑚𝑎𝑥 (𝛽(1 + 𝛿𝜏 1 )𝐴 -𝛿𝐈 𝐍×𝐍 ) ≤ - 𝛾 2 (47)
where 𝜆 𝑚𝑎𝑥 denotes the largest eigenvalue, then for all 0 ≤ 𝑝 𝑖,0 ≤ 1, the DFE of model (31) will converge to zero exponentially, i.e. the DFE is globally exponentially stable (GES).

Proof. We define the following Lyapunov function for 𝑡 ≥ 𝐷:

𝑉 (𝑡) = 1 2 𝑝 𝑇 (𝑡)𝑝(𝑡). (48)
Differentiating (48) with respect to time, we get:

V (𝑡) = 𝑝 𝑇 (𝑡) ( 𝛽(1 + 𝛿𝜏 1 )𝐴 -𝛽𝑃 𝑑 (𝑡)𝐴 -𝛿𝐈 𝐍×𝐍 -𝛽𝛿 ( 𝑋 𝑑 (𝑡) -𝑋 𝑑 (𝑡 -𝐷) ) 𝐴 ) 𝑝(𝑡). ( 49
)
)
𝑖𝑗 ≥ 0 for all 𝑖, 𝑗 = 1, … , 𝑁. Therefore, we can bound the derivative of the Lyapunov function as follows:

V (𝑡) ≤ 𝑝 𝑇 (𝑡) ( 𝛽(1 + 𝛿𝜏 1 )𝐴 -𝛿𝐼 𝑁×𝑁 ) 𝑝(𝑡). ( 50
)
Because the graph is undirected, i.e. 𝐴 = 𝐴 𝑇 , then the matrix 𝛽(1+𝛿𝜏 1 )𝐴-𝛿𝐼 𝑁×𝑁 is symmetric. Now, we use the Rayleigh-Ritz quotient to simplify more (50) and obtain:

V (𝑡) ≤ 𝜆 𝑚𝑎𝑥 ( 𝛽(1 + 𝛿𝜏 1 )𝐴 -𝛿𝐼 𝑁×𝑁 ) ||𝑝|| 2 . ( 51
)
If there exists 𝛾 > 0 such that condition (47) is fulfilled, then by (51), we have that V (𝑡) ≤ -𝛾𝑉 (𝑡) and this gives the exponential decay of the Lyapunov function for 𝑡 ≥ 𝐷 and the proof is complete. < 𝜏 𝑐 = 1 𝜆 𝑚𝑎𝑥 (𝐴) . If 𝜏 1 = 0, we require 𝜆 𝑚𝑎𝑥 (𝐴) < 𝛿 𝛽 and this leads to the stability results in Condition 1 in the article of Ganesh et al. [START_REF] Ganesh | The effect of network topology on the spread of epidemics[END_REF] , Theorem 1 in the article of Paré et al. [START_REF] Paré | Epidemic processes over time-varying networks[END_REF] and Corollary 1 in the article of Wang et al. [START_REF] Wang | Epidemic spreading in real networks: An eigenvalue viewpoint[END_REF] . Equation (47) means that 𝜏 0 does not affect the long-term stability. 𝜏 0 is the incubation time of the virus. Since we assume static links between agents, i.e. fixed topology, we expect that an infected agent will eventually infect its neighbors after the incubation period 𝜏 0 is passed. Therefore, 𝜏 0 must not affect the steady state infection probability and this is what appears in condition (47). However, equation (47) shows that 𝜏 1 is important for stabilizing the epidemics. In practice, it means that having drugs decreasing the time delay 𝜏 1 is decisive in putting an end to epidemics rather than the intrinsic characteristics 𝜏 0 of viruses in the social bubble strategy inducing a graph having a fixed topology. To know how large is the conservativeness of (47), one way is to find if condition (47) is also a necessary condition. This is difficult as we cannot obtain analytically the solution of (42). But we can do a linearization approximation in neglecting the nonlinear terms in (42), the terms in 𝑝 2 (𝑡), as done in Wang et al. [START_REF] Wang | Epidemic spreading in real networks: An eigenvalue viewpoint[END_REF] . This reduces (42) to:

ṗ(𝑡) = (𝛽(1 + 𝛿𝜏 1 )𝐴 -𝛿𝐼 𝑁×𝑁 )𝑝(𝑡) (52)
From (52), one can notice that condition (47) can also be a necessary condition for the infection probability to stabilize at zero. If the nonlinearity do not play an important role, the marge of conservativneness of (47) is not large.

Remark 5. Theorem 2 specifies the stability conditions of the DFE of the 𝑁-intertwined model (31). For the exact model ( 18)-( 20), if the graph is connected, the DFE (𝑝 𝑒𝑞 𝑖 = 0) is always stable even for large delays. In other words, the exact model ( 18)-( 20) will always converge to zero if the delays are not infinite and if the graph is connected. This is because that the disease-free state is an absorbing state for the semi-Markov process 𝑍(𝑡) modeling the SIS infection in Section 3. To clear up the point, let us consider for example a population of three agents. As we have explained in Section 3, an agent can be either infected I or susceptible S. Therefore at any time, this population of three agents is found in one of the following 8 states: 𝑆𝑆𝑆, 𝐼𝑆𝑆, 𝑆𝐼𝑆, 𝑆𝑆𝐼, 𝐼𝐼𝑆, 𝐼𝑆𝐼, 𝑆𝐼𝐼 and 𝐼𝐼𝐼. For example, state 𝐼𝑆𝑆 means that agent 1 is susceptible, agent 2 is susceptible and agent 3 is infected. The flow diagram of this semi-Markov chain is shown on Figure 2. The green arrows represent the recovery paths while the red arrows represent the infection paths. One can directly notice that once the process reaches the state SSS, it cannot go out and it remains there forever. When 𝑡 → +∞, and since their is a probability to arrive at 𝑆𝑆𝑆 (𝛿 ≠ 0), eventually the process will reach 𝑆𝑆𝑆 once and it will remain there. Hence, the state 𝑆𝑆𝑆 is an absorbing state and the exact model ( 18)- (20) will always converge to it at infinity. If the time delays are big, the process will be slow in reaching 𝑆𝑆𝑆. In theory, if for example 𝜏 1 → +∞, the SIS model reduces to the SI model. In an SI model, once an agent becomes infected it will remain infected forever and there is no curing rate 𝛿 (𝛿 = 0). If we retake the 𝑁 = 3 agents example in Figure 2, the state 𝐼𝐼𝐼 becomes an absorbing state for the SI model. But for the SIS model, there is no practical interest in considering large delays because each virus has a finite incubation and a finite recovery period. To recap, the exact SIS model will converge to zero at the steady state, i.e. when 𝑡 ←→ +∞, because zero is the absorbing state of the semi-Markov process. However, for some large enough effective spreading rate 𝜏 and for large number of agents 𝑁, the exact SIS model can settle in the so called meta-stable state. This is the state where the process remains constant for a long period before converging to zero when time goes to infinity. The meta-stable state is not yet clearly defined for the exact SIS processes. In the delay-free case, the authors Van Mieghem et al. [START_REF] Van Mieghem | Epidemics in networks with nodal self-infection and the epidemic threshold[END_REF] have proposed that this meta-stable state of the exact SIS process is the steady state of the so called 𝜖-SIS model. For more details on the 𝜖-SIS model, we orient the readers to the work of Van Mieghem et al. [START_REF] Van Mieghem | Epidemics in networks with nodal self-infection and the epidemic threshold[END_REF] and Li Cong et al. [START_REF] Li | Susceptible-infected-susceptible model: A comparison of N-intertwined and heterogeneous mean-field approximations[END_REF] . For our case with time delays, we will give an insight on this meta-stable state of the exact SIS process using numerical simulations while leaving the theoretical developments for future works. Because the 𝑁-intertwined model approximates the exact SIS model, the non-zero steady state of the 𝑁-intertwined model defined in (45) approximates the meta-stable state of the exact SIS process. Hence, the interest behind studying the stability of the 𝑁-intertwined model in Theorem 2 is because this model approximates the exact model. It has 𝑁 states instead of 2 𝑁 which makes it computationally much lighter. It is extremely important to know the convergence conditions of the 𝑁-intertwined model in order to evaluate the limits of the mean field approximation.

NUMERICAL SIMULATIONS

Delayed SIS model versus SEIRS model

The first scenario that we would like to test through numerical simulations is the advantage of injecting two time delays 𝜏 0 and 𝜏 1 into the SIS model instead of simply replacing them by two extra states 𝐸 and 𝑅 to have the SEIRS model. This is legitimate because one can think of 𝜏 0 as the exposed state E before becoming infected and 𝜏 1 as the recovery state R before returning back susceptible. Theoretically speaking, an ODE state cannot replace the endogenous time delay. The time delay corresponds to an infinite state which is the output of a transport partial differential equation (PDE), see the work of the authors Krstic et al. [START_REF] Krstic | Backstepping boundary control for first-order hyperbolic PDEs and application to systems with actuator and sensor delays[END_REF] . However, the question that remains is: can the delay be approximately compensated by adding an extra state? To answer this

) + 𝛿𝑟 𝑖 (𝑡) ė𝑖 (𝑡) = 𝛽 ( 𝑁 ∑ 𝑗=1 𝑎 𝑖𝑗 𝑒 𝑗 (𝑡) + 𝑎 𝑖𝑗 𝑝 𝑗 (𝑡) ) -𝜎𝑒 𝑖 (𝑡) ṗ𝑖 (𝑡) = 𝜎𝑒 𝑖 (𝑡) -𝛾𝑝 𝑖 (𝑡) ̇𝑟𝑖 (𝑡) = 𝛾𝑝 𝑖 (𝑡) -𝛿𝑟 𝑖 (𝑡) (53)
where 𝑠 𝑖 (𝑡) is the probability that agent 𝑖 is susceptible at time 𝑡, 𝑒 𝑖 (𝑡) is the probability that agent 𝑖 is exposed to infection at time 𝑡, 𝑝 𝑖 (𝑡) is the probability that agent 𝑖 is infected at time 𝑡 and 𝑟 𝑖 (𝑡) is the probability that agent 𝑖 is recovered at time 𝑡. The goal of the evaluation is to check if adding two extra states E and R to the SIS model can compensate the effect of the two time delays 𝜏 0 and 𝜏 1 . We simulate the models in (31) and (53) for the same parameters: 𝛽 = 0.001, 𝛿 = 0.05 and for the same initial conditions. But we vary the value of the incubation delay 𝜏 0 and the minumum time to recovery 𝜏 1 . We take two values of 𝜏 0 = {20, 50} and two values for 𝜏 1 = {50, 100}. The rate of transition from E to I is computed as 𝜎 = 1∕𝜏 0 , and the rate of transition from I to R is computed as 𝛾 = 1∕𝜏 1 . Figure 3 shows the comparison between the delayed SIS model and the SEIRS model. The network is a scale-free network composed of fifty agents 𝑁 = 50, constructed using the Barabasi and Albert algorithm. One can directly observe on the two upper plots of Figure 3 that for relatively small time delays, the SEIRS model (53) approximate the delayed SIS model (31). When the time delays are relatively large on the two bottom plots of Figure 3, the SEIRS fails to approximate the delayed SIS model. Hence, one can deduce that adding a multi-compartmental states to compensate the effect of the delay is reasonable if the delay is small. However, when the time delay is large, this approximation is no more useful. This decreases the predictability of the model. Therefore, it is of great interest to work on the time-delay system directly instead of adding extra states. The multi-agent model that we have presented in the paper is of general setting. It works for small and large delays depending on the intrinsic nature of the disease or of the virus.

Validity of the mean field approximation

The second scenario is dedicated to the analysis of the effect of the time delays on the validity of the approximation of the exact 2 𝑁 -state semi-Markov model by the 𝑁-intertwined model of Section 4. The exact model ( 18)-( 20) of 2 𝑁 states is difficult to simulate directly when the number of agents is large. So we distinguish between two cases: the first case is for small 𝑁 ≤ 9, where we simulate the exact SIS process using equations ( 18)- (20), and the second case is for large 𝑁 where we simulate the SIS exact process by performing realizations/trials. We will explain later the simulation algorithm of these trials. Now, we start the case of small number of agents.

The case of small number of agents 𝑁 ≤ 9

We are able to simulate the exact model ( 18)-( 20) for a maximum number of agents 𝑁 = 9. Recall that when 𝜏 0 = 𝜏 1 = 0, our models in ( 18)-( 20) and (31) coincide with the models of Van Mieghem et al. 7 and hence the same explanations follow, see Section VIII in Van Mieghem et al. 7 . If 𝜏 0 ≠ 0 and 𝜏 1 ≠ 0, we have explained in Remark 4 that if 𝜏 = 𝛽(1+𝛿𝜏 1 ) 𝛿 < 𝜏 𝑐 = 1 𝜆 𝑚𝑎𝑥 (𝐴) , the 𝑁-intertwined model will converge to the DFE. Since 𝜏 0 does not affect the long term stability, we have simulated the models of ( 18)-( 20) and (31) for the same value of 𝜏 0 but different values of 𝜏 1 . We first plot the time response of the exact SIS process and the 𝑁-intertwined model on the complete graph (𝐾 9 ) and the scale-free graph with 9 agents. Starting with the homogeneous case, i.e. the graph is complete and the number of agents is 𝑁 = 9 (𝐾 9 graph). Figure 4 shows the variation of the density of the infected nodes 𝑦(𝑡) = 1 𝑁 ∑ 𝑖=𝑁 𝑖=0 𝑝 𝑖 (𝑡) as a function of time. Note that we have replicated the simulation scenario done in Section VIII in Van Mieghem et al. 7 The evaluation is done for four values of 𝜏 1 ∈ {5, 10, 20, 25}. The rate of infection 𝛽 is set to 𝛽 = 0.01 and the curing rate 𝛿 is set to 𝛿 = 0.5. The first up-left plot of Figure 4 corresponds to 𝜏 1 = 5, and 𝜏 = 0.07 << 𝜏 𝑐 = 0.125. We can observe that the 𝑁-intertwined model is stable, as predicted by Theorem 2, and the approximation of the exact model is accurate. Moving to the bottom-left plot of Figure 4 that corresponds to a value 𝜏 1 = 10 and 𝜏 = 0.12 < 𝜏 𝑐 = 0.125, the 𝑁-intertwined model stays stable, as predicted by Theorem 2, but the approximation of the exact model starts to become less accurate. Increasing more the value of the time delay 𝜏 1 to 𝜏 1 = 20 in the up-right plot of Figure 4, one can notice that the 𝑁-intertwined model differs significantly from the exact model. Stability is not attained by the 𝑁-intertwined model and this confirms the results predicted by Theorem 2 because in this case 𝜏 = 0.22 > 𝜏 𝑐 = 0.125. For large 𝜏 1 = 25, the 𝑁-intertwined model diverges and the aprroximation is not accurate. Now, we consider the setting on a more general heterogeneous graphs. We simulate the exact and the 𝑁-intertwined models on a two scale-free graphs with 9 agents. The first scale-free graph, we call it "scale-free-3" graph, has an average number of connections equal to three while the second has an average number of connections equal to seven and we call it "scale-free-7" graph. Both graphs are constructed using the Barabasi and Albert algorithm. We perform the same simulation scenario as before, with the same parameters 𝛽 = 0.01, 𝛿 = 0.5 and 𝜏 0 = 5 but for different values of 𝜏 1 . Figure 5 shows the evolution The 𝑁-intertwined model is stable and converges to zero since 𝜏 = 0.22 < 𝜏 𝑐 = 0.306. Increasing 𝜏 1 to 𝜏 1 = 50, we see on the right plot of Figure 5 that the 𝑁-interwined model diverge from the exact model and it is not stable since 𝜏 = 0.52 > 𝜏 𝑐 = 0.306. The same analysis apply to "scale-free-7" graph on Figure 6. However, one can observe that the mean field approximation become worst on the "scale-free-7" graph for small values of the delay 𝜏 1 . For 𝜏 1 = 20, the 𝑁-intertwined model is stable on the "scalefree-3" graph and not stable on the "scale-free-7" graph. Also note that in all the cases, the exact SIS model is converging to zero at the steady state. This is expected because as we have explained in Remark 5 that the healthy state is the absorbing state of the exact model, and since the number of agents is very small, the model reaches its equilibrium very fast and the meta-stable state is not observed. The previous analysis was done only for two types of graphs and several values of 𝜏 1 . To conclude the analysis, we now plot the mean error 𝜖(𝜏) as a function of 𝜏 = 𝛽(1

+ 𝛿𝜏 𝜖(𝜏) = 𝑇 0 | 𝑦 𝑒𝑥𝑎𝑐𝑡 (𝑡, 𝜏) -𝑦 𝑁-𝑖𝑛𝑡𝑒𝑟𝑡𝑤𝑖𝑛𝑒𝑑 (𝑡, 𝜏) | | 𝑑𝑡 (54)
where 𝜖(𝜏) is the 𝐿 1 -norm of the mean field approximation error and 𝑇 is the simulation time. 𝜖(𝜏) is the transient error and not the steady state error. So we choose a period 𝑇 = 300 much larger than the time delays. We fix 𝛽 = 0.01, 𝛿 = 0.5 and we increase 𝜏 by increasing the value of 𝜏 1 . 𝜏 0 is kept equal to zero since the stability condition (47) is not related to 𝜏 0 . The evaluation is done on 6 different graphs: complete graph, line graph, star graph, bipartite graph, Erdos-Renyi graph and the Barabsi-Alber graph [START_REF] Kermack | Contributions to the mathematical theory of epidemics. II.-The problem of endemicity[END_REF] . The results are plotted on Figure 7. We can conclude from Figure 7 three important points. The first one is that the mean field approximation is accurate for all types of graphs if 𝜏 << 𝜏 𝑐 . As we can see on Figure 7, the error starts to increase after crossing the red line corresponding to the value of 𝜏 𝑐 for all the graphs. Van Mieghem et al. 7 have also arrived to the same conclusion. The second point is that the error is increasing with 𝜏. We observe for the delay-free case that the approximation error decrease with 𝜏. But for the model with delays and for small number of agents, the error is increasing with 𝜏. The third point is that the time delay 𝜏 1 have significant effect on highly connected topologies. Remember that 𝜏 1 is implicit in the value of 𝜏 = 𝛽(1 + 𝛿𝜏 1 )∕𝛿. We have calculated different values of 𝜏 based on different values of 𝜏 1 while keeping 𝛽 and 𝛿 fixed. It is also helpful to recall that for connected graphs, the complete graph represents maximum connectivity while the line graph represents minimum connectivity. We can see that the approximation on the complete graph deteriorates for 𝜏 = 0.125 (corresponding to 𝜏 1 = 10 ) while it deteriorates for 𝜏 = 0.52 (corresponding to 𝜏 1 = 50) on the line graph. This conclusion is [START_REF] Kermack | Contributions to the mathematical theory of epidemics. II.-The problem of endemicity[END_REF] The bipartite graph is constructed using two sets of nodes. The first set contains 𝑁 3 nodes and the second set contains 2𝑁 3 nodes. The Erdos-Renyi graph is constructed using a link probability of 𝑝 = 2𝑝 𝑐 = 2 𝑙𝑛𝑁 𝑁 to ensure high probability of graph connectivity. The scale free graph is constructed using Barabsi-Albert algorithm and we start by 2 connected nodes (𝑚 = 2) as the seed. quite reasonable since in highly connected graphs, a node has many connections and when infected it will stay for a time 𝜏 1 in the infected state and eventually infect all the neighbors.

The case of large number of agents 𝑁

The goal of this section is to study the effect of time delays on the mean field approximation for graphs with large number of agents. The idea is to plot the transient error 𝜖(𝜏) in (54) as a function of 𝜏 on different graphs and for different number of agents 𝑁. We cannot simulate the exact equations ( 18)-( 20) for large 𝑁 as the number of states grows in the order of 2 𝑁 . However, we can simulate the exact SIS process for large number of agents 𝑁 by performing several realizations/trials of the process, in the same way as for the Figure 1 made by Li Cong et al. [START_REF] Li | Susceptible-infected-susceptible model: A comparison of N-intertwined and heterogeneous mean-field approximations[END_REF] , and then to take the average of these trials to get an approximation of the exact solution. First, we detail the simulation algorithm for exact SIS process with large 𝑁, and after that we plot the mean field approximation error 𝜖(𝜏).

Consider on Figure 1 a sample path of the stochastic process 𝑋 𝑖 (𝑡) attached to agent 𝑖. At any time 𝑡, the agent 𝑖 can be susceptible 𝑋 𝑖 (𝑡) = 𝑆 or infected 𝑋 𝑖 (𝑡) = 𝐼. According to our model, if agent 𝑖 is susceptible, the time for an arrival of an infection is determined by the rate 𝛽 ∑ 𝑗∈𝑁 𝑖 𝑎 𝑖𝑗 𝑋 𝑗 (𝑡) where 𝛽 is the infection rate, 𝑁 𝑖 is the set of neighbors of agent 𝑖 and 𝑎 𝑖𝑗 are the entities of the adjacency matrix 𝐴. Upon arrival of the infection, the agent waits for a time 𝜏 0 in the susceptible state due to the incubation period of the virus. If agent 𝑖 is infected, he has to stay for a time 𝜏 1 in the infected state, after which he can transition with a rate 𝛿. We fix a simulation time 𝑇 , and we simulate the processes 𝑋 𝑖 (𝑡), 𝑖 = 1, 2, … , 𝑁. The simulation is done as follows: we start with only one infected agent. For any agent in the susceptible state, determine if he can transition using the rate 𝛽 ∑ 𝑗∈𝑁 𝑖 𝑎 𝑖𝑗 𝑋 𝑗 (𝑡). If he can, wait for the incubation period 𝜏 0 , then do the transition. For any agent in the infected state, wait a time 𝜏 1 in the state 𝐼. Then determine if he can transition using the rate 𝛿. Perform this algorithm at each small interval of time 𝑑𝑡 between [0, 𝑇 ] and one gets a single realization of the SIS exact process with time delays. We then repeat the algorithm for a big number of trials 𝑁 𝑡𝑟𝑖𝑎𝑙𝑠 and we take the average to have

𝑋 𝑎𝑣𝑔 𝑖 (𝑡) = 1 𝑁 𝑡𝑟𝑖𝑎𝑙𝑠 𝑁 𝑡𝑟𝑖𝑎𝑙𝑠 ∑ 𝑗=1 𝑋 𝑗 𝑖 (𝑡). (55)
Since the trials are independent, 𝑋 𝑎𝑣𝑔 𝑖 (𝑡) will converge to 𝑋 𝑖 (𝑡) as 𝑁 𝑡𝑟𝑖𝑎𝑙𝑠 goes to +∞. Practically and to have a better understanding of the simulation algorithm, consider for example simulating the exact SIS process on the complete graph of 100 agents 𝑁 = 100 on Figure 8. We have chosen the simulation time 𝑇 = 300, 𝛽 = 0.005 and 𝛿 = 0.6. We have run the latter algorithm for 𝑁 𝑡𝑟𝑖𝑎𝑙𝑠 = 200 and plotted the density of infected nodes 𝑦(𝑡) = 1 𝑁 ∑ 𝑁 𝑖 𝑋 𝑎𝑣𝑔 𝑖 (𝑡) as a function of time. The left plot of Figure 8 is 𝑦(𝑡) without time delays 𝜏 0 = 𝜏 1 = 0 while the right plot is 𝑦(𝑡) with time delays 𝜏 0 = 0 and 𝜏 1 = 2. For the left plot, we have 𝜏 = 𝛽 𝛿 (1 + 𝛿 * 0) = 0.0083 < 𝜏 𝑐 = 1 𝑁-1 = 0.01. The exact SIS process converges to zero which is the healthy state. When 𝜏 is increased by adding the time delay 𝜏 1 = 2, corresponding to 𝜏 = 0.015 > 𝜏 𝑐 = 0.01, the exact SIS process converges to its meta-stable state as shown on the right plot of Figure 8. This is expected and confirmed by the authors in Li Cong et al. [START_REF] Li | Susceptible-infected-susceptible model: A comparison of N-intertwined and heterogeneous mean-field approximations[END_REF] because for high enough effective spreading rate 𝜏 and for high number of agents, the exact SIS process is expected to converge to the meta-stable state as it will take a long time to converge to zero in this case. Now, we plot the mean field approximation error as a function of 𝜏 on Figure 9. We follow the same simulation scenario as for the case of small number of agents, i.e. choose a period 𝑇 = 300 much larger than the time delays and fix 𝛽 = 0.01, 𝛿 = 0.5. We increase 𝜏 by increasing the value of 𝜏 1 and 𝜏 0 is kept equal to zero. The evaluation is done on 6 different graphs: complete graph, line graph, star graph, bipartite graph, Erdos-Renyi graph and the Barabsi-Alber graph [START_REF] Pastor-Satorras | Evolution and structure of the Internet: A statistical physics approach[END_REF] . The simulations show that, in the presence of time delays, the approximation becomes generally better as 𝑁 increases and this is the same conclusion as for the delay-free case. The error is lowest on the black lines corresponding to 𝑁 = 300 while it increases when 𝑁 is decreased on the blue and red lines corresponding to 𝑁 = 150 and 𝑁 = 50 respectively. But we note that this behaviour is not so clear on the Barabasi-Albert graph. With respect to 𝜏, if 𝜏 < 𝜏 𝑐 , the approximation is also very accurate [START_REF] Barrat | Dynamical processes on complex networks[END_REF] . If 𝜏 > 𝜏 𝑐 then the error decreases with 𝜏 for the complete graph, star graph, bipartite graph and Erdos-Renyi graph while it increases with 𝜏 on the line graph and on the Barabsi Albert graph. The increase on the line graph is fairly expected and this was the conclusion of the authors in Li Cong et al. [START_REF] Li | Susceptible-infected-susceptible model: A comparison of N-intertwined and heterogeneous mean-field approximations[END_REF] and Van Miegham et al. 7 that the approximations are worst on the line graphs because the node degree does not increase with 𝑁. But surprisingly, the error is also increasing with 𝜏 on the Barabsi-Albert graph. We will pay a special attention to the spread of epidemics with delays on the scale-free graphs in our future works.

Remark 4 .

 4 Condition (47) is equivalent to require 𝜏 = 𝛽(1+𝛿𝜏 1 ) 𝛿

FIGURE 2

 2 FIGURE 2Semi-Markov chain for three agents..

4 FIGURE 3

 43 FIGURE 3 Comparison between the delayed SIS model (31) and the SEIRS model (53). 𝑦 𝑆 (𝑡) = 1 𝑁 ∑ 𝑁 𝑖=1 𝑠 𝑖 (𝑡) is the density of sustepible nodes at time 𝑡 and 𝑦 𝐼 (𝑡) = 1 𝑁 ∑ 𝑁 𝑖=1 𝑝 𝑖 (𝑡) is the density of infected nodes at time 𝑡.

6 FIGURE 4

 64 FIGURE 4 Evolution of the density of infected nodes 𝑦(𝑡) with time for 𝑁 = 9 agents over the complete graph 𝐾 9 for the 2 𝑁 -state semi Markov exact model and the 𝑁-intertwined model. The subplots differ in the values of the time delay 𝜏 1 .

8 FIGURE 5

 85 FIGURE5 Evolution of the density of infected nodes 𝑦(𝑡) with time for 𝑁 = 9 agents over the "scale-free-3" graph for the 2 𝑁 -state semi Markov exact model and the 𝑁-intertwined model.

4 FIGURE 6

 46 FIGURE 6Evolution of the density of infected nodes 𝑦(𝑡) with time for 𝑁 = 9 agents over the "scale-free-7" graph for the 2 𝑁 -state semi Markov exact model and the 𝑁-intertwined model.

FIGURE 7

 7 FIGURE 7Variation of the error 𝜖(𝜏) with respect to 𝜏 over different types of graphs and for 9 agents. The vertical red line is the value of 𝜏 𝑐 = 1∕𝜆 𝑚𝑎𝑥 (𝐴).

FIGURE 8

 8 FIGURE 8Simulation of the exact SIS process on the complete graph. Density of infected nodes 𝑦(𝑡) as a function of time for 100 agents (N=100). Left plot corresponds to time delay 𝜏 1 = 𝜏 0 = 0. Right plot corresponds to 𝜏 0 = 0 and 𝜏 1 = 3.

FIGURE 9

 9 FIGURE 9Variation of the error 𝜖(𝜏, 𝑁) with respect to 𝜏 over different types of graphs and for different number of agents. The vertical green line is the value of 𝜏 𝑐 = 1∕𝜆 𝑚𝑎𝑥 (𝐴).

TABLE 1

 1 States binary notation of the 2 𝑁 -state semi-Markov chain State number 𝑖 𝑥 𝑁 𝑥 𝑁-1

  𝑥 𝑘 is either 0 or 1, 𝑖, 𝑗 ∈ Ω. If one agent in a state 𝑖 ∈ Ω is healthy (𝑥 𝑚 = 0), it will be infected with a rate 𝑞 𝐴 𝑖𝑗 = 𝛽 ∑ 𝑁 𝑘=1 𝑎 𝑖𝑗 𝑥 𝑘 . 𝑞 𝐴 𝑖𝑗 is simply 𝛽 multiplied by the number of the infected neighbors of the healthy agent. Otherwise, if an agent in a state 𝑖 ∈ Ω is infected (𝑥 𝑚 = 1), it will be cured with a rate 𝑞 𝛿 𝑖𝑗 = 𝛿. Note that 𝑄 is written as a sum of three matrices 𝑄 = 𝑄 𝛿 + 𝑄 𝐴 + 𝑄 𝑑 where 𝑄 𝛿 is a lower triangular matrix of elements 𝑞 𝛿 𝑖𝑗 , 𝑄 𝐴 is an upper triangular matrix of elements 𝑞 𝐴 𝑖𝑗 and 𝑄 𝑑 is a diagonal matrix of elements 𝑞 𝑖𝑖 . This formulation is well explained in Van Mieghem et al.

The bipartite graph is constructed using two sets of nodes. The first set contains 𝑁 5 nodes and the second set contains4𝑁 5 nodes.The Erdos-Renyi graph is constructed using a link probability of 𝑝 = 2𝑝 𝑐 = 2 𝑙𝑛𝑁 𝑁 to ensure high probability of graph connectivity. The scale free graph is constructed using Barabsi-Albert algorithm and we start by 4 connected nodes (𝑚 = 4) as the seed.

The interval 𝜏 < 𝜏 𝑐 is not shown on Figure9for the complete and the bipartite graphs because 𝜏 𝑐 is very small in this case and this is not compatible with the scale of the figure. For the other graphs, because the graphs are of different sizes, we have plotted the minimum value of 𝜏 𝑐 for each graph type.

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

SUPPORTING INFORMATION

This paper has supplementary downloadable material provided by the authors. This includes two videos illustrating Theorem 2 and a readme file.