Within the last two decades, network science has become a huge interdisciplinary field [1, 2, 3] recently driven by the significant upswing of online (social) networks and search engines with a burst of works focusing on human mobility and encounter networks [4] . It turned out that random walks in networks are especially powerful to cover spreading and diffusion phenomena widely observed in nature. These diffusion phenomena include so-called 'anomalous diffusion' which have been successfully described by space-time fractional partial differential diffusion equations [5] . On the other hand within the last two decades an impressive amount of scientific work has been devoted to epidemic spreading models. For an introduction in epidemic modeling and state-of-the-art models such as the 'SIR model' (S = susceptible, I = infected, R = recovered) we refer to [6] . It is natural that the present 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 susceptible-infected-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 Section 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 Section 3 we perform numerical simulations employing above mentioned assumptions to generate space-time 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 continuoustime 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 ij ≤ 1, i.e. per construction the transition matrices W (r) are row-stochastic. Further we assume that all walk-ers 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
where for t = 0 we assume here the initial condition P irjr (t)| t=0 = δ ir jr . 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
ij (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 δ ir jr . In order to develop such a model we are interested in the 'state-probabilities', 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
ij (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 (7) indeed is a normalized distribution. We further observe that since (6) is a homogeneous function of total degree Z, namely
As an important case let us consider when all walkers have identical transition matrix W (r) ij = W ij and identical departure node i r = i ∀r = 1, . . . Z. Then with (2) (P ij (t) = P (r) ij (t)) we get for (6) the relation
with the state-probabilities given by the multinomial-coefficients
For illustration let us consider a network of two nodes i = 1, 2 (N = 2) where we have Z independent walkers 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) ij = W ij . Then the collective generating function (6) is given by
where
Hence the stateprobabilities, 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 ij with A ij = 1 if the nodes i, j are connected by an edge and A ij = 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 ij = 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 (non-zero) 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 quantitŷ R(τ 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 
with the stationary distribution W (∞) j = Kj 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 (1)
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 r(0) = (Z − 1)W (∞) j (as P ij (0) = δ ij ). Hence (24) yieldŝ
ij (τ ) 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 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
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
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 neighbours, 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
with 1/µ as the relaxation time of healing. We define particle i as infectious at time n if s (n) i > s 1 and as susceptible if s (n) i ≤ 0. In the range 0 < s (n) i < 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.
≤ 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
is the total number of particles with s (n) i ≤ 0 at time n. As long as R e > 1 the disease spreads and more and more particles get infected. The number of insusceptible particles is given as Z
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. Figs. 2, 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.
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 stay 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 Fig. 6 Cluster formation for small P = 0.2 and h = 1, initial condition of 1000 equally distributed infectious particles (red). After t = 35000 P was increased to P = 0.3 and the clusters grow.
variance
where brackets denote the average over all 100 coarse cells. Fig. 8 shows f H over time for the situation plotted in figs. 6,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 depends 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 nonmarkovian 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.