key: cord-0335958-psx8rxxx
authors: Ferraz, Carlos Handrey Araujo
title: A non-absorbing SIR stochastic lattice gas model on hybrid lattices
date: 2021-10-03
journal: nan
DOI: 10.1016/j.physleta.2021.127871
sha: 1521a967d7b2782cfd4c8127d6c7b5f9b1402a23
doc_id: 335958
cord_uid: psx8rxxx

In this paper, we perform Monte Carlo calculations to study the critical behavior of the spread of infectious diseases through a novel approach to the SIR epidemiological model. A stochastic lattice gas version of the model was applied on hybrid lattices which, in turn, are generated from typical square lattices when inserting a connection probability $p$ that a given lattice site has both first- and second-nearest neighbor interactions. By combining percolation theory and finite-size scaling analysis, we estimate both the critical threshold and leading critical exponent ratios of the non-absorbing SIR model in different cases of hybrid lattices. An analysis of the average size of the percolating cluster and the size distribution of non-percolating clusters of recovered individuals was carried out to determine the universality class of the model.

In the late 1920s, Kermack and McKendrick [1] proposed a system of ordinary differential equations for determining the temporal evolution of the population of individuals who interacted with each other when exposed to an infectious disease. These individuals were initially divided into three classes, namely, susceptible (S), infected (I), and recovered (R) individuals. Taking the acronym for the classes, the model became known as the SIR model. Over the years, variations of this model (e.g., SIRS, SIS, SEIR models) have been proposed considering the different classes of individuals that constitute a given population. The relevance of these approximations cannot be underestimated because they have been successful in modeling the most varied forms of epidemics, including cholera [2, 3] , measles [4, 5] , rubella [6, 7] , hepatitis [8, 9] , influenza [10] , AIDS [11, 12] , COVID-19 [13] [14] [15] , among many others. Epidemiological models encompass both features of collective dynamics [16, 17] and complex systems [18] ; hence, they are important tools for obtaining information regarding the rate of disease spread and for testing protocols adopted by public bodies to contain or mitigate such diseases.

The SIR model is one of the simplest epidemiological models, in which the infected individuals can be removed from the dynamics either by permanent immunity or death. This peculiarity makes the model suitable for simulating the epidemic outbreaks of influenza, SARS, AIDS, etc. The SIR model is in the same universality class as that of dynamic percolation [19] [20] [21] [22] [23] . This study considers a stochastic lattice-gas version of the SIR model with asynchronous site updates. Similar models have been applied to other population dynamics [24, 25] . For both synchronous and asynchronous versions, a phase transition occurs when the model's control parameters are varied. This transition is found to be of second order between two distinct regimes: one, in which the population remains susceptible (inactive or endemic regime), and other, in which the disease spreads throughout the network (active or epidemic regime), where a significant portion of the population becomes infected or eventually recovers (immune or dead). At the transition point, the system becomes critical and corresponds to an epidemic outbreak threshold. It is worth remarking that the connection between the stochastic and deterministic descriptions (via coupled differential equations) can be achieved through meanfield approximation [26, 27] , which results in the socalled Langevin equations associated with the Fokker-Planck equation that describes the temporal evolution of the probability density of the system states.

The lattice gas method [28] used herein has been widely used to address problems involving fluid dynamics [29, 30] , general diffusion processes [31] , damage spreading [32, 33] , and transport phenomena [34] . This method basically consists of discretization both in time and space, where sites in the network are considered as the most likely particle locations. In molecular dynamics, for example, we only have a temporal discretization in the particle state description. The lattice gas method is especially useful when the dynamics of a particle system needs to be described without considering the detailed microscopic aspects of that system. These aspects are often irrelevant to the description of the general behavior of the system. In this study, we will employ this method to describe the spread of the cloud of recovered individuals over time as the epidemic evolves.

