Asynchronous SIR model on Two-Dimensional Quasiperiodic Lattices

TL;DR

This paper investigates the critical behavior of the asynchronous SIR epidemic model on quasiperiodic lattices using a novel algorithm, revealing that it belongs to the 2D percolation universality class and quasiperiodic order is irrelevant.

Contribution

It introduces a Newman-Ziff based algorithm for simulating SIR on non-periodic lattices and analyzes its critical properties.

Findings

System exhibits 2D percolation universality class.

Quasiperiodic order does not affect critical behavior.

Identifies epidemic threshold and cluster properties.

Abstract

We considered the Asynchronous SIR (susceptible-infected-removed) model on Penrose and Ammann-Beenker quasiperiodic lattices, and obtained its critical behavior by using Newman-Ziff algorithm to track cluster propagation by making a tree structure of clusters grown at the dynamics, allowing to simulate SIR model on non-periodic lattices and measure any observable related to percolation. We numerically calculated the order parameter, defined in a geographical fashion by distinguish between an epidemic state, characterized by a spanning cluster formed by the removed nodes and the endemic state, where there is no spanning cluster. We obtained the averaged mean cluster size which plays the role of a susceptibility, and a cumulant ratio defined for percolation to estimate the epidemic threshold. Our numerical results suggest that the system falls into two-dimensional dynamic percolation universality class and the quasiperiodic order is irrelevant, in according to results for classical percolation.

Equations

$$S+I\stackrel{{\mu }_i}{\to }2I$$

$$I\stackrel{{\mu }_c}{\to }R.$$

$$N_r=\sum _{i=1}^n\left|{\psi }_i\right|.$$

$$P_s=s{n}_{\mathrm{cluster}}(s).$$

$$p=P_{\infty }+\sum _{s}s{n}_{\mathrm{cluster}}(s),$$

$$n_{\mathrm{cluster}}(s,{\lambda }_c)\propto s^{-\tau },$$

$$n_{\mathrm{cluster}}(s)=s^{-\tau }\mathrm{F}\left[\left(\lambda -{\lambda }_c\right)s^{\sigma }\right],$$

$$P=⟨P_{\infty }⟩$$

$$S=\frac{1}{N_r}\sum _{s}s{P}_s=\frac{1}{N_r}\sum _s{s}^2{n}_{\mathrm{cluster}}(s),$$

$$\chi =⟨S⟩.$$

$$S^{\prime }=\frac{1}{N_r}\left(s_{\mathrm{perc}}^2+\sum _s{s}^2{n}_{\mathrm{cluster}}(s)\right),$$

$$M^{\prime }=\frac{1}{N_r}\left(s_{\mathrm{perc}}^3+\sum _s{s}^3{n}_{\mathrm{cluster}}(s)\right).$$

$$P$$

$$\chi$$

$$⟨S^{\prime }⟩$$

$$⟨M^{\prime }⟩$$

$$U=P\frac{⟨M^{\prime }⟩}{{⟨S^{\prime }⟩}^2}$$

$$U\approx f_U\left(L^{1/\nu }\left|\lambda -{\lambda }_c\right|\right).$$

$${\lambda }_{\max }={\lambda }_c+b{L}^{1/\nu },$$

Full text

I Introduction

Asynchronous SIR (susceptible-infected-removed) modelBernoulli (1760); Ohtsuki and Keyes (1986); Renshaw (1991); Keeling and Rohani (2007); Henkel et al. (2008); Tomé and Ziff (2010); de Souza et al. (2011); Tomé and de Oliveira (2011, 2015); Ruziska et al. (2017); Pastor-Satorras et al. (2015) is a stochastic model meant to simulate epidemic outbreaks. The essential ingredients to accomplish that are all present in the model: first, it has the needed spatial structure, composed of a lattice with bonds which can represent the contacts between individuals, and second, it has a dynamics composed of stoichiometric equations describing temporal evolution of a disease. Stochasticity itself is the third ingredient, which introduces fluctuations. Fluctuations are an important distinction between stochastic models and deterministic ones to simulate a population of individuals where a disease can spreadRenshaw (1991); Marro and Dickman (1999); Keeling and Rohani (2007); Tomé and de Oliveira (2015).

