Stochastic Compartment Model With Mortality and Its Application to Epidemic Spreading in Complex Networks

1. Introduction

Sudden or recurrent emergence of epidemics has been an everlasting thread to humanity. Highly infectious and fatal diseases such as pestilence, typhus, cholera, and leprosy were among the main causes of death in medieval times in Europe and until 20th century a major scourge of humanity [1]. This permanent challenge has naturally driven interest in protective measures and predictive models.

The systematic mathematical study of epidemic spreading began only a century ago with the seminal work of Kermack and McKendrick [2]. They were the first to introduce what we call nowadays a “compartment model”. In their so called SIR-model the individuals are categorized in compartments susceptible - S, infected - I, recovered (immune) - R characterizing their states of their health. Whereas standard SIR type models are able to capture main features of a certain class of infectious diseases such as mumps, measles and rubella, they fail to describe persistent oscillatory behaviors and spontaneous outbursts which are observed in many epidemics.

A large amount of work still is devoted to compartmental models [3,4] where an impressive field has emerged [5,6] and the interest was again considerably enhanced by the context of COVID-19 pandemics [7]. Beside purely macroscopic models, the study of epidemic dynamics in complex networks has attracted considerable attention [9,10,11]. In these works the importance of the graph topology for spreading phenomena has been highlighted. In particular, Pastor-Satorras and Vespignani showed that for a wide range of scale free networks no critical threshold for epidemic spreading exists [10]. The topological features crucial for epidemic spreading include the small world property (short average network distances1) and a high clustering coefficient measuring the existence of redundant paths between pairs of nodes [25,27].

Further interesting directions are represented by combinations of network science and stochastic compartmental models [16,17,28,29,31]. Such models include Markovian and non-Markovian approaches [23,30,35,36,37] where non-Markovianity is introduced by non-exponentially distributed sojourn times in the compartments [33,34]. In these works explicit formulae for the endemic equilibrium in terms of mean compartmental sojourn times and the basic reproduction number are derived and numerically validated in random walk simulations. A further non-Markovian model appeared recently [12] where non-Markovianity comes into play by introducing an "age of infection" allowing individuals to recover when their infection period exceeds a certain threshold, generalizing the initial idea of Kermack and McKendrick.

Other recent works emphasize the importance of spatial heterogeneity effects of infection patterns in epidemic spreading phenomena [38] and the role of local clusters to generate periodic epidemic outbursts [39] and see [40] for a review of these effects. A cluster model to explain periodic behavior was introduced a long time ago [41]. The role of the complex interplay of retardation (delayed compartmental transitions) and fluctuations for oscillatory behavior has been investigated in one of our recent works [35].

The aim of the present paper is to study the spreading of a certain class of diseases in a population of individuals (random walkers) moving on complex graphs aiming to mimic human mobility patterns in complex environments such as cities, street, and transportation networks. Essential elements in our model are the account for the mortality of infected individuals (random walkers) and an indirect transmission pathway via vectors (detailed below).

The present paper is organized as follows. In Section 2 we establish a stochastic mean field model for the evolution of the compartmental populations. The special case of zero mortality is considered in Section 3 where we obtain an explicit formula for the endemic equilibrium (stationary constant compartmental populations for infinite time). In this way we identify a crucial parameter controlling the stability of the healthy state having the interpretation of the basic reproduction number $R_0$ (Section 4) where the healthy state is stable for $R_0<1$ and unstable for $R_0>1$. A detailed proof of the stability of the endemic state for $R_0>1$ is provided in Appendix A.2. In Section 5 we analyze stability of the healthy state with mortality and derive the basic reproduction number $R_M$ and prove that $R_M<R_0$, i.e., mortality reduces the basic reproduction number. In Section 6 we test robustness of our mean field model under "complex real world conditions" by implementing its assumptions into multiple random random walkers simulations on Barabási-Albert (BA), Erdös-Rényi (ER) and Watts-Strogatz (WS) type graphs ([16,17,24] and Appendix A.3). These graph types have different complexity and connectivity features with impact on the spreading.

2. Compartmental Model with Mortality

The goal of this section is to develop a mean field model for a certain class of diseases with indirect infection transmission via vectors which includes Dengue, Malaria (transmission by mosquitos) or Pestilence (transmission by fleas) and others [8,9]. To that end we consider a population of Z random walkers navigating independently on a connected (ergodic) graph. Each walker performs independent steps from one to another connected node on the network (specified subsequently). We assume that walkers and nodes are in one of the compartments, S (susceptible) and I (infected). In addition, walkers can be in compartment D (dead) whereas nodes never die.

Let $Z_S\left(t\right),Z_I\left(t\right)$ ($N_S\left(t\right),N_I\left(t\right)$) be the number of walkers (nodes) in compartments S and I, and $Z_D\left(t\right)$ the non-decreasing number of walkers (in compartment D) died from the disease up to time t. We consider $Z=Z_I\left(t\right)+Z_S\left(t\right)+Z_D\left(t\right)$ walkers (Z independent of time) and a constant number of nodes $N=N_I\left(t\right)+N_S\left(t\right)$. We assume at instant $t=0$ the first spontaneous occurrence of the disease of a few infected walkers $Z_I\left(0\right)\ll Z$ or nodes $N_I\left(0\right)\ll N$ (and no dead walkers $Z_D\left(0\right)=0$). We introduce the compartmental fractions $S_w\left(t\right)=\frac{Z_S\left(t\right)}{Z}$, $J_w\left(t\right)=\frac{Z_I\left(t\right)}{Z}$, $d_w\left(t\right)=\frac{Z_d\left(t\right)}{Z}$ for the walkers (normalized with respect to Z) with $S_w\left(t\right)+J_w\left(t\right)+d_w\left(t\right)=1$, and $S_n\left(t\right)=\frac{N_S\left(t\right)}{N}$, $J_n\left(t\right)=\frac{N_I\left(t\right)}{N}$ with $S_n\left(t\right)+J_n\left(t\right)=1$. To limit the complexity of our model we do not consider the demographic evolution, i.e., there are no natural birth and dead processes. We denote with ${\mathcal{A}}_w\left(t\right),{\mathcal{A}}_n\left(t\right)$ the infection rates (rates of transitions S → I) of walkers and nodes, respectively. We assume the following simple bi-linear forms

$$\begin{array}{ccc}{\displaystyle {\mathcal{A}}_w\left(t\right)}& {\displaystyle ={\mathcal{A}}_w[S_w\left(t\right),J_n\left(t\right)]={\beta }_w{S}_w\left(t\right)J_n\left(t\right)}\hfill & \\ \\ {\displaystyle {\mathcal{A}}_n\left(t\right)}& {\displaystyle ={\mathcal{A}}_n[S_n\left(t\right),J_w\left(t\right)]={\beta }_n{S}_n\left(t\right)J_w\left(t\right)}\hfill \end{array}$$

(1)

with constant rate parameters ${\beta }_w,{\beta }_n>0$ (independent of time). ${\mathcal{A}}_w\left(t\right)$ indicates the infection rate of walkers where its dependence of $S_w,J_n$ is telling us that susceptible walkers can be infected only by (visiting) infected nodes. ${\mathcal{A}}_n\left(t\right)$ stands for the infection rate of nodes depending on $S_n\left(t\right),J_w\left(t\right)$ indicating that susceptible nodes can only be infected by (visits of) infected walkers. There are no direct transmissions among walkers and among nodes. Infections of walkers (nodes) take place with specific transmission probabilities in a contact of a node and a walker which are captured by (yet not identical with) the transmission rate constants ${\beta }_{w,n}$.

The infection time spans $t_I^{w,n}>0$ without mortality (waiting times in compartment I) of walkers and nodes are assumed to be independent random variables drawn from specific distributions specified hereafter. As only admitted dead process we assume that infected walkers may die within the time span of their infection. To capture this kind of mortality caused by the disease, we introduce a further independent random variable $t_M>0$ which indicates the life span of an infected walker. Both the infection and life time spans $t_I^w,t_M$ are counted from the time instant of the infection. A walker survives the disease if $t_M>t_I^w$ and dies from it for $t_M<t_I^w$. With these assumptions, we give first a stochastic formulation of the evolution equations

$$\begin{array}{ccc}{\displaystyle \frac{d}{dt}{S}_w\left(t\right)}& {\displaystyle =-{\mathcal{A}}_w\left(t\right)+〈{\mathcal{A}}_w(t-t_I^w)\mathsf{\Theta }(t_M-t_I^w)〉+J_w\left(0\right)\langle \delta (t-t_I^w)\mathsf{\Theta }(t_M-t_I^w)\rangle }\hfill & \\ \\ {\displaystyle \frac{d}{dt}{J}_w\left(t\right)}& {\displaystyle ={\mathcal{A}}_w\left(t\right)-〈{\mathcal{A}}_w(t-t_I^w)\mathsf{\Theta }(t_M-t_I^w)〉-J_w\left(0\right)〈\delta (t-t_I^w)\mathsf{\Theta }(t_M-t_I^w)〉-\frac{d}{dt}{d}_w\left(t\right)}\hfill & \\ \\ {\displaystyle \frac{d}{dt}{d}_w\left(t\right)=}& {\displaystyle 〈{\mathcal{A}}_w(t-t_M)\mathsf{\Theta }(t_I^w-t_M)〉+J_w\left(0\right)\langle \delta (t-t_M)\mathsf{\Theta }(t_I^w-t_M)\rangle }\hfill & \\ \\ {\displaystyle \frac{d}{dt}{S}_n\left(t\right)}& {\displaystyle =-{\mathcal{A}}_n\left(t\right)+〈{\mathcal{A}}_n(t-t_I^n)〉+J_n\left(0\right)\langle \delta (t-t_I^n)\rangle }\hfill & \\ \\ {\displaystyle \frac{d}{dt}{J}_n\left(t\right)}& {\displaystyle =-\frac{d}{dt}{S}_n\left(t\right)}\hfill \end{array}$$

(2)

where $\frac{d}{dt}{d}_w\left(t\right)$ indicates the (non-negative) mortality rate of walkers. We indicate with $〈..〉$ average over the contained (set of independent) random variables $t_I^w,t_I^n,t_M$ outlined hereafter and in Appendix A.1. $\mathsf{\Theta }(..)$ stands for the Heaviside function (A2), and $\delta (..)$ for the Dirac’s $\delta$-distribution. An epidemic always starts from “natural” initial conditions $S_w\left(0\right)=1$, $S_n\left(0\right)=1$ (globally healthy state) where at $t=0$ the first infections occur spontaneously and can be “generated” by adding the source terms $J_{w,n}\left(0\right)\delta \left(t\right)$ to the infection rates of walkers and nodes, respectively. Equivalently, we introduce initial conditions $S_{w,n}\left(0\right)=1-J_{w,n}\left(0\right)$ ($d_w\left(0\right)=0$) with $J_{w,n}\left(0\right)>0$ consisting typically of a few infected walkers and/or nodes in a large susceptible population without dead walkers $d_w\left(0\right)=0$.

The interpretation of system (2) is as follows. The instantaneous infection rate ${\mathcal{A}}_w\left(t\right)$ governs the transitions S → I of walkers (due to visits of infected nodes). The term $〈{\mathcal{A}}_w(t-t_I^w)\mathsf{\Theta }(t_M-t_I^w)〉$ describes the rate of walkers recovering at time t and infected at $t-t_I^w$, i.e., their infection time span has elapsed and they survived as $t_M>t_I^w$ (indicated by $\mathsf{\Theta }(t_M-t_I^w)=1$). Then $〈{\mathcal{A}}_w(t-t_M)\mathsf{\Theta }(t_I^w-t_M)〉$ captures the rate of walkers infected at $t-t_M$ dying at at time t during the infection time span (indicated by $\mathsf{\Theta }(t_I^w-t_M)=1$ for $t_I^w>t_M$).

Remark I