In particular, an analysis of the average size of the percolating cluster and the size distribution of nonpercolating clusters of recovered individuals will be performed using the Newman-Ziff algorithm [35] to determine the critical exponent ratios of the model. The SIR stochastic lattice gas model is an absorbing-like model because its active phase is characterized by an infinite number of absorbing configurations, in which the final system state only comprises recovered individuals. However, in this work, the simulations are halted as long as the existence of the percolating cluster (spanning cluster) in the system is verified, thereby setting the nonabsorbing state of the model. This procedure speeds up the analysis of the simulation data without compromising the guarantee of reaching the asymptotic limit of the system. Only a single percolation cluster is generated in the critical regime of the system. By exploiting the analysis of these clusters, we will determine the universality class of the model. The network topology represents an important feature of the system and directly affects the dynamics of the involved processes [32, 33, 36] ; hence, we will study a hybrid network topology in this work.

Hybrid lattices are formed from regular square lattices (N = L×L) by inserting a probability p that a given site has both first-and second-nearest neighbor interactions. Such networks would simulate a more realistic population mainly formed by two types of individuals, i.e., type I with low connectivity and type II with high connectivity. The special cases with p = 0 and p = 1 are also treated here and correspond to pure lattices having all nodes with the connectivity of first neighbors and first and second neighbors, respectively.

This study mainly aims to understand how extended connectivity effects can affect the critical behavior of the non-absorbing SIR model. Since the connection probability p can be understood as a quenched topological disorder, we want to determine if this kind of disorder is relevant to changing the universality class of the model. There are some criteria devised to try to predict if the quenched disorder can change the critical exponents of a given model, such as the Harris criterion [37] and its refinement, the Harris-Barghathi-Vojta criterion (HBV) [38] . The Harris criterion states that a second-order transition in a d-dimensional system, with original correlation length exponent ν, is stable against the quenched spatial disorder if ν > 2/d. Whereas the HBV criterion states that quenched topological disorder is irrelevant with respect to the phase transition stability if the system satisfies ν > 1/a, where a is the disorder decay exponent that measures how fast coordination number fluctuations decay with increasing system length scale. Nevertheless, such criteria are known to fail [39] .

The contents of the article are organized as follows. In section 2, we outline the SIR stochastic lattice gas model. In section 3, we describe details of the Monte Carlo (MC) simulation background and lattice generation. In section 4, we present and discuss the results. Finally, in section 5, we make the conclusions.

The SIR stochastic lattice gas model is defined on a given lattice of N sites in which each site can be occupied by just one individual who can be either a susceptible (state S ), an infected (state I), or a recovered (immunized/dead) individual (state R). The dynamics consists of two subprocesses, namely, an auto-catalytic one, I + S → I + I; and a spontaneous one, I → R. At each time step a site is randomly chosen and then a set of dynamic rules are taken in the following way: i) If the site is in the state S and there is at least one neighboring site in the sate I then the site becomes I with probability proportional to a parameter µ and the number z of neighboring sites, i.e., µm/z, where m is the number of neighboring I sites. ii) If the site is in state I it becomes R spontaneously with probability λ. iii) If the site is R it remains unchanged. At each site i of a 2D lattice we assign a stochastic variable σ i that takes the values 0, 1 or 2, according to whether the site is in the state S , I, or R, respectively. Since the transitions between states in this model are non-equilibrium ones, the allowed transitions of the state i of a site are cyclic, this is, 0 → 1 → 2. The corresponding transition rate is represented by w i (σ) and describes the transition σ → σ ′ in which the whole microscopic configuration (microstate) 

. . , σ N ) differs from σ only by the state of the i-th site. It is given by

where the summation runs over the nearest neighbors of site i and δ(x, y) denotes the Kronecker delta. The parameters µ and λ are related to the subprocesses above described, and are chosen such that µ + λ = 1. The system evolves in time according to a master equation for the probability distribution P(σ, t) described by

where the microstateσ is obtained from σ by an anticyclic permutation of the state of the site i (2 → 1 → 0).