SIR model is one of the simplest models of epidemic spreading by contact, where infected subjects can recover and gain permanent immunity (or die). Another simple model of disease spreading is the SIS (susceptible-infected-susceptible) model (and its variant, the Contact Process model)Pastor-Satorras et al. (2015); Boguná et al. (2013); Harris (1974a); Castellano and Pastor-Satorras (2006); Ferreira et al. (2016); Jensen (1993); Lübeck and Willmann (2002); Lübeck and Heger (2003); Janssen et al. (2007); Henkel et al. (2008); de Oliveira et al. (2008); da Silva et al. (2013); Almeida et al. (2016); Macedo-Filho et al. (2018); Mota et al. (2018); Schrauth et al. (2018). The main difference between SIR and SIS models is the permanent immunity (or removal) present in SIR model. This important distinction is the cause of these two models fall into distinct universality classes: dynamical percolationOhtsuki and Keyes (1986); Pastor-Satorras et al. (2015); Tomé and Ziff (2010); de Souza et al. (2011) and directed percolationJanssen (1981); Grassberger (1982) for SIR and SIS models, respectively. SIR model is appropriate to simulate the spreading diseases like Ebola, SARS, AIDS, Chickenpox, etc., where death/removal or permanent immunity are one of the states an individual can assume.

It is known that SIR model can display a continuous phase transition between an endemic state to an epidemic stateTomé and Ziff (2010); de Souza et al. (2011), where epidemic and endemic phases are defined in a geographical fashion. The epidemic phase is characterized by the presence of an “infinite component”, identified by a percolating cluster formed by the activated sites (eventually removed ones at the end of dynamics). The critical threshold between the epidemic (percolating) and endemic (non-percolating) phases is controlled by the infection and recover rates. The infection rate $\mu_{i}$ is the probability with an infected particle can successfully infect one of its neighbors in a Poisson process and the recover rate $\mu_{c}$ is the probability of an infected particle can gain permanent immunity (or “die”), being excluded from the dynamics. Every path, starting from a “patient zero” (one randomly chosen individual to be the first infected one of the epidemic outbreak), subjecting to both infection and recover rates with permanent immunity, always evolve to an absorbing state where there are not any infected particles. Presence of absorbing states negates the detailed balance condition, therefore SIR model is a non equilibrium one.

In this present paper, we considered SIR model on Penrose and Ammann-Beenker quasiperiodic lattices. Our objectives are two-fold:

• •

First, we propose for the first time a simple way to simulate SIR model on non-periodic lattices by using Newman-Ziff algorithmNewman and Ziff (2001) to track clusters generated by the dynamics and to obtain percolation observables like the mean cluster size;

• •

Second, we investigate if quasiperiodic order can change the SIR model critical behavior. There are some criteria trying to predict if quenched disorder, or quasiperiodic order, can change the critical exponents of a given model, like HarrisHarris (1974b), or Harris-LuckLuck (1993) criteria, respectively. Both Harris and Harris-Luck criteria are superseded by Harris-Barghathi-Vojta criterionVojta et al. (2009); Barghathi and Vojta (2014), however, even Harris-Barghathi-Vojta criterion is known to failSchrauth et al. (2018).

Therefore, to find if quasiperiodic order can change the universality class of SIR model, we obtained the relevant parameters: the order parameter, defined as the average percolating cluster density, the average mean cluster size, and Binder cumulant defined for percolation in order to estimate the epidemic threshold and critical exponents by using standard finite size scaling (FSS) techniques. To correctly determine the universality class, we estimated three independent critical exponent ratios: $1/\nu$, $\beta/\nu$ and $\gamma/\nu$.

This paper is organized as follows: in section II we describe the SIR model and its dynamics, lattice building algorithm and details of calculated parameters by Newman-Ziff algorithm, in section III we discuss our numerical results and in section IV we present our conclusions.