The infection time span of a walker (sojourn time in compartment I) is $min(t_I^w,t_M)$, i.e., $t_I^w$ if $t_M>t_I^w$ (where the walker survives the disease), and is $t_M$ if the walker dies within the infectious time span ($t_M<t_I^w$). $t_I^w$ is the walker’s infection time span without mortality (retrieved for $t_M\to \infty$). The probability of persistence of a walker’s infection at time t, given the infection starts at time 0 is $\langle \mathsf{\Theta }(t_I^w-t)\mathsf{\Theta }(t_M-t)\rangle$ (see (7)). Note that $\mathsf{\Theta }(t_I^w-t)\mathsf{\Theta }(t_M-t)=1$ only if $t<min(t_I^w,t_M)$, i.e., when the walker is in compartment I. As a crucial element of our model, we will analyze the statistics of the walker’s infection time span $min(t_I^w,t_M)$.

The initially infected walkers and nodes are as well subjected to the transition pathways, i.e., walkers either recover (alive) with rate $J_w\left(0\right)\langle \delta (t-t_I^w)\mathsf{\Theta }(t_M-t_I^w)\rangle$ or they die with rate $J_w\left(0\right)\langle \delta (t-t_M)\mathsf{\Theta }(t_I^w-t_M)\rangle$, and nodes always recover with rate $J_n\left(0\right)\langle \delta (t-t_I^n)\rangle$. For $t\to \infty$ these terms are evanescent thus the initial conditions do not affect large time limits (endemic state for zero mortality). The importance of these terms can be seen by setting ${\beta }_{w,n}=0$ (no infections). Without these terms the initially infected walkers and nodes would stay infected forever, inconsistent with our assumptions.

The rate equations for the nodes can be interpreted in the same way as interplay of instantaneous infections and delayed recovery without mortality. We emphasize that the evolution equations of the nodes and walkers are non-linearly coupled by the implicit dependencies of the infection rates (1). In order to derive an explicit representation of system (2), we need to have a closer look on the averaging procedures and the involved distributions related to the independent random variables $T=\{t_I^w,t_I^n,t_M\}>0$ drawn from specific probability density functions (PDFs) which we define by

