key: cord-0884109-2vs57w5o
authors: Bestehorn, Michael; Riascos, Alejandro P.; Michelitsch, Thomas M.; Collet, Bernard A.
title: A Markovian random walk model of epidemic spreading
date: 2021-01-16
journal: Continuum Mech
DOI: 10.1007/s00161-021-00970-z
sha: affeae344860738a6dd4b9e6127dd7246095fce9
doc_id: 884109
cord_uid: 2vs57w5o

We analyze the dynamics of a population of independent random walkers on a graph and develop a simple model of epidemic spreading. We assume that each walker visits independently the nodes of a finite ergodic graph in a discrete-time Markovian walk governed by his specific transition matrix. With this assumption, we first derive an upper bound for the reproduction numbers. Then, we assume that a walker is in one of the states: susceptible, infectious, or recovered. An infectious walker remains infectious during a certain characteristic time. If an infectious walker meets a susceptible one on the same node, there is a certain probability for the susceptible walker to get infected. By implementing this hypothesis in computer simulations, we study the space-time evolution of the emerging infection patterns. Generally, random walk approaches seem to have a large potential to study epidemic spreading and to identify the pertinent parameters in epidemic dynamics.

worldwide pandemic context of COVID-19 is boosting an additional interest to this topic [7, 8] . Epidemic spreading in complex networks was studied by several authors [9] [10] [11] and the epidemic dynamics in scale-free networks was analyzed in [12] . In a recent paper, the effects of quarantine measures to epidemic spreading in activity-driven adaptive temporal networks were studied [13] including percolation effects in epidemic spreading in small-world networks [14, 15] . A renormalization group approach has been employed to model the second COVID-19 wave in Europe [16] , just to quote a few examples.

Strongly driven by the present world-wide COVID-19 spreading there is a huge and urgent need of reliable models that are able to capture essential aspects of the space-time dynamics of infectious diseases allowing to develop preventive strategies. For an overview of the present world-wide COVID-19 situation as far known, we refer to [17] .

Infectious diseases such as measles, mumps, and rubella can be studied in the framework of nonlinear dynamical systems. For the most simple case of spatially homogeneous infection rates, SIR models have been applied successfully in the past [18, 19] . As mentioned, SIR stands for the three compartments susceptibleinfected-recovered into which the individuals are grouped, depending on their state. A susceptible individual (S) can be infected and become ill (I). After a certain time τ 1 , it will recover and be removed from the system (R) in the subsequent computer simulation model. During time τ 1 , it can infect other susceptible individuals. The mathematical description in the SIR model is achieved by an ordinary first-order differential equation for each rate. If spatial effects are taken into account, the rates can be assumed space-dependent and a set of three nonlinear coupled diffusion equations can be derived [6] .

Instead of using partial differential equations, the individuals or particles can be considered as independent random walkers on a discrete network with a given architecture. Our model is based on the following assumptions. The particles perform random jumps from one node to another connected node of the network. If on the same node an infected particle meets a susceptible one, the susceptible walker may be infected with a given probability P. To describe the process of recovery, each particle has an inner variable parametrizing its state. This variable changes in course of time. If 'time' is assumed to be discrete, the whole dynamics on the network and of the inner variable can be formulated as a (nonlinear) mapping from one-time step to the next. The system has no memory, its state is uniquely defined by the positions of the particles and the values of their inner variables at a certain time step (Markov process).

Our paper is organized as follows. In the subsequent Sect. 2, we give a brief general introduction into the dynamics of Z independent Markovian random walkers on finite connected (ergodic) graphs. Without loss of generality, we confine us here to undirected graphs. We utilize the Markovian walk approach to derive an upper bound for the so-called basic reproduction number R 0 which is defined subsequently. In this part, we consider the situation when there is a single infected walker and Z − 1 susceptible walkers in the network. We derive explicit formulae for the expected number of times the infected walker meets a susceptible one which defines an upper bound for R 0 .

In Sect. 3, we perform numerical simulations employing above mentioned assumptions to generate spacetime patterns of the susceptible/infected walkers where we consider Z independent walkers on a finite 2D square lattices with variable adjacency matrices and connectivity. In this way, we explore how the architecture of a network affects the space-time dynamics of the epidemic spreading and identify pertinent parameters governing the space-time patterns in order to establish predictive measures such as confinement and social distance rules.

