Exploring the threshold of epidemic spreading for a stochastic SIR model with local and global contacts

TL;DR

This paper investigates the epidemic threshold in a spatial stochastic SIR model with local and global contacts, deriving an analytical threshold estimate and analyzing outbreak dynamics at the threshold.

Contribution

It provides an analytical expression for the epidemic threshold in a spatial stochastic SIR model with mixed contact types, validated by stochastic simulations.

Findings

The threshold for epidemic spread depends on the ratio of transmission to recovery rates.

The analytical threshold estimate closely matches stochastic simulation results.

Outbreak size distributions at the threshold follow a power-law decay with exponent -3/2.

Abstract

The spread of an epidemic process is considered in the context of a spatial SIR stochastic model that includes a parameter $0\le p\le 1$ that assigns weights $p$ and $1- p$ to global and local infective contacts respectively. The model was previously studied by other authors in different contexts. In this work we characterized the behavior of the system around the threshold for epidemic spreading. We first used a deterministic approximation of the stochastic model and checked the existence of a threshold value of $p$ for exponential epidemic spread. An analytical expression, which defines a function of the quotient $\alpha$ between the transmission and recovery rates, is obtained to approximate this threshold. We then performed different analyses based on intensive stochastic simulations and found that this expression is also a good estimate for a similar threshold value of $p$ obtained in the stochastic model. The dynamics of the average number of infected individuals and the average size of outbreaks show a behavior across the threshold that is well described by the deterministic approximation. The distributions of the outbreak sizes at the threshold present common features for all the cases considered corresponding to different values of $\alpha>1$. These features are otherwise already known to hold for the standard stochastic SIR model at its threshold, $\alpha=1$: (i) the probability of having an outbreak of size $n$ goes asymptotically as $n^{-3/2}$ for an infinite system, (ii) the maximal size of an outbreak scales as $N^{2/3}$ for a finite system of size $N$.

Figures

a and 2b. For $p$ lower than $p_{e}\simeq 0.383$, the $N_{I}^{(d)}(\tau)$-curves fall down after having reached a maximum with a value barely greater than 1 individual (Fig.2a). For $p$=0.383, $N_{I}^{(d)}(\tau)$ remains almost constant at the scale of Fig.2a, for $\tau>5$. For all $p$-values the curves show an exponential behavior after an initial transitory regime. It is remarkable that there is none qualitative change in the $N_{I}^{(d)}(\tau)$ behavior for $p$=0.176 ($R_{0}^{(d)}\simeq 1$). The subsequent dynamics of $N_{I}^{(d)}(\tau)$ is shown in Fig.2b. For $p>p_{e}$ ($r>0$) the classic

c (cases with $p<p_{e}^{\prime}$) but $\langle N_{R}(\tau_{ext})\rangle$ is plotted here versus $-1/r$, where $r$ is determined fitting $\langle N_{I}(\tau)\rangle$ to exp($r\tau$) for $\tau\in(5,30)$ as in Fig.3a. The broken lines are linear fits to the points corresponding to a fixed $\alpha$. The fitted slopes are 1.154, 1.357, and 1.669 for $\alpha$=1.2, 1.4 and 1.6 respectively. The unit slope line (solid) is the expected behavior for the SIR model ($\alpha$=1, $p=1$). (b) The case $\alpha$=1.4 for different values of $p$, and the case $p=0.237$ for different values of $\alpha$ are consid

a-4c). The monotonous increase of $\langle N_{R}(\tau_{ext})\rangle$ with $p$ observed in this $p$-range (Fig.3c) could be associated with the increasing number of samples that last longer. For values of $p$ distinctively larger than $p_{e}^{\prime}$ (Fig.4e and 4f) a change of behavior is observed: the points of the figures are grouped into two clearly differentiated sets. The set that groups the higher values of $N_{R}$ forms a sort of cloud that is clearly separated from the rest of the points that are distributed similarly to the cases corresponding to $p<p_{e}^{\prime}$. The points in the

b). The asympthotic behavior of $G_{n}$ for an infinite system (solid line) is also indicated in the figure. By comparison, it can be inferred that for the finite systems the deviation of $G_{n}$ from the asymptotic behavior is due to finite size effects. These effects are expected because the power law behavior of $G_{n}$ implies that the outbreaks at the threshold have a probability of reaching any size, and this is limited by the size of the system.

a we replot the information of Fig.3c for the cases where $p<p_{e}^{\prime}$ taking $-1/r$ as the independent variable. The figure shows that $\langle N_{R}(\tau_{ext})\rangle$ follows a linear behavior with $-1/r$, as in the case of the SIR model, but with a slightly increasing slope for increasing $\alpha$. In Fig.7b we compare the behavior of $\langle N_{R}(\tau_{ext})\rangle$ when approaching the threshold: $\alpha$=1.4, $p_{e}^{\prime}=0.237$, by decreasing $p_{e}^{\prime}$ (keeping $\alpha$=1.4) or decreasing $\alpha$ (keeping $p_{e}^{\prime}=0.237$). In both cases a linear behavior with

Tables