II Model and Implementation

We considered the Asynchronous SIR model, coupled to the two-dimensional (2D) Penrose and Ammann-Beenker quasiperiodic lattices. SIR model is a compartmental model, where a population, composed of vertexes in a lattice, is divided in three compartments, corresponding to three possible states, namely:

• •

Susceptible state: Individuals in that state are not infected, are not infecting and are not immune/dead;

• •

Infected state: Individuals in that state can infect another individuals, and are not immune/dead;

• •

Removed state: Individuals in that state do not infect and are not infecting because of permanent immunity/being dead. They do not interact with individuals in the two other compartments.

We can associate the lattice state as a vector $\psi_{i}$ which can take the values [math], $1$ and $-1$, corresponding to susceptible, infected and removed individuals, respectively.

SIR model can be expressed as a reaction-diffusion process, where individuals can change its state and go from one compartment to anothervan Kampen (1981). Transitions of the individuals between compartments are represented by the following stoichiometric reactions

[TABLE]

The stoichiometric equations (1) can be mapped in a lattice Markovian dynamics with the following rules:

• (1)

We start with only one infected vertex, the “patient zero”, randomly chosen from the lattice with a population of $N$ individuals, attached to the respective lattice vertexes. Along the entire dynamics, two lists are updated: a list of infected individuals and a list of removed ones. The infected list starts with the “patient zero” and the list of removed sites begins empty.

• (2)

Next, we update the system state by randomly choosing an infected vertex $i$ from the infected list and proceed as follows:

• (a)

A random uniform number $x$ in the interval $[0,1)$ is generated. If $x<\lambda$, the infected site is removed from the infected list and placed in the removed list;

• (b)

If $x>\lambda$, we pick a neighboring vertex randomly to add it on the infected list;

• (3)

The step (2) is repeated several times until the system have no infected sites, i.e., reach any absorbing state. For each pass of the dynamics we can increment the dynamics clock by a time interval $1/N_{i}$ where $N_{i}$ is the number of infected particles. Note that the update is made in an asynchronous way.

From the above kinetic Monte Carlo rules, we can see that the recover and infection rates are given by $\mu_{c}=\lambda$, and $\mu_{i}=1-\lambda$, respectively.

Some observables can be obtained at the absorbing state, for example, the number of removed vertexes $N_{r}$ as function of $\lambda$, expressed as

[TABLE]

However, just $N_{r}$ is not sufficient to fully determine the critical behavior because it is not the order parameter. Following classical percolation theoryStauffer and Aharony (1992); Christensen and Moloney (2005), we need the observables related to the clusters formed by the removed sites and to accomplish that, we have to determine the cluster distribution $n_{\mathrm{cluster}}(s)$, i.e., the number of clusters with $s$ removed vertexes. The cluster distribution can be measured by using Newman-Ziff algorithmNewman and Ziff (2001) with the feature of identifying if there is a percolating cluster as a result of the dynamics. What is a percolating cluster depends on the boundary conditions: in the case of periodic boundaries, it is a wrapping cluster; for non periodic boundaries, a percolating cluster is a spanning clusterNewman and Ziff (2001).

From the cluster distribution, we have the fraction of removed vertexes in the finite (non percolating) clusters with size $s$

[TABLE]

The following identity is then obvious

[TABLE]

where $p$ is the total fraction of removed vertexes and $P_{\infty}$ is the fraction of removed vertexes in the percolating cluster. The epidemic phase is then defined in a geographical way, where we should have a percolating cluster emulating an epidemic spreading reaching the most remote points of a lattice. At the epidemic threshold, size distribution of clusters should be a power lawTomé and Ziff (2010), i.e.

[TABLE]

and, in general, we have the following scaling ansatz in the vicinity of the epidemic threshold

[TABLE]

where $\mathrm{F}$ is a fast decaying functionStauffer and Aharony (1992); Christensen and Moloney (2005).

