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

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 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 states 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 transmission of the disease via vectors (mosquitos). Infected walkers may die during the time span of their infection introducing an additional compartment D of dead walkers. Infected nodes never die and always recover from their infection after a random finite time. This assumption is based on the observation that infectious vectors (mosquitos) 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 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 observe 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, among many others.

Keywords:

Subject: Physical Sciences - Other

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} (t), Z_{I} (t)

(

N_{S} (t), N_{I} (t)

) be the number of walkers (nodes) in compartments S and I, and

Z_{D} (t)

the non-decreasing number of walkers (in compartment D) died from the disease up to time t. We consider

Z = Z_{I} (t) + Z_{S} (t) + Z_{D} (t)

walkers (Z independent of time) and a constant number of nodes

N = N_{I} (t) + N_{S} (t)

. We assume at instant

t = 0

the first spontaneous occurrence of the disease of a few infected walkers

Z_{I} (0) ≪ Z

or nodes

N_{I} (0) ≪ N

(and no dead walkers

Z_{D} (0) = 0

). We introduce the compartmental fractions

S_{w} (t) = \frac{Z_{S} (t)}{Z}

J_{w} (t) = \frac{Z_{I} (t)}{Z}

d_{w} (t) = \frac{Z_{d} (t)}{Z}

for the walkers (normalized with respect to Z) with

S_{w} (t) + J_{w} (t) + d_{w} (t) = 1

, and

S_{n} (t) = \frac{N_{S} (t)}{N}

J_{n} (t) = \frac{N_{I} (t)}{N}

with

S_{n} (t) + J_{n} (t) = 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

A_{w} (t), A_{n} (t)

the infection rates (rates of transitions S → I) of walkers and nodes, respectively. We assume the following simple bi-linear forms

\begin{matrix} A_{w} (t) & = A_{w} [S_{w} (t), J_{n} (t)] = β_{w} S_{w} (t) J_{n} (t) \\ A_{n} (t) & = A_{n} [S_{n} (t), J_{w} (t)] = β_{n} S_{n} (t) J_{w} (t) \end{matrix}

(1)

with constant rate parameters

β_{w}, β_{n} > 0

(independent of time).

A_{w} (t)

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.

A_{n} (t)

stands for the infection rate of nodes depending on

S_{n} (t), J_{w} (t)

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

β_{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{matrix} \frac{d}{d t} S_{w} (t) & = - A_{w} (t) + 〈A_{w} (t - t_{I}^{w}) Θ (t_{M} - t_{I}^{w})〉 + J_{w} (0) 〈 δ (t - t_{I}^{w}) Θ (t_{M} - t_{I}^{w}) 〉 \\ \frac{d}{d t} J_{w} (t) & = A_{w} (t) - 〈A_{w} (t - t_{I}^{w}) Θ (t_{M} - t_{I}^{w})〉 - J_{w} (0) 〈δ (t - t_{I}^{w}) Θ (t_{M} - t_{I}^{w})〉 - \frac{d}{d t} d_{w} (t) \\ \frac{d}{d t} d_{w} (t) = & 〈A_{w} (t - t_{M}) Θ (t_{I}^{w} - t_{M})〉 + J_{w} (0) 〈 δ (t - t_{M}) Θ (t_{I}^{w} - t_{M}) 〉 \\ \frac{d}{d t} S_{n} (t) & = - A_{n} (t) + 〈A_{n} (t - t_{I}^{n})〉 + J_{n} (0) 〈 δ (t - t_{I}^{n}) 〉 \\ \frac{d}{d t} J_{n} (t) & = - \frac{d}{d t} S_{n} (t) \end{matrix}

(2)

where

\frac{d}{d t} d_{w} (t)

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.

Θ (. .)

stands for the Heaviside function (A2), and

δ (. .)

for the Dirac’s

δ

-distribution. An epidemic always starts from “natural” initial conditions

S_{w} (0) = 1

S_{n} (0) = 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} (0) δ (t)

to the infection rates of walkers and nodes, respectively. Equivalently, we introduce initial conditions

S_{w, n} (0) = 1 - J_{w, n} (0)

(

d_{w} (0) = 0

) with

J_{w, n} (0) > 0

consisting typically of a few infected walkers and/or nodes in a large susceptible population without dead walkers

d_{w} (0) = 0

The interpretation of system (2) is as follows. The instantaneous infection rate

A_{w} (t)

governs the transitions S → I of walkers (due to visits of infected nodes). The term

〈A_{w} (t - t_{I}^{w}) Θ (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

Θ (t_{M} - t_{I}^{w}) = 1

). Then

