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

Abstract

We study epidemic spreading in complex networks by a multiple random walker approach. Each walker performs an independent simple Markovian random walk on a complex undirected (ergodic) random graph where we focus on the Barabási–Albert (BA), Erdös–Rényi (ER), and Watts–Strogatz (WS) types. Both walkers and nodes can be either susceptible (S) or infected and infectious (I), representing their state of health. Susceptible nodes may be infected by visits of infected walkers, and susceptible walkers may be infected by visiting infected nodes. No direct transmission of the disease among walkers (or among nodes) is possible. This model mimics a large class of diseases such as Dengue and Malaria with the transmission of the disease via vectors (mosquitoes). Infected walkers may die during the time span of their infection, introducing an additional compartment D of dead walkers. Contrary to the walkers, there is no mortality of infected nodes. Infected nodes always recover from their infection after a random finite time span. This assumption is based on the observation that infectious vectors (mosquitoes) are not ill and do not die from the infection. The infectious time spans of nodes and walkers, and the survival times of infected walkers, are represented by independent random variables. We derive stochastic evolution equations for the mean-field compartmental populations with the mortality of walkers and delayed transitions among the compartments. From linear stability analysis, we derive the basic reproduction numbers $R_M,R_0$ with and without mortality, respectively, and prove that $R_M<R_0$. For $R_M,R_0>1$, the healthy state is unstable, whereas for zero mortality, a stable endemic equilibrium exists (independent of the initial conditions), which we obtained explicitly. We observed that the solutions of the random walk simulations in the considered networks agree well with the mean-field solutions for strongly connected graph topologies, whereas less well for weakly connected structures and for diseases with high mortality. Our model has applications beyond epidemic dynamics, for instance in the kinetics of chemical reactions, the propagation of contaminants, wood fires, and others.

1. Introduction