The epidemic outbreak should result in a percolating cluster at the epidemic phase, i.e. $P_{\infty}\neq 0$. Therefore, the order parameter is the percolating cluster density

[TABLE]

where the average is done over dynamics realizations. Another relevant parameter is the mean cluster size $S$

[TABLE]

The mean cluster size $S$ plays the role of the susceptibility in classical percolation theory, which is defined by the following average

[TABLE]

We can define another two observables: first is the mean cluster size $S$ where the summation now includes the percolating cluster, with size $s_{\mathrm{perc}}$

[TABLE]

and the second is the mean quadratic cluster size

[TABLE]

In the infinite size lattice limit, $s_{\mathrm{perc}}\rightarrow\infty$ and in the endemic phase, $s_{perc}=0$, therefore observables $S^{\prime}$ and $M^{\prime}$ only make sense for finite lattice sizes. They are needed to express a quantity analogous to Binder cumulants for percolationde Souza et al. (2011) as follows.

From the scaling ansatz written in Eq.(6), the relevant observables in Eqs. (7), (9), (10), and (11) should scale as

[TABLE]

according to classical percolation theoryStauffer and Aharony (1992); Christensen and Moloney (2005). Combining Eqs. (12), (14), and (15), we find that the quantity

[TABLE]

should be universal at the epidemic threshold where the scaling dependencies cancel out, being analogous to the Binder cumulant for ferromagnetic spin modelsde Souza et al. (2011). Therefore, $U$ obeys the following scaling

[TABLE]

Binder cumulant curve crossings for different lattice sizes give an estimate for the epidemic threshold $\lambda_{c}$, and scaling relations on Eqs. (12), (13), and (17) can be used to estimate the critical exponents $\beta$, $\gamma$ and $\nu$ by standard data collapses, respectively, to completely determine the universality class.

Concerning the particular case of SIR model, only one cluster of size $s=N_{r}$ is grown for each realization of the dynamics. Typical clusters are shown in Fig.(1). In this way, we have for SIR model at the endemic phase $n_{\mathrm{cluster}}(s)=\delta_{s,N_{r}}$ and $s_{\mathrm{perc}}=0$, where $\delta$ is a Kronecker index, and for the epidemic (percolating) phase we should have $n_{\mathrm{cluster}}(s)=0$ and $s_{\mathrm{perc}}=N_{r}$. Therefore $S^{\prime}=N_{r}$ and $M^{\prime}=N_{r}^{2}$. On the other hand, the order parameter reduces to $P_{\infty}=0$ if there is not a percolating cluster and to $P_{\infty}=1$ for the converse. One should note that this particular simpler cluster structure can be simulated without Newman-Ziff algorithmde Souza et al. (2011), but at doing this, we can lose the power to generalize the model for non-periodic and random lattices. Indeed, Newman-Ziff algorithm can be applied with modifications to identify a spanning cluster in the case of non-periodic boundaries or a wrapping cluster, in the case of periodic boundariesYang (2012).

Regarding the lattice structure, Penrose and Ammann-Beenker lattices were generated by projection algorithmde Bruijn (1981); Naumis and Aragón (2003); Mota et al. (2018), where we can define a generation parameter $g$ for these lattices. The number of lattice vertexes is a function of $g$. We depict Penrose and Ammann-Beenker quasiperiodic lattices for $g=3$ in Figs.(2a) and (2b), respectively. Note that the boundary nodes are all identified, which are needed to determine if there is a spanning cluster.

To correctly identify all the boundary nodes, one should search by four classes of vertexes in the following order

• 1:

Vertexes with only two neighbors;

• 2:

The two neighbors of class 1 vertexes;

• 3:

In the case of two neighboring boundary rhombi sharing a vertex, the shared vertex is a boundary one;

• 4:

In the case of two neighboring boundary rhombi sharing an edge, one of two nodes shared is a boundary one. This node have exactly three neighbors, all contained on boundary rhombi. If the two shared edge vertexes have three neighbors, the vertex located farther from the central lattice node is the boundary node.