$$Prob[T\in [\tau ,\tau +\mathrm{d}\tau ]=K(\tau )\mathrm{d}\tau ,$$

(3)

with their respective PDFs (kernels) $K\left(\tau \right)=\{K_I^{w,n}\left(\tau \right),K_M\left(\tau \right)\}$ which are normalized $Prob[T>0]={\int }_0^{\infty }K\left(\tau \right)\mathrm{d}\tau =1$. Then recall the averaging rule for (suitable) functions $f\left(T\right)$ of the random variable T which we use throughout the paper

$$〈f\left(T\right)〉={\int }_0^{\infty }K\left(\tau \right)f\left(\tau \right)\mathrm{d}\tau ,$$

(4)

see also Appendix A.1. An important case is $\langle \delta (t-T)\rangle =K\left(t\right)$. Then by applying (4) we introduce the persistence probabilities of walker’s (node’s) infection (without mortality)

$${\mathsf{\Phi }}_I^{w,n}\left(t\right)=Prob(t_I^{w,n}>t)=\langle \mathsf{\Theta }(t_I^{w,n}-t)\rangle ={\int }_t^{\infty }{K}_I^{w,n}\left(\tau \right)\mathrm{d}\tau$$

(5)

and the probability of walker’s survival up to time t (given $t_I^w=\infty$)

$${\mathsf{\Phi }}_M\left(t\right)=Prob(t_M>t)=\langle \mathsf{\Theta }(t_M-t)\rangle ={\int }_t^{\infty }{K}_M\left(\tau \right)\mathrm{d}\tau .$$

(6)

The persistence probabilities fulfill the initial condition ${\mathsf{\Phi }}_M\left(0\right)={\mathsf{\Phi }}_I^{w,n}\left(0\right)=1$ corresponding to the normalization of the waiting time PDFs $K\left(\tau \right)=\{K_I^{w,n}\left(\tau \right),K_M\left(\tau \right)\}$ and are vanishing at infinity ${\mathsf{\Phi }}_M(\infty )={\mathsf{\Phi }}_I^{w,n}(\infty )=0$. To evaluate the averages in (2) we will use the following quantities:

$$\begin{array}{ccc}{\displaystyle \langle \delta (t-T)\rangle }& {\displaystyle =K\left(t\right),T=\{t_I^w,t_I^n,t_M\}}\hfill & \\ \\ {\displaystyle \langle \mathsf{\Theta }(t_M-t)\mathsf{\Theta }(t_I^w-t)\rangle }& {\displaystyle =\langle \mathsf{\Theta }(t_M-t)\rangle \langle \mathsf{\Theta }(t_I^w-t)\rangle ={\mathsf{\Phi }}_I^w\left(t\right){\mathsf{\Phi }}_M\left(t\right)}\hfill & \\ \\ {\displaystyle b_d\left(t\right)=\langle \delta (t-t_M)\mathsf{\Theta }(t_I^w-t_M)\rangle }& {\displaystyle =\langle \delta (t-t_M)\rangle \langle \mathsf{\Theta }(t_I^w-t)\rangle =K_M\left(t\right){\mathsf{\Phi }}_I^w\left(t\right)}\hfill & \\ \\ {\displaystyle b_r\left(t\right)=\langle \delta (t-t_I^w)\mathsf{\Theta }(t_M-t_I^w)\rangle }& {\displaystyle =\langle \delta (t-t_I^w)\rangle \langle \mathsf{\Theta }(t_M-t)\rangle =K_I^w\left(t\right){\mathsf{\Phi }}_M\left(t\right)}\hfill & \\ \\ {\displaystyle b_d\left(t\right)+b_r\left(t\right)={\mathcal{K}}_{I,M}^w\left(t\right)}& {\displaystyle =-\frac{d}{dt}\left[\langle \mathsf{\Theta }(t_M-t)\rangle \langle \mathsf{\Theta }(t_I^w-t)\rangle \right]=-\frac{d}{dt}\left[{\mathsf{\Phi }}_I^w\left(t\right){\mathsf{\Phi }}_M\left(t\right)\right]}\hfill & \\ \\ {\displaystyle {\int }_0^{\infty }{\mathcal{K}}_{I,M}^w\left(t\right)\mathrm{d}t}& {\displaystyle =1}\hfill & \\ \\ {\displaystyle \mathcal{R}\left(t\right)=\langle \mathsf{\Theta }(t-t_I^w)\mathsf{\Theta }(t_M-t_I^w)\rangle }& {\displaystyle ={\int }_0^t{b}_r\left(\tau \right)\mathrm{d}\tau ={\int }_0^t{K}_I^w\left(\tau \right){\mathsf{\Phi }}_M\left(\tau \right)\mathrm{d}\tau }\hfill & \\ \\ {\displaystyle \mathcal{D}\left(t\right)=\langle \mathsf{\Theta }(t-t_M)\mathsf{\Theta }(t_I^w-t_M)\rangle }& {\displaystyle ={\int }_0^t{b}_d\left(\tau \right)\mathrm{d}\tau ={\int }_0^t{K}_M\left(\tau \right){\mathsf{\Phi }}_I^w\left(\tau \right)\mathrm{d}\tau }\hfill & \\ \\ {\displaystyle \mathcal{R}\left(t\right)+\mathcal{D}\left(t\right)={\int }_0^t{\mathcal{K}}_{I,M}^w\left(\tau \right)\mathrm{d}\tau }& {\displaystyle ={\int }_0^t[b_d\left(\tau \right)+b_r\left(\tau \right)]\mathrm{d}\tau =1-{\mathsf{\Phi }}_I^w\left(t\right){\mathsf{\Phi }}_M\left(t\right)}\hfill & \\ \\ {\displaystyle \mathcal{D}(\infty )+\mathcal{R}(\infty )=1}& \\ \\ {\displaystyle \langle A(t-t_I)\mathsf{\Theta }(t_M-t_I)\rangle }& {\displaystyle =\langle A(t-t_I){\mathsf{\Phi }}_M\left(t_I\right)\rangle ={\int }_0^{t}A(t-\tau ){\mathsf{\Phi }}_M\left(\tau \right)K_I\left(\tau \right)\mathrm{d}\tau }\hfill & \\ \\ {\displaystyle \langle A(t-t_M\mathsf{\Theta }(t_I^w-t_M)\rangle }& {\displaystyle =\langle A(t-t_M){\mathsf{\Phi }}_I^w\left(t_M\right)\rangle ={\int }_0^{t}A(t-\tau ){\mathsf{\Phi }}_I^w\left(\tau \right)K_M\left(\tau \right)\mathrm{d}\tau }\hfill \end{array}$$

(7)

In these averages we make use of the independence of the waiting times $t_M,t_I^{w,n}$, and of causality of $A\left(\tau \right)$ and the kernels $K\left(\tau \right)$ (i.e., $A\left(\tau \right),K\left(\tau \right)=0$ for $\tau <0$). Of utmost importance are the "defective" PDFs (DPDFs) $b_{d,r}\left(t\right)$ of death and recovery. "Defective" means that $b_{d,r}\left(t\right)$ are no proper PDFs since they are not normalized to one, but rather to $\mathcal{D}(\infty ),\mathcal{R}(\infty )<1$, respectively. Consult [13] for a recent account of defective distributions and related stochastic processes. They have the following interpretation. $b_d\left(t\right)\mathrm{d}t=K_M\left(t\right){\mathsf{\Phi }}_I^w\left(t\right)\mathrm{d}t$ is the probability of transition I → D within $[t,t+\mathrm{d}t]$ of an infected walker (infected at $t^{\prime }=0$). $b_r\left(t\right)\mathrm{d}t=K_I^w\left(t\right){\mathsf{\Phi }}_M\left(t\right)\mathrm{d}t$ is the probability of transition I → S within $[t,t+\mathrm{d}t]$ of a walker infected at $t^{\prime }=0$. Therefore,

$${\mathcal{K}}_{I,M}^w\left(t\right)=b_r\left(t\right)+b_d\left(t\right)=-\frac{d}{dt}〈\mathsf{\Theta }(t_M-t)\mathsf{\Theta }(t_I^w-t)〉=-\frac{d}{dt}\left[{\mathsf{\Phi }}_I^w\left(t\right){\mathsf{\Phi }}_M\left(t\right)\right]$$

(8)

is non-negative (as are $K_I^w=-\frac{d}{dt}{\mathsf{\Phi }}_I^w\ge 0$, $K_M=-\frac{d}{dt}{\mathsf{\Phi }}_M\ge 0$) and is a proper well-normalized PDF of all exits of walkers from compartment I (i.e., I → S + I → D). Without mortality (${\mathsf{\Phi }}_M\left(t\right)=1$) this PDF retrieves ${\mathcal{K}}_{I,M}^w\left(t\right)=K_I^w\left(t\right)$.

The quantities $\mathcal{R}\left(t\right),\mathcal{D}\left(t\right)$ introduced in (7) have the following interpretation. $\mathcal{R}\left(t\right)$ is the probability that a walker infected at instant 0 is at time t in compartment S (i.e., recovered prior or up to time t). $\mathcal{D}\left(t\right)$ is the probability that a walker infected at instant 0 is at time t in compartment D (i.e., died prior and up to time t). Important are the infinite time limits: $\mathcal{R}(\infty )$ has the interpretation of the overall probability that an infected walker survives the infection, and $\mathcal{D}(\infty )$ is the overall probability for an infected walker to die from the disease. We refer $\mathcal{D}(\infty )$ also to as "overall mortality". It must not be confused with the infinite time limit of dead walkers fraction $d_w(\infty )$ which is different from $\mathcal{D}(\infty )$ as we will see in details subsequently. A small value $\mathcal{D}(\infty )$ may cause a high value of $d_w(\infty )$ for instance for short infectious periods where walkers may be repeatedly infected.

In most cases not all infected walkers die from their disease (in an infinite observation time), hence $\mathcal{D}(\infty )<1$ (as $b_d$ is defective). $\mathcal{D}(\infty )\to 1$ represents the limit of a fatal disease, and $\mathcal{D}(\infty )\to 0$ a disease without mortality. $\mathcal{R}(\infty )<1$ (as $b_r$ is defective) is the complementary probability with $\mathcal{D}(\infty )+\mathcal{R}(\infty )=1$.

With these remarks, system (2) reads

$$\begin{array}{ccc}{\displaystyle \frac{d}{dt}{S}_w\left(t\right)}& {\displaystyle =-{\mathcal{A}}_w\left(t\right)+{\int }_0^t{\mathcal{A}}_w(t-\tau )K_I^w\left(\tau \right){\mathsf{\Phi }}_M\left(\tau \right)\mathrm{d}\tau +J_w\left(0\right)K_I^w\left(t\right){\mathsf{\Phi }}_M\left(t\right)}\hfill & \\ \\ {\displaystyle \frac{d}{dt}{J}_w\left(t\right)}& {\displaystyle =\frac{d}{dt}{\int }_0^t{\mathcal{A}}_w\left(\tau \right){\mathsf{\Phi }}_M(t-\tau ){\mathsf{\Phi }}_I^w(t-\tau )\mathrm{d}\tau -J_w\left(0\right)\left[K_I^w\left(t\right){\mathsf{\Phi }}_M\left(t\right)+K_M\left(t\right){\mathsf{\Phi }}_I^w\left(t\right)\right]}\hfill & \\ \\ {\displaystyle \frac{d}{dt}{S}_n\left(t\right)}& {\displaystyle =-{\mathcal{A}}_n\left(t\right)+{\int }_0^t{\mathcal{A}}_n(t-\tau )K_I^n\left(\tau \right)\mathrm{d}\tau +J_n\left(0\right)K_I^n\left(t\right)}\hfill & \\ \\ {\displaystyle \frac{d}{dt}{J}_n\left(t\right)}& {\displaystyle =-\frac{d}{dt}{S}_n\left(t\right).}\hfill & \\ & \end{array}$$

(9)

The PDF (8) that a walker leaves compartment I (either by recovery or by death) allows to rewrite the second equation of (9) as

$$\frac{d}{dt}{J}_w\left(t\right)={\mathcal{A}}_w\left(t\right)-{\int }_0^t{\mathcal{A}}_w(t-\tau ){\mathcal{K}}_{I,M}^w\left(\tau \right)\mathrm{d}\tau -J_w\left(0\right){\mathcal{K}}_{I,M}^w\left(t\right).$$

(10)

Worthy of closer consideration is the mortality rate of the infected walkers (representing the total mortality – entry rate into the D compartment)

$$\begin{array}{ccc}{\displaystyle \frac{d}{dt}{d}_w\left(t\right)=-\frac{d}{dt}(S_w\left(t\right)+J_w\left(t\right))}& {\displaystyle =〈{\mathcal{A}}_w(t-t_M)\mathsf{\Theta }(t_I^w-t_M)〉+J_w\left(0\right)\langle \delta (t-t_M)\mathsf{\Theta }(t_I^w-t_M)\rangle }\hfill & \\ \\ & {\displaystyle ={\int }_0^t{\mathcal{A}}_w(t-\tau )K_M\left(\tau \right){\mathsf{\Phi }}_I^w\left(\tau \right)\mathrm{d}\tau +J_w\left(0\right)K_M\left(t\right){\mathsf{\Phi }}_I^w\left(t\right)}\hfill \end{array}$$

(11)

where clearly $\frac{d}{dt}{d}_w\left(t\right)\ge 0$. Integrating this relation yields the fraction $d_w\left(t\right)$ of dead walkers

$$\begin{array}{ccc}{\displaystyle d_w\left(t\right)}& {\displaystyle =1-S_w\left(t\right)-J_w\left(t\right)}\hfill & \\ \\ & {\displaystyle ={\int }_0^t{\mathcal{A}}_w(t-\tau )\langle \mathsf{\Theta }(\tau -t_M)\mathsf{\Theta }(t_I^w-t_M)\rangle \mathrm{d}\tau +J_w\left(0\right)\langle \mathsf{\Theta }(t-t_M)\mathsf{\Theta }(t_I^w-t_M)\rangle }\hfill & \\ \\ & {\displaystyle ={\int }_0^t{\mathcal{A}}_w(t-\tau )\mathcal{D}\left(\tau \right)\mathrm{d}\tau +J_w\left(0\right)\mathcal{D}\left(t\right).}\hfill \end{array}$$

(12)

An interesting quantity is the cumulative recovery rate of walkers (integrated entry rates of walkers into the S compartment, see first equation in (2))

$$\begin{array}{ccc}{\displaystyle r_w\left(t\right)}& {\displaystyle ={\int }_0^t{\mathcal{A}}_w(t-\tau )〈\mathsf{\Theta }(\tau -t_I^w)\mathsf{\Theta }(t_M-t_I^w)〉\mathrm{d}\tau +J_w\left(0\right)〈\mathsf{\Theta }(t-t_I^w)\mathsf{\Theta }(t_M-t_I^w)〉}\hfill & \\ \\ & {\displaystyle ={\int }_0^t{\mathcal{A}}_w(t-\tau )\mathcal{R}\left(\tau \right)\mathrm{d}\tau +J_w\left(0\right)\mathcal{R}\left(t\right).}\hfill \end{array}$$

(13)

The quantity $r_w\left(t\right)$ records all recovery events of walkers up to time t, where individual walkers may suffer repeated infections and recoveries. We observe that (see (7))

$$r_w\left(t\right)+d_w\left(t\right)={\int }_0^t{\mathcal{A}}_w(t-\tau )[1-{\mathsf{\Phi }}_I^w\left(\tau \right){\mathsf{\Phi }}_M\left(\tau \right)]\mathrm{d}\tau +J_w\left(0\right)\left(1-{\mathsf{\Phi }}_I^w\left(t\right){\mathsf{\Phi }}_M\left(t\right)\right].$$

(14)

Relation (12) records all dead events of walkers up to time t. Since each walker may die only once, it follows indeed that $d_w\left(t\right)\in [0,1]$. Contrarily, the quantity $r_w\left(t\right)$ is not restricted to this interval as walkers may be repeatedly infected and recovered but due to mortality eventually only a finite number of times ($r_w(\infty )<\infty$, see (18)). Mortality renders the chain of infection and recovery events transient (due to the defective feature of $b_r=K_I^w{\mathsf{\Phi }}_M$). To shed more light on the behavior of $r_w\left(t\right)$ consider for a moment zero mortality ($\mathcal{R}(\infty )=1$) and $t\to \infty$: We then have ${\mathcal{A}}_w(\infty )={\beta }_w{S}_w^e{J}_n^e>0$ (shown in Section 3) thus $r_w(\infty )=\infty$ coming along with an infinite chain of recurrent infection and recovery events (as $b_r\left(t\right)$ turns into the proper non-defective PDF $b_r=K_I^w$).

Using (7) we can rewrite (2) in equivalent integral form

$$\begin{array}{ccc}{\displaystyle S_w\left(t\right)}& {\displaystyle =1-J_w\left(0\right)\left[{\mathsf{\Phi }}_M\left(t\right){\mathsf{\Phi }}_I^w\left(t\right)+\mathcal{D}\left(t\right)\right]-{\int }_0^t{\mathcal{A}}_w\left(\tau \right)[{\mathsf{\Phi }}_M(t-\tau ){\mathsf{\Phi }}_I^w(t-\tau )+\mathcal{D}(t-\tau )]\mathrm{d}\tau }\hfill & \\ \\ {\displaystyle J_w\left(t\right)}& {\displaystyle =J_w\left(0\right){\mathsf{\Phi }}_M\left(t\right){\mathsf{\Phi }}_I^w\left(t\right)+{\int }_0^t{\mathcal{A}}_w\left(\tau \right){\mathsf{\Phi }}_M(t-\tau ){\mathsf{\Phi }}_I^w(t-\tau )\mathrm{d}\tau }\hfill & \\ \\ {\displaystyle S_n\left(t\right)}& {\displaystyle =1-J_n\left(t\right)}\hfill & \\ \\ {\displaystyle J_n\left(t\right)}& {\displaystyle =J_n\left(0\right){\mathsf{\Phi }}_I^n\left(t\right)+{\int }_0^t{\mathcal{A}}_n\left(\tau \right){\mathsf{\Phi }}_I^n(t-\tau )\mathrm{d}\tau }\hfill \end{array}$$

(15)

and with (redundant) Equation (12) for the fraction of dead walkers. (15) is a self-consistent system since the infection rates are implicit functions of susceptible and infected population fractions ${\mathcal{A}}_w\left(t\right)=A_w[S_w\left(t\right),J_n\left(t\right)]$, ${\mathcal{A}}_n\left(t\right)=A_w[S_n\left(t\right),J_w\left(t\right)]$ (see (1)). Explore now the infinite time limit of (15) where it is convenient to consider the Laplace transformed equations. We introduce the Laplace transform (LT, denoted with a hat) of a function $g\left(t\right)$ as

$$\widehat{g}\left(\lambda \right)={\int }_0^{\infty }g\left(t\right)e^{-\lambda t}\mathrm{d}t$$

where $\lambda$ denotes the (suitably chosen) Laplace variable. In order to retrieve infinite time limits we use the asymptotic feature

$$g(\infty )=\underset{\lambda \to 0}{lim}\lambda \widehat{g}\left(\lambda \right)\left(=\underset{\lambda \to 0}{lim}{\int }_0^{\infty }g\left(\frac{\tau }{\lambda }\right)e^{-\tau }\mathrm{d}\tau \to g(\infty ){\int }_0^{\infty }{e}^{-\tau }\mathrm{d}\tau \right).$$

(16)

Observing that the LT of ${\mathsf{\Phi }}_I^w\left(t\right){\mathsf{\Phi }}_M\left(t\right)$ is ${\lambda }^{-1}[1-{\widehat{\mathcal{K}}}_{I,M}\left(\lambda \right)]$ and $\widehat{D}\left(\lambda \right)={\lambda }^{-1}{\widehat{b}}_d\left(\lambda \right)$ where ${\widehat{b}}_d\left(0\right)=\mathcal{D}(\infty )$ (see (7)), we arrive at

$$\begin{array}{ccc}{\displaystyle J_w(\infty )}& {\displaystyle =\underset{\lambda \to 0}{lim}\lambda {\widehat{J}}_w\left(\lambda \right)=[1-{\widehat{\mathcal{K}}}_{I,M}\left(0\right)][J_w\left(0\right)+{\widehat{A}}_w\left(0\right)]=0}\hfill & \\ \\ {\displaystyle J_n(\infty )}& {\displaystyle =\underset{\lambda \to 0}{lim}\lambda {\widehat{J}}_n\left(\lambda \right)=[1-{\widehat{K}}_n\left(0\right)][J_n\left(0\right)+{\widehat{A}}_n\left(0\right)]=0}\hfill \\ \\ {\displaystyle d_w(\infty )}& {\displaystyle =\underset{\lambda \to 0}{lim}\lambda {\widehat{d}}_w\left(\lambda \right)=\mathcal{D}(\infty )[J_w\left(0\right)+{\widehat{A}}_w\left(0\right)],}\hfill & \hfill {\widehat{A}}_w\left(0\right)={\beta }_w{\int }_0^{\infty }{S}_w\left(\tau \right)J_n\left(\tau \right)\mathrm{d}\tau \\ \\ {\displaystyle S_w(\infty )}& {\displaystyle =1-d_w(\infty )}\hfill & \\ \\ {\displaystyle S_n(\infty )}& {\displaystyle =1}\hfill \end{array}$$

(17)

where ${\widehat{A}}_{w,n}\left(0\right)={\int }_0^{\infty }{\mathcal{A}}_{w,n}\left(t\right)\mathrm{d}t<\infty$. In the same way one obtains

$$r_w(\infty )=\mathcal{R}(\infty )(J_w\left(0\right)+{\widehat{A}}_w\left(0\right)).$$

(18)

Since $\mathcal{D}(\infty )$ is non-zero, the asymptotic values $S_w(\infty ),d_w(\infty )$ depend on the initial condition $J_w\left(0\right)$ and the infection (rate) history. This is not any more true for zero mortality ($\mathcal{D}(\infty )=0$) and considered in Section 3. We define the overall probability $P_D$ that a walker dies ($P_R$ survives) the disease

$$P_D=\frac{d_w(\infty )}{r_w(\infty )+d_w(\infty )}=\mathcal{D}(\infty ),P_R=1-P_D=\frac{r_w(\infty )}{r_w(\infty )+d_w(\infty )}=\mathcal{R}(\infty )$$

(19)

consistent with our previous interpretation of $\mathcal{D}(\infty ),\mathcal{R}(\infty )$, and the ratio

$$\frac{d_w(\infty )}{r_w(\infty )}=\frac{\mathcal{D}(\infty )}{\mathcal{R}(\infty )}.$$

(20)

The quantities (19) and (20) depend only on the stochastic features of $t_I^w$ and $t_M$. They are independent of the infection rates and initial conditions and therefore of the topological properties of the network. In addition, they also do not depend on the stochastic features of the node’s infection time span $t_I^n$.

Markovian (memoryless) case

Generally the system (9) contains the history of the process (memory) which makes the process non-Markovian. The only exception is when all waiting times are exponentially distributed, namely ${\mathsf{\Phi }}_I^w\left(t\right)=e^{-{\xi }_I^{w}t}$, ${\mathsf{\Phi }}_M\left(t\right)=e^{-{\xi }_{M}t}$, ${\mathsf{\Phi }}_I^n\left(t\right)=e^{-{\xi }_I^{n}t}$ ($\langle t_I^{w,n}\rangle ={\xi }_I^{w,n}{)}^{-1}$, $\langle t_M\rangle ={\left({\xi }_M\right)}^{-1}$). Then (9) takes with (15) the memoryless form

$$\begin{array}{ccc}{\displaystyle \frac{d}{dt}{S}_w\left(t\right)}& {\displaystyle =-{\beta }_w{S}_w\left(t\right)J_n\left(t\right)+{\xi }_I^w{J}_w\left(t\right)}\hfill & \\ \\ {\displaystyle \frac{d}{dt}{J}_w\left(t\right)}& {\displaystyle ={\beta }_w{S}_w\left(t\right)J_n\left(t\right)-({\xi }_I^w+{\xi }_M)J_w\left(t\right)}\hfill & \\ \\ {\displaystyle \frac{d}{dt}{d}_w\left(t\right)}& {\displaystyle ={\xi }_M{J}_w\left(t\right)}\hfill & \\ \\ {\displaystyle \frac{d}{dt}{S}_n\left(t\right)}& {\displaystyle =-{\beta }_n{S}_n\left(t\right)J_w\left(t\right)+{\xi }_I^n{J}_n\left(t\right)}\hfill & \\ \\ {\displaystyle \frac{d}{dt}{J}_n\left(t\right)}& {\displaystyle ={\beta }_n{S}_n\left(t\right)J_w\left(t\right)-{\xi }_I^n{J}_n\left(t\right).}\hfill \end{array}$$

(21)

Putting the left-hand sides to zero yields the stationary state

$$\begin{array}{ccc}{\displaystyle J_w(\infty )=J_n(\infty )={\mathcal{A}}_w(\infty )={\mathcal{A}}_n(\infty )=0}& & \\ \\ {\displaystyle S_w(\infty )=1-d_w(\infty ),}& {\displaystyle d_w(\infty )={\xi }_M{\int }_0^{\infty }{J}_w\left(\tau \right)\mathrm{d}\tau }\hfill & \\ \\ {\displaystyle S_n(\infty )=1.}& \end{array}$$

(22)

Let us check whether this result is consistent with (17). To this end, we integrate the second equation in (21) knowing that $J_w(\infty )=0$ leading to

$$0=J_w\left(0\right)+{\int }_0^{\infty }{\mathcal{A}}_w\left(t\right)\mathrm{d}t-({\xi }_I^w+{\xi }_M){\int }_0^{\infty }{J}_w\left(t\right)\mathrm{d}t$$

(23)

thus ${\int }_0^{\infty }{J}_w\left(t\right)\mathrm{d}t=\frac{1}{{\xi }_M+{\xi }_I^w}(J_w\left(0\right)+{\widehat{\mathcal{A}}}_w\left(0\right))$. Plugging this relation into (22) and accounting for $\mathcal{D}(\infty )=\frac{{\xi }_M}{{\xi }_I^w+{\xi }_M}$ recovers indeed the representation of expression (17).

For zero mortality ${\xi }_M=0$ one can straight-forwardly obtain the constant endemic equilibrium values $J_w^e,J_n^e$ by setting the left-hand side of (21) to zero leading to subsequent Equation (32) derived in Section 3 for general waiting time distributions with finite means.

A few more words on waiting time distributions

In our simulations we assume that the time spans $t_I^w,t_I^n,t_M$ are independent random variables drawn from specific Gamma distributions. The advantage to use Gamma distributions is that they may realize a large variety of shapes, see Figure 1 for a few examples. To generate Gamma distributed random numbers we employ the PYTHON random number generator (library numpy.random). Recall the Gamma distribution

$$K_{\alpha ,\xi }\left(t\right)=\frac{{\xi }^{\alpha }{t}^{\alpha -1}}{\Gamma \left(\alpha \right)}{e}^{-\xi t},\xi ,\alpha >0$$

(24)

where $\alpha$ is the so called `shape parameter’ and $\xi$ the rate parameter (often is used the term ’scale parameter’ $\theta ={\xi }^{-1}$) and $\Gamma \left(\alpha \right)$ stands for the Gamma function. We also will subsequently use the LT of the Gamma PDF

$${\widehat{K}}_{\alpha ,\xi }\left(\lambda \right)={\int }_0^{\infty }{K}_{\alpha ,\xi }\left(t\right)e^{-\lambda t}\mathrm{d}t=\frac{{\xi }^{\alpha }}{{(\lambda +\xi )}^{\alpha }}$$

(25)

The Gamma PDF has finite mean ${\langle t\rangle }_{\alpha ,\xi }={\int }_0^{\infty }t{K}_{\alpha ,\xi }\left(t\right)\mathrm{d}t=\frac{\alpha }{\xi }$ and for $\alpha <1$ the Gamma-PDF is weakly singular at $t=0$ and $\alpha =1$ recovers exponential PDFs. For $\alpha \le 1$ the Gamma PDF is completely monotonic (CM) (see Appendix, (A9) for a definition). For the range $\alpha >1$ the Gamma-PDF has a maximum at $t_{max}=\frac{\alpha -1}{\xi }$ and becomes narrower the larger $\xi$ (while keeping its mean $\alpha /\xi$ fixed), especially we can generate sharp waiting times using the feature

$$\underset{\xi \to \infty }{lim}{K}_{\alpha =\xi T_0,\xi }\left(t\right)=\delta (t-T_0).$$

(26)

We also will use subsequently the persistence probability of the Gamma distribution (see right frame of Figure 1)

$${\mathsf{\Phi }}_{\alpha ,\xi }\left(t\right)={\int }_t^{\infty }\frac{{\xi }^{\alpha }{t}^{\alpha -1}}{\Gamma \left(\alpha \right)}{e}^{-\xi t}\mathrm{d}t=\frac{\Gamma (\alpha ,\xi t)}{\Gamma \left(\alpha \right)}$$

(27)

where $\Gamma (\alpha ,x)$ indicates the upper incomplete Gamma function with $\Gamma (\alpha ,0)=\Gamma \left(\alpha \right)$. (27) fulfills necessarily the initial condition ${\mathsf{\Phi }}_{\alpha ,\xi }\left(0\right)=1$ and is vanishing at infinity ${\mathsf{\Phi }}_{\alpha ,\xi }(\infty )=0$.

3. Endemic Equilibrium for Zero Mortality

Here we consider the large time asymptotics of the compartment populations without mortality ($S_w\left(t\right)+J_w\left(t\right)=1$) where the self-consistent system (15) reads

$$\begin{array}{ccc}{\displaystyle J_w\left(t\right)}& {\displaystyle =J_w\left(0\right){\mathsf{\Phi }}_I^w\left(t\right)+{\int }_0^t{\mathcal{A}}_w(t-\tau ){\mathsf{\Phi }}_I^w\left(\tau \right)\mathrm{d}\tau }\hfill & \\ \\ {\displaystyle J_n\left(t\right)}& {\displaystyle =J_n\left(0\right){\mathsf{\Phi }}_I^n\left(t\right)+{\int }_0^t{\mathcal{A}}_n(t-\tau ){\mathsf{\Phi }}_I^n\left(\tau \right)\mathrm{d}\tau }\hfill \end{array}$$

(28)

The endemic state emerging in the large time asymptotics does not depend on the initial conditions $J_{w,n}\left(0\right)$ as ${\mathsf{\Phi }}_I^{w,n}\left(t\right)\to 0$. For what follows it is convenient to consider the LTs of (28) which read

$$\begin{array}{ccc}{\displaystyle {\widehat{J}}_w\left(\lambda \right)}& {\displaystyle =\left[J_w\left(0\right)+{\widehat{\mathcal{A}}}_w\left(\lambda \right)\right]\frac{1-{\widehat{K}}_I^w\left(\lambda \right)}{\lambda }}\hfill & \\ \\ {\displaystyle {\widehat{J}}_n\left(\lambda \right)}& {\displaystyle =\left[J_n\left(0\right)+{\widehat{\mathcal{A}}}_n\left(\lambda \right)\right]\frac{1-{\widehat{K}}_I^n\left(\lambda \right)}{\lambda }}\hfill \end{array}$$

(29)

where ${\mathsf{\Phi }}_I^{w,n}\left(\lambda \right)=\frac{1-{\widehat{K}}_I^{w,n}\left(\lambda \right)}{\lambda }$ are the LTs of the persistence distributions, and ${\widehat{S}}_{w,n}\left(\lambda \right)+{\widehat{J}}_{w,n}\left(\lambda \right)=\frac{1}{\lambda }$ reflecting constant populations of walkers and nodes. In order to determine the endemic equilibrium, we assume that the mean infection time spans for the nodes and walkers exist

$$\langle t_I^{w,n}\rangle =\underset{\lambda \to 0}{lim}\frac{1-{\widehat{K}}_I^{w,n}\left(\lambda \right)}{\lambda }=-\frac{d}{d\lambda }{\widehat{K}}_I^{w,n}\left(\lambda \right){|}_{\lambda =0}={\int }_0^{\infty }{\mathsf{\Phi }}_I^{w,n}\left(t\right)\mathrm{d}t={\int }_0^{\infty }\tau K_I^{w,n}\left(\tau \right)\mathrm{d}\tau <\infty$$

(30)

thus the admissible range of the Laplace variable is $\lambda \ge 0$ (if chosen real). Now using (16) we obtain the endemic equilibrium from $J_{w,n}(\infty )={lim}_{\lambda \to 0}\lambda {\widehat{J}}_{w,n}\left(\lambda \right)$ where the initial conditions are wiped out at infinity as ${\widehat{K}}_{w,n}\left(\lambda \right){|}_{\lambda =0}=1$. Assuming that the infection rates are constant in the endemic equilibrium we have ${\mathcal{A}}_{w,n}\left(\lambda \right)\sim {\mathcal{A}}_{w,n}(\infty )/\lambda$, ($\lambda \to 0$) and arrive at

$$\begin{array}{ccc}{\displaystyle J_w(\infty )}& {\displaystyle ={\mathcal{A}}_w(\infty )〈t_I^w〉,({\mathcal{A}}_w(\infty )={\beta }_w{S}_w(\infty )J_n(\infty ))}\hfill & \\ \\ {\displaystyle J_n(\infty )}& {\displaystyle ={\mathcal{A}}_n(\infty )〈t_I^n〉,({\mathcal{A}}_n(\infty )={\beta }_n{S}_n(\infty )J_w(\infty ))}\hfill & \\ & \end{array}$$

(31)

thus

$$\begin{array}{cc}{\displaystyle \frac{J_w(\infty )}{1-J_w(\infty )}-{\beta }_w\langle t_I^w\rangle J_n(\infty )=0}& \\ \\ {\displaystyle \frac{J_n(\infty )}{1-J_n(\infty )}-{\beta }_n\langle t_I^n\rangle J_w(\infty )=0.}\end{array}$$

(32)

One can see that the globally healthy state $J_{w,n}\left(0\right)=0$ is also a solution of this equation. Beside that, only solutions $J_n(\infty ),J_w(\infty )\in (0,1)$ correspond to an endemic equilibrium. One gets

$$\begin{array}{c}{\displaystyle J_w(\infty )=J_w^e=\frac{{\beta }_w{\beta }_n\langle t_I^w\rangle \langle t_I^n\rangle -1}{{\beta }_n\langle t_I^n\rangle [1+{\beta }_w\langle t_I^w\rangle ]}=\frac{R_0-1}{R_0}\frac{{\beta }_w\langle t_I^w\rangle }{1+{\beta }_w\langle t_I^w\rangle }}\\ \\ {\displaystyle J_n(\infty )=J_n^e=\frac{{\beta }_w{\beta }_n\langle t_I^w\rangle \langle t_I^n\rangle -1}{{\beta }_w\langle t_I^w\rangle [1+{\beta }_n\langle t_I^n\rangle ]}=\frac{R_0-1}{R_0}\frac{{\beta }_n\langle t_I^n\rangle }{1+{\beta }_n\langle t_I^n\rangle }}\end{array}$$

(33)

for the endemic equilibrium which is independent of the initial conditions $J_{w,n}\left(0\right)$. It depends only on the phenomenological rate constants ${\beta }_{w,n}$ and the mean infection time spans $\langle t_I^{w,n}\rangle$. We point out that the endemic equilibrium (33) has exchange symmetry $w\leftrightarrow n$ between walkers and nodes reflecting this feature in the system (2) of evolution equations without mortality. The endemic values $J_{w,n}^e$ are within $(0,1)$, i.e., existing only if $R_0={\beta }_w{\beta }_n\langle t_I^w\rangle \langle t_I^n\rangle >1$. We interpret $R_0$ as the basic reproduction number (average number of new infections at $t=0$ – nodes or walkers – due to one infected node or walker at $t=0$). That this is really the appropriate interpretation can be seen by the following somewhat rough consideration of the infection rates at $t=0$. Assume we have initially one single infected node $J_n\left(0\right)=1/N$ and no infected walkers $S_w\left(0\right)=1$. The expected number of walkers infected by this first infected node during its infectious period $t_I^n$ is

$$\langle Z_I\left(t_I^n\right)\rangle \sim Z{\mathcal{A}}_w\left(0\right)\langle t_I^n\rangle =Z\langle t_I^n\rangle {\beta }_w/N\sim Z\langle J_w\left(t_I^w\right)\rangle .$$

The number of nodes getting infected by these $\langle Z_I\left(t_I^n\right)\rangle$ infected walkers during their infectious time $t_I^w$ is

$$N_I(\langle t_I^n\rangle +\langle t_I^w\rangle )\sim N\langle {\mathcal{A}}_n\left(t_I^n\right)\rangle \langle t_I^w\rangle \sim N{\beta }_n\langle J_w\left(t_I^w\right)\rangle ={\beta }_n{\beta }_w\langle t_I^n\rangle \langle t_I^w\rangle =R_0.$$

Hence $R_0$ is the average number of infected nodes caused by the first infected node (with zero initially infected walkers). Due to the exchange symmetry (nodes ↔ walkers) of the infection rates, this argumentation remains true when we start with one infected walker and no infected nodes.

We infer that the globally healthy state is unstable for $R_0>1$ where the endemic equilibrium (33) is a unique stable fixed point and attractor for all initial conditions $J_{w,n}\left(0\right)$. We will confirm this in the next section by a linear stability analysis of the globally healthy state. The stability of the endemic state is demonstrated in the next section with Appendix A.2.

Remarkable is the limit ${\beta }_w\langle t_I^w\rangle \to \infty$ (while ${\beta }_n\langle t_I^n\rangle$ are kept constant) where all walkers become infected $J_w^e\to 1$ but not all nodes $J_n^e\to \frac{{\beta }_n\langle t_I^n\rangle }{1+{\beta }_n\langle t_I^n}<1$ and vice versa. We plot the endemic state in Figure 2 versus $R_0$ where positive values for $J_{w,n}(\infty )$ occur only for $R_0>1$ which correspond to endemic states. Next we perform a linear stability analysis of the endemic and healthy state where we will indeed identify $R_0$ as crucial control parameter.

4. Stability Analysis of Endemic and Healthy State without Mortality

Here we are interested in the condition of spreading for zero mortality, or equivalently in the condition for which the globally healthy state (endemic state) is unstable (stable). To check stability of the endemic fixed point $S_e^w=1-J_e^w,J_e^w$, $S_e^n=1-J_e^n,J_e^n$ we set

$$\begin{array}{ccc}{S}_w\left(t\right)& =S_w^e+u_w{e}^{\mu t},J_w\left(t\right)=J_w^e-u_w{e}^{\mu t}\hfill & \\ \\ S_n\left(t\right)& =S_n^e+u_n{e}^{\mu t},J_n\left(t\right)=J_n^e-u_n{e}^{\mu t}\hfill \end{array}$$

(34)

where $u_w,u_n$ are `small’ constant amplitudes. This ansatz accounts for the constant populations of nodes and walkers. Then we have for the infection rates up to linear orders in the amplitudes

$$\begin{array}{ccc}{\mathcal{A}}_w\left(t\right)& ={\beta }_w{S}_w\left(t\right)J_n\left(t\right)={\beta }_w{S}_w^e{J}_n^e+{\beta }_w(u_w{J}_n^e-u_n{S}_w^e)e^{\mu t}\hfill & \\ \\ {\mathcal{A}}_n\left(t\right)& ={\beta }_n{S}_n\left(t\right)J_w\left(t\right)={\beta }_n{S}_n^e{J}_w^e+{\beta }_n(u_n{J}_w^e-u_w{S}_n^e)e^{\mu t}\hfill \end{array}$$

(35)

Plugging these relations in our evolution equations (2) without mortality, omitting two redundant equations leads to the system

$$\left[\begin{array}{cc}\mu +{\beta }_w{J}_n^e[1-{\widehat{K}}_I^w\left(\mu \right)];\hfill & -{\beta }_w{S}_w^e[1-{\widehat{K}}_I^w\left(\mu \right)]\hfill \\ \\ -{\beta }_n{S}_n^e[1-{\widehat{K}}_I^n\left(\mu \right)];\hfill & \mu +{\beta }_n{J}_w^e[1-{\widehat{K}}_I^n\left(\mu \right)]\hfill \end{array}\right]·\left(\begin{array}{c}{u}_w\hfill \\ \\ u_n\hfill \end{array}\right)=\left(\begin{array}{c}0\hfill \\ \\ 0\hfill \end{array}\right)$$

(36)

where we have used $\langle e^{-\mu t_I^{w,n}}\rangle ={\widehat{K}}_I^{w,n}\left(\mu \right)$ and the cases of $\delta$-distributed $t_I^{w,n}$ are contained for ${\widehat{K}}_I^{w,n}\left(\mu \right)=e^{-\mu t_I^{w,n}}$. We point out that in ansatz (35) we relax causality i.e., we admit ${\mathcal{A}}_{w,n}(t-\tau )\ne 0$ for $t-\tau <0$ thus

$$\langle e^{\mu (t-t_I^{w,n})}\rangle =e^{\mu t}\langle e^{-\mu t_I^{w,n}}\rangle =e^{\mu t}{\widehat{K}}_I^{w,n}\left(\mu \right).$$

(37)

The solvability of this matrix equation requires the determinant to vanish leading to a transcendental characteristic equation for $\mu$

$${\mu }^2+\mu \left({\beta }_w{J}_n^e[1-{\widehat{K}}_I^w\left(\mu \right)]+{\beta }_n{J}_w^e[1-{\widehat{K}}_I^n\left(\mu \right)]\right)+{\beta }_w{\beta }_n[1-{\widehat{K}}_I^w\left(\mu \right)][1-{\widehat{K}}_I^n\left(\mu \right)](J_n^e{J}_w^e-S_n^e{S}_w^e)=0.$$

(38)

Generally, a fixed point is unstable if solutions with positive real part of $\mu$ exist. Consider this first for the globally healthy state $J_n=0,J_w=0$ for which Equation (38) reads

$$G\left(\mu \right)=1-{\beta }_w{\beta }_n\frac{[1-{\widehat{K}}_I^w\left(\mu \right)]}{\mu }\frac{[1-{\widehat{K}}_I^n\left(\mu \right)]}{\mu }=1-{\beta }_w{\beta }_n{\widehat{\mathsf{\Phi }}}_I^w\left(\mu \right){\widehat{\mathsf{\Phi }}}_I^n\left(\mu \right)=0$$

(39)

where we notice that $\frac{[1-{\widehat{K}}_I^{w,n}\left(\mu \right)]}{\mu }={\widehat{\mathsf{\Phi }}}_I^{w,n}\left(\mu \right)$ are the LTs of the persistence probabilities of the infection time spans. Consider this equation for $\mu \to 0$ and take into account (30) we arrive at

$$G\left(0\right)=1-{\beta }_w{\beta }_n\langle t_I^w\rangle \langle t_I^n\rangle .$$

(40)

We observe that $G\left(0\right)<0$ for $R_0={\beta }_n{\beta }_w\langle t_I^w\rangle \langle t_I^n\rangle >1$. On the other hand, we have for $\mu \to \infty$ that ${\widehat{\mathsf{\Phi }}}_I^{w,n}\left(\mu \right)\to 0$ and hence

$$G(\infty )=1.$$

(41)

One can hence infer from complete monotony of ${\widehat{\mathsf{\Phi }}}_I^{w,n}\left(\mu \right)$ and therefore of ${\widehat{\mathsf{\Phi }}}_I^w\left(\mu \right){\widehat{\mathsf{\Phi }}}_I^n\left(\mu \right)$ (see Appendix A.2, Equation (A9) for a precise definition), that $\frac{d}{d\mu }G\left(\mu \right)>0$ ($\mu \ge 0$) thus $G\left(\mu \right)=0$ has one single positive zero only if $G\left(0\right)<0$, i.e., for $R_0>1$ which therefore is the condition of instability of the healthy state (spreading of the disease). Conversely, for $R_0<1$ the healthy state turns into a stable fixed point where there is no spreading of the disease. In particular, the healthy state is always unstable ($R_0=\infty$) if at least one of mean infection time spans $\langle t_I^{w,n}\rangle =\infty$. This is true for fat-tailed kernels scaling as $K_I^{w,n}\left(t\right)\propto t^{-\alpha -1}$ ($\alpha \in (0,1)$) for $t\to \infty$. We consider such a distribution briefly in subsequent section. We plot function $G\left(\mu \right)$ versus $\mu$ for different $R_0$ for Gamma-distributed $t_I^{w,n}$ in Figure 3.

Now we consider the stability of the endemic state with $G_e\left(\mu \right)=0$ where from (38) this function reads

$$G_e\left(\mu \right)=1-{\beta }_w{\beta }_n{\widehat{\mathsf{\Phi }}}_I^w\left(\mu \right){\widehat{\mathsf{\Phi }}}_I^n\left(\mu \right)+{\beta }_w{J}_n^e{\widehat{\mathsf{\Phi }}}_I^w\left(\mu \right)+{\beta }_n{J}_w^e{\widehat{\mathsf{\Phi }}}_I^n\left(\mu \right)+{\beta }_w{\beta }_n(J_w^e+J_n^e){\widehat{\mathsf{\Phi }}}_I^w\left(\mu \right){\widehat{\mathsf{\Phi }}}_I^n\left(\mu \right)$$

(42)

with

$$\begin{array}{ccc}{\displaystyle G_e\left(0\right)}& {\displaystyle =1-R_0+{\beta }_w\langle t_I^w\rangle J_n^e+{\beta }_n\langle t_I^n\rangle J_w^e+(J_w^e+J_n^e)R_0}\hfill & \\ \\ & {\displaystyle =R_0-1.}\hfill \end{array}$$

(43)

On the other hand, $G_e(\infty )=1$ (as ${\widehat{\mathsf{\Phi }}}_I^{w,n}(\infty )=0$) and from monotony of $G_e\left(\mu \right)$ follows that there is no positive solution of $G_e\left(\mu \right)=0$ for $R_0>1$. We plot $G_e\left(\mu \right)$ in Figure 4 for different values of $R_0$ and Gamma distributed $t_I^{w,n}$. In Appendix A.2 we complete the analytical proof that $G_e\left(\mu \right)>0$ for $R_0>1$.

5. Stability Analysis of the Healthy State with Mortality

An important question is, how mortality does modify the instability of the healthy state and the basic reproduction number. To shed light on this question we perform a linear stability analysis of the healthy state $S_{w,n}=1$ where we set

$$S_w\left(t\right)=1+a{e}^{\mu t},J_w\left(t\right)=b{e}^{\mu t},d_w\left(t\right)=-(a+b)e^{\mu t},S_n\left(t\right)=1-c{e}^{\mu t},J_n\left(t\right)=c{e}^{\mu t}$$

(44)

with ${\mathcal{A}}_w\left(t\right)={\beta }_{w}c{e}^{\mu t}$ and ${\mathcal{A}}_n\left(t\right)={\beta }_{n}b{e}^{\mu t}$. Plugging this ansatz for $\mu \ge 0$ into three independent equation of (2), say the first, third and fourth one, and performing the averages (relaxing causality as previously) we arrive at

$$\left[\begin{array}{ccc}\mu ;\hfill & 0;\hfill & {\beta }_w[1-{\widehat{b}}_r\left(\mu \right)]]\hfill \\ \\ \mu ;\hfill & \mu ;\hfill & {\beta }_w{\widehat{b}}_d\left(\mu \right)\hfill \\ \\ 0\hfill & -{\beta }_n[1-{\overline{K}}_I^n\left(\mu \right)];\hfill & \mu \hfill \end{array}\right]·\left(\begin{array}{c}a\hfill \\ \\ b\hfill \\ \\ c\hfill \end{array}\right)=\left(\begin{array}{c}0\hfill \\ \\ 0\hfill \\ \\ 0\hfill \end{array}\right).$$

