key: cord-0269530-w3dm3g89
authors: Blavatska, V.; Holovatch, Yu.
title: Spreading processes in"post-epidemic"environments
date: 2020-10-19
journal: nan
DOI: 10.1016/j.physa.2021.125980
sha: 9067ad3494e327944bf8e9718ec2848965623d98
doc_id: 269530
cord_uid: w3dm3g89

We analyze infection spreading processes in a system where only a fraction $p$ of individuals can be affected by disease, while remaining $1-p$ individuals are immune. Such a picture can emerge as a natural consequence of previously terminated epidemic process or arise in formerly vaccinated population. To this end, we apply the synchronous cellular automata algorithm studying stationary states and spatial patterning in SI, SIS and SIR models on a square lattice with the fraction $p$ of active sites. A concept of"safety patterns"of susceptible agents surrounded by immune individuals naturally arises in a proposed system, which plays an important role in the course of epidemic processes under consideration. Detailed analysis of distribution of such patterns is given, which in turn determine the fraction of infected agents in a stationary state $I^*(p)$. Estimates for the threshold values of the basic reproduction number $R_0^c$ as a function of active agents fraction $p$ are obtained as well. In particular, our results allow to predict the optimal fraction of individuals, needed to be vaccinated in advance in order to get the maximal values of unaffected agents in a course of epidemic process with a given curing rate.

The frequency and magnitude of epidemic processes is observed to increase in the last decades mainly due to globalization, growth of large centers of population in industrialized society and thus the significant growth of social contacts. Mathematical modeling in epidemiology allows to explain the dynamics patterns of epidemic process course and to propose mechanisms of disease control by providing a deeper understanding of factors controlling disease persistence [1, 2] .

The basic concepts of epidemics modeling were established in pioneering work of Kermack and McKendrick [4] . The population is assumed to be homogeneously mixed and contains three main classes -compartments -of individuals: susceptible S (which are healthy and can be infected), infected I (able to transmit the disease via air or physical contact), and removed R (either by recovery and acquiring immunity or by death). The way each particular infection is spread and cured is described by simple differential equations governed by parameters β and γ [3] . The infection rate β is dependent on the mechanism of infection transmission (air or direct physical contact) and may in general be tuned by changing the averaged frequency of contacts of individuals (e.g. it is sensitive to seasonal fluctuations [5, 6] or age of individuals [7] ). The curing rate γ is basically intrinsic for the given disease and its progression. Their ratio R 0 = β/γ called basic reproductive number is the expected number of secondary cases produced by an infectious individual in a population [5, 8] and is directly measured for a number of diseases.

Three principally different scenarios, directly observed in epidemiology, are perfectly captured within the frames of three corresponding simplest models called SI, SIS and SIR. Within the frames of SI model, individuals without proper treatment stay infected and infectious throughout their life and remain in contact with the susceptible population. This model matches the behavior of diseases like cytomegalovirus (CMV) or herpes. The SIS model describes the disease with no immunity, where the cured individuals can be infected repeatedly. This model is appropriate for diseases like the common cold (rhinoviruses). It consequently describes the emergence of so-called endemic state of the disease, where starting with R 0 > 1 the infection spreading process continues in time without termination. The SIR model describes the opposite situation when recovered person becomes immune, as in the case of influenza or measles. This mechanism leads to possibility of epidemic outbreak (at R 0 > 1). Further refinements have been done by taking into account the latency period (with its duration dependent on the type of the disease, the dose of infectious agent, the host immune response etc.) by introducing the additional compartment for so-called exposed individuals (E). This leads to so-called SEIS and SEIR models. In particular, SIR and SEIR models have been used to describe and predict epidemic processes of respiratory infections (influenza A, H1N1 etc.) [13, 9, 10, 11, 12] as well as COVID-19 [14, 15, 16, 17, 18, 19, 20, 21, 22, 23] .

In the classical models of infection spreading it is assumed, that at initial stage the population is homogeneous and individuals are equally susceptible to the disease. However, heterogeneity in susceptibility and infection is an important feature of many infectious illnesses. The impact of variable infection rates, which are heterogeneously distributed in population, on epidemic growth have been analyzed in Refs. [24, 25, 26, 27] .