Table 1. Table 1: Threshold value for the global contact parameter, p 𝑝 p , above which exponential epidemic spread is expected in the DA ( p e subscript 𝑝 𝑒 p_{e} ) and the SM ( p e ′ superscript subscript 𝑝 𝑒 ′ p_{e}^{\prime} ) for different α 𝛼 \alpha  values. The value of p e subscript 𝑝 𝑒 p_{e} is given by ( 10 ), while the value of p e ′ superscript subscript 𝑝 𝑒 ′ p_{e}^{\prime} is determined numerically (see text). R 0 ​ ( p e ′ ) subscript 𝑅 0 superscript subscript 𝑝 𝑒 ′ R_{0}(p_{e}^{\prime}) is the basic reproductive ratio for the SM at p = p e ′ 𝑝 superscript subscript 𝑝 𝑒 ′ p=p_{e}^{\prime} .

$\alpha$	$p_e$	$p_e^{\prime }$	$R_0\left(p_e^{\prime }\right)$
1.2	0.383	0.401	1.08
1.4	0.202	0.237	1.17
1.6	0.103	0.151	1.25
2.0	0.000	0.0672	1.39

Equations

$$W_{inf}^j=\left[p\beta \frac{N_I}{N}+(1-p)\beta \frac{1}{4}\sum _{j^{\prime }\in {\nu }_j}{\delta }_{E_{j^{\prime }},I}\right]{\delta }_{E_j,S}$$

$$W_{rec}^j=\gamma {\delta }_{E_j,I}$$

$$\Gamma (t_1)\to \Gamma (t_2)\to ...\Gamma (t_{ext})$$

$$\frac{d{X}_s}{dt}=$$

$$\frac{d{X}_i}{dt}=$$

$$\frac{d\text{SS}}{dt}=$$

$$\frac{d{X}_{si}}{dt}=$$

$$\displaystyle+2p\beta X_{i}\SS-\left.{\bigg{(}}\frac{(1-p)\beta}{4}+\gamma\right.{\bigg{)}}X_{si}$$

$$\frac{d{X}_{ii}}{dt}=$$

$$R_{0}^{(d)}=\int_{0}^{\infty}\beta\bigg{(}\frac{(1-p)X_{si}(t)}{4X_{i}(t)}+pX_{s}(t)\bigg{)}e^{-\gamma t}dt.$$

$${\int }_0^{\infty }\beta X_s(t)e^{-\gamma t}dt=\frac{\beta }{\gamma }=\alpha .$$

$$X_i(\tau )\simeq X_i(0)e^{(\alpha -1)\tau }$$

$$r=\frac{1}{4}\left(1+p+\sqrt{1+10p-7p^2}\right)\alpha -1$$

$$\alpha >\frac{4}{1+p+\sqrt{1+10p-7p^2}}$$

$$p_e=\frac{\alpha +1-\sqrt{(\alpha +7)(\alpha -1)}}{2\alpha }$$

$$r(p,\alpha )=0$$

$$⟨n⟩=\frac{1}{1-\alpha }$$

$$⟨N_R({\tau }_{ext})⟩=\frac{1}{1-\alpha }=-\frac{1}{r}$$

$$⟨N_R({\tau }_{ext})⟩=\frac{a(\alpha )}{r}+b(\alpha )=\frac{c(p)}{r}+d(p)$$

$$\frac{U_n(N)}{U_n(\infty )}=f\left(n/N^{2/3}\right)$$

Full text

1 Introduction

The SIR model is probably the most widely used model in mathematical epidemiology [1, 2, 3]. Since W. Kermack and A. McKendrick employed it to describe the development of a plague epidemic in Bombay [4], it has been extensively used for many purposes. Sometimes as the kernel of more complex epidemiological models in the field of infectious disease transmission [2, 5] and other times to study spreading phenomena in other fileds such as rumor propagation [6], computer viruses [7], information diffusion in Web forums [8] or investors’ behavior in stock markets. [9].

In the deterministic version of the model, the three epidemiological classes are represented by three continuous variables (S: susceptible, I: infected, and R: recovered) that evolve in time according to a set of ordinary coupled differential equations. It is easy to treat mathematically, has a straightforward interpretation and leads to simple and important predictions [4]. However, from the very beginning it was discovered that for some applications, it was necessary to consider the fluctuations and the discrete nature of the population and of the processes involved (infection and recovery) [1, 10]. For example, using a stochastic version of the SIR model, M. S. Bartlett showed that stochasticity was an essential aspect to be considered to explain the persistence of measles as a function of the city size in several cities of England and Wales [10, 11].

The threshold for epidemic spreading that in both versions of the SIR model (deterministic and stochastic) occurs when the transmission rate ($\beta$) equals the recovery rate ($\gamma$) and so, the quotient $\alpha\equiv\beta/\gamma=1$ is of particular interest. The properties of the stochastic SIR model at the threshold have been extensively studied lately by rigorous mathematical theory, more empirical treatments, and computer simulations [12, 13, 14, 15].

In the present work we study a stochastic SIR model on the lattice with local and global contacts where the weight of the global contacts is taken into account through a single parameter $p$. It is inspired in Watts-Strogatz model [16] and has already been used by other authors to study transmission of diseases of high $\alpha$ values [17, 18, 19, 20]. The purpose of using these simple models in this field is to gain insight into the qualitative trends observed when the global contacts are lowered. To obtain an accurate description of the transmission process, a more realistic network well suited to represent human interactions should be used. On the other hand, SIR type models that combine local contacts in a square lattice with some kind of global connection have been proposed to study the spread of infectious diseases in plants and animals [21, 22, 23]. In these cases, where the lattice could eventually be a good enough representation of the “spatiality” of the problem, predictions could also have a quantitative character.