(45)

Putting the determinant of the matrix to zero yields the condition

$${\mu }^2-{\beta }_n{\beta }_w[1-{\overline{K}}_I^n\left(\mu \right)][1-{\widehat{b}}_r\left(\mu \right)-{\widehat{b}}_d\left(\mu \right)]=0$$

(46)

where the LTs ${\widehat{b}}_r\left(\mu \right),{\widehat{b}}_d\left(\mu \right)$ of the DPDFs $b_{r,d}\left(t\right)$ defined in (7) come into play. We are interested under which conditions there is a positive solution (instability of healthy state) of (46). Since ${\widehat{b}}_r\left(0\right)=\mathcal{R}(\infty )$ and ${\widehat{b}}_d\left(0\right)=\mathcal{D}(\infty )$ with $\mathcal{R}(\infty )+\mathcal{D}(\infty )=1$ we see that $\mu =0$ is solution of (46). Recall from (7) that $b_d\left(t\right)+b_r\left(t\right)={\mathcal{K}}_{I,M}^w\left(t\right)$ is the (properly normalized) PDF (8). Condition (46) then reads

$$G_M\left(\mu \right)=1-{\beta }_n{\beta }_w\frac{[1-{\overline{K}}_I^n\left(\mu \right)]}{\mu }\frac{[1-{\widehat{\mathcal{K}}}_{I,M}^w\left(\mu \right)]}{\mu }=0$$

