Generalized Diffusive Epidemic Process with Permanent Immunity in Two Dimensions

This work considers the diffusive epidemic process (DEP) with permanent immunity. The standard DEP has the feature of a metapopulation model [1], which the lattice is a substrate to a fluctuating infected population [2, 3, 4, 5, 6, 7, 8, 9, 10, 11], unlike sedentary models where the individuals are the lattice nodes on typical stochastic lattice processes [12]. In the standard DEP, individuals can diffuse by hopping to neighboring sites and interact only when they divide the same lattice node. The DEP on lattices has a rich critical behavior, whose exponents should differ according to three regimes: $D_{S}>D_{I}$, $D_{S}=D_{I}$, and $D_{S}<D_{I}$ where $D_{S}$ and $D_{I}$ are the diffusion rates of the susceptible and infected individuals, respectively [13, 2, 3, 8, 11]. The critical behavior of standard DEP in any diffusion regime also differs from the paradigmatic directed percolation universality class [14, 8].

In general, one can simulate the DEP as a cellular automata model where reaction-diffusion dynamics can be separated into two synchronous stages [3, 7, 11]. In diffusion, susceptible and infected individuals jump to neighboring nodes according to the $D_{S}$ and $D_{I}$ rates, respectively. All lattice nodes should update their populations after the diffusion stage. Eventually, the reaction stage takes place according to the channels

where Eq. (1) means that susceptible individuals would be contaminated if they share the same lattice node with at least one infected individual with rate $\mu_{c}$ and infected individuals can spontaneously recover with rate $\mu_{r}$. According to the standard definition of DEP [3], individuals become susceptible before they recover (compartment $R$ is just the same as $S$), which means the susceptible individuals can be infected again. We can turn recovered individuals in compartment $R$ permanently immune to the disease to prevent them from interacting with the active particles of compartments $S$ and $I$. Therefore, the population of active particles will not be conserved after the modification.

In sedentary stochastic epidemic models, perfect immunity is introduced in the generalized epidemic process (GEP) [15, 16]. The GEP model describes the spreading of a non-conserved population and falls into the dynamic percolation universality class [17, 16, 14]. Permanent immunity generally makes a stochastic process fall into the dynamic percolation universality class [18]. Permanent immunity is also essential to introduce the effects of vaccination protocols [19]. Recently, some reported results of a vaccination epidemic model using Monte Carlo (MC) simulations, finite-size scaling (FSS), and mean-field approximations [20] presented continuous phase transitions for square lattices with critical exponents compatible with the directed percolation universality class.

An interesting question is the role of permanent immunity when added to the general DEP dynamics and metapopulation models. Therefore, we analyze the effect of permanent immunity in the DEP. We define the analogous model for diffusive individuals, called the generalized diffusive epidemic process (GDEP), and study its properties. By activating lattice nodes with successful contaminations, we address the problem of mapping the epidemic dynamics into a growing process and exploring it as a percolation problem. In addition, the final removed population and typical related observables are analyzed.

The definition of GDEP presented in this manuscript is close to one extension of the susceptible-infected-removed (SIR) model [21] where mobile point particles are free to move in a square space with periodic boundary conditions and interact according to the reaction channels in Eq. (1) if they share a single square cell, a subdivision of the square space. The main difference is the diffusive movement of the individuals is restricted to a two-dimensional lattice in the GDEP, like in transportation networks [22]. Limiting particle movement to a lattice allows one to see the dependence of the critical exponents on the lattice dimension [23, 24] and also analyze the phase transition as controlled by the population size. One can expect the percolation transition to happen when one increases the particle concentration.

We organized this paper as follows. In Sec. II is presented the Markovian dynamics of the GDEP model. In Sec. III are discussed our numerical results. Finally, we give our final considerations in Sec. IV.

This section presents the GDEP at a square lattice with $N=L^{2}$ nodes. The reaction-diffusion dynamics is formulated as a cellular automata model and separated into two synchronous stages where all lattice nodes are updated simultaneously. The dynamics are summarized in the following kinetic Monte-Carlo rules [25, 26]:

1. 1.

   We begin by placing the walker population in the lattice nodes. The walker population size $N_{w}$ is related to the concentration $\rho$ as

   $$N_{w}=\rho N,$$

   (2)

   Increasing the concentration $\rho$, a continuous phase transition from a non-percolating phase to a percolating one will happen, in a critical threshold $\rho_{c}$. Dynamics start from a state with one infected individual and $N_{w}-1$ susceptibles and the stochastic arrays $\textbf{S}=\{S_{1},S_{2},...,S_{N}\}$ and $\textbf{I}=\{I_{1},I_{2},...,I_{N}\}$ store the susceptible and infected populations per node, respectively;