In this work we study the model for low $\alpha$ values since we are interested in the behavior around the threshold for epidemic spreading. In models where the hypothesis of uniform mixing holds, the threshold is given by the condition $R_{0}$=1, where $R_{0}$ is the basic reproductive number, defined as the average number of secondary infections produced by one infected individual in a completely susceptible population [2]. This is the case of the classic SIR model where $R_{0}=\alpha$. In models where uniform mixing is broken, as in the one studied in the present work, it is known that the concept of $R_{0}$ is meaningless [2, 24, 25], and so, the threshold condition is something to be explored. We here show that it is possible to define a threshold for exponential spreading that in this model (with two free parameters) is not a point as in the classical SIR model but a curve in the $(p,\alpha)$ plane. We found that for any $\alpha>1$ considered, by sufficiently lowering the value of $p$, the threshold can be crossed producing a drastic reduction of the probability for a major outbreak. We managed to characterize several model properties around the threshold through the analysys of stochastic simulations and by using a deterministic approximation developed in a previous work [20].

The work is organized as follows: we first introduce the stochastic model (SM) and the deterministic approximation (DA). Then, using the DA, we can build a phase diagram in the parameter space detecting the region in which it is possible to prevent the exponential spread. We perform intensive stochastic simulations and find that average magnitudes follow the same trends predicted by the DA. In particular the dynamical behavior of the average number of infected individuals makes it possible to unambiguously define the threshold for exponential spreading in the SM for all the values of $\alpha$ considered. We then analyze the distributions of the outbreak sizes at the epidemic threshold and obtain several features that are known to hold for the classical SIR model. A summary of our findings is finally given in the conclusions.

2 The stochastic model and the deterministic approximation

2.1 The stochastic model

We consider the stochastic model studied in [19], but since in the present work we focus on the epidemic spread only, we do not consider the births and deaths. The population consists of $N$ individuals identified with the sites of an $L\times L$ square lattice with periodic boundary conditions. They may be in one of the three epidemiological states: $S$, $I$ or $R$ (susceptible, infected or recovered). The dynamics of the model are described by a stochastic Markovian process in which an individual may experience one of these two changes in its state: $S\rightarrow I$ (infection) or $I\rightarrow R$ (recovery). Infections occur through infective contacts among susceptible and infected individuals. We define an infective contact as a contact between two individuals such that if one individual is susceptible and the other infected, the former becomes infected. We assume that an individual at a given site has an infective contact with a randomly chosen individual on the lattice with transition rate $p\beta$, and with one of its four nearest neighbors (also randomly chosen) with transition rate $(1-p)\beta$. By changing $p$, we may change the relative weight of the global and local contacts in the system. The case $p$=1 corresponds to the classical SIR model (uniform mixing) where an individual may contact any other individual in the system. On the other hand, the case $p$=0 corresponds to the square lattice where an individual may only contact one of its four nearest neighbors. Recovery from infection in this model is the same for every site and occurs at a transition rate $\gamma$.

The state of the system, $\Gamma=(E_{1},E_{2},...,E_{N})$, is defined by specifying $E_{j}$ (the state of site $j$) for the $N$ sites of the lattice. Three variables of interest are the number of individuals in the system that are in state $S$, $I$ and $R$ that we call: $N_{S}$, $N_{I}$, and $N_{R}$ respectively.

The probability transition rates for infection and recovery processes at site “$j$” are

[TABLE]

where $\delta_{AB}$ is one if states $A$ and $B$ are the same, and zero if not. Index $j^{\prime}$ in the sum runs over the $4$ neighbors of site $j$ (we call this set of sites $\nu_{j}$).

Stochastic simulations are performed using the Gillespie algorithm [26]. Each simulation begins in the same initial state where $N-1$ individuals are susceptible and one individual is infected. A Markov chain

[TABLE]

with a set of exponentially distributed times $t_{1}$, $t_{2}$, …$t_{ext}$ is generated, where $t_{ext}$ (the time at extinction) is the first time with $N_{I}=0$. $\Gamma(t_{ext})$ is an absorbing state and $N_{R}(t_{ext})$ (the number of recovered individuals) is the total number of individuals that have experienced the infection during the dynamic evolution of the disease from the initial state to its extinction. We also refer to $N_{R}(t_{ext})$ as the ”outbreak size”.

For most of the calculations in this work we take $L=800$, $N=N_{0}=640,000$, but in Section 3.2.3 we also consider $N=2N_{0}$ and $N=4N_{0}$ to study size effects.

2.2 Deterministic approximation

In the present work we use a deterministic approximation of the stochastic model in [19]. In this approximation, developed in [20], the local infective contacts are treated using a pair approximation scheme with second moment closure.

In the deterministic context we denote $N_{S}^{(d)}$, $N_{I}^{(d)}$, $N_{R}^{(d)}$, the time functions whose values are the number of individuals of each type, and $N_{SS}^{(d)}$, $N_{SI}^{(d)}$, $N_{SR}^{(d)}$, $N_{II}^{(d)}$, $N_{IR}^{(d)}$, $N_{RR}^{(d)}$, the similar ones for the number of pairs formed by two neighboring individuals whose type corresponds to the subscripts; for example, $N_{SR}^{(d)}$ represents the number of pairs formed by a susceptible and a recovered neighboring individual.