〈A_{w} (t - t_{M}) Θ (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

Θ (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

m i n (t_{I}^{w}, t_{M})

, i.e.,

t_{I}^{w}

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

〈 Θ (t_{I}^{w} - t) Θ (t_{M} - t) 〉

(see (7)). Note that

Θ (t_{I}^{w} - t) Θ (t_{M} - t) = 1

only if

t < m i n (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

m i n (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} (0) 〈 δ (t - t_{I}^{w}) Θ (t_{M} - t_{I}^{w}) 〉

or they die with rate

J_{w} (0) 〈 δ (t - t_{M}) Θ (t_{I}^{w} - t_{M}) 〉

, and nodes always recover with rate

J_{n} (0) 〈 δ (t - t_{I}^{n}) 〉

. 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

β_{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

P r o b [T \in [τ, τ + d τ] = K (τ) d τ,

(3)

with their respective PDFs (kernels)

K (τ) = {K_{I}^{w, n} (τ), K_{M} (τ)}

which are normalized

P r o b [T > 0] = \int_{0}^{\infty} K (τ) d τ = 1

. Then recall the averaging rule for (suitable) functions

f (T)

of the random variable T which we use throughout the paper

〈f (T)〉 = \int_{0}^{\infty} K (τ) f (τ) d τ,

(4)

see also Appendix A.1. An important case is

〈 δ (t - T) 〉 = K (t)

. Then by applying (4) we introduce the persistence probabilities of walker’s (node’s) infection (without mortality)

Φ_{I}^{w, n} (t) = P r o b (t_{I}^{w, n} > t) = 〈 Θ (t_{I}^{w, n} - t) 〉 = \int_{t}^{\infty} K_{I}^{w, n} (τ) d τ

(5)

and the probability of walker’s survival up to time t (given

t_{I}^{w} = \infty

)

Φ_{M} (t) = P r o b (t_{M} > t) = 〈 Θ (t_{M} - t) 〉 = \int_{t}^{\infty} K_{M} (τ) d τ .

(6)

The persistence probabilities fulfill the initial condition

Φ_{M} (0) = Φ_{I}^{w, n} (0) = 1

corresponding to the normalization of the waiting time PDFs

K (τ) = {K_{I}^{w, n} (τ), K_{M} (τ)}

and are vanishing at infinity

Φ_{M} (\infty) = Φ_{I}^{w, n} (\infty) = 0

. To evaluate the averages in (2) we will use the following quantities:

\begin{matrix} 〈 δ (t - T) 〉 & = K (t), T = {t_{I}^{w}, t_{I}^{n}, t_{M}} \\ 〈 Θ (t_{M} - t) Θ (t_{I}^{w} - t) 〉 & = 〈 Θ (t_{M} - t) 〉 〈 Θ (t_{I}^{w} - t) 〉 = Φ_{I}^{w} (t) Φ_{M} (t) \\ b_{d} (t) = 〈 δ (t - t_{M}) Θ (t_{I}^{w} - t_{M}) 〉 & = 〈 δ (t - t_{M}) 〉 〈 Θ (t_{I}^{w} - t) 〉 = K_{M} (t) Φ_{I}^{w} (t) \\ b_{r} (t) = 〈 δ (t - t_{I}^{w}) Θ (t_{M} - t_{I}^{w}) 〉 & = 〈 δ (t - t_{I}^{w}) 〉 〈 Θ (t_{M} - t) 〉 = K_{I}^{w} (t) Φ_{M} (t) \\ b_{d} (t) + b_{r} (t) = K_{I, M}^{w} (t) & = - \frac{d}{d t} [〈 Θ (t_{M} - t) 〉 〈 Θ (t_{I}^{w} - t) 〉] = - \frac{d}{d t} [Φ_{I}^{w} (t) Φ_{M} (t)] \\ \int_{0}^{\infty} K_{I, M}^{w} (t) d t & = 1 \\ R (t) = 〈 Θ (t - t_{I}^{w}) Θ (t_{M} - t_{I}^{w}) 〉 & = \int_{0}^{t} b_{r} (τ) d τ = \int_{0}^{t} K_{I}^{w} (τ) Φ_{M} (τ) d τ \\ D (t) = 〈 Θ (t - t_{M}) Θ (t_{I}^{w} - t_{M}) 〉 & = \int_{0}^{t} b_{d} (τ) d τ = \int_{0}^{t} K_{M} (τ) Φ_{I}^{w} (τ) d τ \\ R (t) + D (t) = \int_{0}^{t} K_{I, M}^{w} (τ) d τ & = \int_{0}^{t} [b_{d} (τ) + b_{r} (τ)] d τ = 1 - Φ_{I}^{w} (t) Φ_{M} (t) \\ D (\infty) + R (\infty) = 1 \\ 〈 A (t - t_{I}) Θ (t_{M} - t_{I}) 〉 & = 〈 A (t - t_{I}) Φ_{M} (t_{I}) 〉 = \int_{0}^{t} A (t - τ) Φ_{M} (τ) K_{I} (τ) d τ \\ 〈 A (t - t_{M} Θ (t_{I}^{w} - t_{M}) 〉 & = 〈 A (t - t_{M}) Φ_{I}^{w} (t_{M}) 〉 = \int_{0}^{t} A (t - τ) Φ_{I}^{w} (τ) K_{M} (τ) d τ \end{matrix}

(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 (τ)

and the kernels

K (τ)

(i.e.,

A (τ), K (τ) = 0

for

τ < 0

). Of utmost importance are the "defective" PDFs (DPDFs)

b_{d, r} (t)

of death and recovery. "Defective" means that

b_{d, r} (t)

are no proper PDFs since they are not normalized to one, but rather to

D (\infty), R (\infty) < 1

, respectively. Consult [13] for a recent account of defective distributions and related stochastic processes. They have the following interpretation.

b_{d} (t) d t = K_{M} (t) Φ_{I}^{w} (t) d t

is the probability of transition I → D within

[t, t + d t]

of an infected walker (infected at

t^{'} = 0

b_{r} (t) d t = K_{I}^{w} (t) Φ_{M} (t) d t

is the probability of transition I → S within

[t, t + d t]

of a walker infected at

t^{'} = 0

. Therefore,

K_{I, M}^{w} (t) = b_{r} (t) + b_{d} (t) = - \frac{d}{d t} 〈Θ (t_{M} - t) Θ (t_{I}^{w} - t)〉 = - \frac{d}{d t} [Φ_{I}^{w} (t) Φ_{M} (t)]

(8)

is non-negative (as are

K_{I}^{w} = - \frac{d}{d t} Φ_{I}^{w} \geq 0

K_{M} = - \frac{d}{d t} Φ_{M} \geq 0

) and is a proper well-normalized PDF of all exits of walkers from compartment I (i.e., I → S + I → D). Without mortality (

Φ_{M} (t) = 1

) this PDF retrieves

K_{I, M}^{w} (t) = K_{I}^{w} (t)

The quantities

R (t), D (t)

introduced in (7) have the following interpretation.

R (t)

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

D (t)

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:

R (\infty)

has the interpretation of the overall probability that an infected walker survives the infection, and

D (\infty)

is the overall probability for an infected walker to die from the disease. We refer

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

D (\infty)

as we will see in details subsequently. A small value

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

D (\infty) < 1

(as

b_{d}

is defective).

D (\infty) \to 1

represents the limit of a fatal disease, and

D (\infty) \to 0

a disease without mortality.

R (\infty) < 1

(as

b_{r}

is defective) is the complementary probability with

D (\infty) + R (\infty) = 1

With these remarks, system (2) reads

\begin{matrix} \frac{d}{d t} S_{w} (t) & = - A_{w} (t) + \int_{0}^{t} A_{w} (t - τ) K_{I}^{w} (τ) Φ_{M} (τ) d τ + J_{w} (0) K_{I}^{w} (t) Φ_{M} (t) \\ \frac{d}{d t} J_{w} (t) & = \frac{d}{d t} \int_{0}^{t} A_{w} (τ) Φ_{M} (t - τ) Φ_{I}^{w} (t - τ) d τ - J_{w} (0) [K_{I}^{w} (t) Φ_{M} (t) + K_{M} (t) Φ_{I}^{w} (t)] \\ \frac{d}{d t} S_{n} (t) & = - A_{n} (t) + \int_{0}^{t} A_{n} (t - τ) K_{I}^{n} (τ) d τ + J_{n} (0) K_{I}^{n} (t) \\ \frac{d}{d t} J_{n} (t) & = - \frac{d}{d t} S_{n} (t) . \end{matrix}

(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}{d t} J_{w} (t) = A_{w} (t) - \int_{0}^{t} A_{w} (t - τ) K_{I, M}^{w} (τ) d τ - J_{w} (0) K_{I, M}^{w} (t) .

(10)

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

\begin{matrix} \frac{d}{d t} d_{w} (t) = - \frac{d}{d t} (S_{w} (t) + J_{w} (t)) & = 〈A_{w} (t - t_{M}) Θ (t_{I}^{w} - t_{M})〉 + J_{w} (0) 〈 δ (t - t_{M}) Θ (t_{I}^{w} - t_{M}) 〉 \\ = \int_{0}^{t} A_{w} (t - τ) K_{M} (τ) Φ_{I}^{w} (τ) d τ + J_{w} (0) K_{M} (t) Φ_{I}^{w} (t) \end{matrix}

(11)

where clearly

\frac{d}{d t} d_{w} (t) \geq 0

. Integrating this relation yields the fraction

d_{w} (t)

of dead walkers

\begin{matrix} d_{w} (t) & = 1 - S_{w} (t) - J_{w} (t) \\ = \int_{0}^{t} A_{w} (t - τ) 〈 Θ (τ - t_{M}) Θ (t_{I}^{w} - t_{M}) 〉 d τ + J_{w} (0) 〈 Θ (t - t_{M}) Θ (t_{I}^{w} - t_{M}) 〉 \\ = \int_{0}^{t} A_{w} (t - τ) D (τ) d τ + J_{w} (0) D (t) . \end{matrix}

(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{matrix} r_{w} (t) & = \int_{0}^{t} A_{w} (t - τ) 〈Θ (τ - t_{I}^{w}) Θ (t_{M} - t_{I}^{w})〉 d τ + J_{w} (0) 〈Θ (t - t_{I}^{w}) Θ (t_{M} - t_{I}^{w})〉 \\ = \int_{0}^{t} A_{w} (t - τ) R (τ) d τ + J_{w} (0) R (t) . \end{matrix}

(13)

The quantity

r_{w} (t)

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} (t) + d_{w} (t) = \int_{0}^{t} A_{w} (t - τ) [1 - Φ_{I}^{w} (τ) Φ_{M} (τ)] d τ + J_{w} (0) (1 - Φ_{I}^{w} (t) Φ_{M} (t)] .

(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} (t) \in [0, 1]

. Contrarily, the quantity

r_{w} (t)

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} Φ_{M}

). To shed more light on the behavior of

r_{w} (t)

consider for a moment zero mortality (

R (\infty) = 1

) and

t \to \infty

: We then have

A_{w} (\infty) = β_{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} (t)

turns into the proper non-defective PDF

b_{r} = K_{I}^{w}

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

\begin{matrix} S_{w} (t) & = 1 - J_{w} (0) [Φ_{M} (t) Φ_{I}^{w} (t) + D (t)] - \int_{0}^{t} A_{w} (τ) [Φ_{M} (t - τ) Φ_{I}^{w} (t - τ) + D (t - τ)] d τ \\ J_{w} (t) & = J_{w} (0) Φ_{M} (t) Φ_{I}^{w} (t) + \int_{0}^{t} A_{w} (τ) Φ_{M} (t - τ) Φ_{I}^{w} (t - τ) d τ \\ S_{n} (t) & = 1 - J_{n} (t) \\ J_{n} (t) & = J_{n} (0) Φ_{I}^{n} (t) + \int_{0}^{t} A_{n} (τ) Φ_{I}^{n} (t - τ) d τ \end{matrix}

(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

A_{w} (t) = A_{w} [S_{w} (t), J_{n} (t)]

A_{n} (t) = A_{w} [S_{n} (t), J_{w} (t)]

(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 (t)

\hat{g} (λ) = \int_{0}^{\infty} g (t) e^{- λ t} d t

where

λ

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

g (\infty) = lim_{λ \to 0} λ \hat{g} (λ) (= lim_{λ \to 0} \int_{0}^{\infty} g (\frac{τ}{λ}) e^{- τ} d τ \to g (\infty) \int_{0}^{\infty} e^{- τ} d τ) .

(16)

Observing that the LT of

Φ_{I}^{w} (t) Φ_{M} (t)

λ^{- 1} [1 - {\hat{K}}_{I, M} (λ)]

and

\hat{D} (λ) = λ^{- 1} {\hat{b}}_{d} (λ)

where

{\hat{b}}_{d} (0) = D (\infty)

(see (7)), we arrive at

\begin{matrix} J_{w} (\infty) & = lim_{λ \to 0} λ {\hat{J}}_{w} (λ) = [1 - {\hat{K}}_{I, M} (0)] [J_{w} (0) + {\hat{A}}_{w} (0)] = 0 \\ J_{n} (\infty) & = lim_{λ \to 0} λ {\hat{J}}_{n} (λ) = [1 - {\hat{K}}_{n} (0)] [J_{n} (0) + {\hat{A}}_{n} (0)] = 0 \\ d_{w} (\infty) & = lim_{λ \to 0} λ {\hat{d}}_{w} (λ) = D (\infty) [J_{w} (0) + {\hat{A}}_{w} (0)], & {\hat{A}}_{w} (0) = β_{w} \int_{0}^{\infty} S_{w} (τ) J_{n} (τ) d τ \\ S_{w} (\infty) & = 1 - d_{w} (\infty) \\ S_{n} (\infty) & = 1 \end{matrix}

(17)

where

{\hat{A}}_{w, n} (0) = \int_{0}^{\infty} A_{w, n} (t) d t < \infty

. In the same way one obtains

r_{w} (\infty) = R (\infty) (J_{w} (0) + {\hat{A}}_{w} (0)) .

(18)

Since

D (\infty)

is non-zero, the asymptotic values

S_{w} (\infty), d_{w} (\infty)

depend on the initial condition

J_{w} (0)

and the infection (rate) history. This is not any more true for zero mortality (

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)} = D (\infty), P_{R} = 1 - P_{D} = \frac{r_{w} (\infty)}{r_{w} (\infty) + d_{w} (\infty)} = R (\infty)

(19)

consistent with our previous interpretation of

D (\infty), R (\infty)

, and the ratio

\frac{d_{w} (\infty)}{r_{w} (\infty)} = \frac{D (\infty)}{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

Φ_{I}^{w} (t) = e^{- ξ_{I}^{w} t}

Φ_{M} (t) = e^{- ξ_{M} t}

Φ_{I}^{n} (t) = e^{- ξ_{I}^{n} t}

(

〈 t_{I}^{w, n} 〉 = ξ_{I}^{w, n})^{- 1}

〈 t_{M} 〉 = {(ξ_{M})}^{- 1}

). Then (9) takes with (15) the memoryless form

\begin{matrix} \frac{d}{d t} S_{w} (t) & = - β_{w} S_{w} (t) J_{n} (t) + ξ_{I}^{w} J_{w} (t) \\ \frac{d}{d t} J_{w} (t) & = β_{w} S_{w} (t) J_{n} (t) - (ξ_{I}^{w} + ξ_{M}) J_{w} (t) \\ \frac{d}{d t} d_{w} (t) & = ξ_{M} J_{w} (t) \\ \frac{d}{d t} S_{n} (t) & = - β_{n} S_{n} (t) J_{w} (t) + ξ_{I}^{n} J_{n} (t) \\ \frac{d}{d t} J_{n} (t) & = β_{n} S_{n} (t) J_{w} (t) - ξ_{I}^{n} J_{n} (t) . \end{matrix}

(21)

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

\begin{matrix} J_{w} (\infty) = J_{n} (\infty) = A_{w} (\infty) = A_{n} (\infty) = 0 \\ S_{w} (\infty) = 1 - d_{w} (\infty), & d_{w} (\infty) = ξ_{M} \int_{0}^{\infty} J_{w} (τ) d τ \\ S_{n} (\infty) = 1 . \end{matrix}

(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} (0) + \int_{0}^{\infty} A_{w} (t) d t - (ξ_{I}^{w} + ξ_{M}) \int_{0}^{\infty} J_{w} (t) d t

(23)

thus

\int_{0}^{\infty} J_{w} (t) d t = \frac{1}{ξ_{M} + ξ_{I}^{w}} (J_{w} (0) + {\hat{A}}_{w} (0))

. Plugging this relation into (22) and accounting for

D (\infty) = \frac{ξ_{M}}{ξ_{I}^{w} + ξ_{M}}

recovers indeed the representation of expression (17).

For zero mortality

ξ_{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_{α, ξ} (t) = \frac{ξ^{α} t^{α - 1}}{Γ (α)} e^{- ξ t}, ξ, α > 0

(24)

where

α

is the so called `shape parameter’ and

ξ

the rate parameter (often is used the term ’scale parameter’

θ = ξ^{- 1}

) and

Γ (α)

stands for the Gamma function. We also will subsequently use the LT of the Gamma PDF

{\hat{K}}_{α, ξ} (λ) = \int_{0}^{\infty} K_{α, ξ} (t) e^{- λ t} d t = \frac{ξ^{α}}{{(λ + ξ)}^{α}}

(25)

The Gamma PDF has finite mean

{〈 t 〉}_{α, ξ} = \int_{0}^{\infty} t K_{α, ξ} (t) d t = \frac{α}{ξ}

and for

α < 1

the Gamma-PDF is weakly singular at

t = 0

and

α = 1

recovers exponential PDFs. For

α \leq 1

the Gamma PDF is completely monotonic (CM) (see Appendix, (A9) for a definition). For the range

α > 1

the Gamma-PDF has a maximum at

t_{m a x} = \frac{α - 1}{ξ}

and becomes narrower the larger

ξ

(while keeping its mean

α / ξ

fixed), especially we can generate sharp waiting times using the feature

lim_{ξ \to \infty} K_{α = ξ T_{0}, ξ} (t) = δ (t - T_{0}) .

(26)

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

Φ_{α, ξ} (t) = \int_{t}^{\infty} \frac{ξ^{α} t^{α - 1}}{Γ (α)} e^{- ξ t} d t = \frac{Γ (α, ξ t)}{Γ (α)}

(27)

where

Γ (α, x)

indicates the upper incomplete Gamma function with

Γ (α, 0) = Γ (α)

. (27) fulfills necessarily the initial condition

Φ_{α, ξ} (0) = 1

and is vanishing at infinity

Φ_{α, ξ} (\infty) = 0

3. Endemic Equilibrium for Zero Mortality

Here we consider the large time asymptotics of the compartment populations without mortality (

S_{w} (t) + J_{w} (t) = 1

) where the self-consistent system (15) reads

\begin{matrix} J_{w} (t) & = J_{w} (0) Φ_{I}^{w} (t) + \int_{0}^{t} A_{w} (t - τ) Φ_{I}^{w} (τ) d τ \\ J_{n} (t) & = J_{n} (0) Φ_{I}^{n} (t) + \int_{0}^{t} A_{n} (t - τ) Φ_{I}^{n} (τ) d τ \end{matrix}

(28)

The endemic state emerging in the large time asymptotics does not depend on the initial conditions

J_{w, n} (0)

Φ_{I}^{w, n} (t) \to 0

. For what follows it is convenient to consider the LTs of (28) which read

\begin{matrix} {\hat{J}}_{w} (λ) & = [J_{w} (0) + {\hat{A}}_{w} (λ)] \frac{1 - {\hat{K}}_{I}^{w} (λ)}{λ} \\ {\hat{J}}_{n} (λ) & = [J_{n} (0) + {\hat{A}}_{n} (λ)] \frac{1 - {\hat{K}}_{I}^{n} (λ)}{λ} \end{matrix}

(29)

where

Φ_{I}^{w, n} (λ) = \frac{1 - {\hat{K}}_{I}^{w, n} (λ)}{λ}

are the LTs of the persistence distributions, and

{\hat{S}}_{w, n} (λ) + {\hat{J}}_{w, n} (λ) = \frac{1}{λ}

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

〈 t_{I}^{w, n} 〉 = lim_{λ \to 0} \frac{1 - {\hat{K}}_{I}^{w, n} (λ)}{λ} = - \frac{d}{d λ} {\hat{K}}_{I}^{w, n} (λ) |_{λ = 0} = \int_{0}^{\infty} Φ_{I}^{w, n} (t) d t = \int_{0}^{\infty} τ K_{I}^{w, n} (τ) d τ < \infty

(30)

thus the admissible range of the Laplace variable is

λ \geq 0

(if chosen real). Now using (16) we obtain the endemic equilibrium from

J_{w, n} (\infty) = {lim}_{λ \to 0} λ {\hat{J}}_{w, n} (λ)

where the initial conditions are wiped out at infinity as

{\hat{K}}_{w, n} (λ) |_{λ = 0} = 1

. Assuming that the infection rates are constant in the endemic equilibrium we have

A_{w, n} (λ) \sim A_{w, n} (\infty) / λ

, (

λ \to 0

) and arrive at

\begin{matrix} J_{w} (\infty) & = A_{w} (\infty) 〈t_{I}^{w}〉, (A_{w} (\infty) = β_{w} S_{w} (\infty) J_{n} (\infty)) \\ J_{n} (\infty) & = A_{n} (\infty) 〈t_{I}^{n}〉, (A_{n} (\infty) = β_{n} S_{n} (\infty) J_{w} (\infty)) \end{matrix}

(31)

thus

\begin{matrix} \frac{J_{w} (\infty)}{1 - J_{w} (\infty)} - β_{w} 〈 t_{I}^{w} 〉 J_{n} (\infty) = 0 \\ \frac{J_{n} (\infty)}{1 - J_{n} (\infty)} - β_{n} 〈 t_{I}^{n} 〉 J_{w} (\infty) = 0 . \end{matrix}

(32)

One can see that the globally healthy state

J_{w, n} (0) = 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{matrix} J_{w} (\infty) = J_{w}^{e} = \frac{β_{w} β_{n} 〈 t_{I}^{w} 〉 〈 t_{I}^{n} 〉 - 1}{β_{n} 〈 t_{I}^{n} 〉 [1 + β_{w} 〈 t_{I}^{w} 〉]} = \frac{R_{0} - 1}{R_{0}} \frac{β_{w} 〈 t_{I}^{w} 〉}{1 + β_{w} 〈 t_{I}^{w} 〉} \\ J_{n} (\infty) = J_{n}^{e} = \frac{β_{w} β_{n} 〈 t_{I}^{w} 〉 〈 t_{I}^{n} 〉 - 1}{β_{w} 〈 t_{I}^{w} 〉 [1 + β_{n} 〈 t_{I}^{n} 〉]} = \frac{R_{0} - 1}{R_{0}} \frac{β_{n} 〈 t_{I}^{n} 〉}{1 + β_{n} 〈 t_{I}^{n} 〉} \end{matrix}

(33)

for the endemic equilibrium which is independent of the initial conditions

J_{w, n} (0)

. It depends only on the phenomenological rate constants

β_{w, n}

and the mean infection time spans

〈 t_{I}^{w, n} 〉

. 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} = β_{w} β_{n} 〈 t_{I}^{w} 〉 〈 t_{I}^{n} 〉 > 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} (0) = 1 / N

and no infected walkers

S_{w} (0) = 1

. The expected number of walkers infected by this first infected node during its infectious period

t_{I}^{n}

〈 Z_{I} (t_{I}^{n}) 〉 \sim Z A_{w} (0) 〈 t_{I}^{n} 〉 = Z 〈 t_{I}^{n} 〉 β_{w} / N \sim Z 〈 J_{w} (t_{I}^{w}) 〉 .

The number of nodes getting infected by these

〈 Z_{I} (t_{I}^{n}) 〉

infected walkers during their infectious time

t_{I}^{w}

N_{I} (〈 t_{I}^{n} 〉 + 〈 t_{I}^{w} 〉) \sim N 〈 A_{n} (t_{I}^{n}) 〉 〈 t_{I}^{w} 〉 \sim N β_{n} 〈 J_{w} (t_{I}^{w}) 〉 = β_{n} β_{w} 〈 t_{I}^{n} 〉 〈 t_{I}^{w} 〉 = 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} (0)

. 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

β_{w} 〈 t_{I}^{w} 〉 \to \infty

(while

β_{n} 〈 t_{I}^{n} 〉

are kept constant) where all walkers become infected

J_{w}^{e} \to 1

but not all nodes

J_{n}^{e} \to \frac{β_{n} 〈 t_{I}^{n} 〉}{1 + β_{n} 〈 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{matrix} S_{w} (t) & = S_{w}^{e} + u_{w} e^{μ t}, J_{w} (t) = J_{w}^{e} - u_{w} e^{μ t} \\ S_{n} (t) & = S_{n}^{e} + u_{n} e^{μ t}, J_{n} (t) = J_{n}^{e} - u_{n} e^{μ t} \end{matrix}

(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{matrix} A_{w} (t) & = β_{w} S_{w} (t) J_{n} (t) = β_{w} S_{w}^{e} J_{n}^{e} + β_{w} (u_{w} J_{n}^{e} - u_{n} S_{w}^{e}) e^{μ t} \\ A_{n} (t) & = β_{n} S_{n} (t) J_{w} (t) = β_{n} S_{n}^{e} J_{w}^{e} + β_{n} (u_{n} J_{w}^{e} - u_{w} S_{n}^{e}) e^{μ t} \end{matrix}

(35)

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

[\begin{matrix} μ + β_{w} J_{n}^{e} [1 - {\hat{K}}_{I}^{w} (μ)]; & - β_{w} S_{w}^{e} [1 - {\hat{K}}_{I}^{w} (μ)] \\ - β_{n} S_{n}^{e} [1 - {\hat{K}}_{I}^{n} (μ)]; & μ + β_{n} J_{w}^{e} [1 - {\hat{K}}_{I}^{n} (μ)] \end{matrix}] \cdot (\begin{matrix} u_{w} \\ u_{n} \end{matrix}) = (\begin{matrix} 0 \\ 0 \end{matrix})

(36)

where we have used

〈 e^{- μ t_{I}^{w, n}} 〉 = {\hat{K}}_{I}^{w, n} (μ)

and the cases of

δ

-distributed

t_{I}^{w, n}

are contained for

{\hat{K}}_{I}^{w, n} (μ) = e^{- μ t_{I}^{w, n}}

. We point out that in ansatz (35) we relax causality i.e., we admit

A_{w, n} (t - τ) \neq 0

for

t - τ < 0

thus

〈 e^{μ (t - t_{I}^{w, n})} 〉 = e^{μ t} 〈 e^{- μ t_{I}^{w, n}} 〉 = e^{μ t} {\hat{K}}_{I}^{w, n} (μ) .

(37)

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

μ

μ^{2} + μ (β_{w} J_{n}^{e} [1 - {\hat{K}}_{I}^{w} (μ)] + β_{n} J_{w}^{e} [1 - {\hat{K}}_{I}^{n} (μ)]) + β_{w} β_{n} [1 - {\hat{K}}_{I}^{w} (μ)] [1 - {\hat{K}}_{I}^{n} (μ)] (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

μ

exist. Consider this first for the globally healthy state

J_{n} = 0, J_{w} = 0

for which Equation (38) reads

G (μ) = 1 - β_{w} β_{n} \frac{[1 - {\hat{K}}_{I}^{w} (μ)]}{μ} \frac{[1 - {\hat{K}}_{I}^{n} (μ)]}{μ} = 1 - β_{w} β_{n} {\hat{Φ}}_{I}^{w} (μ) {\hat{Φ}}_{I}^{n} (μ) = 0

(39)

where we notice that

\frac{[1 - {\hat{K}}_{I}^{w, n} (μ)]}{μ} = {\hat{Φ}}_{I}^{w, n} (μ)

are the LTs of the persistence probabilities of the infection time spans. Consider this equation for

μ \to 0

and take into account (30) we arrive at

G (0) = 1 - β_{w} β_{n} 〈 t_{I}^{w} 〉 〈 t_{I}^{n} 〉 .

(40)

We observe that

G (0) < 0

for

R_{0} = β_{n} β_{w} 〈 t_{I}^{w} 〉 〈 t_{I}^{n} 〉 > 1

. On the other hand, we have for

μ \to \infty

that

{\hat{Φ}}_{I}^{w, n} (μ) \to 0

and hence

G (\infty) = 1 .

(41)

One can hence infer from complete monotony of

{\hat{Φ}}_{I}^{w, n} (μ)

and therefore of

{\hat{Φ}}_{I}^{w} (μ) {\hat{Φ}}_{I}^{n} (μ)

(see Appendix A.2, Equation (A9) for a precise definition), that

\frac{d}{d μ} G (μ) > 0

(

μ \geq 0

) thus

G (μ) = 0

has one single positive zero only if

G (0) < 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

〈 t_{I}^{w, n} 〉 = \infty

. This is true for fat-tailed kernels scaling as

K_{I}^{w, n} (t) \propto t^{- α - 1}

(

α \in (0, 1)

) for

t \to \infty

. We consider such a distribution briefly in subsequent section. We plot function

G (μ)

versus

μ

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} (μ) = 0

where from (38) this function reads

G_{e} (μ) = 1 - β_{w} β_{n} {\hat{Φ}}_{I}^{w} (μ) {\hat{Φ}}_{I}^{n} (μ) + β_{w} J_{n}^{e} {\hat{Φ}}_{I}^{w} (μ) + β_{n} J_{w}^{e} {\hat{Φ}}_{I}^{n} (μ) + β_{w} β_{n} (J_{w}^{e} + J_{n}^{e}) {\hat{Φ}}_{I}^{w} (μ) {\hat{Φ}}_{I}^{n} (μ)

(42)

with

\begin{matrix} G_{e} (0) & = 1 - R_{0} + β_{w} 〈 t_{I}^{w} 〉 J_{n}^{e} + β_{n} 〈 t_{I}^{n} 〉 J_{w}^{e} + (J_{w}^{e} + J_{n}^{e}) R_{0} \\ = R_{0} - 1 . \end{matrix}

(43)

On the other hand,

G_{e} (\infty) = 1

(as

{\hat{Φ}}_{I}^{w, n} (\infty) = 0

) and from monotony of

G_{e} (μ)

follows that there is no positive solution of

G_{e} (μ) = 0

for

R_{0} > 1

. We plot

G_{e} (μ)

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} (μ) > 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} (t) = 1 + a e^{μ t}, J_{w} (t) = b e^{μ t}, d_{w} (t) = - (a + b) e^{μ t}, S_{n} (t) = 1 - c e^{μ t}, J_{n} (t) = c e^{μ t}

(44)

with

A_{w} (t) = β_{w} c e^{μ t}

and

A_{n} (t) = β_{n} b e^{μ t}

. Plugging this ansatz for

μ \geq 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

[\begin{matrix} μ; & 0; & β_{w} [1 - {\hat{b}}_{r} (μ)]] \\ μ; & μ; & β_{w} {\hat{b}}_{d} (μ) \\ 0 & - β_{n} [1 - {\bar{K}}_{I}^{n} (μ)]; & μ \end{matrix}] \cdot (\begin{matrix} a \\ b \\ c \end{matrix}) = (\begin{matrix} 0 \\ 0 \\ 0 \end{matrix}) .

(45)

Putting the determinant of the matrix to zero yields the condition

μ^{2} - β_{n} β_{w} [1 - {\bar{K}}_{I}^{n} (μ)] [1 - {\hat{b}}_{r} (μ) - {\hat{b}}_{d} (μ)] = 0

(46)

where the LTs

{\hat{b}}_{r} (μ), {\hat{b}}_{d} (μ)

of the DPDFs

b_{r, d} (t)

defined in (7) come into play. We are interested under which conditions there is a positive solution (instability of healthy state) of (46). Since

{\hat{b}}_{r} (0) = R (\infty)

and

{\hat{b}}_{d} (0) = D (\infty)

with

R (\infty) + D (\infty) = 1

we see that

μ = 0

is solution of (46). Recall from (7) that

b_{d} (t) + b_{r} (t) = K_{I, M}^{w} (t)

is the (properly normalized) PDF (8). Condition (46) then reads

G_{M} (μ) = 1 - β_{n} β_{w} \frac{[1 - {\bar{K}}_{I}^{n} (μ)]}{μ} \frac{[1 - {\hat{K}}_{I, M}^{w} (μ)]}{μ} = 0

(47)

where

\frac{[1 - {\hat{K}}_{I, M}^{w} (μ)]}{μ}

is the LT of the persistence probability

Φ_{M} (t) Φ_{I}^{w} (t)

of the walker’s infection, i.e., the probability that

t < m i n (t_{I}^{w}, t_{M})

(see Remark I). For zero mortality we have

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

〈 m i n (t_{I}^{w}, t_{M}) 〉 = 〈 t_{I M}^{w} 〉 = \frac{[1 - {\hat{K}}_{I, M}^{w} (μ)]}{μ} |_{μ = 0} = \int_{0}^{\infty} t K_{I, M}^{w} (t) d t = \int_{0}^{\infty} Φ_{M} (t) Φ_{I}^{w} (t) d t \leq \int_{0}^{\infty} Φ_{I}^{w} (t) d t = 〈 t_{I}^{w} 〉

(48)

where we arrive at

G_{M} (0) = 1 - β_{n} β_{w} 〈 t_{I}^{n} 〉 〈 t_{I M}^{w} 〉 .

(49)

Relation (48) shows that

〈 t_{I M}^{w} 〉 \leq 〈 t_{I}^{w} 〉

(equality only for zero mortality). On the other hand we have

G_{M} (\infty) = 1

, so there is a positive solution of

G_{M} (μ) = 0

only if

R_{M} = β_{n} β_{w} 〈 t_{I}^{n} 〉 〈 t_{I M}^{w} 〉 > 1

(50)

where

R_{M}

is the basic reproduction number modified by mortality with

R_{M} \leq R_{0}

(equality only for zero mortality). To visualize the effect of mortality on the instability of the healthy state we plot

G_{M} (μ)

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

ξ_{M}

and plot the result in Figure 6 where one can see that

R_{M}

is monotonous decreasing with mortality rate

ξ_{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

ξ_{M}

approaching zero for

ξ_{M} \to \infty

. Diseases with high mortality stabilize the healthy state even for

〈 t_{I}^{w} 〉 \to \infty

Consider briefly the case where

t_{M}

is exponentially distributed (i.e.,

α_{M} = 1

in the Gamma distribution of

t_{M}

) with

Φ_{M} (t) = e^{- ξ_{M} t}

. Then we have

〈 t_{I M}^{w} 〉 = {\hat{Φ}}_{I}^{w} (ξ_{M})

thus

R_{M} = β_{w} β_{n} 〈 t_{I}^{n} 〉 {\hat{Φ}}_{I}^{w} (ξ_{M}) .

(51)

The zero mortality case is recovered for

ξ_{M} = 0

with

{\hat{Φ}}_{I}^{w} (0) = 〈 t_{I}^{w} 〉

. For Gamma distributed

t_{I}^{w, n}

this yields

R_{M} = β_{w} β_{n} \frac{α_{I}^{n}}{ξ_{I}^{n} ξ_{M}} (1 - \frac{{(ξ_{I}^{w})}^{α_{I}^{w}}}{{(ξ_{M} + ξ_{I}^{w})}^{α_{I}^{w}}}) .

(52)

where

α_{I}^{w, n}

ξ_{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

α_{I}^{w, n} = 1

and yields

R_{M} = \frac{β_{w} β_{n}}{ξ_{I}^{n} (ξ_{I}^{w} + ξ_{M})}

(53)

containing the zero mortality case for

ξ_{M} = 0

Fat-tailed distributed

t_{I}^{w}

Finally, an interesting case emerges if

t_{I}^{w}

follows a fat-tailed distribution, i.e.,

Φ_{I}^{w} (t) \propto t^{- α}

for t large (

α \in (0, 1)

) and

〈 t_{I}^{w} 〉 = \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

Φ_{I}^{w} (t) = E_{α} (- ξ_{I}^{w} t^{α})

where

E_{α} (τ)

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

α = 1

(

E_{1} (- τ) = e^{- τ}

). Assuming exponential mortality

Φ_{M} (t) = e^{- ξ_{M} t}

one obtains with (51)

R_{M} = β_{w} β_{n} 〈 t_{I}^{n} 〉 \frac{{(ξ_{M})}^{α - 1}}{ξ_{I}^{w} + {(ξ_{M})}^{α}}, α \in (0, 1)

(54)

containing the LT of the ML persistence probability distribution

{\hat{Φ}}_{I}^{w} (λ) = λ^{α - 1} / (ξ_{I}^{w} + λ^{α})

. The essential feature here is that

R_{M}

is weakly singular at

ξ_{M} = 0

with a monotonously decreasing

ξ_{M}^{α - 1}

scaling law, where the healthy state becomes stable for mortality parameters larger as

ξ_{M} \approx 1

. We depict this case in Figure 6 for

α = ξ = 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} (t) + Z_{I} (t) + Z_{D} (t)

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

β_{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

(A_{i j})

where

A_{i j} = 1

if the pair of nodes

i, j

is connected by an edge, and

A_{i j} = 0

if the pair is disconnected. Further, we assume

A_{i i} = 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 = Δ t, 2 Δ 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

Π (i \to j) = \frac{A_{i j}}{k_{i}}, z = 1, \dots, Z, i, j = 1, \dots, N

(55)

which is normalized

\sum_{j = 1}^{N} Π (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

lim_{n \to \infty} \frac{1}{n} \sum_{k = 1}^{n} A (t - T_{k}) = \int_{0}^{t} A (t - τ) K (τ) d τ

(56)

for random variables T drawn from PDFs

K (τ)

. We perform this average for any time increment

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

β_{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

Δ t

of walker’s steps. The so determined rate parameters

β_{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} (0) = 10 ≪ N

(

N_{I} (0) \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

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

ξ_{M}

and a much smaller overall mortality

D (\infty) \approx 1 %

. The different network connectivity leads to different values of

β_{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

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

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

〈 t_{I}^{w} 〉 = 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

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} \leq 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

β_{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.

Acknowledgments

T.G. gratefully acknowledges to have been hosted at the Institut Jean le Rond d’Alembert in the framework of an internship (stage M1 Physique

n^{o} 28281

) for the development of the PYTHON simulation codes and participation in the present study.

Appendix A

Appendix A.1. Some Basic Notions

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

A (t) = A (t) Θ (t)

(A1)

where

Θ (t)

indicates the Heaviside step function defined by

Θ (ζ) = \{\begin{matrix} 1, ζ \geq 0 \\ 0, ζ < 0 \end{matrix}

(A2)

with

Θ (ζ) + Θ (- ζ) = 1

and in or definition

Θ (0) = 1

. Its derivative yields the Dirac

δ

-distribution

\frac{d}{d ζ} Θ (ζ) = δ (ζ)

. We use throughout the paper mutually independent strictly positive random variables

T_{1}, \dots, T_{n} \in R_{+}

The random variables

T_{j}

are assumed to follow their specific PDFs

P r o b [T_{j} \in [u, u + d u]] = 〈δ (u - T_{j})〉 = K_{j} (u) 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} d t_{1} \dots d t_{n} f (t; t_{1}, t_{2}, \dots, t_{n}) K_{1} (t_{1}) \dots K_{n} (t_{n})

(A4)

for suitable functions f. For

f (t; T_{1}, T_{2}, \dots, T_{n}) = g_{1} (T_{1}) \dots g_{n} (T_{n})

using independence of the

T_{j}

this yields

〈g_{1} (T_{1}) \dots g_{n} (T_{n})〉 = 〈g_{1} (T_{1}),〉 \dots 〈g_{n} (T_{n})〉 .

Important cases emerge by applying (A4) to exponentials

〈e^{- λ (T_{1} + \dots + T_{n})}〉 = \int_{0}^{\infty} e^{- λ t} 〈δ (t - T_{1} - \dots - T_{n})〉 d t = {\hat{K}}_{1} (λ) \dots {\hat{K}}_{n} (λ), ℜ {λ} \geq 0

(A5)

In this relation the LTs of the PDFs come into play

{\hat{K}}_{j} (λ) = \int_{0}^{\infty} e^{- λ t} K_{j} (t) d t

(A6)

where

{\hat{K}}_{j} (λ) |_{λ = 0} = 1

reflects the normalization of PDFs (A3). A further observation is

〈δ (t - T_{1} - \dots - T_{n})〉 = (K_{1} ★ \dots ★ K_{n}) (t)

(A7)

where ★ stands for convolution

(K_{1} ★ K_{2}) (t) = \int_{0}^{t} K_{1} (τ) K_{2} (t - τ) d τ

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

μ > 0

with

R_{0} - 1 > 0

. First, we observe that

{\hat{Φ}}_{I}^{w, n} (μ) \leq 〈 t_{I}^{w, n} 〉

with

{\hat{Φ}}_{I}^{w, n} (0) = 〈 t_{I}^{w, n} 〉

and

{\hat{Φ}}_{I}^{w, n} (\infty) = 0

. For our convenience, we introduce the functions

λ_{w, n} (μ) = \frac{{\hat{Φ}}_{I}^{w, n} (μ)}{〈 t_{I}^{w, n} 〉} \in (0, 1]

(A8)

which the LTs of the normalized

Φ_{I}^{w, n} (t) / 〈 t_{I}^{w, n} 〉

and which are by virtue of Bernstein’s theorem [44] completely monotonic (CM) with respect to

μ

, i.e.,

{(- 1)}^{n} \frac{d^{n}}{d μ^{n}} λ_{w, n} (μ) \geq 0, μ \in (0, \infty)

(A9)

inheriting this feature from the exponential

e^{- μ τ}

(

t > 0

). Therefore,

\frac{d}{d μ} λ_{w, n} (μ) = - \frac{1}{〈 t_{I}^{w, n} 〉} \int_{0}^{\infty} e^{- μ t} t Φ_{I}^{w, n} (t) d t < 0

(as

Φ_{I}^{n, w} (t) \in (0, 1]

) exists thus

λ_{w, n} (μ)

is monotonously decreasing with

μ

with

λ_{w, n} (0) = 1 \geq λ_{w, n} (μ) > 0

. Further we observe in Eqs. (33) that

J_{n}^{e} β_{w} 〈 t_{I}^{w} 〉 = \frac{R_{0} - 1}{1 + β_{n} 〈 t_{I}^{n} 〉}

J_{w}^{e} β_{n} 〈 t_{I}^{n} 〉 = \frac{R_{0} - 1}{1 + β_{w} 〈 t_{I}^{w} 〉}

thus

R_{0} (J_{n}^{e} + J_{w}^{e}) = (R_{0} - 1) (\frac{β_{n} 〈 t_{I}^{n} 〉}{1 + β_{n} 〈 t_{I}^{n} 〉} + \frac{β_{w} 〈 t_{I}^{w} 〉}{1 + β_{w} 〈 t_{I}^{w} 〉})

Then

G_{e} (μ)

reads

\begin{matrix} G_{e} (μ) & = 1 + λ_{w} (μ) λ_{n} (μ) \{- R_{0} + (R_{0} - 1) (\frac{β_{n} 〈 t_{I}^{n} 〉}{1 + β_{n} 〈 t_{I}^{n} 〉} + \frac{β_{w} 〈 t_{I}^{w} 〉}{1 + β_{w} 〈 t_{I}^{w} 〉})\} \\ + (R_{0} - 1) (\frac{λ_{w} (μ)}{1 + β_{n} 〈 t_{I}^{n} 〉} + \frac{λ_{n} (μ)}{1 + β_{w} 〈 t_{I}^{w} 〉}) \end{matrix}

(A10)

Now we observe that

λ_{w} (μ) λ_{n} (μ) \leq λ_{w, n} (μ)

thus a lower bound function

H (μ) \leq G_{e} (μ)

is generated by replacing

λ_{w, n} (μ) \to λ_{w} (μ) λ_{n} (μ)

in the second line. Now it is sufficient to prove that

0 < H (μ)

. We hence get for this lower bound function

\begin{matrix} H (μ) & = 1 + λ_{w} (μ) λ_{n} (μ) \{- R_{0} + (R_{0} - 1) (\frac{β_{n} 〈 t_{I}^{n} 〉}{1 + β_{n} 〈 t_{I}^{n} 〉} + \frac{β_{w} 〈 t_{I}^{w} 〉}{1 + β_{w} 〈 t_{I}^{w} 〉}) + (R_{0} - 1) (\frac{1}{1 + β_{n} 〈 t_{I}^{n} 〉} + \frac{1}{1 + β_{w} 〈 t_{I}^{w} 〉})\} \\ = 1 + λ_{w} (μ) λ_{n} (μ) (R_{0} - 2) \\ = 1 - λ_{w} (μ) λ_{n} (μ) + (R_{0} - 1) λ_{w} (μ) λ_{n} (μ) \end{matrix}

(A11)

and with

1 - λ_{w} (μ) λ_{n} (μ) \geq 0

and

(R_{0} - 1) λ_{w} (μ) λ_{n} (μ) > 0

it follows that

0 < H (μ) < G_{e} (μ)

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} (μ)

G_{e} (μ) = 1 + σ \hat{g} (μ) \geq 1 - \hat{g} (μ)

(A12)

where

σ = \pm 1

and

\hat{g} (μ)

is a non-negative CM function (see Figure 4) with maximum value

\hat{g} (0) = | R_{0} - 2 |

. Then the following two cases may occur.

Case (i) $σ = - 1$ ; $0 < G_{e} (0) = R_{0} - 1 < 1$ ( $1 < R_{0} < 2$ ):

Then

\hat{g} (μ)

can be represented as LT of a non-negative function

g (t)

and consider now

μ = μ_{1} + i μ_{2}

with

μ_{1} \geq 0

thus the real part of

\hat{g} (μ_{1} + i μ_{2})

can be written as

ℜ g (μ_{1} + i μ_{2}) = \int_{0}^{\infty} g (t) e^{- μ_{1} t} cos (μ_{2} t) d t, μ_{1} \geq 0

(A13)

with

\hat{g} (0) = 2 - R_{0}

and clearly

- \hat{g} (μ_{1}) \leq ℜ g (μ_{1} + i μ_{2}) \leq \hat{g} (μ_{1})

. Therefore,

ℜ G_{e} (μ_{1} + i μ_{2}) = 1 - ℜ g (μ_{1} + i μ_{2}) \geq 1 - \hat{g} (μ_{1}) \geq 1 - \hat{g} (0) = 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) $σ = + 1$ ; $R_{0} - 1 > 1$ :

Here we have two pertinent ranges of

R_{0}

. The first is the range (a)

\hat{g} (0) = R_{0} - 2 < 1

(i.e.,

R_{0} < 3

) and the second one (b) is

\hat{g} (0) = R_{0} - 2 > 1

. Clearly in the range (a) (A14) remains true and

ℜ G_{e} (μ_{1} + i μ_{2})

strictly positive. Hence for

1 < R_{0} < 3

no Hopf instability is possible.

This changes in the range (b) since

\hat{g} (0) = R_{0} - 2 > 1

thus

ℜ G_{e} (μ_{1} + i μ_{2}) = 1 + ℜ g (μ_{1} + i μ_{2})

may become negative. Therefore, for

R_{0} > 3

a Hopf instability of the endemic state becomes possible. However, the possibility that

ℜ G_{e} (μ_{1} + i μ_{2}) = 0

is only necessary but not sufficient for a Hopf instability. One also needs simultaneously that the imaginary part

ℑ G_{e} (μ_{1} + i μ_{2}) = 0

is vanishing for the same

μ = μ_{1} + i μ_{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} (k)

that a node has

0 \leq k \leq N - 1

connections is given by a binomial distribution

P_{N} (k) = (\begin{matrix} N - 1 \\ k \end{matrix}) p^{k} {(1 - p)}^{N - 1 - k} \to \frac{{〈 k 〉}^{k}}{k!} e^{- 〈 k 〉}

where

〈 k 〉 = (N - 1) p \sim N p

denotes the average degree. For

N \to \infty

(while

N p

is kept constant) the degree distribution

P_{N} (k)

converges to a Poisson distribution representing the infinite graph limit of the ER

G (N, p)

-model. Therefore,

P_{N} (k)

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{l o g N}{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

2 m

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

2 m

for all nodes. In the limit

p = 1

an ER graph is emerging with probability

2 m / (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 \leq 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 (k) \propto k^{- 2 - μ}, μ \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.

References

Rhodes, P.; Bryant, J.H. Public Health. Encycl. Br. 2024. Available online: https://www.britannica.com/topic/public-health.
Kermack, W.O.; McKendrick, A.G. A contribution to the mathematical theory of epidemics. Proc. Roy. Soc. A 1927, 115, 700–721. [Google Scholar]
Liu, W.M.; Levin, H.W.H.H.W.S.A. Dynamical behavior of epidemiological models with non-linear incidence rate. J. Math. Biol. 1987, 25, 359–380. [Google Scholar] [CrossRef]
Li, M.Y.; Graef, J.R.; Wang, L.; Karsai, J. Global dynamics of a SEIR model with varying total population size. Math. Biosci. 1999, 160, 191–213. [Google Scholar] [CrossRef]
Anderson, R.M.; May, R.M. Infectious Diseases in Humans; Oxford University Press: Oxford, 1992. [Google Scholar]
Martcheva, M. An Introduction to Mathematical Epidemiology; Springer, 2015. [Google Scholar]
Harris, J.E. Population-Based Model of the Fraction of Incidental COVID-19 Hospitalizations during the Omicron BA.1 Wave in the United States. COVID 2023, 3, 728–743. [Google Scholar] [CrossRef]
Whitehead, S.S.; Blaney, J.E.; Durbin, A.P.; Murphy, B.R. Prospects for a dengue virus vaccine. Nat Rev Microbiol 2017, 5, 518–528. [Google Scholar] [CrossRef] [PubMed]
Pastor-Satorras, R.; Castellano, C.; Mieghem, P.V.; Vespignani, A. Epidemic processes in complex networks. Rev. Mod. Phys. 2015, 87, 925–979. [Google Scholar] [CrossRef]
Pastor-Satorras, R.; Vespignani, A. Epidemic dynamics and endemic states in complex networks. Phys. Rev. E 2001, 63, 066117. [Google Scholar] [CrossRef]
Shudo, Y.O.Y.A. Microscopic Numerical Simulations of Epidemic Models on Networks. Mathematics 2021, 9, 932. [Google Scholar]
Basnarkov, L.; Tomovski, I.; Sandev, T.; Kocarev, L. on-Markovian SIR epidemic spreading model of COVID-19. Chaos, Solitons and Fractals 2022, 160, 112286. [Google Scholar] [CrossRef]
d’Onofrio, G.; Michelitsch, T.M.; Polito, F.; Riascos, A.P. On discrete-time arrival processes and related random motions. arXiv 2024, arXiv:2403.06821. [Google Scholar]
Metzler, R.; Klafter, J. The random walk’s guide to anomalous diffusion : A fractional dynamics approach. Phys. Rep. 2000, 339, 1–77. [Google Scholar] [CrossRef]
Mainardi, F.; Gorenflo, R.; Scalas, E. A fractional generalization of the Poisson processes. Vietnam J. Math. 2004, 32, 53–64. [Google Scholar]
Barabási, A.-L. Network Science; Cambridge University Press: Cambridge, 2016. [Google Scholar]
Barabási, A.-L.; Albert, R. Emergence of Scaling in Random Networks. Science 1999, 286, 509. [Google Scholar] [CrossRef] [PubMed]
Barabási; Albert, R. ; Jeong, H. Mean-field theory for scale-free random networks. Physica A 1999, 272, 173–187. [Google Scholar] [CrossRef]
Jeong, H.; Tombor, B.; Albert, R.; Oltvai, Z.N.; Barabási, A.L. The large-scale organization of metabolic networks. Nature 2000, 407, 651–4. [Google Scholar] [CrossRef]
Erdös, P.; Rényi, A. On Random Graphs I. Publ. Math. 1959, 6, 290–297. [Google Scholar] [CrossRef]
Erdös, P.; Rényi, P. On the evolution of random graphs. Publ. Math. Inst. Hung. Acad. Sci 1960, 5, 17. [Google Scholar]
Gilbert, E.N. Random Graphs. Annals Math. Sci. 1959, 30, 1141–1144. [Google Scholar] [CrossRef]
Ross, S.M. Stochastic Processes; John Wiley & Sons: New York, 1996. [Google Scholar]
Newman, M.E.J. Networks: An Introduction; Oxford University Press: Oxford, 2010. [Google Scholar]
Newman, M.E.J.; Watts, D.J.; Strogatz, S.H. Random graph models of social networks. PNAS 2002, 99, 2566–2572. [Google Scholar] [CrossRef]
Watts, D.J.; Strogatz, S.H. Collective dynamics of ‘small-world’ networks. Nature 1998, 393, 440. [Google Scholar] [CrossRef]
Eraso-Hernandez, L.K.; Riascos, A.P.; Michelitsch, T.M.; Wang-Michelitsch, J. Evolution of transport under cumulative damage in metro systems. Int. J. Mod. Phys. C 2023, 2450037. [Google Scholar] [CrossRef]
Barrat, A.; Barthélemy, M.; Vespignani, A. Epidemic spreading in population networks. In Dynamic Processes on Complex Networks; Cambridge University Press, 2008; pp. 180–215. [Google Scholar] [CrossRef]
Riascos, A.P.; Sanders, D.P. Mean encounter times for multiple random walkers on networks. Phys. Rev. E 2021, 103, 042312. [Google Scholar] [CrossRef] [PubMed]
Bestehorn, M.; Riascos, A.P.; Michelitsch, T.M.; Collet, B.A. A Markovian random walk model of epidemic spreading. Continuum Mech. Thermodyn. 2021, 33, 1207–1221. [Google Scholar] [CrossRef] [PubMed]
Michelitsch, T.M.; Riascos, A.P.; Collet, B.A.; Nowakowski, A.F.; Nicolleau, F.C.G.A. Fractional Dynamics on Networks and Lattices; ISTE/Wiley: London, 2019. [Google Scholar]
Riascos, A.P.; Mateos, J.L. Long-range navigation on complex networks using Lévy random walks. Phys. Rev. E 2012, 86, 056110. [Google Scholar] [CrossRef] [PubMed]
Bestehorn, M.; Michelitsch, T.M.; Collet, B.A.; Riascos, A.P.; Nowakowski, A.A.F. Simple model of epidemic dynamics with memory effects. Phys. Rev. E 2022, 105, 024205. [Google Scholar] [CrossRef] [PubMed]
Granger, T.; Michelitsch, T.M.; Bestehorn, M.; Riascos, A.P.; Collet, B.A. Four-compartment epidemic model with retarded transition rates. Phys. Rev. E 2023, 107, 044207. [Google Scholar] [CrossRef]
Bestehorn, M.; Michelitsch, T.M. Oscillating Behavior of a Compartmental Model with Retarded Noisy Dynamic Infection Rate. Int. J. Bifurcation Chaos 2023, 33, 2350056. [Google Scholar] [CrossRef]
van Kampen, N.G. Stochastic processes in chemistry and physics; North Holland: Amsterdam, 1981. [Google Scholar]
Mieghem, P.V. Exact Markovian SIR and SIS epidemics on networks and an upper bound for the epidemic threshold. arXiv 2014, arXiv:1402.1731. [Google Scholar]
Zhu, Y.; Shen, R.; Dong, H.; Wang, W. Spatial heterogeneity and infection patterns on epidemic transmission disclosed by a combined contact-dependent dynamics and compartmental model. PLoS ONE 2023, 18, e0286558. [Google Scholar] [CrossRef]
Gostiaux, L.; Bos, W.J.T.; Bertoglio, J.-P. Periodic epidemic outbursts explained by local saturation of clusters. Phys. Rev. E 2023, 107, L012201. [Google Scholar] [CrossRef]
Peyrard, M. What can we learn from the dynamics of the Covid-19 epidemic? arXiv 2023, arXiv:2308.14090. [Google Scholar] [CrossRef] [PubMed]
Soper, H.E. The interpretation of periodicity in disease prevalence. J. Royal Statistical Society 1929, 92, 34–61. [Google Scholar] [CrossRef]
Noh, J.-D.; Rieger, H. Random walks on complex networks. Phys. Rev. Lett. 2004, 92. [Google Scholar] [CrossRef] [PubMed]
Michelitsch, T.; Riascos, A.P.; Collet, B.A.; Nowakowski, A.; Nicolleau, F. Fractional Dynamics on Networks and Lattices; ISTE-Wiley: London, 2019. [Google Scholar]
Schilling, R.; Song, R.; Vondraček, Z. Bernstein functions. In Theory and Applications, Studies in Mathematics, 37; de Gruyter: Berlin, 2010. [Google Scholar]
Simon, C.M. The SIR dynamic model of infectious disease transmission and its analogy with chemical kinetics. PeerJ Physical Chemistry 2020. [Google Scholar] [CrossRef]
Supplementary materials: PYTHON codes (© Téo Granger 2023) and animated films. Available online: https://sites.google.com/view/scirs-model-supplementaries/accueil.

1	The network `distance’ of two nodes is the number of edges of the shortest path connecting them.
2	Free to download and non-commercial use.
3	A similar consideration of function $G (μ)$ of (39) shows as well that the healthy state for $R_{0} < 1$ does not exhibit an oscillatory instability.

Figure 1. Left frame: Gamma distribution for four different cases: Weakly singular (at

t = 0

) [

〈 t 〉 = 0.5

ξ = 0.7

], exponential [

〈 t 〉 = 2

ξ = 0.5

], broad [

〈 t 〉 = 1.5

ξ = 4

], and narrow [

〈 t 〉 = 1.5

ξ = 30

]. Right frame: Their Persistence (survival) probability distributions of Equation (27) where the same color code is used.

Figure 1. Left frame: Gamma distribution for four different cases: Weakly singular (at

t = 0

) [

〈 t 〉 = 0.5

ξ = 0.7

], exponential [

〈 t 〉 = 2

ξ = 0.5

], broad [

〈 t 〉 = 1.5

ξ = 4

], and narrow [

〈 t 〉 = 1.5

ξ = 30

]. Right frame: Their Persistence (survival) probability distributions of Equation (27) where the same color code is used.

Figure 2. Endemic states of infected walkers/nodes

J_{w, n} (\infty) = (R_{0} - 1) / (R_{0} + r)

versus

R_{0}

for various values of parameter r which has to be read

r = β_{n} 〈 t_{I}^{n} 〉

(

r = β_{w} 〈 t_{I}^{w} 〉

) for the walker’s (node’s) endemic states.

Figure 2. Endemic states of infected walkers/nodes

J_{w, n} (\infty) = (R_{0} - 1) / (R_{0} + r)

versus

R_{0}

for various values of parameter r which has to be read

r = β_{n} 〈 t_{I}^{n} 〉

(

r = β_{w} 〈 t_{I}^{w} 〉

) for the walker’s (node’s) endemic states.

Figure 3.

G (μ)

of (39) for a some Gamma distributed

t_{I}^{w, n}

. Positive zeros of

G (μ)

exist only for

R_{0} > 1

(instability of globally healthy state).

Figure 3.

G (μ)

of (39) for a some Gamma distributed

t_{I}^{w, n}

. Positive zeros of

G (μ)

exist only for

R_{0} > 1

(instability of globally healthy state).

Figure 4.

G_{e} (μ)

of (42) for different values of

R_{0}

where

G_{e} (μ) > 0

for

R_{0} > 1

(stability of the endemic state).

Figure 4.

G_{e} (μ)

of (42) for different values of

R_{0}

where

G_{e} (μ) > 0

for

R_{0} > 1

(stability of the endemic state).

Figure 5. We depict function

G_{M}

of Equation (47) for a few values of

R_{M}

for exponentially distributed

t_{I}^{w, n}, t_{M}

. The basic reproduction number

R_{M}

is monotonously decreasing with increasing mortality parameter

ξ_{M}

(see Figure 6). The parameters are

β_{w, n} = 1

α_{I}^{w} = 1

ξ_{I}^{w} = 1

α_{n} = 1

ξ_{I}^{n} = 0.5

with

R_{0} = 2

where here

〈 t_{I, M} 〉 = R_{0} / (1 + ξ_{M})

Figure 5. We depict function

G_{M}

of Equation (47) for a few values of

R_{M}

for exponentially distributed

t_{I}^{w, n}, t_{M}

. The basic reproduction number

R_{M}

is monotonously decreasing with increasing mortality parameter

ξ_{M}

(see Figure 6). The parameters are

β_{w, n} = 1

α_{I}^{w} = 1

ξ_{I}^{w} = 1

α_{n} = 1

ξ_{I}^{n} = 0.5

with

R_{0} = 2

where here

〈 t_{I, M} 〉 = R_{0} / (1 + ξ_{M})

Figure 6. Basic reproduction number

R_{M}

of Equation (50) versus mortality rate parameter

ξ_{M}

for Gamma distributed

t_{I}^{w, n}, t_{M}

for various

α_{M}

where we have set

β_{n} = β_{w} = 〈 t_{I}^{w} 〉 = 〈 t_{I}^{n} 〉 = 1

, (

α_{I}^{w} = ξ_{I}^{w} = 0.3

) and

α_{M} = 1

α_{w} = ξ_{I}^{w} = 1

for the Markovian case which is recovered by Equation (53).

Figure 6. Basic reproduction number

R_{M}

of Equation (50) versus mortality rate parameter

ξ_{M}

for Gamma distributed

t_{I}^{w, n}, t_{M}

for various

α_{M}

where we have set

β_{n} = β_{w} = 〈 t_{I}^{w} 〉 = 〈 t_{I}^{n} 〉 = 1

, (

α_{I}^{w} = ξ_{I}^{w} = 0.3

) and

α_{M} = 1

α_{w} = ξ_{I}^{w} = 1

for the Markovian case which is recovered by Equation (53).

Figure 7. Barabási-Albert, Erdös-Renyi and Watts-Strogatz types with 300 nodes and connectivity parameters used in some of the simulations. The WS graph for connectivity parameter

m = 2

lacks the small world property resembling a complex real world structure. The ER network has a broad degree distribution and the small world property. The BA graph is for

N \to \infty

asymptotically scale-free with a power law degree distribution and the small world feature where a large number of nodes have small degrees and a few (hub) nodes with very large degrees. Almost all nodes are only a few links away from hub nodes.

Figure 7. Barabási-Albert, Erdös-Renyi and Watts-Strogatz types with 300 nodes and connectivity parameters used in some of the simulations. The WS graph for connectivity parameter

m = 2

lacks the small world property resembling a complex real world structure. The ER network has a broad degree distribution and the small world property. The BA graph is for

N \to \infty

Figure 8. Snapshots of spreading in a WS graph (

Z = 2000

walkers,

N = 2000

nodes, connectivity parameter

m = 2

) and mortality parameter

ξ_{M} = 0.4

with

D (\infty) \approx 16 %

. The remaining parameters are the same as in Figure 10. One observes

d_{w} (\infty) \approx 0.99

and

S_{w} (\infty) \approx 1 %

with only about 20 surviving walkers after extinction of the disease. S walkers are in cyan color, I walkers red, D walkers invisible and nodes without walkers are represented in black. Consult here an animated video of this process.

Figure 8. Snapshots of spreading in a WS graph (

Z = 2000

walkers,

N = 2000

nodes, connectivity parameter

m = 2

) and mortality parameter

ξ_{M} = 0.4

with

D (\infty) \approx 16 %

. The remaining parameters are the same as in Figure 10. One observes

d_{w} (\infty) \approx 0.99

and

S_{w} (\infty) \approx 1 %

Figure 9. The plots show the evolution on WS graph with

Z = 1000

walkers for connectivity parameter

m = 8

and rewiring probability

p = 0.7

(

n x . c o n n e c t e d_w a t t s_s t r o g a t z_g r a p h (N = 1000, m = 8, p = 0.7, s e e d = s e e d)

) without mortality (left frame) and with mortality (right frame).

t_{I}^{w, n}, t_{M}

are Gamma distributed with the parameters

〈 t_{M} 〉 = 14

ξ_{M} = 2

〈 t_{I}^{w} 〉 = 8

ξ_{I}^{w} = 10

, and

〈 t_{I}^{n} 〉 = 15

ξ_{n} = 10^{5}

, see histogram.

D (\infty) \approx 0.01

and is determined by numerical integration of (7).

Figure 9. The plots show the evolution on WS graph with

Z = 1000

walkers for connectivity parameter

m = 8

and rewiring probability

p = 0.7

(

n x . c o n n e c t e d_w a t t s_s t r o g a t z_g r a p h (N = 1000, m = 8, p = 0.7, s e e d = s e e d)

) without mortality (left frame) and with mortality (right frame).

t_{I}^{w, n}, t_{M}

are Gamma distributed with the parameters

〈 t_{M} 〉 = 14

ξ_{M} = 2

〈 t_{I}^{w} 〉 = 8

ξ_{I}^{w} = 10

, and

〈 t_{I}^{n} 〉 = 15

ξ_{n} = 10^{5}

, see histogram.

D (\infty) \approx 0.01

and is determined by numerical integration of (7).

Figure 10. Evolution on WS graph with

Z = 1000

walkers and

N = 1000

nodes for the same parameters as in Figure 9 but with reduced connectivity parameter

m = 2

. The upper frame shows a snapshot (

t = 15

) of the spreading in one random walk realization (susceptible walkers green, susceptible nodes black, infected nodes red).

Figure 10. Evolution on WS graph with

Z = 1000

walkers and

N = 1000

nodes for the same parameters as in Figure 9 but with reduced connectivity parameter

m = 2

. The upper frame shows a snapshot (

t = 15

) of the spreading in one random walk realization (susceptible walkers green, susceptible nodes black, infected nodes red).

Figure 11. Evolution on ER graph (

n x . e r d o s_r e n y i_g r a p h (N = 1000, p = 0.1, s e e d = s e e d)

) with

Z = 1000

walkers and small rewiring probability

p = 0.1

(above the percolation limit

p_{c} = 0.01

to ensure a connected structure). The parameters are

〈 t_{I}^{n} 〉 = 5

ξ_{I}^{n} = 10

〈 t_{I}^{w} 〉 = 10

ξ_{I}^{w} = 0.05

〈 t_{M} 〉 = 65

ξ_{M} = 1

. The left upper frame shows a snapshot of the evolution (same color code as in Figure 10).

Figure 11. Evolution on ER graph (

n x . e r d o s_r e n y i_g r a p h (N = 1000, p = 0.1, s e e d = s e e d)

) with

Z = 1000

walkers and small rewiring probability

p = 0.1

(above the percolation limit

p_{c} = 0.01

to ensure a connected structure). The parameters are

〈 t_{I}^{n} 〉 = 5

ξ_{I}^{n} = 10

〈 t_{I}^{w} 〉 = 10

ξ_{I}^{w} = 0.05

〈 t_{M} 〉 = 65

ξ_{M} = 1

. The left upper frame shows a snapshot of the evolution (same color code as in Figure 10).

Figure 12. Evolution on BA graph with

Z = 50

walkers and

N = 5000

nodes (

n x . b a r a b a s i_a l b e r t_g r a p h (N = 5000, m = 5, s e e d = s e e d)

) with parameters

〈 t_{I}^{n} 〉 = 32

ξ_{I}^{n} = 10^{4}

〈 t_{I}^{w} 〉 = 8

ξ_{I}^{w} = 10^{4}

(sharp

t_{I}^{w, n}

t_{M} = 500

ξ_{M} = 10^{- 3}

. The basic reproduction number

R_{M}

is here only slightly smaller than

R_{0}

without mortality. The left upper frame shows a snapshot of the evolution (same color code as in Figure 11).

Figure 12. Evolution on BA graph with

Z = 50

walkers and

N = 5000

nodes (

n x . b a r a b a s i_a l b e r t_g r a p h (N = 5000, m = 5, s e e d = s e e d)

) with parameters

〈 t_{I}^{n} 〉 = 32

ξ_{I}^{n} = 10^{4}

〈 t_{I}^{w} 〉 = 8

ξ_{I}^{w} = 10^{4}

(sharp

t_{I}^{w, n}

t_{M} = 500

ξ_{M} = 10^{- 3}

. The basic reproduction number

R_{M}

is here only slightly smaller than

R_{0}

without mortality. The left upper frame shows a snapshot of the evolution (same color code as in Figure 11).

Figure 13. Evolution with the same parameters and number of walkers (

Z = 50

) as in Figure 12 but less nodes (

N = 2100

) for one random walk realization. We interpret the increase of

R_{M}

and

R_{0}

due to more frequent passages of susceptible walkers on infected nodes (higher infection rates).

Figure 13. Evolution with the same parameters and number of walkers (

Z = 50

) as in Figure 12 but less nodes (

N = 2100

) for one random walk realization. We interpret the increase of

R_{M}

and

R_{0}

due to more frequent passages of susceptible walkers on infected nodes (higher infection rates).

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

MDPI Initiatives

Important Links

Choose an area of interest and we will send you notifications of new preprints at your preferred frequency.

Disclaimer

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

Abstract

1. Introduction

2. Compartmental Model with Mortality

Remark I

Markovian (memoryless) case

A few more words on waiting time distributions

3. Endemic Equilibrium for Zero Mortality

4. Stability Analysis of Endemic and Healthy State without Mortality

5. Stability Analysis of the Healthy State with Mortality

6. Random Walk Simulations

Case study and discussion

7. Conclusions

Acknowledgments

Appendix A

Appendix A.1. Some Basic Notions

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

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

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

References

MDPI Initiatives

Important Links

Subscribe