(47)

where $\frac{[1-{\widehat{\mathcal{K}}}_{I,M}^w\left(\mu \right)]}{\mu }$ is the LT of the persistence probability ${\mathsf{\Phi }}_M\left(t\right){\mathsf{\Phi }}_I^w\left(t\right)$ of the walker’s infection, i.e., the probability that $t<min(t_I^w,t_M)$ (see Remark I). For zero mortality we have ${\mathcal{K}}_{I,M}^w=K_I^w$, ($b_r=K_I^w$ and $b_d=0$) retrieving condition (39). The mean sojourn time in compartment I with mortality yields

$$\langle min(t_I^w,t_M)\rangle =\langle t_{IM}^w\rangle =\frac{[1-{\widehat{\mathcal{K}}}_{I,M}^w\left(\mu \right)]}{\mu }{|}_{\mu =0}={\int }_0^{\infty }t{\mathcal{K}}_{I,M}^w\left(t\right)\mathrm{d}t={\int }_0^{\infty }{\mathsf{\Phi }}_M\left(t\right){\mathsf{\Phi }}_I^w\left(t\right)\mathrm{d}t\le {\int }_0^{\infty }{\mathsf{\Phi }}_I^w\left(t\right)\mathrm{d}t=\langle t_I^w\rangle$$