Then we define the nine functions $X_{i}=N_{I}^{(d)}/N$, $X_{s}=N_{S}^{(d)}/N$, $X_{r}=N_{R}^{(d)}/N$, $\SS=N_{SS}^{(d)}/N$, $X_{si}=N_{SI}^{(d)}/N$, $X_{sr}=N_{SR}^{(d)}/N$,   $X_{ii}=N_{II}^{(d)}/N$, $X_{ir}=N_{IR}^{(d)}/N$, $X_{rr}=N_{RR}^{(d)}/N$.

In [20] we construct a system of nine differential equations having these nine unknowns. Then, using several relationships between the unknowns, the system is reduced to five equations and the unknowns $X_{s}$, $X_{i}$ ,$\SS$, $X_{si}$, $X_{ii}$. In the present case, where birth and death are not considered, the five equations of the DA are:

[TABLE]

We take $X_{s}(0)=(N-1)/N,\ X_{i}(0)=1/N,\ \SS(0)=(2N-4)/N,\ X_{si}(0)=4/N,\ X_{ii}(0)=0$ as initial values, which correspond to the presence of only one infected individual, and $N-1$ susceptible ones.

Now we define $R_{0}^{(d)}$, which we will call “$R_{0}$ of the deterministic approximation”. Consider a fixed time $t$. By equation (2) the amount of new infective cases in the (infinitesimal) time interval between $t$ and $t+dt$ is N\beta\big{(}{(1-p)X_{si}(t)}/{4}+pX_{s}(t)X_{i}(t)\big{)}dt. If the only initial infected individual remains in this state at time $t$, we estimate the number of new infections generated between times $t$ and $t+dt$ by this individual as the quotient between that amount and $N_{I}^{(d)}(t)=NX_{i}(t)=$ number of infected individuals at time $t$. To take into account the recovery possibility of this individual, we multiply by the “probabilistic” factor $e^{-\gamma t}$, thus arriving at

[TABLE]

As expected, for $p=1$ we obtain the value $\alpha$ corresponding to the SIR model (see [27]):

[TABLE]

When necessary, we will calculate $R_{0}^{(d)}$ numerically, for different values of $p$, simultaneously with the resolution of the system, which is performed using the Euler’s Method with a time step of 0.01 day.

3 Results and Discussion

Henceforth, we take $1/\gamma$ (the mean duration of infection) as the unit of time. This is equivalent to making the change of variables: $\tau=\gamma t$. Doing this in equations (1) to (7) is equivalent to substituting $\gamma$ by 1, $\beta$ by $\beta/\gamma$ and $t$ by $\tau$. Since $\beta/\gamma$=$\alpha$, there are only two free parameters in our study: $\alpha$ and $p$. In this work we focus on the case $0\leq p\leq 1$ and $1\leq$ \alpha$$\leq 2.

3.1 Deterministic approximation predictions

In this section we explore the dynamical behavior of the DA system for low $\alpha$ values.

3.1.1 Condition for exponential epidemic spread

For $p=1$ the DA reduces to the SIR model where it is well known that $X_{i}(\tau)$ behaves as

[TABLE]

as long as $X_{s}(\tau)$ does not fall appreciably from $1$. Then the condition \alpha$$>1 determines the exponential spread of the epidemic in this case.

For $p<1$, in our previous work [20], we checked that after an initial transient time, for $\alpha$=7 and 17, $X_{i}(\tau)$ grows exponentially but with a reduced exponent $r$ well approximated by:

[TABLE]

Assuming for the moment that $X_{i}(\tau)$ presents this behavior also for low $\alpha$ values, the condition for exponential growth $r>0$ leads to:

[TABLE]

Since the right side of (9) is a strictly decreasing function of $p$ that takes all the values of the interval [1,2], the condition \alpha$$>1 does not guarantee that inequality (9) holds.

In Fig. 1 we plot the curves $r=0$ and $R_{0}^{(d)}=1$ in the $(p,\alpha)$ plane. While the curve for $r=0$ is explicitly obtained from (8), the curve for $R_{0}^{(d)}=1$ is computed numerically from (6) by sweeping the $(p,\alpha)$ plane. These curves divide the plane into three regions. In principle, one would expect exponential epidemic growth only in region III. We denote $p_{e}$ the value of $p$ that satisfies $r=0$ for a given $\alpha$. From (8), $p_{e}$ is well approximated by the value:

[TABLE]

3.1.2 Numerical results

We first show that after an initial transient time, $X_{i}(\tau)$ behaves as exp($-r\tau$) with $r$ given by expression (8) also for low $\alpha$ values, and that therefore Fig.1 contains useful information. For this purpose we first keep $\alpha$=1.2 fixed and compute $N_{I}^{(d)}(\tau)$ for different values of $p$ crossing the different regions of the phase diagram of Fig.1.