In the present section, we recall some basic features of random walks with independent multiple walkers on the network [20] (See also [4, 21, 22] for outlines and analysis of the emergent space-temporal dynamics which we employ in our model). We focus on unbiased Markovian walks, however, this approach can be generalized to biased walks on directed graphs and also to continuous-time random walks (CTRWs). We consider Z independent random walkers r = 1, . . . Z on a connected undirected network of p = 1, . . . , N nodes. Despite the results of the present section can be derived in a simpler way, the approach recalled here allows to be applied to elaborate more sophisticated models such as for instance when the walkers perform independent CTRWs.

We assume that the walkers move independently through the network where the jumps of a walker r are governed by his own N × N one-step transition matrix W (r ) with the elements W (r ) i j (i, j = 1, . . . , N , r = 1, . . . , Z ) indicating the probability of the transition between the nodes i → j in one jump with N j=1 W (r ) i j = 1 and 0 ≤ W (r ) i j ≤ 1, i.e., per construction the transition matrices W (r ) are row-stochastic. Further, we assume that all walkers jump synchronously at integer times t = 1, 2, . . . ∈ N and occupy at t = 0 their respective departure nodes. Performing its nth jump at t = n, each walker remains during t ∈ [n, n +1) on the node he has reached at t = n ∈ N 0 . For convenience, we employ Dirac's bra|−|ket notation with |i = |i 1 , i 2 , . . . , i Z = |i 1 |i 2 . . . |i Z where i r indicates the node occupied by walker r . We refer |i to as 'state-vector' containing the positions of the walkers in the network. The collective dynamics of the Z independent walkers is then characterized by the collective one-step transition matrix 1

Assuming walker r starts at t = 0 at node i r , then the probability to find the walker r on node j r at time t is given by P

where for t = 0 we assume here the initial condition P (r ) i r j r (t)| t=0 = δ i r j r . In this relation we assume that each walker r = 1, . . . , Z moves independently through the graph in a Markovian walk governed by the master equation

thus P (r )

i j (t) = i|W (r ) ) t | j indicates the probability of walker r to reach node j in t jumps when departing at t = 0 from node i. When all walkers hop synchronously at t ∈ N 0 the probability to find the Z walkers in the state |j = | j 1 , j 2 , . . . , j Z at time t becomes

with P (r ) (i, j, t)| t=0 = δ i, j = Z r =1 δ i r j r . In order to develop such a model we are interested in the 'stateprobabilities,' i.e., the probabilities that the nodes j = 1, . . . , N are occupied by s 1 , . . . , s N ( N j=1 s j = Z ) walkers. For our convenience, we introduce the following generating functions

with G r i r (u 1 = 1, . . . , u N = 1, t) = N j=1 P (r ) i j (t) = 1 reflecting normalization. Now consider the collective generating function

which is a multinomial of total degree Z . The coefficients (s i = 0, 1, . . . , Z )

indicate the state-probabilities, i.e., the probabilities that the nodes 1, 2, . . . , N at time t and with the given initial condition are occupied by s 1 , s 2 , . . . , s N walkers (where s 1 +s 2 +· · ·+s N = Z recovers the total number of walkers). We observe that G i (u 1 , . . . , u N , t) u 1 =...=u N =1 = 1, i.e., (7) indeed is a normalized distribution. We further observe that since (6) is a homogeneous function of total degree Z ,

As an important case let us consider when all walkers have identical transition matrix W

with the state-probabilities given by the multinomial-coefficients

where s 1 + s 2 + · · · + s N = Z and s j ∈ [0, Z ]. Case: N = 2 For illustration, let us consider a network of two nodes i = 1, 2 (N = 2) where we have Z independent walkers and let us assume the initial condition i r = 1 for all Z walkers. Let us assume all walkers have the same transition matrix W (r ) i j = W i j . Then, the collective generating function (6) is given by

where

Hence, the state-probabilities, i.e., probabilities that (with the given initial condition) at time t node 1 is occupied by s walkers and node 2 by Z − s walkers are obtained as

The normalization of the state-probability distribution again is easily verified

Now we need to relate the architecture of the graph with its random walk features. The information of the topology of an undirected graph is contained in the one-step transition matrix [1, 24] 

where we assume that each walker undertakes jumps on the graph governed by the same one-step transition matrix. In (14), we introduced the N × N Laplacian matrix

which contains the adjacency matrix A i j with A i j = 1 if the nodes i, j are connected by an edge and A i j = 0 else. Further, we do not allow self-connections which is expressed by A ii = 0. In undirected networks the edges do not have a direction, i.e., the adjacency matrix and the Laplacian matrix are symmetric. Further important is the degree K i of a node i which counts the number of nodes connected with i, namely