(48)

where we arrive at

$$G_M\left(0\right)=1-{\beta }_n{\beta }_w\langle t_I^n\rangle \langle t_{IM}^w\rangle .$$

(49)

Relation (48) shows that $\langle t_{IM}^w\rangle \le \langle t_I^w\rangle$ (equality only for zero mortality). On the other hand we have $G_M(\infty )=1$, so there is a positive solution of $G_M\left(\mu \right)=0$ only if

$$R_M={\beta }_n{\beta }_w\langle t_I^n\rangle \langle t_{IM}^w\rangle >1$$

(50)

where $R_M$ is the basic reproduction number modified by mortality with $R_M\le R_0$ (equality only for zero mortality). To visualize the effect of mortality on the instability of the healthy state we plot $G_M\left(\mu \right)$ for a few values of $R_M$ in Figure 5. Increasing mortality turns an unstable healthy state into a stable one.

In the random walk simulations we deal with Gamma distributed $t_I^{w,n},t_M$ where the persistence probabilities are then normalized incomplete Gamma functions (27). To explore the effect of mortality for such cases, we determine $R_M$ by numerical integration of (48) as a function of the mortality rate parameter ${\xi }_M$ and plot the result in Figure 6 where one can see that $R_M$ is monotonous decreasing with mortality rate ${\xi }_M$. We also include a case of a fat-tailed (Mittag-Leffler) distributed $t_I^w$ which we discuss hereafter. The parameters in Figure 6 are such that the zero mortality case occurs with $R_0=1$ as the upper bound. The essential feature is that $R_M$ decays monotonically with increasing mortality rate parameter ${\xi }_M$ approaching zero for ${\xi }_M\to \infty$. Diseases with high mortality stabilize the healthy state even for $\langle t_I^w\rangle \to \infty$.

Consider briefly the case where $t_M$ is exponentially distributed (i.e., ${\alpha }_M=1$ in the Gamma distribution of $t_M$) with ${\mathsf{\Phi }}_M\left(t\right)=e^{-{\xi }_{M}t}$. Then we have $\langle t_{IM}^w\rangle ={\widehat{\mathsf{\Phi }}}_I^w\left({\xi }_M\right)$ thus

$$R_M={\beta }_w{\beta }_n\langle t_I^n\rangle {\widehat{\mathsf{\Phi }}}_I^w\left({\xi }_M\right).$$

(51)

The zero mortality case is recovered for ${\xi }_M=0$ with ${\widehat{\mathsf{\Phi }}}_I^w\left(0\right)=\langle t_I^w\rangle$. For Gamma distributed $t_I^{w,n}$ this yields

$$R_M={\beta }_w{\beta }_n\frac{{\alpha }_I^n}{{\xi }_I^n{\xi }_M}\left(1-\frac{{\left({\xi }_I^w\right)}^{{\alpha }_I^w}}{{({\xi }_M+{\xi }_I^w)}^{{\alpha }_I^w}}\right).$$

(52)

where ${\alpha }_I^{w,n}$, ${\xi }_I^{w,n}$ are the parameters of the respective Gamma distributions of the infection times of nodes and walkers. The Markovian case where all waiting times are exponentially distributed is covered for ${\alpha }_I^{w,n}=1$ and yields

$$R_M=\frac{{\beta }_w{\beta }_n}{{\xi }_I^n({\xi }_I^w+{\xi }_M)}$$

(53)

containing the zero mortality case for ${\xi }_M=0$.

Fat-tailed distributed $t_I^w$:

Finally, an interesting case emerges if $t_I^w$ follows a fat-tailed distribution, i.e., ${\mathsf{\Phi }}_I^w\left(t\right)\propto t^{-\alpha }$ for t large ($\alpha \in (0,1)$) and $\langle t_I^w\rangle =\infty$, $R_0=\infty$. Let us have a look, how mortality is affecting this situation. Fat tailed $t_I^w$ distributions describe diseases where the infectious periods are very long and the healthy state without mortality is extremely unstable ($R_0=\infty$). Infected walkers can infect many nodes during their long infection time spans. An important case of this class is constituted by the Mittag-Leffler (ML) distribution ${\mathsf{\Phi }}_I^w\left(t\right)=E_{\alpha }(-{\xi }_I^w{t}^{\alpha })$ where $E_{\alpha }\left(\tau \right)$ indicates the Mittag-Leffler function, see [14,15] and references therein for representations and connections with fractional calculus. The ML function recovers the exponential for $\alpha =1$ ($E_1(-\tau )=e^{-\tau }$). Assuming exponential mortality ${\mathsf{\Phi }}_M\left(t\right)=e^{-{\xi }_{M}t}$ one obtains with (51)

$$R_M={\beta }_w{\beta }_n\langle t_I^n\rangle \frac{{\left({\xi }_M\right)}^{\alpha -1}}{{\xi }_I^w+{\left({\xi }_M\right)}^{\alpha }},\alpha \in (0,1)$$

(54)

containing the LT of the ML persistence probability distribution ${\widehat{\mathsf{\Phi }}}_I^w\left(\lambda \right)={\lambda }^{\alpha -1}/({\xi }_I^w+{\lambda }^{\alpha })$. The essential feature here is that $R_M$ is weakly singular at ${\xi }_M=0$ with a monotonously decreasing ${\xi }_M^{\alpha -1}$ scaling law, where the healthy state becomes stable for mortality parameters larger as ${\xi }_M\approx 1$. We depict this case in Figure 6 for $\alpha =\xi =0.3$ (violet curve).

6. Random Walk Simulations

The remainder of our paper is devoted to test the mean field model under "real world conditions" which we mimic by $Z=Z_S\left(t\right)+Z_I\left(t\right)+Z_D\left(t\right)$ random walkers navigating independently on an undirected connected (ergodic) graph. In our simulations we focused on Barabási-Albert (BA), Erdös-Rényi (ER) and Watts-Strogatz (WS) graphs [29,42,43] (see Appendix A.3 for a brief recap) and implemented the compartments and transmission pathway for walkers and nodes outlined in Section 2. A susceptible walker gets infected with probability $p_w$ by visiting an infected node, and a susceptible node gets infected with probability $p_n$ at a visit of an infectious walker. We assume that the infection probabilities $p_{n,w}$ are constant for all nodes and walkers, respectively. They are related yet not identical with the macroscopic rate constants ${\beta }_{w,n}$. A critical issue is whether the simple bi-linear forms for the mean field infection rates (1) still capture well the complexity of the spreading in such "real world" networks. One goal of the subsequent case study is to explore this question.

We characterize the network topology by $i=1,\dots N$ nodes with the $N\times N$ adjacency matrix $\left(A_{ij}\right)$ where $A_{ij}=1$ if the pair of nodes $i,j$ is connected by an edge, and $A_{ij}=0$ if the pair is disconnected. Further, we assume $A_{ii}=0$ to avoid self-connections of nodes. We confine us to undirected networks where the edges have no direction and the adjacency matrix is symmetric. The degree $k_i$ of a node i counts the number of neighbor nodes (edges) of this node. Each walker $z=1,..,Z$ performs simultaneous independent random steps at discrete time instants $t=\Delta t,2\Delta t,\dots$ from one to another connected node. The steps from a node i to one of the neighbor nodes are chosen with probability $1/k_i$, following for all walkers the same transition matrix

$$\mathsf{\Pi }(i\to j)=\frac{A_{ij}}{k_i},z=1,\dots ,Z,i,j=1,\dots ,N$$

(55)

which is normalized ${\sum }_{j=1}^N\mathsf{\Pi }(i\to j)=1$. This is a common way to connect the network topology with simple Markovian random walks [24,42]. In the simulations the departure nodes at $t=0$ of the walkers are randomly chosen. The path of each walker is independent and not affected by contacts with other walkers or by transition events from one to another compartment.

Case study and discussion