The results are shown in Fig. 2a and 2b. For $p$ lower than $p_{e}\simeq 0.383$, the $N_{I}^{(d)}(\tau)$-curves fall down after having reached a maximum with a value barely greater than 1 individual (Fig.2a). For $p$=0.383, $N_{I}^{(d)}(\tau)$ remains almost constant at the scale of Fig.2a, for $\tau>5$. For all $p$-values the curves show an exponential behavior after an initial transitory regime. It is remarkable that there is none qualitative change in the $N_{I}^{(d)}(\tau)$ behavior for $p$=0.176 ($R_{0}^{(d)}\simeq 1$). The subsequent dynamics of $N_{I}^{(d)}(\tau)$ is shown in Fig.2b. For $p>p_{e}$ ($r>0$) the classical epidemic behavior is observed. It is remarkable that for $r=0$, $N_{I}^{(d)}(\tau)$ decreases slightly from 1 individual after 800 infectious periods. It means that the DA predicts that even when an epidemic does not develop, the disease may persist for a very long time. The total number of individuals that get infected in the population is the asymptotic number of recovered individuals, $N_{R}^{(d)}(\infty)$. In Fig.2c it is plotted as a function of $p$ for $\alpha$=1.2, 1.4 and 1.6. The curves show a sharp increase for $p=p_{e}$, giving support to the idea that in the phase diagram of Fig.1 region III corresponds to a system behavior qualitatively different than in regions I and II, not only for the dynamic evolution of $N_{I}^{(d)}(\tau)$, but also for the total number of infections produced in the population.

3.2 Stochastic model predictions

3.2.1 Averaged magnitudes. Comparison between deterministic

and stochastic approaches

It is not obvious that the DA equations are useful to approach the dynamics of the stochastic model. For high values of $\alpha$, in our previous work [20] we checked that the number of infected individuals averaged over several samples was well approximated by the DA during the epidemic spread. Moreover, the fraction of infected individuals for each single epidemic sample was almost identical to each other except for a shift in time that was caused by the different instants at which the epidemic was triggered (see Fig.1 in [20]). Of course, there was always a probability that extinction would occur before the infected individual can infect anyone, and a fraction of the samples did not lead to an epidemic.

For low values of $\alpha$, the situation for individual samples is quite different. Almost half of the samples (48%) become extinct within the infectious period ($\tau=1$), but the samples that spread for longer times present a very different dynamical behavior. Then, in this case, with the DA we could at best obtain a good approximation for the average behavior of the SM.

We now present some results showing that the DA captures the essential behavior of the averaged magnitudes obtained with the SM. In Fig.3a and 3b the dynamical behavior of the averaged number of infected individuals for $\alpha$=1.2 and different values of $p$ is shown for the SM. As can be observed in the figure, when $p$ is changed, the time evolution of $\langle N_{I}\rangle$ in the SM (Fig.3a and 3b) presents a similar behavior to that obtained with the DA for $N_{I}^{(d)}$ (Fig.2a and 2b). In particular, it can be seen that there should be a value, $p_{e}^{\prime}$, for $p$ (which corresponds to the $p_{e}$ of the DA) such that $\langle N_{I}\rangle$ grows exponentially with time for $p>p_{e}^{\prime}$. To determine $p_{e}^{\prime}$ we analyze $\langle N_{I}\rangle$ for $\tau\in(5,30)$ and check whether it is approximately constant in this interval. For $\alpha$=1.2 a value $p_{e}^{\prime}\simeq 0.401\gtrsim p_{e}\simeq 0.383$ is obtained (see Fig.3a). The magnitude of $\langle N_{I}\rangle$ is, however, strongly overestimated by $N_{I}^{(d)}(\tau)$ computed with the DA because of extinctions.

The values of $p_{e}^{\prime}$ for other $\alpha$ values (determined as described above) are compared with the deterministic values $p_{e}$ in Table 1. It can be seen that $p_{e}^{\prime}\gtrsim p_{e}$ where the difference $p_{e}^{\prime}-p_{e}$ increases with $\alpha$. The value $R_{0}(p_{e}^{\prime})$ is also listed in the table to visualize how much $R_{0}$ deviates from the value of 1 when exponential growth has already been prevented. $R_{0}$ is very well approximated by the DA in the whole range of $p$ values for $\alpha$=1.2. For larger $\alpha$ the agreement worsens but even for $\alpha$=1.6 the relative agreement between SM and DA computations of $R_{0}$ is within 3%.

Fig.3c shows the total number of infected cases until extinction averaged over several samples in the SM for different values of $p$ and $\alpha$. As can be observed in the figure, a sudden increase of this magnitude occurs for $p\sim p_{e}^{\prime}$, as was previously observed in Fig.2c for the DA results. In the SM, however, the transition is not as sharp as the one observed for the DA. $N_{R}^{(d)}(\infty)$ also strongly overestimates $\langle N_{R}(\tau_{ext})\rangle$. For $\alpha$=1.2 and $p=1$, for example, $N_{R}^{(d)}(\infty)\simeq 200,770$ while $\langle N_{R}(\tau_{ext})\rangle\simeq 33,326$. We may conclude that the DA predicts the correct trend for $\langle N_{R}(\tau_{ext})\rangle$ and its qualitative change of behavior when the proportion of global contacts, $p$, is increased. The origin of the quantitative discrepancy between $N_{R}^{(d)}(\infty)$ and $\langle N_{R}(\tau_{ext})\rangle$ will be studied in the following section considering the contribution of each single sample to the last average.

3.2.2 Final size of epidemics.

The stochastic model has an inherent probabilistic nature and, therefore, so do its predictions. When the infected individual enters the fully susceptible population, an epidemic (of a given size) may be unleashed, or there may be no epidemic. Average magnitudes are useful analysis tools but they do not have an epidemiological correlate. In order to study the behavior of individual samples, and how it changes as a function of the global contact parameter, $p$, in Fig.4 we plot the value of the final size of each epidemic as a function of the time it lasts for a set of 100,000 samples.