We can implement an asynchronous, non-absorbing SIR model on computer by following the kinetic Monte Carlo rules below:

1. First, we start with a single central infected site (seed) and the remaining ones being all susceptible on a two-dimensional lattice in which each individual of the population is attached to its respective lattice site. In order to speed up the simulation we create two lists that are updated at each algorithm step: a list of infected individuals (infected list) and a list of recovered individuals (recovered list), which begins empty. 2. Next, we update the system state by randomly choosing an available infected site from the infected list and proceed as follows: (a) Generate a random number x in the interval (0, 1). If x ≤ λ, the infected site is removed from the infected list and placed in the recovered list; (b) Otherwise (if x > λ), pick randomly one nearest neighbor of the infected site and make it also infected provided that it is susceptible, adding it to infected list. 3. Repeat asynchronously the step (2) several times until either there is no infected sites (endemic phase) or there is a percolating cluster of recovered sites (non-absorbing epidemic phase).

One could determine the MC time t by incrementing t by δt = 1/n I , where n I is the current number of infected sites, each time an infected site is pick from the list. However, we do not keep track of time here since we are more interested in static quantities such as the fraction and the mean cluster size of R sites.

Remarkably, it was shown that the SIR model on square lattices belongs to the same universality class then dynamic percolation (DP). This allows us to investigate the phase transition which takes place in the present non-absorbing SIR model by making use of the percolation theory. Thus we can define the epidemic phase of the model when it is formed a percolating cluster of recovery sites in the system and the endemic phase when it is not. Such a graphical analogy has been used for other compartmental models as well. Following the classical percolation theory it is important first to determine the cluster distribution of recovery sites, i.e., the number of clusters with s recovery sites n p (s). That can be accomplished by using the Newmann-Ziff algorithm, which possess also the built-in feature of identifying whether a percolating cluster was formed or not, depending on the considered λ value. Notice that for systems with non-periodic boundaries like the ones concerned here, the percolating cluster is actually a spanning cluster [35, 40] . From the cluster size distribution, we have the fraction of recovery sites in the finite (non-percolating) cluster with s size

where n R is the total number of recovery sites and n p (s) is the number of clusters with s recovery sites. Furthermore, the fraction of recovery sites in the percolating cluster P ∞ can be obtained by

such that the above summation excludes the percolation cluster. Now we can define the order parameter from Eq. (4) as

where < x > means an average taken over different dynamic realizations. The epidemic phase of the model is reached when P 0, that is, when the percolating cluster density is non-zero; while the endemic phase happens when P = 0. Other important quantities are the mean cluster size

which plays the rule of the susceptibility in classical percolation theory [41] [42] [43] when taking the average over different runs, i.e.,

the overall mean cluster size

and the mean quadratic cluster size

where the primed summations above also include the percolating cluster. It is worth remarking that the last two quantities S ′ and M ′ only make sense for finite lattice as in the asymptotic limit (N → ∞), the percolating cluster size diverges. At criticality, the cluster size distribution should obey a power-law scaling [22, 42, 44] as

where λ c is the epidemic threshold and F is a scaling function. Similarly, the remaining quantities also obey scaling relations in accordance to classical percolation theory given by

The reciprocal correlation-length exponents 1/ν can be obtained by calculating the modulus of the logarithmic derivative of P at the critical threshold point λ C

where the derivative of the function f ≡ ln(P) was evaluated numerically by using a finite central difference scheme in the form

which has an truncation error of the order of O(h) 2 . The error function δ f ′ of d f /dλ was obtained via error propagation from the uncertainties in the values of f (δ f ), being expressed by

In our computations h was taken equal to 2.0 × 10 −3 . Close to λ c , the quantity φ obeys a power-law scaling as