In order to compare the epidemic dynamics of the mean field model and random walk simulations we integrate the stochastic evolution Eqs. (2) numerically as follows. We average the increments of the compartmental fractions in a generalized Monte-Carlo sense converging towards the convolutions of the right hand side of (9) where we use the Monte-Carlo convergence feature

$$\underset{n\to \infty }{lim}\frac{1}{n}\sum _{k=1}^{n}A(t-T_k)={\int }_0^{t}A(t-\tau )K\left(\tau \right)\mathrm{d}\tau$$

(56)

for random variables T drawn from PDFs $K\left(\tau \right)$. We perform this average for any time increment $\mathrm{d}t$ with respect to all involved independent random time spans $t_I^{w,n},t_M$ (see Appendix A.1) and integrate the averaged compartmental increments in a fourth order Runge-Kutta scheme (RK4). We use in the random walk simulations and the Monte-Carlo (mean field) integration exactly the same (Gamma distributed) random values (PYTHON seeds) for the $t_I^{w,n},t_M$. The values of the infection rate parameters ${\beta }_{w,n}$ used in the mean field integration are determined from Equation (32) by plugging in the large time asymptotic values of the random walk simulation with identical parameters (without mortality). The compartmental fractions in the random walk simulations are determined by simply counting the compartmental populations at each time increment $\Delta t$ of walker’s steps. The so determined rate parameters ${\beta }_{w,n}$ plugged into the mean field integration depend in a complex manner on the infection probabilities $p_{w,n}$ and topology of the network. In this way this information is also contained in the basic reproduction numbers with and without mortality.

We explore the spreading in random graphs of different complexity such as represented in Figure 7. The BA graph is small world with power law distributed degree (Appendix A.3) which means that there are many nodes having a few connections, and a few (hub) nodes with a huge number of connections. The average distance between nodes becomes small, as it is sufficient that almost every node is only a few links away from a hub node. The ER graph is small world due to a broad degree distribution. The WS graph with the choice of connectivity parameter $m=2$ in Figure 7 has long average distances and is large world. Intuitively, one infers that a small world structure is favorable for spreading processes, a feature which was already demonstrated in the literature [9,10]. In our simulations spreading in network architectures with increased connectivity comes along with increased values of $R_0$ and $R_M$, respectively.

We identify the starting time instant ($t=0$) of the evolution in the mean field model with the time instant of the first infection of a walker in the random walk simulations. In all cases we start with a small number of randomly chosen initially infected nodes $N_I\left(0\right)=10\ll N$ ($N_I\left(0\right)\approx 10$) and no infected or dead walkers. To reduce the numbers of parameters and to concentrate on topological effects we have put in all simulations the transmission probabilities $p_w=p_n=1$. We refer to [46] for the PYTHON codes2 and animated simulation videos related to the present study. In order to visualize a typical spreading process, we depict in Figure 8 a few snapshots in a Watt-Strogatz graph with rather high overall mortality probability of $\mathcal{D}(\infty )\approx 16\%$. In this case a single infection wave emerges where a large part of walkers gets repeatedly infected increasing their probability to die. This leads to a very high fraction of eventually dead walkers $d_w(\infty )\approx 99\%$ and small fraction $S_w(\infty )\approx 1\%$ of surviving walkers corresponding to the stationary state (17) which is taken as soon as the disease gets extinct $J_w=J_n=0$. Figure 8 shows that first the infection gains large parts of the network consistent to the large value of $R_M$ observed in this case. After the first wave the disease gets extinct by the high mortality of the walkers. A disease with a similar high mortality characteristics is for instance Pestilence. The process of Figure 8 is visualized in an animated video. Figure 9 and Figure 10 show the evolution in WS graphs with identical parameters and Gamma distributions of $t_I^{w,n},t_M$ as in Figure 8 but with different mortality rate parameter ${\xi }_M$ and a much smaller overall mortality $\mathcal{D}(\infty )\approx 1\%$. The different network connectivity leads to different values of ${\beta }_{w,n}$ and mean field solutions in Figure 9 and Figure 10. In addition, the networks of Figure 9 and Figure 10 have different connectivity features. The graph of Figure 9 is small world (highly connected) whereas the WS graph in Figure 10 is weakly connected and large world. One observes in Figure 9 that the infection numbers exhibit strong and immediate increases followed by attenuated oscillations around the endemic equilibrium (for zero mortality) with high values $J_w^e\approx 0.9$ and $J_n^e\approx 0.95$. The basic reproduction numbers with mortality are in both graphs only slightly smaller as $R_0$. This is due to a rather small overall mortality of $\mathcal{D}(\infty )\approx 0.01$. This effect can also be seen in the small overlap of the Gamma distributions of $t_I^w$ and $t_M$ in the histogram. Recall that a small value of $\mathcal{D}(\infty )$ does not necessarily mean small $d_w(\infty )$ as this quantity depends also on the infection rates and network topology (see (17)) and is sensitive to repeated infections. repeated infections may indeed play an important role here as $\langle t_I^w\rangle =8$ is rather small.

In Figure 10 the infections of the random walk simulations are increasing slower (red curves) compared to Figure 9. The structure with higher connectivity Figure 9 shows excellent quantitative agreement of random walk and mean field solutions for the walkers and nodes capturing well the attenuated oscillations, especially for zero mortality. In the network with smaller connectivity of Figure 10 the increase of the infections is delayed compared to the mean field. On the other hand, for non-zero mortality the mean field and random walk dynamics for the walkers diverge slightly with time. We infer that mortality may deviate the infection rates from (1).

The comparison of the spreading in Figure 9 and Figure 10 shows clearly the role of the connectivity: The mean field model captures better the spreading in networks with higher connectivity (short average distances between nodes) and with low mortality. The following cases give further evidence for these observations.

Next we explore the spreading on an ER graph in Figure 11. The agreement of random walk simulations and mean field model is impressive where this holds for both with and without mortality. One can see by the degree distribution in Figure 7 that for these connectivity parameters the graph is well connected and small world giving strong evidence that the mean field approach is here well capturing the spreading dynamics.

Finally we explore in Figure 12 the dynamics on a BA network. In the right frame we have high overall mortality of $\mathcal{D}(\infty )\approx 10\%$ probability for a walker to die from an infection. In this example the disease is starting to spread as $R_M\approx 3.48>1$ where only a single infection wave emerges which is extinct by the high mortality. Recall that that $R_M>1$ is only telling us that the healthy state is unstable, i.e., that the disease is starting to spread. It does not contain the information whether the spreading is persistent or whether the disease is eventually extinct. To explore the role of topological features such as average distances between nodes we perform the same simulation experiment with identical parameters and less ($N=2100$) nodes, i.e., higher density of walkers (Figure 13).

The accordance of mean field model and random walk simulation is also in Figure 13 indeed excellent. We explain this by the fact that the BA network is a strongly connected structure with pronounced small world property. The higher density of walkers lead to increased $R_M$ and $R_0$ compared to Figure 12. There is also only a single infection wave occurring with a higher maximum value compared to Figure 12. In both cases (Figure 12 and Figure 13, right frames) the infection waves are extinct by the high mortality of walkers where stationary states (17) with $d_w(\infty )\approx 80\%$ of dead walkers are taken. When we switch off mortality (left frames), stable endemic states emerge more rapidly in Figure 13 (case with higher density of walkers).

Further simulation experiments (not shown here) reveal that the mean field model and random walk simulations exhibit excellent accordance when we further increase the attachment parameters m or the density of walkers with otherwise identical parameters. For higher mortality the agreement becomes less well and diverges with increasing observation time. This observation suggests that mortality modifies the infection rates in the network for larger observation times. We leave this issue for future research.

Our overall finding from this case study is that the mean field approach (with infection rates (1)) is particularly well suited to mimic spreading in strongly connected environments with pronounced small world feature, but is less well for higher mortality.

7. Conclusions

We studied epidemic spreading in complex graphs where we focused on transmission pathway via vectors mimicking the spreading of a certain class of diseases such as Dengue, Malaria or Pestilence and others. We developed a stochastic compartment model for the walkers and nodes with mortality for the walkers. For zero mortality we obtained the endemic equilibrium in explicit form (Equations (33)). Stability analysis of the endemic and healthy states reveal the crucial control parameter for spreading, the basic reproduction number. We obtained the basic reproduction numbers $R_M$ and $R_0$ with and without mortality, respectively, where we proved that $R_M\le R_0$, (relations (50) with (48)). For $R_M,R_0>1$ the healthy state is an unstable fixed point where the endemic equilibrium exists for zero mortality as a unique stable fixed point independent of the initial conditions. The basic reproduction numbers depend on the means of the compartmental sojourn times in compartment I of the nodes and walkers and on the topology of the network captured by the mean field rate constants ${\beta }_{w,n}$.

Our model has applications beyond epidemic dynamics, for instance in chemical reaction models [45]. An interesting question to explore in a follow-up project is whether our class of compartment models with indirect transmission pathways may exhibit (for zero mortality) persistent oscillations (Hopf instabilities) (See our brief discussion at the end of Appendix A.2.). An open problem also remains how large world network topology may be included into such a mean field model, modifying the infection rates. Further promising directions include an account for immune and incubation compartments, effects of mitigation measures, vaccination among many others.

Appendix A.1. Some Basic Notions

Here we recall briefly some basic notions used in the paper. The infection rates are a causal functions

$$\mathcal{A}\left(t\right)=A\left(t\right)\mathsf{\Theta }\left(t\right)$$

(A1)

where $\mathsf{\Theta }\left(t\right)$ indicates the Heaviside step function defined by