The figure shows that for $p\lesssim p_{e}^{\prime}$ there is a strong correlation between the time elapsed until extinction and the number of infections produced during the epidemic spread (Fig. 4a-4c). The monotonous increase of $\langle N_{R}(\tau_{ext})\rangle$ with $p$ observed in this $p$-range (Fig.3c) could be associated with the increasing number of samples that last longer. For values of $p$ distinctively larger than $p_{e}^{\prime}$ (Fig.4e and 4f) a change of behavior is observed: the points of the figures are grouped into two clearly differentiated sets. The set that groups the higher values of $N_{R}$ forms a sort of cloud that is clearly separated from the rest of the points that are distributed similarly to the cases corresponding to $p<p_{e}^{\prime}$. The points in the cloud have to be identified with the samples that experienced exponential epidemic spread, while the other set includes the ones that became extinct without developing a major outbreak. For the case $p=1$, for example, if $\langle N_{R}(\tau_{ext})\rangle$ is computed only considering the points in the cloud, a value of around 200,799 individuals is obtained, which is very close to that predicted by the DA (200,770). This confirms that the main failure of the DA is not accounting for extinctions. Finally, for $p\gtrsim p_{e}^{\prime}$ (Fig.4c and 4d) there is a transition region where an intermediate behavior is observed between those described for low and high values of $p$. In particular, for $p=0.45$ (Fig.4d), the grouping of the points for high $N_{R}$ begins to be noticeable, but the gap in $N_{R}$-values observed in Fig.4e and 4f is absent yet. In Fig.5a we present the histograms with the probability $G_{n}$ of obtaining an outbreak of size $N_{R}(\tau_{ext})=n$ for the cases considered in Fig.4. The curves quantify the changes in the distribution of outbreak sizes around $p_{e}^{\prime}$, observed qualitatively in Fig.4. In particular, for $p=p_{e}^{\prime}$, a power law behavior for $G_{n}\propto n^{-3/2}$ is observed for $10\lesssim n\lesssim 2000$. The same behavior has been reported for the $G_{n}$ corresponding to the classical stochastic SIR model when $\alpha$=1 [13, 15].

3.2.3 System behavior at and near the threshold. Comparison with the SIR model.

The stochastic SIR model presents a threshold for epidemic spreading at $\alpha\equiv\beta/\gamma$=1 that has been studied by other authors [12, 13, 15, 14]. For the SM we defined the threshold from the equation

[TABLE]

where $r$ is the exponent obtained by fitting $\langle N_{I}\rangle$ to $k$.exp($r\tau$) for $\tau\in(5,30)$. Thus, for each value of $\alpha$ there is a value of $p=p_{e}^{\prime}$, that defines the threshold: $(p_{e}^{\prime},\alpha)$. We computed $G_{n}$ for the SM, for different ($p_{e}^{\prime}$, $\alpha$), and obtained the same behavior in all cases (Fig.5b). The asympthotic behavior of $G_{n}$ for an infinite system (solid line) is also indicated in the figure. By comparison, it can be inferred that for the finite systems the deviation of $G_{n}$ from the asymptotic behavior is due to finite size effects. These effects are expected because the power law behavior of $G_{n}$ implies that the outbreaks at the threshold have a probability of reaching any size, and this is limited by the size of the system.

In ref. [13, 14, 15] the authors studied the size effects for the SIR model at the threshold and argue that the maximal size of an outbreak should scale as $N^{2/3}$. To check their hypothesis they computed the probability $U_{n}(N)$ of having an outbreak of size greater than $n$ in a system of size $N$, and found that the curves $U_{n}(N)/U_{n}(\infty)$ vs $n/N^{2/3}$ collapse for different values of $N$. To see whether the SM verifies the same scaling law we performed simulations for the cases $\alpha$=1.4, $p=0.237$ and $\alpha$=1, $p=1$ (SIR-model) considering three system sizes: $N=N_{0}=640,000$, $N=2N_{0}$, and $N=4N_{0}$. The collapse of the curves, shown in Fig.6, indicate that the SM and the SIR model obey the same scaling law. However, the collapse curve for $\alpha$=1.4 is different from that for $\alpha$=1.

In the case of the SIR model, the asymptotic behavior of $G_{n}$ as $n^{-3/2}$ implies that, for an infinite system, the average size of an outbreak, $\langle n\rangle=\sum nG_{n}$ diverges. However, below the threshold, it is well known [13, 15, 14] that

[TABLE]

and a finite number of average infected individuals is expected for $\alpha$ $<1$. For a finite system of size $N$, if $1/(1-\alpha)\ll N$, one would expect eqn.(12) to hold even for a finite system [14].

In the case of the SM studied in the present work we do not have an expression equivalent to eqn.(12). In order to explore the behavior of $\langle N_{R}(\tau_{ext})\rangle$ with the ”distance” to the threshold, in Fig.7a we replot the information of Fig.3c for the cases where $p<p_{e}^{\prime}$ taking $-1/r$ as the independent variable. The figure shows that $\langle N_{R}(\tau_{ext})\rangle$ follows a linear behavior with $-1/r$, as in the case of the SIR model, but with a slightly increasing slope for increasing $\alpha$. In Fig.7b we compare the behavior of $\langle N_{R}(\tau_{ext})\rangle$ when approaching the threshold: $\alpha$=1.4, $p_{e}^{\prime}=0.237$, by decreasing $p_{e}^{\prime}$ (keeping $\alpha$=1.4) or decreasing $\alpha$ (keeping $p_{e}^{\prime}=0.237$). In both cases a linear behavior with similar values for the slope of the straight lines was obtained. Of course, if the threshold is approached such that $-1/r$ increases beyond the scale of the figure, the points depart from the linear behavior due to finite size effects.