Boundary rhombi can be identified by looking for extremal grid points in grid lines when calculate the multigrid crossingsMota et al. (2018). Spanning clusters are defined in a different way for Penrose and Ammann-Beenker lattices because of rotational symmetry:

• •

For Penrose lattice, a spanning cluster is a cluster who contains any pair of vertexes symmetrical over a lattice rotation of $\pm 4\pi/5$ about the central node;

• •

For Ammann-Beenker lattice, a spanning cluster is a cluster who contains any pair of vertexes symmetrical over a lattice rotation of $\pi$ about the central node.

Note that we used non periodic boundaries and our results should be prone to finite size boundary effects, but they can be avoided by just using bigger lattice sizes as expected for amenable graphs. Generations from $g=4$ to $g=10$, used in our numerical results, give good data collapses for both lattices, indicating that these sizes are enough to avoid boundary effects.

We repeated SIR model dynamics $10^{6}$ times for every control parameter $\lambda$, starting from “patient zero” until the system evolved to an absorbing state, growing $10^{6}$ clusters. We calculated an ensemble composed of $10^{6}$ measurements of order parameter, mean cluster size and Binder cumulant, one for each cluster generated, to take the relevant ensemble averages. We calculated statistical errors by using the “jackknife” resamplingTukey (1958). In the following section, we discuss our obtained numerical results.

III Results and Discussion

In this section, we show our numerical results of Asynchronous SIR model coupled to Penrose and Ammann-Beenker quasiperiodic lattices. First, we show results for the Binder cumulant in Figs.(3a) and (3b). From Binder cumulant crossings, we estimated the epidemic thresholds at $\lambda_{c}=0.1713(2)$ for Penrose lattice and $\lambda_{c}=0.1732(5)$ for Ammann-Beenker lattice. Binder cumulant expressed in Eq.(16) should not depend on the lattice size at the critical point, but it should depend on the dimensionality and boundary conditions as seen from comparing Figs.(3a) and (3b). Note that the Binder cumulant values at the critical threshold are not the same for Penrose and Ammann-Beenker lattices because of different boundary conditions we adopted to identify spanning clusters.

Following Binder cumulant discussion, we show order parameter results in Figs.(3c) and (3d), i.e, results for the average spanning cluster density. We see that the curves for both Penrose and Ammann-Beenker lattices have the typical sigmoidal shape, indicating a continuous phase transition from the epidemic phase to the endemic phase, where the order parameter vanishes, i.e., where we do not have any spanning cluster. By using Newman-Ziff algorithm, we were able to numerically calculate the averaged mean cluster size. In Figs.(3e) and (3f), we show the susceptibilities for Penrose and Ammann-Beenker lattices, respectively. They diverge at the epidemic threshold $\lambda_{c}$ in the infinite size lattice limit, according to the finite size scaling relation

[TABLE]

which means that the averaged mean cluster size maxima $\lambda_{\mathrm{max}}$ gets closer to the epidemic threshold $\lambda_{c}$ when increasing $N$, as expected for a continuous phase transition and as observed from our data.

Now, we discuss the critical exponent values. We estimated the critical exponent ratios $1/\nu$, $\beta/\nu$, and $\gamma/\nu$ by using finite size scaling regressions of $\left|\frac{\mathrm{d}\left(\ln{P}\right)}{\mathrm{d}\lambda}\right|$, $\ln P$, and $\ln\chi$, evaluated at the critical threshold $\lambda_{c}$, respectively. To extract $\left|\frac{\mathrm{d}\left(\ln{P}\right)}{\mathrm{d}\lambda}\right|$ data, we made linear regressions of the logarithm of order parameter in the vicinity of the epidemic threshold. Error bars are extracted in this case from the least squares method. In Figs.(4a) and Figs.(4b) we show data for $\left|\frac{\mathrm{d}\left(\ln{P}\right)}{\mathrm{d}\lambda}\right|$ as function of $\left(\ln N\right)/2$, allowing us to estimate $1/\nu\approx 0.7883\pm 0.0267$ and $1/\nu\approx 0.7788\pm 0.0260$ for Penrose and Ammann-Beenker lattices, respectively. These numerical values deviate $5\%$ and $4\%$, respectively, from the exact critical exponent ratio $1/\nu=3/4$ of 2D dynamical percolation.