Another approach of mathematical modeling of disease spreading is based on discrete lattice models, motivated by cellular automaton implementation of synchronous-update Markovian process on a lattice [28, 29] . In particular, such approach allows to take into account the local characteristics of the system, e.g. to include variable susceptibility for different individuals, and to analyze the patterns (clusters) formation in spreading process [30, 27, 31, 32] . The evolution of system is discrete in time and is governed by a specific update algorithm which takes into account the state of each agent [29, 33, 34, 13] . In particular, the special case when β = 1 − γ in SIS cellular automata model (with R 0 = (1 − γ)/γ), can be directly connected to dynamics of the so-called contact process [35, 36] , which exhibits a non-equilibrium phase transition with critical behavior of the universality class of directed percolation [37] . The critical value of creation rate λ c of this model thus is directly related with the threshold value of reproduction number R c 0 in disease spreading process. In contrast to continuous SIS model, in this case the threshold value of reproduction number which separates the non-spreading and endemic states is R c 0 = 1.64874(4) [36] , and thus is larger than the result obtained within homogeneous mixing assumption R c 0 = 1. Corresponding critical value of curing rate γ c = 0.3775. In the SIR cellular automaton, it was argued [38] that a continuous phase transition between non-spreading and spreading regimes is of the bond percolation universality class. A deep connection between the epidemy spreading and percolation theory is discussed in Ref. [39] .

Again for the special case when β = 1 − γ one obtains the critical curing rate γ c = 0.1765005(10) and corresponding threshold value for reproduction number R c 0 = 4.665706 [40] . It thus follows that the epidemic spreading for the SIS model occurs at a smaller value of the reproduction number comparing with the SIR model. This is intuitively clear since the susceptible individuals in the SIS model can be infected multiple times.

In the present work, we consider the epidemic spreading on the lattice model of population with the fraction p of individuals susceptible for a particular disease, whereas the rest 1 − p randomly chosen agents are considered to pertain immunity to it. On the one hand, such situation can describe an existence of the congenital immunity or vaccination performed against the given disease. On the other hand, it may correspond to the starting point of the "second wave" of an epidemic process, when some fraction of population already get removed (or obtained lifelong immunity) after termination of the first wave. The proposed model is related to the contact processes on heterogeneous lattices, where the spreading is restricted only to the fraction p of available sites [41] . There, the threshold value of the reproduction number R c 0 was shown to increase when p decreases. The related problems of contact processes on heterogeneous surfaces in the shape of deterministic fractals like Sierpinski gasket [42, 43, 44, 45] or two-dimensional percolation cluster [43, 46] have also been studied. We will use the cellular automata implementation of the three basic epidemic models SI, SIS and SIR on a two-dimensional regular lattice with variable ratio p (0 ≤ p ≤ 1) of sites containing agents involved in epidemic process.

The layout of the paper is as follows. In the next Section, we give the scheme of cellular automata algorithm applied for modeling the spreading processes. In Section 3, the peculiarities of cluster size distributions on disordered lattices are discussed. Section 4 contains our main results for description of SI, SIS and SIR spreading scenarios on spatially-inhomogeneous lattices. We end up with Conclusions in Section 5.

Let S(t), I(t) and R(t) be fractions at time t of healthy (susceptible), infected and recovered individuals, and S(t) + I(t) + R(t) = 1. Parameters β and γ are correspondingly infection and curing rate. We consider the population of individuals (agents), located on the sites of a regular two-dimensional square lattice of size L × L. The links between each site and its four nearest neighbours denote the possible local contacts of a corresponding individual. At time t, the kth site of the lattice is in a state σ k (t), with σ k (t) taking values from a set {0, 1, 2}, corresponding to S, I, or R respectively. Thus, the global fractions S, I, and R are given by:

where δ is the Kronecker delta and N = L × L is total number of agents. System evolution is governed by an update algorithm which changes the state of a particular individual according to the states of its neighbors. We will consider the synchronous process, where one time step implies a sweep throughout the whole lattice [28] . Time update t−1 → t consecutively considers all k = 1, . . . , N lattice sites and changes each state σ k (t − 1) → σ k (t) according to the spreading scenario rules. The updated state of each lattice site is stored in a separate array, and after the sweep is finished, the old state is replaced by the updated one. Update rules for the three spreading scenarios considered here read:

• If σ k (t − 1) = 1, then one of its randomly chosen susceptible neighbors j (with σ j (t − 1) = 0) changes its state to 1 with probability β, so that σ j (t) = 1.

(ii) SIS model

• If σ k (t − 1) = 1, then with probability γ it is cured (so that σ k (t) = 0), otherwise one of its randomly chosen susceptible neighbors j changes its state to 1 with probability 1 − γ, so that σ j (t) = 1.

(iii) SIR model

• If σ k (t− 1) = 1, than with probability γ it is recovered (so that σ k (t) = 2), otherwise one of its randomly chosen susceptible neighbors j changes its state to 1 with probability 1 − γ, so that σ j (t) = 1.

We consider a square lattice of size L = 300, the total number of individuals thus being N = 90000. The periodic boundary conditions are applied in all directions. As already mentioned above, we will analyze spreading processes on a lattice, where only a randomly selected part of sites (agents) is susceptible to disease. To this end, we consider each randomly chosen site k of the lattice to be either susceptible to disease with probability p or immune with probability 1 − p (the state of corresponding sites σ k = 0 always, and these sites are considered as non-active in updating algorithm). We perform averaging over 5000 replicas (different realizations of random distributions of active and non-active agents). The maximum number of time steps is taken t = 400 (typical times for reaching the stationary state for the cases studied below are found in general at t < 100).

The problem under consideration is closely related to the well known site percolation problem [53] . Indeed, randomly choosing the fraction p of lattice sites, one can identify the subgraphs of s linked susceptible individuals as clusters of size s. The probability P s (p) to find a cluster containing s susceptible sites is analyzed next. It is straightforward to obtain the probability to find a single susceptible site surrounded by immune. As long as these events are independent, this is obtained as a probability to find a susceptible site p times probabilities of four neighbor sites to be immune, (1 − p) 4 :

In turn

and so on. One has: P s (p) = sn s (p) with n s (p) being the number of s-cluster per one site. Exact values for n s (p) with s up to 17 were obtained in Ref. [50] for the problem of percolation on a simple square lattice. In general, one can write n s (p) = p s D s . In a low density limit of small values of p (concentration of susceptible sites well below the percolation threshold p c = 0.592746 [52] ), only clusters of relatively small size can be found in a system. At the percolation threshold, when the spanning cluster appears, one obtains the scaling behaviour

with the Fisher exponent Θ = 187/91 [53] . Note that for a finite system, the power law is restricted to small s only and naturally breaks down for cluster sizes comparable with the size of a system [54] .

To extract clusters of different sizes numerically, we apply an algorithm developed by Hoshen and Kopelman [55] . This algorithm is successfully applied in studies of percolation phenomena in disordered environments [56, 57, 58, 59, 60, 61, 62] . As the first step of the algorithm, all susceptible sites of the lattice are labeled (numbered in an increasing order). At the second step, for each of the labeled sites (say, for the site with the label n), we check whether its nearest neighbors are also susceptible. If yes, two possibilities appear. If the label of the neighbor is larger than n, we change the label of the neighbor to n. If the label of the neighbor is smaller than n, we change the label of the site n to that of the neighbor. This procedure is applied until no more changes of site labels is needed. As a result, we obtain groups of clusters of susceptible sites of different sizes, where all the sites in a given group have the same label. Our numerical data for n s (p) are given in Fig. 1 in a double logarithmic scale. On the upper panel of Fig. 1 we present results for p values increasing from p = 0.1 up to the percolation threshold. Whereas only small clusters are present in system at small p, the probability to find larger clusters increases with increasing p. The scaling behavior in vicinity of the percolation threshold p ≃ 0.6 is observed, corresponding to a straight line in a log-log plot. One notices the pronounced fluctuations in sizes of large clusters emerging in a system at the threshold point. Also, as noted above, the violation of the scaling law at large s is caused by finite system size. The lower panel of Fig. 1 shows results for p values above the percolation threshold. The probability of observing small clusters is decreasing with increasing p, whereas the fraction of sites in spanning clusters increases. Let us estimate also the "diversity" of clusters, i.e. the number of all possible clusters, that can be found on a lattice with given fraction p of susceptible sites. The lattice of size L × L contains n s (p) × L 2 different clusters of size s, and thus s n s (p) × L 2 different clusters of any size. Effective "diversity" of clusters per one site n(p) = s n s (p) can be estimated, using available analytical values for n s (p) from Ref. [50] mentioned above. These estimates are compared with our numerical data in Fig. 2 . As expected, n(p) increases with p, when p values are relatively small, and there is a large amount of small clusters. At larger p, small clusters tend to segregate into larger ones, so the number of different clusters decreases. Nice coincidence of analytical results that take into account only small clusters up to s = 17 and numerical data is due to the fact, that though at large p the prevailing number of sites is located in very large clusters, the number of such clusters is very small, so that the main contribution into n(p) is given by numerous small clusters.