where the condition K i > 0 tells us that there are no isolated disconnected nodes. With (15) , the transition matrix (14) can also be written as

where we directly verify row-stochasticity N j=1 W i j = 1. Per construction, we have W ii = 0 thus the walkers at any time step have to move and change the node. The transition matrix is non-symmetric if there are nodes with variable degree K i = K j . For later use, we introduce the canonical representation (For a detailed spectral analysis of spectral properties, see [24] )

where |Φ s and Φ s | denote the right-and left eigenvectors of W, respectively, and we assume an aperiodic ergodic (connected) network with the eigenvalue structure |λ s | ≤ 1 with real eigenvalues λ s ∈ R where the largest unique (Frobenius-) eigenvalue is λ 1 = 1 and −1 < λ m < 1 for m = 2, . . . N . We thus have the unique stationary distribution

where K is called the total degree and K denotes the average degree of the network. It is important to notice that in (aperiodic) ergodic (i.e., connected) networks the stationary distribution has uniquely (nonzero) positive elements W

> 0 and is given by the normalized degrees independent of the departure node i. The stationary transition matrix is a matrix consisting of identical rows (see, e.g., [24] for an analysis of the related spectral properties of the transition matrix in ergodic graphs). Having recalled these general features, we can now use these properties to derive estimates for the reproduction numbers which are key quantities in epidemic models.

We now consider the situation of Z independent walkers where one walker is infectious in the time interval 0 ≤ t ≤ τ 1 . We denote the infectious walker by r = 1 and Z − 1 walkers (denoted by r = 2, . . . Z ) are susceptible. For later use let us introduce the 'effective reproduction number' R e (τ 1 ) as the number of infections an infectious walker causes up to time τ 1 while he is infectious. Apart of this quantity the so-called basic reproduction number R 0 (τ 1 ) is of interest. R 0 (τ 1 ) indicates the number of newly infected walkers (up to time τ 1 ) by one infected walker under the assumption the infected walker meets only susceptible walkers. In fact R e also depends on time by the time-dependence of the number of susceptible walkers. In the present part, in order to derive an upper bound, we ignore this time-dependence. On the other hand, the quantity R 0 ignores the fact that an infectious walker does not only meet susceptible ones, but also infected and recovered walkers. Therefore, R 0 ≥ R e , i.e., the basic reproduction number overestimates the 'real' effective reproduction number R e . For R e > 1, the number of infected walkers is increasing. If R e > 1 is persisting over longer times, then we are in the regime of (exponential) epidemic spreading. For R e = 1, the number of infected walkers remains stable, and for R e < 1, the number of infected walkers is decreasing and when persisting over longer times then the epidemics dies out. Now for the sake of simplicity in the formulas to be derived, we assume for the susceptible walkers random initial conditions and stationary distributions, namely

independent of time. In order to get an upper bound for the basic reproduction number, we are now interested in the expected number of timesR(τ 1 ) the infectious walker meets another walker (no matter whether or not susceptible) during the time τ 1 of his infection. ClearlyR(τ 1 ) ≥ R 0 (τ 1 ), i.e.,R(τ 1 ) represents an upper bound for the basic reproduction number R 0 (τ 1 ). The quantityR(τ 1 ) ignores also the fact that the infectious walker may multiply meet the same susceptible walker. We come back to the issue of variable 'susceptibility' with a probability P of infection as a crucial parameter later on. For the susceptible walkers in the stationary state, the generating function (10) becomes independent of their initial nodes and of time (as we ignore transitions from susceptible to the infectious state) and takes the form

The 'state-probabilities' that s j susceptible walkers are on node j ( j = 1, . . . N ) with j s j = Z − 1 then are obtained as