The sudden or recurrent emergence of epidemics has been an everlasting threat 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 the 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 currently a “compartment model”. In their so-called SIR model, the individuals are categorized into the compartments susceptible (S), infected (I), and recovered (immune) (R), characterizing the states of their health. While standard SIR-type models are able to capture the 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 the COVID-19 pandemic [7]. In addition to purely macroscopic models, the study of epidemic dynamics in complex networks has attracted considerable attention [8,9,10]. 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 [9]. The topological features crucial for epidemic spreading include the small world property (short average network distances (the network “distance” of two nodes is the number of edges of the shortest path connecting them) and a high clustering coefficient, measuring the existence of redundant paths between pairs of nodes [11,12].

Further interesting directions are represented by combinations of network science and stochastic compartmental models [13,14,15,16,17]. Such models include Markovian and non-Markovian approaches [18,19,20,21,22], where non-Markovianity is introduced by non-exponentially distributed sojourn times in the compartments [23,24]. 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 [25], 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 works emphasize the importance of the spatial heterogeneity effects of infection patterns in epidemic spreading phenomena [26]. The role of local clusters in generating periodic epidemic outbursts has been highlighted in the references [27,28]. A cluster model to explain periodic behavior was introduced a long time ago [29]. 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 [22].

The aim of the present paper is to study the spreading of a certain class of vector-transmitted diseases in a population of individuals (random walkers) moving on complex graphs aiming to mimic human mobility patterns in complex environments such as cities, streets, and transportation networks. Essential elements in our model are the accounting of the mortality of infected individuals (random walkers) and an indirect transmission pathway via vectors (nodes of the network). The random walkers are mimicking individuals that are navigating on the network. The nodes of the network are assumed to mimic the vectors (for instance, the mosquitoes in the case of Malaria). The reason why we chose the nodes of the network to mimic the vectors is based on the observation that, in real-world situations, the vectors live in stationary well-defined areas such as swamps and others. The nodes are immortal, since we assume that the vectors (mosquitoes) are not falling ill during their infection time span.

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 the stability of the healthy state with mortality, 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 the robustness of our mean-field model under “complex real-world conditions” by implementing its assumptions in multiple random random walker simulations on Barabási–Albert (BA)-, Erdös–Rényi (ER)-, and Watts–Strogatz (WS)-type graphs (see [14,15,30] and Appendix A.4). These graph types have different complexity and connectivity features with an 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 mosquitoes), pestilence (transmission by fleas), and others [8,31]. 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) that 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 death 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 on $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 the contact of a node and a walker, which are captured by (yet not identical to) 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 the only admitted death 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 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 first give 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)\mathrm{\Theta }(t_M-t_I^w)〉+J_w\left(0\right)\langle \delta (t-t_I^w)\mathrm{\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)\mathrm{\Theta }(t_M-t_I^w)〉-J_w\left(0\right)〈\delta (t-t_I^w)\mathrm{\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)\mathrm{\Theta }(t_I^w-t_M)〉+J_w\left(0\right)\langle \delta (t-t_M)\mathrm{\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 $〈..〉$ the average over the contained (set of independent) random variables $t_I^w,t_I^n,t_M$ outlined hereafter and in Appendix A.1. $\mathrm{\Theta }(..)$ stands for the Heaviside function (A2), and $\delta (..)$ for 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 the system (2) is as follows. The instantaneous infection rate ${\mathcal{A}}_w\left(t\right)$ governs the transitions S → I of walkers (due to the visits of infected nodes). The term $〈{\mathcal{A}}_w(t-t_I^w)\mathrm{\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 $\mathrm{\Theta }(t_M-t_I^w)=1$). Then, $〈{\mathcal{A}}_w(t-t_M)\mathrm{\Theta }(t_I^w-t_M)〉$ captures the rate of walkers infected at $t-t_M$ dying at time t during the infection time span (indicated by $\mathrm{\Theta }(t_I^w-t_M)=1$ for $t_I^w>t_M$).

Remark 1.

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 the persistence of a walker’s infection at time t, given the infection starts at time 0, is $\langle \mathrm{\Theta }(t_I^w-t)\mathrm{\Theta }(t_M-t)\rangle$ (see (7)). Note that $\mathrm{\Theta }(t_I^w-t)\mathrm{\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)\mathrm{\Theta }(t_M-t_I^w)\rangle$ or they die with rate $J_w\left(0\right)\langle \delta (t-t_M)\mathrm{\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 the 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 the system (2), we need to take a closer look at 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 the (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 the walker’s (node’s) infection (without mortality):

$${\mathrm{\Phi }}_I^{w,n}\left(t\right)=Prob(t_I^{w,n}>t)=\langle \mathrm{\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 the walker’s survival up to time t (given $t_I^w=\infty$):

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

(6)

The persistence probabilities fulfill the initial condition ${\mathrm{\Phi }}_M\left(0\right)={\mathrm{\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 ${\mathrm{\Phi }}_M(\infty )={\mathrm{\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 \mathrm{\Theta }(t_M-t)\mathrm{\Theta }(t_I^w-t)\rangle }& {\displaystyle =\langle \mathrm{\Theta }(t_M-t)\rangle \langle \mathrm{\Theta }(t_I^w-t)\rangle ={\mathrm{\Phi }}_I^w\left(t\right){\mathrm{\Phi }}_M\left(t\right)}\hfill & \\ \\ {\displaystyle b_d\left(t\right)=\langle \delta (t-t_M)\mathrm{\Theta }(t_I^w-t_M)\rangle }& {\displaystyle =\langle \delta (t-t_M)\rangle \langle \mathrm{\Theta }(t_I^w-t)\rangle =K_M\left(t\right){\mathrm{\Phi }}_I^w\left(t\right)}\hfill & \\ \\ {\displaystyle b_r\left(t\right)=\langle \delta (t-t_I^w)\mathrm{\Theta }(t_M-t_I^w)\rangle }& {\displaystyle =\langle \delta (t-t_I^w)\rangle \langle \mathrm{\Theta }(t_M-t)\rangle =K_I^w\left(t\right){\mathrm{\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 \mathrm{\Theta }(t_M-t)\rangle \langle \mathrm{\Theta }(t_I^w-t)\rangle \right]=-\frac{d}{dt}\left[{\mathrm{\Phi }}_I^w\left(t\right){\mathrm{\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 \mathrm{\Theta }(t-t_I^w)\mathrm{\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){\mathrm{\Phi }}_M\left(\tau \right)\mathrm{d}\tau }\hfill & \\ \\ {\displaystyle \mathcal{D}\left(t\right)=\langle \mathrm{\Theta }(t-t_M)\mathrm{\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){\mathrm{\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-{\mathrm{\Phi }}_I^w\left(t\right){\mathrm{\Phi }}_M\left(t\right)}\hfill & \\ \\ {\displaystyle \mathcal{D}(\infty )+\mathcal{R}(\infty )=1}& \\ \\ {\displaystyle \langle A(t-t_I)\mathrm{\Theta }(t_M-t_I)\rangle }& {\displaystyle =\langle A(t-t_I){\mathrm{\Phi }}_M\left(t_I\right)\rangle ={\int }_0^{t}A(t-\tau ){\mathrm{\Phi }}_M\left(\tau \right)K_I\left(\tau \right)\mathrm{d}\tau }\hfill & \\ \\ {\displaystyle \langle A(t-t_M)\mathrm{\Theta }(t_I^w-t_M)\rangle }& {\displaystyle =\langle A(t-t_M){\mathrm{\Phi }}_I^w\left(t_M\right)\rangle ={\int }_0^{t}A(t-\tau ){\mathrm{\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 the 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 not proper PDFs since they are not normalized to one, but rather to $\mathcal{D}(\infty ),\mathcal{R}(\infty )<1$, respectively. Consult [32] 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){\mathrm{\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){\mathrm{\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}〈\mathrm{\Theta }(t_M-t)\mathrm{\Theta }(t_I^w-t)〉=-\frac{d}{dt}\left[{\mathrm{\Phi }}_I^w\left(t\right){\mathrm{\Phi }}_M\left(t\right)\right]$$

(8)

is non-negative (as are $K_I^w=-\frac{d}{dt}{\mathrm{\Phi }}_I^w\ge 0$, $K_M=-\frac{d}{dt}{\mathrm{\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 (${\mathrm{\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). The infinite time limits are important: $\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 to $\mathcal{D}(\infty )$ also as “overall mortality”. It must not be confused with the infinite time limit of the dead walker’s fraction $d_w(\infty )$, which is different from $\mathcal{D}(\infty )$, as we will see in detail 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, the 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){\mathrm{\Phi }}_M\left(\tau \right)\mathrm{d}\tau +J_w\left(0\right)K_I^w\left(t\right){\mathrm{\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){\mathrm{\Phi }}_M(t-\tau ){\mathrm{\Phi }}_I^w(t-\tau )\mathrm{d}\tau -J_w\left(0\right)\left[K_I^w\left(t\right){\mathrm{\Phi }}_M\left(t\right)+K_M\left(t\right){\mathrm{\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 & \\ & \hfill \end{array}$$

(9)

The PDF (8) for which a walker leaves compartment I (either by recovery or by death) allows rewriting 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)\mathrm{\Theta }(t_I^w-t_M)〉+J_w\left(0\right)\langle \delta (t-t_M)\mathrm{\Theta }(t_I^w-t_M)\rangle }\hfill & \\ \\ & {\displaystyle ={\int }_0^t{\mathcal{A}}_w(t-\tau )K_M\left(\tau \right){\mathrm{\Phi }}_I^w\left(\tau \right)\mathrm{d}\tau +J_w\left(0\right)K_M\left(t\right){\mathrm{\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 \mathrm{\Theta }(\tau -t_M)\mathrm{\Theta }(t_I^w-t_M)\rangle \mathrm{d}\tau +J_w\left(0\right)\langle \mathrm{\Theta }(t-t_M)\mathrm{\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 the first equation in (2)):

$$\begin{array}{ccc}{\displaystyle r_w\left(t\right)}& {\displaystyle ={\int }_0^t{\mathcal{A}}_w(t-\tau )〈\mathrm{\Theta }(\tau -t_I^w)\mathrm{\Theta }(t_M-t_I^w)〉\mathrm{d}\tau +J_w\left(0\right)〈\mathrm{\Theta }(t-t_I^w)\mathrm{\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-{\mathrm{\Phi }}_I^w\left(\tau \right){\mathrm{\Phi }}_M\left(\tau \right)]\mathrm{d}\tau +J_w\left(0\right)\left[1-{\mathrm{\Phi }}_I^w\left(t\right){\mathrm{\Phi }}_M\left(t\right)\right].$$

(14)

Relation (12) records all death 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{\mathrm{\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[{\mathrm{\Phi }}_M\left(t\right){\mathrm{\Phi }}_I^w\left(t\right)+\mathcal{D}\left(t\right)\right]-{\int }_0^t{\mathcal{A}}_w\left(\tau \right)[{\mathrm{\Phi }}_M(t-\tau ){\mathrm{\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){\mathrm{\Phi }}_M\left(t\right){\mathrm{\Phi }}_I^w\left(t\right)+{\int }_0^t{\mathcal{A}}_w\left(\tau \right){\mathrm{\Phi }}_M(t-\tau ){\mathrm{\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){\mathrm{\Phi }}_I^n\left(t\right)+{\int }_0^t{\mathcal{A}}_n\left(\tau \right){\mathrm{\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 the unknown susceptible and infected population fractions, i.e., ${\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 ${\mathrm{\Phi }}_I^w\left(t\right){\mathrm{\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 no longer true for zero mortality ($\mathcal{D}(\infty )=0$) and considered in Section 3. We define the overall probability $P_D$ that a walker dies from ($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$.

2.1. 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 ${\mathrm{\Phi }}_I^{w,n}\left(t\right)=e^{-{\xi }_I^{w,n}t}$, ${\mathrm{\Phi }}_M\left(t\right)=e^{-{\xi }_{M}t}$, where $\langle t_I^{w,n}\rangle ={\left({\xi }_I^{w,n}\right)}^{-1}$ and $\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}& \hfill & \\ \\ {\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.}& \hfill \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}$, we recover indeed the representation of the expression (17).

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

2.2. A Few More Words on Waiting Time Distributions

In our simulations, we assumed that the time spans $t_I^w,t_I^n,t_M$ are independent random variables drawn from specific Gamma distributions. The advantage of using 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 employed the Python 3.10.8 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, the term “scale parameter” is used, $\theta ={\xi }^{-1}$) and $\Gamma \left(\alpha \right)$ stands for the Gamma function. It is worthy of mention that the Gamma distribution for $\alpha \in \mathbb{N}$ (also referred to as the Erlang distribution) is the PDF of the sum of independent and identically exponentially distributed random variables. 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); specifically, 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 subsequently use the persistence probability of the Gamma distribution (see the right frame of Figure 1):

$${\mathrm{\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) necessarily fulfills the initial condition ${\mathrm{\Phi }}_{\alpha ,\xi }\left(0\right)=1$ and is vanishing at infinity ${\mathrm{\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){\mathrm{\Phi }}_I^w\left(t\right)+{\int }_0^t{\mathcal{A}}_w(t-\tau ){\mathrm{\Phi }}_I^w\left(\tau \right)\mathrm{d}\tau }\hfill & \\ \\ {\displaystyle J_n\left(t\right)}& {\displaystyle =J_n\left(0\right){\mathrm{\Phi }}_I^n\left(t\right)+{\int }_0^t{\mathcal{A}}_n(t-\tau ){\mathrm{\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 ${\mathrm{\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 ${\mathrm{\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 }{\mathrm{\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 & \\ & \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. Besides that, only solutions $J_n(\infty ),J_w(\infty )\in (0,1)$ correspond to an endemic equilibrium. One obtains

$$\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 further mention the useful relation $R_0{S}_w(\infty )S_n(\infty )=1$ following straightforwardly from (32), connecting $R_0$ directly with the endemic state. 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 together with Appendix A.2.

The limit ${\beta }_w\langle t_I^w\rangle \to \infty$ (while ${\beta }_n\langle t_I^n\rangle$ are kept constant) is remarkable, 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 a 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 the stability of the endemic fixed point $S_w^e=1-J_w^e,J_e^w$, $S_n^e=1-J_n^e,J_n^e$, 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 Equation (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 the 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 the 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{\mathrm{\Phi }}}_I^w\left(\mu \right){\widehat{\mathrm{\Phi }}}_I^n\left(\mu \right)=0$$

(39)

where we notice that $\frac{[1-{\widehat{K}}_I^{w,n}\left(\mu \right)]}{\mu }={\widehat{\mathrm{\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{\mathrm{\Phi }}}_I^{w,n}\left(\mu \right)\to 0$, and hence,

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

(41)

One can, hence, infer from the complete monotony of ${\widehat{\mathrm{\Phi }}}_I^{w,n}\left(\mu \right)$ and, therefore, of ${\widehat{\mathrm{\Phi }}}_I^w\left(\mu \right){\widehat{\mathrm{\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 the 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 the 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 the 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{\mathrm{\Phi }}}_I^w\left(\mu \right){\widehat{\mathrm{\Phi }}}_I^n\left(\mu \right)+{\beta }_w{J}_n^e{\widehat{\mathrm{\Phi }}}_I^w\left(\mu \right)+{\beta }_n{J}_w^e{\widehat{\mathrm{\Phi }}}_I^n\left(\mu \right)+{\beta }_w{\beta }_n(J_w^e+J_n^e){\widehat{\mathrm{\Phi }}}_I^w\left(\mu \right){\widehat{\mathrm{\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{\mathrm{\Phi }}}_I^{w,n}(\infty )=0$), and from the monotony of $G_e\left(\mu \right)$, it 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 modifies 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 the three independent equations 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 in under which conditions there is a positive solution (instability of the 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 a 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 ${\mathrm{\Phi }}_M\left(t\right){\mathrm{\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 1). For zero mortality, we have ${\mathcal{K}}_{I,M}^w=K_I^w$, ($b_r=K_I^w$ and $b_d=0$), retrieving the 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 }{\mathrm{\Phi }}_M\left(t\right){\mathrm{\Phi }}_I^w\left(t\right)\mathrm{d}t\le {\int }_0^{\infty }{\mathrm{\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 monotonously 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 for the Gamma-distributed cases. 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 ${\mathrm{\Phi }}_M\left(t\right)=e^{-{\xi }_{M}t}$. Then, we have $\langle t_{IM}^w\rangle ={\widehat{\mathrm{\Phi }}}_I^w\left({\xi }_M\right)$; thus,

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

(51)

The zero mortality case is recovered for ${\xi }_M=0$ with ${\widehat{\mathrm{\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., ${\mathrm{\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 take a look at 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 ${\mathrm{\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 [33,34] and the 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 ${\mathrm{\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{\mathrm{\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 testing 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 [16,17,35] (see Appendix A.4 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 assumed that the infection probabilities $p_{n,w}$ are constant for all nodes and walkers, respectively. They are related, yet not identical to 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 the complexity of the spreading in such “real world” networks well. 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 assumed $A_{ii}=0$ to avoid self-connections of nodes. We confined ourselves 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,\dots ,Z$ performs simultaneous independent random steps at discrete time instants $t=\mathsf{\Delta }t,2\mathsf{\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:

$$\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\Pi (i\to j)=1$. This is a common way to connect the network topology with simple Markovian random walks [30,35]. 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 integrated the stochastic evolution Equation (2) numerically as follows. We averaged 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 performed 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 integrated the averaged compartmental increments in a fourth-order Runge–Kutta scheme (RK4). We used in the random walk simulations and the Monte Carlo (mean-field) integration exactly the same (Gamma-distributed) random values (Python 3.10.8 seeds) for $t_I^{w,n},t_M$. The values of the infection rate parameters ${\beta }_{w,n}$ used in the mean-field integration were 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 were determined by simply counting the compartmental populations at each time increment $\mathsf{\Delta }t$ of the 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. Indeed, the importance of the structural features (topology) of the network is crucial for the spreading of the disease. We refer to the recent review article [36] exploring a large variety of these effects. In that work, a detailed study was performed on the role of structural elements such as the average distance of the nodes and the effective network size (among others) on the epidemic spreading. Naturally, pure mean-field (compartment) approaches ignore the topological features of the environments where the diseases are spreading. In our mean-field model, the structural network properties are captured by the infection rate parameters ${\beta }_{w,n}$ as outlined above.

We explore now the spreading in random graphs of different complexity such as represented in Figure 7. The BA graph is small-world with a power law-distributed degree (Appendix A.4), 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 that was already demonstrated in the literature [8,9]. In our simulations, spreading in network architectures with increased connectivity comes along with increased values of $R_0$ and $R_M$, respectively.

We identified 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 started 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 number 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 [37] for the Python codes (free to download and non-commercial use) 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 a rather high overall mortality probability of $\mathcal{D}(\infty )\approx 16\%$. In this case, a single infection wave emerges where a large part of the 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 a small fraction $S_w(\infty )\approx 1\%$ of surviving walkers corresponding to the stationary state (17), which is taken as soon as the disease becomes extinct $J_w=J_n=0$. Figure 8 shows that, first, the infection gains large parts of the network, consistent with the large value of $R_M$ observed in this case. After the first wave, the disease becomes extinct by the high mortality of the walkers. A disease with similar high mortality characteristics is, for instance, pestilence. The process of Figure 8 is visualized in an animated video https://drive.google.com/file/d/1-fhroAsoAVDKGR5H9yWtqjD7A1ZU5pQt/view (accessed on 22 April 2024).

Figure 9 and Figure 10 show the evolution in the WS graphs with identical parameters and Gamma distributions of $t_I^{w,n},t_M$ as in Figure 8, but with a 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 WS graph of Figure 9 is small-world (strongly connected), whereas the WS graph in Figure 10 is weakly connected and large-world. Comparing Figure 9 and Figure 10 shows that the network with higher connectivity of Figure 9 has higher values of $R_0$ or $R_M$, respectively. The graph in Figure 9 is small-world with larger connectivity parameter $m=\langle k\rangle =8$ coinciding with the average degree, i.e., has shorter average distances between nodes as the graph of Figure 10, which has connectivity parameter $m=\langle k\rangle =2$ with clearly larger average distances between the nodes. The graph of Figure 9 with higher connectivity has $R_0\approx 183.3$ ($R_M\approx 183.17$), whereas the graph of Figure 10 with lower connectivity has much smaller $R_0\approx 108$ ($R_M\approx$ 107.92). Hence, the spreading (at $t=0$) in the graph with higher connectivity is much stronger as in the less well-connected structure. The effect of the smaller connectivity on the spreading can also be clearly seen by means of the delayed increase of the infection numbers produced by the random walk simulations (red curves) in Figure 10. One further 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 )$ since this quantity depends also on the infection rates and network topology (see (17)) and is sensitive to repeated infections, which 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 more slowly (red curves) compared to Figure 9. The structure with higher connectivity in Figure 9 shows excellent quantitative agreement of the random walk and mean-field solutions for the walkers and nodes, well capturing 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 better captures 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 the mean-field model is impressive, where this holds for both with and without mortality. One can see by means of the high average degree $\langle k\rangle =100$ and 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 capturing the spreading dynamics well.

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.46>1$, where only a single infection wave emerges, which is made extinct by the high mortality. Recall 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 about whether the spreading is persistent or whether the disease is eventually made extinct. To explore the role of topological features such as the average distances between nodes, we performed the same simulation experiment with identical parameters and less ($N=2100$) nodes, i.e., a higher density of walkers (Figure 13).

The accordance of the mean-field model and random walk simulation is indeed also, as seen in Figure 13, excellent. We explain this by the fact that the BA network is a strongly connected structure with a pronounced small-world property. The higher density of walkers leads 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 made 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 (the case with a higher density of walkers).

Further simulation experiments (not shown here) revealed that the mean-field model and random walk simulations exhibited excellent accordance when we further increased the attachment parameters m or the density of walkers with otherwise identical parameters. For higher mortality, the agreement becomes less good 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 a pronounced small-world feature, but is less good for higher mortality.

7. Conclusions

We studied epidemic spreading in complex graphs, where we focused on the transmission pathway via vectors mimicking the spreading of a certain class of diseases such as Dengue, Malaria, 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 (Equation (33)). The stability analysis of the endemic and healthy states revealed 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$ (the 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}$. For future research, it would be desirable to extend the mean-field model in such a way that it connects the rate parameters ${\beta }_{w,n}$ in explicit form with the topological features of the network such as the mean distance of the nodes, average connectivity, and others. We refer to [38] for a recent study connecting the spreading of COVID-19 with such topological features, where it has to be emphasized that COVID-19 is not a vector-transmitted disease and, hence, does not refer to the class of diseases analyzed here.

Our model has applications beyond epidemic dynamics, for instance in chemical reaction models [39], and can be extended in various directions. For instance, the inclusion of additional compartments and natural birth and death processes based on real-world observations of Malaria transmission dynamics is of interest [40]. 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 a large-world network topology may be included in such a mean-field model, modifying the infection rates. Further promising directions include an accounting of immune and incubation compartments, the effects of mitigation measures, and vaccinations, among many others.