The above distributions will play an important role in further analysis of epidemic processes in such system. Thus, the considered population has a "patchy" structure and contains small and large groups of susceptible individuals, contacting either between themselves or with surrounding immune agents. Related patchy system, when infection is allowed to spread only within separated populated areas was studied recently in Ref. [47] .

We apply the cellular automaton mechanism to study the epidemic process on disordered lattices, described in the previous Section. We take different values of p in a range 0.1 to 1.0 with step 0.1. and perform the averaging of all observables of interest over an ensemble of 5000 replicas of random realizations of susceptible sites configurations.

At time t = 0, a small fraction i 0 of randomly chosen susceptible individuals is supposed to get infected. The disease starts to spread from infected to susceptible agents. Taking into account the patchy structure of population under consideration, we immediately conclude, that if at t = 0 in some cluster of susceptible agents no one site gets infected, this cluster will not be touched by epidemic process in any way and remains safe till the epidemic terminates. The probability P s (p, i 0 ), that cluster of size s is not affected by infection is given by:

and thus

gives the total fraction of individuals, which are lucky to be in such safe patterns, provided that i 0 of initial p susceptible sites have been infected. The sum in (6) spans all clusters of susceptible sites. We can get analytical estimation of this value by taking available exact values for P s (p) up to s = 17 as given in Ref. [50] and substituting them into (5) . Resulting curves at several values of i 0 are presented in Fig. 3 . It is important to note the limitation of the analytical estimates, since they do not take into account the contribution of larger clusters. However, at any p there is a larger amount of small clusters than the large ones, as discussed in the previous Section. The "diversity" (number of different clusters) n(p) increases with p, until it reaches some critical value, as shown in Fig. 2 . This can help us to understand the general tendency of P safe (p, i 0 ) curves. At fixed i 0 , number of unaffected clusters grows with increasing p, since the number of different clusters grows, and thus the probability for one of them to get infected gets lower. On the other hand, above the threshold value of p the number of various clusters decreases (small clusters start to segregate into larger structures), which makes it easier to infect one of them. Let us compare the analytical estimates with our numerical results for P safe (p, i 0 ), as shown in Fig. 3 . Discrepancy at small i 0 is explained by the fact, that also the larger clusters (despite the small amount of them) which are not taken into account in analytical prediction, with higher probability remain unaffected by infection. At larger i 0 , there is perfect coincidence with analytical prediction. This tells us, that mainly the clusters of small size remain safe in such processes.