with the stationary distribution W (∞) j = K j N K . Now we assume that the duration of the infection is τ 1 ∈ N and that each walker performs jumps exactly at integer times t ∈ N. Accounting for the fact that the infectious walker remains on his departure node during the time-interval [0, 1) and performs its first jump at t = 1, then it follows that the infectious walker during his infection, i.e., within the time interval [0, τ 1 ) performs τ 1 − 1 jumps where at each jump he meets susceptible walkers in the stationary distribution (23) . In our calculation, we ignore the transitions of susceptible walkers to the infectious state and assume the number of susceptible walkers remains constant Z −1. The expected number of times R(τ 1 ) the infectious walker meets a susceptible one within the time-interval [0, τ 1 ) then is obtained as (where we assume the infectious walker has departure node i and transition probabilities at time t: P

where

The quantity r (t) indicates the expected number of susceptible walkers met by the infectious one in the time increment Δt = 1 following to his tth jump and we observe that

where T (1) i j (τ ) = τ −1 t=0 P (1) i j (τ ) indicates the expected sojourn time of the infectious walker (with departure node i) on node j in a walk of τ − 1 time steps (i.e., in a walk of duration [0, τ )). For a detailed analysis of this issue, consult [24] . For τ 1 = 0, we have with P

which are at t = 0 the exact values for the effective and basic reproduction numbers since per construction at t = 0 the infectious walker meets on his departure node r (0) = (Z − 1)W (∞) i susceptible walkers. We also can define a global value by averaging (26) over all departure nodes of the infectious walker, namelŷ

It is worthy to consider above result for regular networks, i.e., networks with constant degree K j = K = K (i = 1, . . . N ). Then, we get for (26) which coincides then with (27) the simple expression

where we have used W (∞) j = 1 N and normalization N j=1 P (1) i j (t) = 1 where ρ s = Z −1 N denotes the density of the susceptible walkers.

In the previous section, we ignored the transitions between the states susceptible, infectious, and recovered. In the present section, we present numerical simulations of space-time patterns of infectious/susceptible walkers where we account for transitions between them.

We consider again Z independent random walkers (particles) performing independent jumps at integer times on a two-dimensional undirected graph with N = L 2 nodes where (x (n) i , y (n) i ) indicate the position of walker ('particle') i at time n, namely

where x i , y i , L are integer numbers. Let the walkers jump according to

Here, ξ x,y are equally distributed random integer numbers ξ in [−h, h]. In our simple network, each node has d = (2h + 1) 2 accessible neighbors, where d is the degree of a node. The velocity of each walker (mean distance in one step) is given asv

Let s i be the 'grade of infection' of walker i. Due to recovery, we assume a simple linear decrease < s 1 , we define particle i to be immune. For infection, the following rule applies. If two particles i, j meet on the same node, i.e.,

j ≤ 0 then particle i infects particle j with a given probability P. If particle j gets infected at time-step n, we set

Thus, we may identify three regions (Fig. 1 ):

1. s 1 ≤ s i ≤ 1: particle i is infectious and infects particle j with probability P (duration of infectibility τ 1 ). 2. 0 < s j < s 1 : particle j is immune and cannot be infected by particle i (duration of immunity τ 2 − τ 1 ). 3. s j ≤ 0: particle j is healthy (again) and can be (re)-infected.

From Fig. 1 , the relations

follow. Here, τ 1 is the time while a particle can infect another one, τ 2 denotes the time where a particle is not susceptible after infection (time of infectibility plus time of immunity after recovering). The period of immunity after recovering is τ 2 − τ 1 ≥ 0. In the present model, we assume the characteristic times τ 1,2 to be the same for all infected and immune particles, respectively. After the time τ 2 , a particle is again susceptible and can be re-infected. If τ 2 → ∞, particles stay immune forever after recovering.

The basic reproduction number R 0 as mentioned above is defined as the number of particles that are infected by one particle under the assumption that all other particles are healthy and susceptible. The probability for a particle to meet another one during one time increment Δt = 1 is equal to the density (where we assume Z , N 1, see relation (28) for τ 1 = 1)

To find R 0 , this quantity must be multiplied with the time of infectivity τ 1 and with the probability of infection P to obtain

This simple relation indeed is consistent with expression (28) of the previous section by introducing the probability of infection P. Given τ 1 and P, R 0 is a constant. However, in real life due to hygiene measures P may vary considerably in time but also in space, leading to an inhomogeneously distributed R 0 . Distance rules or lockdowns may rather restrict the mobility of the particles and can be considered by changing the velocity (31) or the connectivity of the network. The effective reproduction number is found by replacing the particle number in (35) by the number of those particles which are not infected or not immune 

s , they can be either ill or immune. In course of time, Z s and therefore R s decreases. If R e = 1, herd immunity is reached and from (36) one finds

From the Z Up to here, we assumed an average (stationary) particle distribution over the nodes. However, if clusters of infected particles are formed, Z s may vary strongly in space thus this assumption does not any more hold true. For an isolated cluster in an elsewhere healthy environment, R e may be locally around one and the number of ill particles saturates due to herd immunity, where in the healthy regions R e can be much larger than one.

Spatial patterns are expected if the initial distribution of infected walkers is localized (clusters). Let us assume that K particles form a cluster in the central node of the layer and that all K particles are infected:

with ξ randomly distributed in [s 1 , 1]. The other N − K particles are healthy and randomly distributed over all nodes:

and η x , η y as random integers in [1, L] .

We present numerical solutions of the system with the fixed parameters N = 30000, L = 1500, N = 2.25 · 10 6 , h = 4, τ 1 = 600, τ 2 = 2400. Figures 2 and 3 show the situation for P = 0.4, leading to a basic R-number of R 0 = 3.2. The thin black line denotes herd immunity. For Fig. 4 , P was much higher, P = 1 and the virus dies out after some sweeps (Fig. 5) .

Depending on P, but also on the mean particle velocityv, different pattern scenarios can be obtained. For the case of smallv = 1.15, corresponding to h = 1 and small P = 0.2 clusters are formed independently from the initial condition (Fig. 6) . The clusters do not connect and large areas of the domain remain healthy. As a consequence, the average number of infected walkers stays relatively low (Fig. 7) . If P or h is increased, the cluster size increases and the clusters connect (percolation point). Then, the number of infected particles increases also strongly.

To characterize cluster formation, we define an inhomogeneity factor f H that is zero if a pattern is completely homogeneous (constant) in space and that becomes large if clusters are formed. Therefore, we introduce a coarse mesh over the domain with 10 × 10 cells and count the number of ill particles laying in each cell with X i , where i = 1, . . . , 100. Then, we compute the normalized variance

where brackets denote the average over all 100 coarse cells. Figure 8 shows f H over time for the situation plotted in Figs. 6 and 7. If clusters are formed, f H increases, but after P is increased, the clusters grow and f H tends to small values, showing that the pattern becomes more and more homogeneous.

In the present paper, we first have developed a simple Markovian random walker model of epidemic spreading in undirected graphs. We derived an upper bound for the reproduction numbers R 0 (τ 1 ) and R e (τ 1 ) in a multiple random walker model where among Z independent random walkers one is infectious and Z −1 are susceptible. We derived the expected number of times the infectious walker meets another (susceptible) walker (relations (26) and (27)) where this quantity constitutes an upper bound for the basic reproduction number. Further, we performed computer simulations of the space-time evolution patterns on a 2D network. We showed that these space-time patterns depend sensitively on the infection probability P but also crucially depend on the characteristic times of infectivity τ 1 and duration of immunity τ 2 − τ 1 after recovering.

Despite its considerable simplicity, the present model allows predictions on the effect of lockdowns and distance rules. For future research, it would be interesting to see what happens in the space-time epidemic dynamics when τ 1,2 become random variables drawn from waiting-time densities such as for instance exponential or Mittag-Leffler with heavy power-law tails and non-Markovian long memory features. An exponential decay in the distribution of τ 2 − τ 1 describes the situation of short-time immunity whereas distributions with heavy power-law tails correspond to long-time immunity. In this way, effects of 'genetic stability' of a virus and its mutation activity could be taken into account. Such models could be important to obtain scenarios for the efficiency of vaccinations. Another interesting feature is introduced by the space-time fractional dynamics of the walkers on biased networks such as analyzed in recent papers [23, 25] . Although epidemic spreading has been widely addressed in many works, there are still many open questions such as effects of social distancing, lockdowns, and others calling for further thorough analysis.

Networks: An Introduction

Network Science

Random Walks and Random Environments

Emergence of encounter networks due to human mobility

The random Walk's guide to anomalous diffusion : a fractional dynamics approach

An Introduction to Mathematical Epidemiology

Recurrent host mobility in spatial epidemics: beyond reaction-diffusion

Epidemic spreading in heterogeneous networks with recurrent mobility patterns

Epidemic dynamics and endemic states in complex networks

Epidemic spreading in scale-free networks

Epidemic processes in complex networks

Epidemics and immunization in scale-free networks

Active and inactive quarantine in epidemic spreading on adaptive activity-driven networks

Epidemics and percolation in small-world networks

Scaling and percolation in the small-world network model

Second wave COVID-19 pandemics in Europe: a temporal playbook

A contribution to the mathematical theory of epidemics

Population biology of infectious diseases: Part I

Mean encounter times for multiple random walkers on networks

Temporal networks

temporal network theory: a colloquium

Non-local biased random walks and fractional transport on directed networks

Fractional Dynamics on Networks and Lattices

Biased continuous-time random walks with Mittag-Leffler jumps. Fractal Fract

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations

Acknowledgements M.B. gratefully acknowledges to have been hosted at the Institut Jean le Rond d'Alembert (Paris) during September 2020 for the aim of the present study.

Conflicts of interest The authors declare that they have no conflict of interest.