$$\mathsf{\Theta }\left(\zeta \right)=\left\{\begin{array}{c}1,\zeta \ge 0\hfill \\ \\ 0,\zeta <0\hfill \end{array}\right.$$

(A2)

with $\mathsf{\Theta }\left(\zeta \right)+\mathsf{\Theta }(-\zeta )=1$ and in or definition $\mathsf{\Theta }\left(0\right)=1$. Its derivative yields the Dirac $\delta$-distribution $\frac{d}{d\zeta }\mathsf{\Theta }\left(\zeta \right)=\delta \left(\zeta \right)$. We use throughout the paper mutually independent strictly positive random variables $T_1,\dots ,T_n\in {\mathbb{R}}_{+}$ The random variables $T_j$ are assumed to follow their specific PDFs

$$Prob[T_j\in [u,u+\mathrm{d}u]]=〈\delta (u-T_j)〉=K_j\left(u\right)\mathrm{d}u$$

(A3)

where the PDFs $K_j$ are causal functions as a consequence of the positiveness of the $T_j$. Then applies the averaging rule

$$〈f(t;T_1,T_2,\dots ,T_n)〉={\int }_0^{\infty }\dots {\int }_0^{\infty }\mathrm{d}{t}_1\dots \mathrm{d}{t}_{n}f(t;t_1,t_2,\dots ,t_n)K_1\left(t_1\right)\dots K_n\left(t_n\right)$$

(A4)

for suitable functions f. For $f(t;T_1,T_2,\dots ,T_n)=g_1\left(T_1\right)\dots g_n\left(T_n\right)$ using independence of the $T_j$ this yields

$$〈g_1\left(T_1\right)\dots g_n\left(T_n\right)〉=〈g_1\left(T_1\right),〉\dots 〈g_n\left(T_n\right)〉.$$

Important cases emerge by applying (A4) to exponentials

$$〈e^{-\lambda (T_1+\dots +T_n)}〉={\int }_0^{\infty }{e}^{-\lambda t}〈\delta (t-T_1-\dots -T_n)〉\mathrm{d}t={\widehat{K}}_1\left(\lambda \right)\dots {\widehat{K}}_n\left(\lambda \right),\Re \left\{\lambda \right\}\ge 0$$

(A5)

In this relation the LTs of the PDFs come into play

$${\widehat{K}}_j\left(\lambda \right)={\int }_0^{\infty }{e}^{-\lambda t}{K}_j\left(t\right)\mathrm{d}t$$

(A6)

where ${\widehat{K}}_j\left(\lambda \right){|}_{\lambda =0}=1$ reflects the normalization of PDFs (A3). A further observation is

$$〈\delta (t-T_1-\dots -T_n)〉=(K_1★\dots ★K_n)\left(t\right)$$

(A7)

where ★ stands for convolution $(K_1★K_2)\left(t\right)={\int }_0^t{K}_1\left(\tau \right)K_2(t-\tau )\mathrm{d}\tau$ of the causal PDFs.

Appendix A.2. Proof of Stability of the Endemic Equilibrium

Here we develop the rest of the proof stability of the endemic equilibrium, i.e., we show that the function (42) is strictly positive for $\mu >0$ with $R_0-1>0$. First, we observe that ${\widehat{\mathsf{\Phi }}}_I^{w,n}\left(\mu \right)\le \langle t_I^{w,n}\rangle$ with ${\widehat{\mathsf{\Phi }}}_I^{w,n}\left(0\right)=\langle t_I^{w,n}\rangle$ and ${\widehat{\mathsf{\Phi }}}_I^{w,n}(\infty )=0$. For our convenience, we introduce the functions

$${\lambda }_{w,n}\left(\mu \right)=\frac{{\widehat{\mathsf{\Phi }}}_I^{w,n}\left(\mu \right)}{\langle t_I^{w,n}\rangle }\in (0,1]$$

(A8)

which the LTs of the normalized ${\mathsf{\Phi }}_I^{w,n}\left(t\right)/\langle t_I^{w,n}\rangle$ and which are by virtue of Bernstein’s theorem [44] completely monotonic (CM) with respect to $\mu$, i.e.,

$${(-1)}^n\frac{d^n}{d{\mu }^n}{\lambda }_{w,n}\left(\mu \right)\ge 0,\mu \in (0,\infty )$$

(A9)

inheriting this feature from the exponential $e^{-\mu \tau }$ ($t>0$). Therefore,

$$\frac{d}{d\mu }{\lambda }_{w,n}\left(\mu \right)=-\frac{1}{\langle t_I^{w,n}\rangle }{\int }_0^{\infty }{e}^{-\mu t}t{\mathsf{\Phi }}_I^{w,n}\left(t\right)\mathrm{d}t<0$$

(as ${\mathsf{\Phi }}_I^{n,w}\left(t\right)\in (0,1]$) exists thus ${\lambda }_{w,n}\left(\mu \right)$ is monotonously decreasing with $\mu$ with ${\lambda }_{w,n}\left(0\right)=1\ge {\lambda }_{w,n}\left(\mu \right)>0$. Further we observe in Eqs. (33) that $J_n^e{\beta }_w\langle t_I^w\rangle =\frac{R_0-1}{1+{\beta }_n\langle t_I^n\rangle }$, $J_w^e{\beta }_n\langle t_I^n\rangle =\frac{R_0-1}{1+{\beta }_w\langle t_I^w\rangle }$ thus

$$R_0(J_n^e+J_w^e)=(R_0-1)\left(\frac{{\beta }_n\langle t_I^n\rangle }{1+{\beta }_n\langle t_I^n\rangle }+\frac{{\beta }_w\langle t_I^w\rangle }{1+{\beta }_w\langle t_I^w\rangle }\right)$$

Then $G_e\left(\mu \right)$ reads

$$\begin{array}{ccc}{\displaystyle G_e\left(\mu \right)}& {\displaystyle =1+{\lambda }_w\left(\mu \right){\lambda }_n\left(\mu \right)\left\{-R_0+(R_0-1)\left(\frac{{\beta }_n\langle t_I^n\rangle }{1+{\beta }_n\langle t_I^n\rangle }+\frac{{\beta }_w\langle t_I^w\rangle }{1+{\beta }_w\langle t_I^w\rangle }\right)\right\}}\hfill & \\ \\ & {\displaystyle +(R_0-1)\left(\frac{{\lambda }_w\left(\mu \right)}{1+{\beta }_n\langle t_I^n\rangle }+\frac{{\lambda }_n\left(\mu \right)}{1+{\beta }_w\langle t_I^w\rangle }\right)}\hfill \end{array}$$

(A10)

Now we observe that ${\lambda }_w\left(\mu \right){\lambda }_n\left(\mu \right)\le {\lambda }_{w,n}\left(\mu \right)$ thus a lower bound function $H\left(\mu \right)\le G_e\left(\mu \right)$ is generated by replacing ${\lambda }_{w,n}\left(\mu \right)\to {\lambda }_w\left(\mu \right){\lambda }_n\left(\mu \right)$ in the second line. Now it is sufficient to prove that $0<H\left(\mu \right)$. We hence get for this lower bound function

$$\begin{array}{ccc}H\left(\mu \right)& =1+{\lambda }_w\left(\mu \right){\lambda }_n\left(\mu \right)\left\{-R_0+(R_0-1)\left(\frac{{\beta }_n\langle t_I^n\rangle }{1+{\beta }_n\langle t_I^n\rangle }+\frac{{\beta }_w\langle t_I^w\rangle }{1+{\beta }_w\langle t_I^w\rangle }\right)+(R_0-1)\left(\frac{1}{1+{\beta }_n\langle t_I^n\rangle }+\frac{1}{1+{\beta }_w\langle t_I^w\rangle }\right)\right\}\hfill & \\ \\ & {\displaystyle =1+{\lambda }_w\left(\mu \right){\lambda }_n\left(\mu \right)(R_0-2)}\hfill & \\ \\ & {\displaystyle =1-{\lambda }_w\left(\mu \right){\lambda }_n\left(\mu \right)+(R_0-1){\lambda }_w\left(\mu \right){\lambda }_n\left(\mu \right)}\hfill \end{array}$$

(A11)

and with $1-{\lambda }_w\left(\mu \right){\lambda }_n\left(\mu \right)\ge 0$ and $(R_0-1){\lambda }_w\left(\mu \right){\lambda }_n\left(\mu \right)>0$ it follows that $0<H\left(\mu \right)<G_e\left(\mu \right)$ concluding the proof of stability of the endemic equilibrium.

A few remarks on the possibility of oscillatory (Hopf) instabilities of the endemic equilibrium

Let us briefly explore whether an oscillatory (Hopf) instability of the endemic equilibrium is possible. To that end we write $G_e\left(\mu \right)$ as

$$G_e\left(\mu \right)=1+\sigma \widehat{g}\left(\mu \right)\ge 1-\widehat{g}\left(\mu \right)$$

(A12)

where $\sigma =\pm 1$ and $\widehat{g}\left(\mu \right)$ is a non-negative CM function (see Figure 4) with maximum value $\widehat{g}\left(0\right)=|R_0-2|$. Then the following two cases may occur.

• Case (i)$\sigma =-1$; $0<G_e\left(0\right)=R_0-1<1$ ($1<R_0<2$):

Then $\widehat{g}\left(\mu \right)$ can be represented as LT of a non-negative function $g\left(t\right)$ and consider now $\mu ={\mu }_1+i{\mu }_2$ with ${\mu }_1\ge 0$ thus the real part of $\widehat{g}({\mu }_1+i{\mu }_2)$ can be written as

$$\Re g({\mu }_1+i{\mu }_2)={\int }_0^{\infty }g\left(t\right)e^{-{\mu }_{1}t}cos\left({\mu }_{2}t\right)\mathrm{d}t,{\mu }_1\ge 0$$

(A13)

with $\widehat{g}\left(0\right)=2-R_0$ and clearly $-\widehat{g}\left({\mu }_1\right)\le \Re g({\mu }_1+i{\mu }_2)\le \widehat{g}\left({\mu }_1\right)$. Therefore,

$$\Re G_e({\mu }_1+i{\mu }_2)=1-\Re g({\mu }_1+i{\mu }_2)\ge 1-\widehat{g}\left({\mu }_1\right)\ge 1-\widehat{g}\left(0\right)=R_0-1>0.$$

(A14)

Hence in the range of case (i) $1<R_0<2$ there is no oscillatory (Hopf) instability of the endemic state possible3.

• Case (ii)$\sigma =+1$; $R_0-1>1$:

Here we have two pertinent ranges of $R_0$. The first is the range (a) $\widehat{g}\left(0\right)=R_0-2<1$ (i.e., $R_0<3$) and the second one (b) is $\widehat{g}\left(0\right)=R_0-2>1$. Clearly in the range (a) (A14) remains true and $\Re G_e({\mu }_1+i{\mu }_2)$ strictly positive. Hence for $1<R_0<3$ no Hopf instability is possible.

This changes in the range (b) since $\widehat{g}\left(0\right)=R_0-2>1$ thus $\Re G_e({\mu }_1+i{\mu }_2)=1+\Re g({\mu }_1+i{\mu }_2)$ may become negative. Therefore, for $R_0>3$ a Hopf instability of the endemic state becomes possible. However, the possibility that $\Re G_e({\mu }_1+i{\mu }_2)=0$ is only necessary but not sufficient for a Hopf instability. One also needs simultaneously that the imaginary part $\Im G_e({\mu }_1+i{\mu }_2)=0$ is vanishing for the same $\mu ={\mu }_1+i{\mu }_2$. In the simulations performed for this paper, we did not observe persistent oscillations. We leave the exploration of this issue for future research in a follow-up project.

Appendix A.3. A Very Brief Recap of Random Graphs

Here we recall briefly some essential features of the three classes of random graphs, which we use in the random walk simulations. For an extended outline, consult e.g., [18]. The three classes of random graph models depicted herafter are motivated by the observation that complex random network structures are encountered ubiquitously and crucially determine human and animal mobility patterns including epidemic propagation.

(i)

Erdös and Rényi (ER) graph

The ER graph is one of the most basic variants of a random graph, which was introduced in 1959 by Erdös and Rényi [20]. We use is the so-called $G(N,p)$ variant of random ER graph model (which is actually due to Gilbert [22]) which is generated as follows [16,25]. Given are N labeled nodes. Any pair of nodes is connected independently by an edge with uniform probability p. The probability $P_N\left(k\right)$ that a node has $0\le k\le N-1$ connections is given by a binomial distribution

$$P_N\left(k\right)=\left(\begin{array}{c}N-1\\ k\end{array}\right)p^k{(1-p)}^{N-1-k}\to \frac{{\langle k\rangle }^k}{k!}{e}^{-\langle k\rangle }$$

where $\langle k\rangle =(N-1)p\sim Np$ denotes the average degree. For $N\to \infty$ (while $Np$ is kept constant) the degree distribution $P_N\left(k\right)$ converges to a Poisson distribution representing the infinite graph limit of the ER $G(N,p)$-model. Therefore, $P_N\left(k\right)$ is rapidly decaying with degree k, so the number of nodes with a high number of connections is very small. In order to obtain in the $G(N,p)$-model a connected graph, it is necessary that $p>p_c=\frac{logN}{N}$ is above the percolation limit [20,24].

(ii)

Watts-Strogatz (WS) network

The WS graph model [21,25,26] starts with a ring of N nodes where each node is connected symmetrically with a number $m<<N$ to left and right neighbor nodes by an edge such that each node has $2m$ connections. In the second step, each of the connections $i,j$ of node i is replaced with probability p by a randomly chosen connection $i,k$ uniformly among other nodes avoiding self-connections and link duplication, so that each connection $i,k$ is chosen only once. There are two noteworthy limits for a WS graph. For $p=0$ (no rewiring of links) we have a regular ring with constant degree $2m$ for all nodes. In the limit $p=1$ an ER graph is emerging with probability $2m/(N-1)$ for a link. The WS graph has the small-world property (short average distances between pairs of nodes) and a high tendency to develop clusters of nodes [26].

(iii)

Barabási-Albert (BA) graph

The BA graph is generated by a preferential attachment mechanism for newly added nodes [16,17,18,19]. One starts with $m_0$ nodes and adds new nodes. Any newly added node is connected with $m\le m_0$ existing nodes (m is referred to as the attachment parameter) where most likely with nodes of high degrees. In this way nodes with high degree receive further links. This leads to an asymptotically scale-free network with a power law degree-distribution

$$P\left(k\right)\propto k^{-2-\mu },\mu \approx 1.$$

As the decrease in this power law is relatively slow, there might exist quite a few nodes with many links (hub nodes) and many nodes with few links. BA graphs are believed to mimic a large class of real-world networks including the world wide web, citation-, social-, and metabolic networks.

Realizations of these three types of random graph types used in our multiple random walkers simulations are shown in Figure 7.