One can estimate the critical value of i 0 for each cluster size s, above which there are more infected clusters than safe at given p. This value i crit 0 can be obtained from equality condition for probabilities to have the safe and infected clusters P s (p) ( 

, and is thus independent on p and can be easily calculated analytically. The results are presented in Fig. 4 .

Another interesting question is, whether there exist a possible minimum value of i 0 at which all the existing susceptible patterns in a system get infected. This is possible in the case, when each cluster of susceptible sites contains at least one infected agent. Thus this minimum value of i 0 is strictly defined by the number of different clusters at given p, as presented in Fig. 2 . Note however, that in the case of random distribution of infected sites probability of such event is negligibly small. So in general, introducing any fraction of infected individuals i 0 < 1 into the patchy lattice one cannot end up with completely infected system and some fraction of individuals P safe (p, i 0 ) > 0 will remain locked in the safe 

The above results concern any spreading process on a disordered lattice. Let us now discuss the outcomes of specific spreading scenarios.

Within this model, the initially infected agents are spreading infection by contact with their neighbors within the clusters, to which they belong, until the equilibrium state is reached. Snapshots of configurations of infected and susceptible agents at the initial time and after reaching the equilibrium are shown in Fig. 5 . The initial concentration of infected sites of the whole lattice I(0) is related to i 0 via: I(0) = p × i 0 . We took the value for infection ratio β = 0.5 to make the evolution process quite fast, though the equilibrium values of ratios of infected and susceptible individuals I * (p) and S * (p) do not depend on β.

Note that the estimates of I * (p, i 0 ) and S * (p, i 0 ) can be obtained straightforwardly by making use of conclusions from the previous Subsection. Indeed, the total fraction P safe (p, i 0 ) of safe patterns, surrounded by immune individuals, gives the fraction of susceptible individuals in equilibrium: S * (p, i 0 ) = P safe (p, i 0 ) and thus we get:

Again, recalling the definition (6), making use of available exact values for P s (p) up to s = 17 as given in Ref. [50] and substituting them into (7), we obtain analytical estimates for I * (p, i 0 ), presented in Fig. 6 at several values for i 0 . These results should be compared with our numerical data for I * (p, i 0 ). Again, at larger values of i 0 , where mainly the smaller clusters remain untouched by infection, there is perfect coincidence with the analytical predictions.

Within this model, the initially infected individuals are spreading infection by contact with their neighbors with infection rate β or get cured with curing rate γ = 1−β until the equilibrium state is reached. Snapshots of configurations of infected and susceptible agents at the initial time and after reaching the equilibrium are shown in Fig. 7 for the same initial conditions as those of Fig.  5 for the SI scenario. As expected, one notes essentially smaller number of infected sites in the equilibrium I * (p). An essential feature here is that any cluster of size s, which gets infected initially, can become completely "cured" in a course of the SIS process and once it happens, all sites in this cluster remain susceptible. Thus at the endemic state of SIS process, when the fraction of infected sites reaches its stationary value I * (p, i 0 ), one observes in general the larger amount of patterns of susceptible agents, than in the case of SI process. Indeed, in addition to the set of "safe" patterns P safe (p, i 0 ), which were not affected by the infection from the very beginning (given by Eq. (6)), also the set of "cured" patterns P cured (p, i 0 ) contributes.

In the case of disease with repeated infection the only condition for a given individual to stay healthy (once it is cured) is the absence of infected individuals in its neighborhood. Thus, the probability for cluster of size s to become completely "cured" is proportional to q s 0 , where q 0 is the probability, that in the stationary state no one of four nearest neighbors of a given susceptible site is infected. The value of q 0 is easily estimated numerically as the averaged number of susceptible nearest neighbors of a given susceptible agent; corresponding numbers are presented in Fig. 8 at i 0 = 0.1. As expected, the higher is the curing rate γ, the larger is the probability for a given site to be surrounded by only non-infectious neighbors. Note that these values are dependent also on the initial concentration i 0 slightly decreasing with an increase of i 0 .

Thus, the total fraction of susceptible clusters which can be found in a system after the equilibration of SIS process, is given by:

The second term in this equation contains the fraction of individuals in clusters, which were originally infected but completely cured in SIS process. In Fig. 9 we give the estimates based on our numerical data. In a course of time evolution, the ratio of infected individuals I(t) reaches the equilibrium value I * (p, γ). Corresponding numerical data are given in Fig.  10 . As expected, for fixed γ the final equilibrium fraction of infected individuals decreases with the decrease of p. The data shown in Fig. 10 demonstrate a continuous transition from the state when I * ≃ 0 (no endemic occurs) and I * > 0 at certain critical value of curing rate γ c . To estimate γ c at various p, we apply the fitting to the power form

where β ′ is the critical exponent for the order parameter with a mean-field value β ′ = 0.5 [48] . In turn, we can estimate the threshold values of reproduction number according to R c 0 = (1 − γ c )/γ c . Results of fitting of our data are given in Fig. 11 .

At p = 1, our estimates are in agreement with R c 0 = 1.646 obtained in [31] as well as the estimates of λ c = 1.64872(3) of contact processes [36] . Note that at p < 1 our results could have been compared with those obtained for contact processes on disordered lattices where only fraction p of sites are available for spreading [41] . However, two different types of ensemble averaging are Figure 11 : The threshold value of the reproductive number R c 0 (p) as a function of p for two initial concentrations of infected sites i 0 performed in our work and in Ref. [41] . Whereas we put the fraction i 0 of initially infected sites randomly and apply averaging over all clusters of initially susceptible sites in the systems, in the mentioned study only averaging over the spanning percolation cluster was performed and the process started with single active site on this cluster. As a result, the values for λ c at various p obtained in that study are larger than our corresponding R c 0 (p). Note also, that at p < 1 the results for R c 0 are sensitive to the value of initial concentration i 0 . Plots of Fig. 10 reflect behavior of a finite-size system and therefore manifest both universal and non-universal features. Indeed, for large values of p > p c ≃ 0.592746 [52] there is a non-zero probability to find a spanning cluster of active sites. Therefore, an infinite system possesses at γ = γ c (p) a well-defined transition to an endemic state I * (p, γ) = 0 for any non-zero value of i 0 . This is reflected by upper curves in Fig. 10 . On opposite, for p < p c probability to have a macroscopic equilibrium value I * depends of an initial value of i 0 , as shown in Fig.  11 .

Within this model, the initially infected individuals are spreading infection by contact with their neighbors with infection rate β or get removed with curing rate γ = 1 − β until the equilibrium state is reached. Typical configurations of infected and susceptible agents at the beginning and after reaching the equilibrium are shown in Fig. 12 for the same initial conditions as those of Figs. 5, 7 for SI and SIS scenarios, correspondingly.

An obvious difference is that since the repeatedly infection is not allowed in course of SIR process, it terminates with some equilibrium fraction R * (p, i 0 , γ) of individuals becoming removed (or with the life-long immunity). The rest of individuals remain susceptible and thus can be affected in the case of the subsequent wave of epidemic. In the case of SIR scenario, the fraction of individuals in susceptible patterns is given by:

Again we expect, that the total fraction P SIR (p, i 0 , γ) of susceptible individuals after the termination of SIR process is larger than the fraction in nonaffected patterns P safe (p, i 0 ). Indeed, besides P safe (p, i 0 ) it includes also contributions from those infected clusters, where the infection spreading has been terminated keeping some agents in the S state. This is further demonstrated by Fig. 13 where we present our simulation data for P SIR (p, i 0 , γ) in comparison with P safe (p, i 0 ).

Note that for each value of the curing rate γ there exists a corresponding value p max (γ) at which P SIR (p, i 0 , γ) reaches the maximum. After performing a more refined simulation with p changed in a range [0 . . . 1] with a step 0.001 we obtain the data for p max (γ) as a function of γ, see Fig. 14. Recalling that 1 − p can be treated as the fraction of individuals that get vaccinated before the spreading process starts, these results can be used for prediction of efficiency of vaccination for disease with given curing rate γ. The larger is the value of γ (and thus the weaker is the infection), the larger is the critical value p max (γ) and reversely the smaller is the fraction of agents 1 − p max (γ) who need to be vaccinated in order to have maximal values of unaffected individuals in the equilibrium state. For the case of diseases with sufficiently large values of γ, the infection is so weak that practically no vaccination is needed: p max (γ) tends to 1, therefore the number of vaccinated persons, needed to obtain the maximum of unaffected individuals, tends to zero.

Infection spreading throughout the lattice results in an increase of clusters of affected sites (infected and eventually removed). The curing rate critical value γ c separates a state with no epidemic (corresponding fraction of removed sites R * ∼ 0) and a post-epidemic state with considerably large value of R * . The value of γ c can be estimated by analyzing the removed cluster distributions from the condition of emergence of the spanning percolation cluster [40, 51] . In the case when spreading starts with a single infected site, the process is closely related to standard percolation growth, when an active (infective) site spreads to susceptible nearest neighbor with the given bond probability, corresponding in this case to the infection rate β. This approach however cannot be applied straightforwardly in our problem, where we put the fraction i 0 of initially infected sites randomly on various clusters of susceptible sites available at p < 1.

In Fig. 15 we present our results obtained for R * (p, γ) as function of γ for various p. The very rough estimates for γ c can be obtained from the condition of pronounced decrease of R * (p, γ) after the plateau regime. Let us remind, that at p = 1 the critical value γ c = 0.1765005(10) [40] , and it is obviously decreasing with the decrease of p. As expected, the smaller is the fraction of susceptible agents in population, the smaller is the critical value of curing rate (and thus the larger is the corresponding value of the basic reproduction number R c 0 ), above which the infection spreading process terminates without causing a global epidemic process.

In this paper, we have presented quantitative analysis of a spreading process in a non-homogeneous environment. To this end, we have considered three archetypal spreading models, SI, SIS and SIR on structurally-disordered lattices, when only a part of lattice sites (agents) is active and takes part in spreading phenomenon. One of possible interpretations of such a problem statement is infection spreading in population where only a part of individuals is susceptible to disease. The rest is non-sensitive to infection, as may happen after the first wave of epidemic in over and a part of population gets immune. Another interpretation may be that some randomly chosen part of agents gets immune due to the formerly performed vaccination of a part of population. In our analysis, we have considered the case, when the probability of arbitrary chosen lattice site to be susceptible is uniquely determined by their concentration p. In turn, this leads to emergence of different size clusters of active susceptible sites. As it becomes evident from the analysis, spreading process in such "patchy" population is characterized by a number of distinctive and sometimes unexpected features. In particular, absence of connetivity between finite-size clusters of susceptible sites causes existence of the so-called "safety patterns": if no infected sites fall into such pattern, its state remains unchanged untill the stationary state is reached. Detailed analysis of distribution of the "safety patterns" en-abled us to determine the fraction of infected agents in a stationary state I * (p). Subsequently, this lead to the estimates for the threshold value of the basic reproduction number R c 0 as a function of active agents fraction p for different spreading scenarios.

Within the homogeneous mixing hypothesis [3] , consideration of a part of population as being susceptible is trivial and leads to renormalization of an active fraction. However, this is not the case for the lattice model considered here. In this respect it is instructive to compare processes of spreading and ordering on an inhomogeneous lattice. To give an example, inhomogeneities in lattice structure caused by dilution of a magnet by a non-magnetic component, may cause severe impact on a magnetic phase transition [48, 66] . However, for a short-range interaction, the spontaneous magnetization can be non-zero only if the concentration of magnetic sites is above the percolation threshold, p > p perc . No ordering is expected for p < p perc : indeed, in this region only the finite-size clusters of magnetic sites exist, and their contribution to overall magnetization vanishes in the thermodynamic limit N → ∞. As we have shown above, this is not the case for the spreading phenomenon on disordered lattices: there, the spreading is non-trivial for any concentration of active sites p, also at p < p perc . The reason is that although in this region the contribution from each finite-size cluster vanishes, their steady-state value may be non-zero also in the limit N → ∞, depending on the basic reproduction number R c 0 . It is obvious that results obtained in our analysis crucially depend on the distribution of susceptible sites. Here, we have considered the simplest case, when correlations in distribution are absent. It may be interesting to re-consider the problem for the case when such correlations are present, taking them from available real-world data. Returning back to the above discussed comparison of spreading and ordering on inhomogeneous lattices, this example may find its counterpart in and impact of long-range-correlated impurities on magnetic phase transitions [67] and scaling [68] . Another work is progress is to apply similar methodology to analyze spreading processes on complex networks [69] .

Mathematical Epidemiology

The Mathematical Theory of Infectious Diseases (Macmillan

Mathematical models in medicine

Advanced Computing and Communication Technologies

Infectious Disease modeling

Infectious Disease modeling

The Critical Point

Order, Disorder and Criticality. Advanced Problems of Phase Transition Theory, vols

Introduction to Percolation Theory

Phase Transitions and Critical Phenomena

This work was supported in part by National Research Foundation of Ukraine, project "Science for Safety of Human and Society" No. 2020.01/0338 (VB) and by and by the National Academy of Sciences of Ukraine, project KP-KBK 6541230 (YuH).. It is our pleasure to acknowledge useful discussions with Jaroslav Ilnytskyi and Mar'jana Krasnytska.