whereφ is a scaling function, b is correlation amplitude and ω is the non-universal correction-to-scaling exponent. We have inserted a correction-to-scaling term in Eq. (18) to improve the fit quality of the data such that the values of b and ω are chosen in order to minimize the reduced chi-squared (χ 2 ) of the fits. As we will see, however, only for the case p = 1, the fit quality is slightly improved by this correction term. The optimal values of b and ω are given in the figure captions for each p case. In addition, we can define a universal quantity U in which the scaling dependencies cancel out by combining Eqs. (11) , (13) , and (14) in the following way [22] 

being analogous to the Binder cumulant for ferromagnetic spin model [45, 46] , and obeying also a scaling relation

The crossing point of the U curves for different lattice sizes allow us to estimate the epidemic threshold λ c , whereas a finite-size scaling analysis of the observables P, χ and φ by using Eqs. (11) , (12) and (18) yields the according critical exponent ratios β/ν, γ/ν and 1/ν. The leading critical exponents β, γ and ν define the universality class of the system. The hybrid lattices used in the present study were constructed starting from regular square lattices with free boundary conditions. First, we begin from a regular square lattice consisting of nodes linked to their four first nearest neighbors by both outgoing and incoming links. Then, with probability p, we connect a chosen site also to their second nearest neighbors. After repeating this procedure for every site, a new lattice is constructed with a density p of nodes with both first-and second-nearest neighbor connections. Such networks would mimic a hypothetical population formed basically by two types of individuals: type I with low connectivity, and type II with high connectivity. The special cases with p = 0 and p = 1 correspond to pure lattices with all nodes having the connectivity of first neighbors and first and second neighbors, respectively. Figs. 1(a) and 1(b) display typical spanning clusters along with few smaller clusters formed close to λ c for the cases with p = 0 and p = 1/2, respectively. These clusters arise from the SIR model dynamics. The spanning cluster is formed when the cloud of recovered individuals reaches any two opposing edges of the lattice, in other words, when such a cloud spans the lattice from one side to the other. We grew more than 10 5 spanning clusters for every considered λ value to take reliable averages of the above quantities. On average it took about 26 ms to grow a single spanning cluster like those shown in Fig. 1 on an 3.70 GHz Intel Xeon workstation.

In this section, we show our numerical results of the non-absorbing SIR model coupled to hybrid lattices. In order to determine both the critical region and the order of the phase transition in this model on hybrid lattices, we calculated the order parameter P, Binder cumulant U, and susceptibility χ for each p case in a wide range of the parameter λ. These quantities were averaged over at least 10 5 different dynamic realizations of the SIR model. Furthermore, for the case p = 1/2 we consider also different lattice configurations upon taking the averages. We deal with several lattice sizes, ranging from N = 400 up to N = 48400. Let us first discuss the case p = 0.

In Figs. 2(a), 2(b) and 2(c) are shown the order parameter, Binder cumulant, and susceptibility as a function of the recovery rate λ for the case p = 0, respectively. As one can see from Fig. 2(a) , a typical secondorder phase transition takes place for the case p = 0. As we shall see, the same conclusion can be drawn for the remaining cases, where a typical sigmoid-shaped curve also occurs. From Binder cumulant crossings, we can estimate the corresponding epidemic thresholds for each case. In the inset of Fig. 2(b) is shown a refinement of the calculations for U inside the critical region. The critical thresholds were estimated with five significant figures. For the case p = 0, we obtained λ c = 0.176 (6) .

As expected, it is approximately equal to that observed in the absorbing version of the SIR model. By taking the slope of the log-log plot of the quantity φ versus the linear size of the system L, we can get an estimate for 1/ν. Fig. 2(d) shows the best linear fit to Eq. (18) for the case p = 0. We obtained 1/ν = 0.719 ± 0.011. Er- Eqs. (11) and (12) yield, respectively, the critical exponent ratios β/ν and γ/ν. Figs. 2(e) and 2(e) show the log-log plot of P and χ (both calculated at λ c ) against L, respectively. The red straight lines in those figures are the best linear fit to Eqs. (11) and (12), respectively. We obtained β/ν = 0.109 ± 0.003 and γ/ν = 1.778 ± 0.005. These values are also in very good agreement with the exact critical exponent ratios of 2D dynamical percola-tion, namely, β/ν = 5/48 and γ/ν = 43/24. These estimates of the critical exponents ratios and critical threshold λ c for case p = 0 are summarized and compared with the corresponding exact values form 2D dynamics percolation in Table 1 .