Finally, concerning the goodness of the DA to approach the SM behavior, it is worth mentioning that the case $p=0$ (where only local contacts are present) was considered by Souza et al. [28] both for the SM and the DA. The authors found that the SM presents a percolation-like transition for $\alpha$ $\simeq$ 4.6657, while their pair-wise approximation (our DA) predicts a transition at $\alpha$ =2. It is worth noting that the consideration of global contacts drastically improves the performance of the DA as can be observed in Table 1.

4 Conclusions

In this work we studied a stochastic epidemiological SIR model with local and global contacts where the weight of global contacts is given by a parameter $p$. Taking appropriate units of time, there is only another free parameter in this model, the quotient between transmission and recovery rates: $\alpha$.

By using a deterministic approximation of the model, we were able to construct a phase diagram in the ($p$, $\alpha$) plane where we identified a region in which exponential epidemic spread is prevented even when the basic reproductive ratio of the model, $R_{0}$, remains above 1. We obtained an analytical expression $p_{e}$($\alpha$) that approximates, for each $\alpha$, the threshold value of $p$ for the exponential growth of the number of infected individuals.

We found that these predictions of the DA are closely linked to the behavior of the average number of infected individuals in the SM, $\langle N_{I}\rangle$. In this case, we could define a threshold value $p_{e}^{\prime}$ (which is slightly larger than $p_{e}$) such that $\langle N_{I}\rangle$ does not grow exponentially with time for $p<p_{e}^{\prime}$. The absolute values of the average number of infected individuals and the average epidemic size do not match the corresponding magnitudes obtained by the DA because of extinctions. For $p>>p_{e}^{\prime}$ the agreement could be recovered if the average is performed using the samples that lead to epidemics exceeding a threshold that, for high $p$ values, is well defined.

We found that the SM behavior around the threshold is closely related to the behavior of the classic stochastic SIR model around its threshold. We summarize below the similarities and differences between both behaviors that emerge from the comparison of our results for the SM and what is known from the literature for the standard SIR model [13, 15, 12, 14].

In the SIR model, the threshold is given by $\beta/\gamma\equiv\alpha$=1. As $R_{0}$=$\alpha$ and the exponent for epidemic spreading is $r=\alpha-1$, the threshold may be expressed unambiguously as $R_{0}=1$ or $r=0$.

In the SM both $R_{0}$ and $r$ are functions of $\alpha$ and $p$. As is well known (and was verified in the present work) if the population is not well mixed $R_{0}$ is not a useful concept [2, 24, 25]. The threshold is obtained from the condition $r=0$, which in this case does not give a parameter value, but a relation between them: $r(p,\alpha)=0$ that leads to $p=p_{e}^{\prime}(\alpha)$.

• -

Below the threshold, the average size of an outbreak in the SIR model verifies:

[TABLE]

For the SM we empirically obtained that $\langle N_{R}(\tau_{ext})\rangle$ grows approximately linearly with $-1/r$. However, if $p$ or $\alpha$ is kept fixed, the behavior is strictly linear:

[TABLE]

with $a$, $b$, $c$ and $d$ smooth functions. In both cases the above relationships are no longer valid in the vicinity of the threshold due to finite size effects.

• -

Strictly at the threshold, the probability of having an outbreak of size $n$ in the SIR model is $G_{n}=kn^{-3/2}$ for $n\in(n_{1},n_{2})$, where $n_{2}$ increases with the size of the system, $N$.

In the SM the same behavior is observed for all the $(p_{e}^{\prime},\alpha)$ cases considered.

• -

The probability $U_{n}(N)$ that the outbreak size is at least $n$ in a population of size $N$ shows the following scaling law for the SIR model [13, 15].

[TABLE]

The same scaling law was observed in the present work for the SM, where the scaling function $f$ depends on the threshold point $\left(p_{e}^{\prime}(\alpha),\alpha\right)$, $f$ being an increasing function of $\alpha$.

It would be interesting to know if those properties at the threshold, which the SM shares with the stochastic SIR model, hold for other networks with global and local contacts

Our study highlights the importance of keeping the global contacts as low as possible as a key measure to prevent large epidemics and points out that a substantial improvement of the epidemiological status (where exponential epidemic spread is prevented) could be accompanied by an insignificant reduction of $R_{0}$, which remains with values well above one. Even though the SM model is very simple, from the epidemiological perspective and in the treatment of the spatial structure, we believe that our conclusions could be taken as a basis for exploration by more complex models in specific contexts.

5 Acknowledgments

This work was supported by Agencia Nacional de Promoción Científica y Tecnológica-ANCPyT grant PICT2010-0707 and Universidad Nacional de La Plata grant X805(2018-2019). G.F. is member of the Scientific Career of Consejo Nacional de Investigaciones Científicas y Tecnológicas-CONICET (Argentina).

6 Figure captions

Figure 1