2. 2.

   In a single step, the arrays S and I are updated at the end of each stage:

   • •

     Diffusion: susceptible and infected individuals in all nodes jump with probabilities $D_{S}$ and $D_{I}$, respectively, along one lattice edge;

   • •

     Reactions: if susceptible individuals share a lattice node with at least one infected individual, they will be infected at a rate $\mu_{c}$. In addition, infected sites are removed with probability $\mu_{r}$. In all cases, we used $\mu_{c}=\mu_{r}=0.5$;

3. 3.

   Step 2 is iterated until the system has no infected individual. The system will always reach an absorbing state with a lattice filled only with removed and susceptible individuals.

The Monte Carlo time unit is the time-lapse of updating all lattice nodes.

The Markovian rules above present the following mean-field behavior [3]

where $\rho_{\text{S,I,R}}$ are the particle densities in each compartment $S$, $I$, and $R$, which obey the following conservation law

One can see the striking similarity with the mean-field equations for the susceptible-infected-removed (SIR) model [27, 28, 29, 30], which motivates us to map the model into a growing model where we can use union-find algorithms to study the system as a percolation problem. The mapping is done by activating nodes with successful contaminations in a GDEP realization. One of these union-find algorithms is the Newman-Ziff algorithm [31, 32, 33], which is suitable for finding the cluster-size distribution $n_{\mathrm{cluster}}(s)$. In the following, we call the nodes with at least one successful contamination as activated ones.

We can interpret the activated nodes as infection sources in a lattice of connected places. Given the Brownian nature of the epidemic spreading by the random walkers, the activity is concentrated at the cluster boundaries. We also expect only one giant component in every realization of GDEP, which comprises all the epidemic spreading. The order parameter of the percolation transition is the wrapping probability [27, 28, 29, 30]

The $\left<\right>$ indicates an average over repetitions of GDEP dynamics starting from only one seed. The averaged order parameter is the wrapping probability of a cluster formed by the activated nodes.

From the cluster size distribution $n_{\mathrm{cluster}}(s)$, we can also obtain the cluster size moments which are given by

where $N_{r}$ is the number of activated nodes. Taking the sum on Eq. (6), excluding the wrapping cluster to calculate the first moment, we have the mean cluster size or the susceptibility [29, 30]

Now, taking the sum including the wrapping cluster [27], we can define a ratio analogous to Binder cumulants for spin models

The cluster moments present the following finite-size behavior [34, 35]

from we can obtain the following finite-size relations

where $f_{U,P,\chi}$ are scaling functions.

Also, one can collect the average removed population and its moments to investigate the critical behavior of GDEP, which should present the same critical behavior as the cluster observables. As already stated, the GDEP dynamics form only one giant component because of the Brownian diffusion of the individuals. Therefore, the removed population can be interpreted as the cluster mass in a way that removed population moments should obey the following finite-size scaling relation

We can define the following observables

In Eq. (12), $r$ is the average removed concentration. From Eq. (11), we have the following finite-size behavior of the observables in Eq. (12)

Moreover, we should have a dynamical critical behavior associated with dynamical exponents. One can investigate dynamical exponents and the dynamical behavior by the following set of observables of GDEP dynamics [36], which scales as

at the critical threshold $\rho_{c}$ where $P_{s}(t)$ is the survival probability, i.e., the fraction of GDEP realizations still active after a Monte-Carlo time. The second one is the total number of removed individuals as a function of the time. The survival probability decays as a power-law at the critical threshold $\rho_{c}$. In addition, the number of infected individuals increases as a power-law at the critical threshold $\rho_{c}$.