An analogous analysis can be done for the cases p = 1/2 and p = 1.0. Figs. 3(a) and 3(a) display the order parameter P as a function of the recovery rate λ for the cases p = 1/2 and p = 1, respectively. As already remarked, both systems undergo also a secondorder transition with their respective P curves exhibiting a typical sigmoidal shape. While Figs. 3(a) and 3(a) show the Binder cumulant for p = 1/2 and p = 1, respectively. From these figures, one see that the crossing points are located at λ c = 0.228(4) for p = 1/2 and λ c = 0.275(0) for p = 1. Likewise, we took the slope of the log-log of φ defined by Eq. (15) versus L to estimate the exponent ratio 1/ν for p = 1/2 and p = 1 cases. Figs. 3(d) and 4(d) display the linear curve fitting to Eq. (18) for p = 1/2 and p = 1 cases, respectively. Error bars were calculated by using Eq. (17) . We got 1/ν = 0.728 ± 0.013 for p = 1/2 (with a least reduced χ 2 = 1.14 and a goodness-of-fit probability Q = 32.2%) and 1/ν = 0.716 ± 0.014 for p = 1 (with a least reduced χ 2 = 0.68 and a goodness-of-fit probability Q = 75.5%). These estimates deviate, respectively, only 3% and 4.5% from the exact value of 1/ν. Such results strongly suggest that both cases are in the same universality class of 2D dynamic percolation. Moreover, these results are fairly close to each other, as one falls within only one standard deviation from the other. Similarly, we took the slopes of the log-log plots of P and χ versus L for both cases. Again, the red straight lines are linear regressions to the corresponding data. From Figs. 3(e) and 4(e), we obtained β/ν = 0.108 ± 0.003 for p = 1/2 and β/ν = 0.112 ± 0.003 for p = 1, respectively. While from Figs. 3(f) and 4(f), we got γ/ν = 1.79 ± 0.01 for p = 1/2 and γ/ν = 1.776 ± 0.009 for p = 1, respectively. As can be seen, these estimates are rather close to each other. These results clearly suggest that the nonabsorbing SIR model on both hybrid lattices (p = 1/2) and square lattices with first-and second-nearest neighbor interactions (case p = 1) belong also to the same universality class as that of 2D dynamic percolation. Nevertheless, the critical threshold in the non-absorbing SIR model increases with increasing p-connection disorder. The estimates of the critical exponent ratios and the critical threshold λ c for all treated p cases are summarized in Table 1 .

We performed Monte Carlo simulations of the nonabsorbing SIR stochastic lattice gas model on hybrid lattices to study the critical behavior presented by these systems. Both the critical threshold and leading critical exponent ratios were estimated for different cases. Our numerical analysis has revealed that the quenched p-connection disorder is irrelevant to changing the critical exponents of the model, irrespective of the considered p value, strongly suggesting that the present model belongs to the same universality class as that of twodimensional dynamical percolation. However, it was found that the critical threshold in the non-absorbing SIR model increases with increasing p-connection disorder.

This study has wide applications to many issues, including not only the spread of infectious diseases, but also in general diffusion processes, damage propagation in random networks, and performance optimization in multi-core architectures. Finally, we expect that the results presented in this paper can be helpful to understand how topological disorders can affect the critical properties of other related complex systems.

Lattice Boltzmann Modeling

Simulation of Diffusion in Lattice Gases and Related Kinetic Phenomena

Introduction to Percolation Theory