Continuing the discussion of critical exponent ratios, in Figs.(4c) and (4d) we show the order parameter logarithm at the epidemic threshold as function of $\left(\ln N\right)/2$. In this case, finite size regressions should give estimates of $\beta/\nu$ critical exponent ratio, in a way we obtained $\beta/\nu\approx 0.1012\pm 0.0067$ and $\beta/\nu\approx 0.1034\pm 0.0078$ for Penrose and Ammann-Beenker lattices, respectively. Our numerical values for $\beta/\nu$ deviate $3\%$ and $1\%$ from exact $\beta/\nu$ value of 2D dynamical percolation universality class. We show the same for averaged mean cluster size in Figs.(4e) and (4f) at the epidemic threshold as function of $\left(\ln N\right)/2$. In an analogous way of the order parameter, we obtained $\gamma/\nu\approx 1.8094\pm 0.0079$ and $\gamma/\nu\approx 1.7845\pm 0.0264$ deviating $1\%$ and $0.5\%$ from the exact critical exponent ratio $\gamma/\nu=43/24$ for Penrose and Ammann-Beenker lattices, respectively.

By noting that the numerical values obtained from the finite size scaling regressions are all close from exact 2D dynamical percolation values, we can expect that the system falls into 2D percolation universality class. To confirm that, we collapsed all our numerical data shown in Fig.(3) by using the finite size scaling relations written in Eqs. (12), (13), and (17) with the known exact values of 2D dynamic percolation exponent ratios $1/\nu=3/4$, $\beta/\nu=5/48$, and $\gamma/\nu=43/24$. All data collapse panels shown in Fig.(5) are respective to the panels shown in Fig.(3). The data collapses for both Penrose and Ammann-Beenker lattices are a strong evidence that the system falls into dynamic percolation universality class in the same way of 2D square latticeTomé and Ziff (2010); de Souza et al. (2011). Therefore the quasiperiodic order is irrelevant, in agreement with results for classical percolation in the same quasiperiodic latticesSakamoto et al. (1989); Babalievski (1995); Ziff and Babalievski (1999).

IV Conclusions

We proposed a new way to simulate the Asynchronous SIR model on non periodic lattices by using Newman-Ziff algorithm to numerically determine the cluster size distribution and calculate any observable related to percolation for the first time. We applied the Newman-Ziff algorithm to identify the epidemic phase, characterized by a spanning cluster and used this technique to investigate SIR model coupled to Penrose and Ammann-Beenker lattices, and our numerical results suggest that the quasiperiodic order is irrelevant and the system falls into the 2D dynamical percolation universality class in the same way of the Asyncronous SIR model on periodic lattices, in agreement with classical percolation results. In order to characterize the critical behavior, we obtained the order parameter, defined as the spanning cluster density, the averaged mean cluster size and the Binder cumulant ratio defined for percolation. We estimated the epidemic thresholds: $\lambda_{c}=0.1713(2)$ and $\lambda_{c}=0.1732(5)$ for Penrose and Ammann-Beenker lattices, respectively.

V Acknowledgments

We would like to thank CNPq (Conselho Nacional de Desenvolvimento Científico e tecnológico), FUNCAP (Fundação Cearense de Apoio ao Desenvolvimento Científico e Tecnológico) and FAPEPI (Fundação de Amparo a Pesquisa do Estado do Piauí) for the financial support. We acknowledge the use of Dietrich Stauffer Computational Physics Lab, Teresina, Brazil, and Laboratório de Física Teórica e Modelagem Computacional - LFTMC, Piripiri, Brazil, where the numerical simulations were performed.