Phase diagram for epidemic growth in the deterministic approximation. The regions above and below the dashed line correspond to the points ($p$, $\alpha$) with $R_{0}^{(d)}$ greater and lower than 1 respectively. The regions above and below the solid line correspond to the rate of epidemic growth, $r$, greater and lower than 0 respectively.

Figure 2

DA results. (a) Number of infected individuals ($N_{I}^{(d)}$) as a function of time ($\tau$) for $\alpha$=1.2 and $p$=0.10, 0.176, 0.30, 0.35, 0.383, 0.40 and 0.45. The values $p$=0.176 and 0.383 correspond to $R_{0}^{(d)}\simeq 1$ and $r\simeq 0$ respectively. (b) $N_{I}^{(d)}$ for $\alpha$=1.2 on an expanded time scale and $p$=0.10, 0.30, 0.35, 0.383, 0.40, 0.45, 0.50, 0.60 and $1.0$. The solid line corresponds to $p=0.383\simeq p_{e}$. (c) Asymptotic number of infected individuals $N_{R}^{(d)}(\infty)$ as a function of $p$ for three $\alpha$ values. Vertical lines indicate the $p_{e}$ values: 0.383, 0.202, 0.103 corresponding to $\alpha$=1.2, 1.4, and 1.6 respectively.

Figure 3

Average magnitudes in the SM. (a) Averaged number of infected individuals, $\langle N_{I}\rangle$, for $\alpha$=1.2 and different values of $p$. The thin horizontal line has been drawn to show that $\langle N_{I}(\tau)\rangle$ for $p=0.401$ remains almost constant for $\tau\in(5,30)$. (b) $\langle N_{I}\rangle$ for $\alpha$=1.2 on an expanded time scale and for $p=0.10,0.30,0.35,0.401,0.45,0.50,0.60$ and $1.0$. The solid line corresponds to $p=0.401\simeq p_{e}^{\prime}$. Each curve is the average of $m$ independent stochastic simulations, where $m=100,000$ for $p\geq 0.50$ and $m=1,000,000$ for $p<0.50$. (c) Averaged number of individuals that have experienced the infection until extinction for different values of $\alpha$ and $p$. Each point corresponds to an average over 100,000 samples when $\langle N_{R}(\tau_{ext})\rangle>5000$, and to an average over 1,000,000 samples when $\langle N_{R}(\tau_{ext})\rangle<5000$. The crosses correspond to $p=0.401\simeq p_{e}^{\prime}(1.2)$, $p=0.237\simeq p_{e}^{\prime}(1.4)$, and $p=0.151\simeq p_{e}^{\prime}(1.6)$. The broken lines have been drawn to guide the eye. The vertical solid lines indicate the $p_{e}$ values corresponding to the beginning of exponential spread in the DA, as in Fig.2c.

Figure 4

Number of individuals that have experienced the infection in a given stochastic simulation, $N_{R}(\tau_{ext})$, as a function of the duration of the corresponding simulation, $\tau_{ext}$, for $\alpha$=1.2. Different panels correspond to simulations performed for different $p$ values. Each panel contains 100,000 points, each one corresponding to an independent stochastic simulation with identical initial conditions.

Figure 5

Probability $G_{n}$ of obtaining an outbreak of size $N_{R}(\tau_{ext})=n$. (a) $G_{n}$ for $\alpha$=1.2 and different $p$ values. The different curves represent the distributions of points in the $y$-axis of the different pannels of Fig.4. (b) $G_{n}$ at the threshold for different cases. Black line: $n^{-3/2}/\sqrt{4\pi}$ (asymptotic behavior for an infinite system in the case $p=1$, $\alpha=1$ [15]).

Figure 6

Scaling behavior of the normalized cumulative distribution $U_{n}(N)/U_{n}(\infty)$ at the threshold. Six different systems, corresponding to two sets of parameters at the threshold ($p(\alpha),\alpha$), and three system sizes ($N$), were considered. When plotted as a function of $n/N^{2/3}$ the distributions corresponding to $\alpha=1.4$ or $\alpha=1$ colapse in the upper or lower curves, respectively. $N_{0}=640,000$.

Figure 7

Averaged number of individuals that have experienced the infection until extinction below the threshold. (a) The cases $\alpha$=1.2, 1.4, and 1.6 for different values of $p$ are considered. The points correspond to the same simulations of Fig.3c (cases with $p<p_{e}^{\prime}$) but $\langle N_{R}(\tau_{ext})\rangle$ is plotted here versus $-1/r$, where $r$ is determined fitting $\langle N_{I}(\tau)\rangle$ to exp($r\tau$) for $\tau\in(5,30)$ as in Fig.3a. The broken lines are linear fits to the points corresponding to a fixed $\alpha$. The fitted slopes are 1.154, 1.357, and 1.669 for $\alpha$=1.2, 1.4 and 1.6 respectively. The unit slope line (solid) is the expected behavior for the SIR model ($\alpha$=1, $p=1$). (b) The case $\alpha$=1.4 for different values of $p$, and the case $p=0.237$ for different values of $\alpha$ are considered. Broken lines are linear fits corresponding to each case with slopes 1.357 ($\alpha=1.4$) and 1.382 ($p=0.237$). Note that the two cases correspond to two perpendicular directions on approaching the threshold ($p=0.237$, $\alpha$=1.4) in the ($p,\alpha$)-plane.