On finite lattices, one can expect the following finite-size critical behavior of the dynamical properties in Eq. (14

where $g_{P_{s},r}$ are scaling functions. The dynamical critical exponents $z$, $\delta$, and $\eta$ in Eqs. (14) and (15) are related with 2d percolation exponents by the scaling relations $\delta z=\beta/\nu$ and $\eta z=\gamma/\nu$.

We performed Monte Carlo simulations of the GDEP on square lattices. In all realizations, we begin from one seed. The dynamics would end when the system falls into one of the possible absorbing states, with no remaining infected individual. For each seed, the GDEP grows one cluster formed by the activated nodes. All curves are given as functions of node concentration $\rho$.

Snapshots of one realization of the growing process are shown in Fig. 1. We choose the case $D_{S}=D_{I}=0.5$ with concentration $\rho=3.682$ (close to the critical threshold $\rho=3.682(5)$, see Tab. 1). One can see that the cluster activity is concentrated at the cluster boundaries as expected. The growing process generates one cluster, which wraps the lattice between Monte Carlo times $t=400$ and $t=450$ for this particular GDEP realization.

The system displays a continuous phase transition if the basic reproduction number $R_{0}=\frac{\mu_{c}}{\mu_{r}}$ satisfies $R_{0}\geq 1$. Moreover, the critical threshold $\rho_{c}$ increases with the relative diffusion ratio $D_{S}/D_{I}$ and vice versa. We considered in detail the cases (1) $D_{S}=0.25$, $D_{I}=0.75$, (2) $D_{S}=0.5$, $D_{I}=0.5$, and (3) $D_{S}=0.75$, $D_{I}=0.25$, with critical thresholds given in Tab. 1. The critical threshold is consistent with a minor concentration needed for a percolating phase if infected individuals can diffuse faster.

The finite-size scaling behavior given in Eqs. (10) and (13) should be compatible with the exact 2d percolation critical exponents $\beta=5/36$, $\gamma=43/18$, and $\nu=4/3$ on the square lattice. Also, the dynamic critical behavior is compatible with the 2d dynamical exponents $z=1.1309(1)$, $\delta=0.0921(3)$, and $\eta=1.5846(2)$. In the following, we show results of the asymptotic behavior $t\to\infty$ of cluster and removed concentration observables. In addition, we also show in the following the GDEP properties as a function of the time at the critical threshold $\rho_{c}$ that separates the percolating and non-percolating phases.

We show results in the asymptotic limit $t\to\infty$ for the cluster observables in Fig. 2 for the case (2) with $D_{S}=D_{I}=0.5$. The system presents a continuous phase transition in the dynamic percolation universality class where the system goes from a non-percolating phase to a percolating one by increasing the concentration $\rho$. The same behavior is found for the cases (1) $D_{S}=0.25$, $D_{I}=0.75$, and (3) $D_{S}=0.75$, $D_{I}=0.25$ (data not shown). We can conclude that permanent immunity is the fundamental ingredient that lets the system fall into the dynamic percolation universality class. A remaining question is how the system goes from one of the standard DEP universality classes to the dynamical percolation criticality if one introduces a reinfection rate [16], which can stimulate further studies on the partial immunity of diffusive metapopulation models.

The continuous phase transition can also be analyzed by looking at the asymptotic removed concentration observables, shown in Fig. 3. These observables can be particularly useful when the system has no proper boundary conditions to identify a wrapping cluster, which should be valid for complex and scale-free networks [37]. Of course, the critical behavior is also compatible with dynamical percolation exponents. The estimated critical concentration $\rho_{c}$ is the same as shown in Fig. 2.

Regarding dynamical exponents, we show results for the survival probability $P_{s}(t)$ and the removed concentration $r(t)$ in Fig. 4. From the scaling relations written in Eq. (15), we can see that the data has the expected behavior of Eq. (14) where the survival probability decays algebraically with the exponent $\delta$ while the removed concentration grows as a power-law with the critical exponent $\eta$. The removed population at $\rho_{c}$ presents subexponential (power-law) growth, a feature shared by other epidemic models with short-range interactions [21], which explains the flattening of the COVID-19 growing curve [38, 39] as a consequence of containment measures.

We considered a non-sedentary metapopulation model, namely the generalized diffusive epidemic process. This process modifies the standard diffusive epidemic, including permanent immunity. Moreover, the process does not conserve the active population. The Markovian dynamics rules were formulated as a cellular automaton with synchronous updates of all lattice nodes. The reaction-diffusion dynamics is separated into two stages: (i) the Brownian diffusion, where the particles can jump to neighboring nodes, and (ii) contamination and recovery reactions. We also map the metapopulation model as a growing process by activating nodes with successful contaminations, interpreted as infection sources.

The generalized diffusive epidemic process was studied using extensive Monte-Carlo simulations. Our results indicate a continuous phase transition in the dynamic percolation universality class controlled by the population size, irrespective of relative diffusion rates $D_{S}$ and $D_{I}$. Conversely, the standard DEP (without immunity) should fall in three different universality classes [8], one for each diffusion regime. Therefore, the critical behavior of the generalized epidemic process differs from the standard diffusive epidemic process, where the universality class depends on the relative ratio of the diffusive rates. Moreover, the removed population at $\rho_{c}$ presents subexponential (power-law) growth related to containment measures of an epidemic spreading.

We would like to thank CAPES (Coordenação de Aperfeiçoamento de Pessoal de Nível Superior), CNPq (Conselho Nacional de Desenvolvimento Científico e tecnológico), and FAPEPI (Fundação de Amparo à Pesquisa do Estado do Piauí) for the financial support. We acknowledge the Dietrich Stauffer Computational Physics Lab., Teresina, Brazil, and Laboratório de Física Teórica e Modelagem Computacional - LFTMC, Teresina, Brazil, where we performed the numerical simulations.