in

A single-agent extension of the SIR model describes the impact of mobility restrictions on the COVID-19 epidemic

Combining agent mobility patterns and SIR model

To take into account agent mobility19 in a scenario compatible with a SIR model, we developed the model pictorially illustrated in Fig. 1. As explained in details in the Methods Section, the agents can move on a lattice through jumps processes, modelled using a Lévy walk of jump parameter (beta)36,37,38. When (beta) becomes large, i.e., for (beta rightarrow 2), agents tend to perform a Brownian random walk with very short jumps. As (beta rightarrow 1), agents can travel long distances in just one step. There are no constraints on the number of agents that can occupy a single cell. In each cells, agents can be infected by neighbours according to the SIR rules. Thus, the parameters that control the model are the jump parameter (beta) plus the standard SIR parameters, infection rate (alpha) and removal rate (gamma). The agent-based lattice model considered here reduces to a standard SIR model when the well-mixed population condition is satisfied, i. e. when large jumps dominate the dynamics (Fig. 2).

Figure 1

Agent-based SIR model on a lattice. (a) Agents of different colors, representing the SIR states, move on a lattice. White cells represent empty sites. Green cells are occupied by susceptible (S) agents, blue cells contain only removed (R) agents. Red cells contain only infected (I) agents. Shaded cells contain agents in a mixture of states. Agents can move among cells performing jumps (black arrows) whose length follows Lévy statistics. The letters i and j, with (i=1,..,N_b) and (j=1,…,N_b) define the location of the cell (ij). (b,c) Agents in the same cell undergo a SIR dynamics: (b) S become I at a rate (alpha); (c) I become R at rate (gamma). (d) The jump dynamics allows an agent to move from the cell (ij) to ((i+k,j+l)). The probability to perform a large/small jump is controlled by the parameter (beta in [1.0,1.99]). Large (beta) values correspond to small jumps, i. e., a random walk that gives rise to Brownian motion. Small (beta) values correspond to large jumps.

Full size image

For reproducing the kinetics of real data we made the following assumptions:

  • In the absence of containing strategies, the infection is characterized by a high infection rate (we take (alpha =0.9)) and a low removal rate ((gamma =0.025) or 0.05). Using as a unit of time the update of all agent positions (see Methods for details), the removal rate introduce a time scale (tau _I = gamma ^{-1}=40) or (20). This characteristic time scale represents the average time an agent remains infected and can thus spread the infection. This condition ensures that we are in an epidemic regime, i. e., the mean-field value is (R_t gg 1). We stress that, since the SIR dynamics with only three sub-populations is a simplification of the real chain of epidemic transmission, the parameters we choose for the epidemic spreading are not strictly related to those of Covid-19. Because we are interested in the effect of mobility restriction on epidemic spreading, we fix the epidemic parameters in a way that, without mobility restrictions, we are sure to stay in the worst-case scenario with an exponentially fast spreading of the infection.

  • The parameter (beta in [1,1.99]) tunes the intensity of mobility restrictions. The higher its value, the stricter the limitations. (beta) is one of the fitting parameters.

  • Other interventions that mitigate the epidemic spreading tend to increase the removal rate (gamma). We thus assume that (gamma) is another fitting parameter. This is because typical measures, for instance, quarantine, remove infected agents from the system. In this way, we reabsorb the presence of many hidden sub-populations into an effective value of (gamma).

  • We define the parameter (delta), i. e., the fraction of infected agents at the epidemic peak with respect to the entire population, that provides a quantitative measure of the reduction of the epidemic peak. In other words, the parameter (delta) represents the efficiency of a given containing strategy compared to the uncontrolled situation where all the agents turn out to contract the infection (which is the case of our model for (gamma ll alpha), (alpha =0.9), and (beta =1)).

To detail how mobility restrictions induce deviations from the SIR model, we calculate, via numerical simulations, the epidemic curves as a function of time for different values of (beta) as illustrated in Fig. 2a. Here, the SIR parameters are (alpha =0.9) and (gamma =0.025), i. e., the corresponding SIR model is in the fully blown epidemic regime. For small (beta) the epidemic growth is well captured by the exponential function, indicating that we are in the epidemic regime. As (beta) increases the curve turns out to be flattened and the peak reduces to (80%). Moreover, the growth of the epidemic for the largest (beta) examined is well described by the power law (I(t) sim t^{2}). The value of the exponent is comparable with those measured in different countries during the COVID(-19) epidemic wave23. The model considered here suggests that the crossover from exponential growth to power-law might be related to changes of the mobility patterns that, in our picture, shift from being dominated by large jumps to small ones. This finding is consistent with the observation that a sub-exponential growth in the number of infected people is a consequence of containing strategies23. Moreover, in the microscopic description adopted here, the crossover in the kinetics of I(t) is driven by just one parameter.

Figure 2

Agent dynamics impacts the epidemic spreading process. (a) The graph shows the dependency of the epidemic curves on (beta =1.20,1.50,1.75,1.80,1.85,1.87,1.90,1.92,1.95,1.97,1.99) (increasing values of (beta) from yellow to violet). As (beta) decreases, the epidemic grows exponentially fast (dotted black curve) and approaches the evolution of SIR model in well-mixed population (dashed red curve). The dash-dot blue curve is a power law (sim t^2). The parameters of the SIR reactions are (alpha =0.9) and (gamma =0.025). (bg) Typical configurations taken at the same fraction of infected agents (I/N sim 0.25) for increasing values of (beta =1.0,1.2,1.4,1.6,1.8,1.9) (red are infected sites, green the susceptible ones, we keep white the sites populated by removed agents). (h) The probability distribution function of the local density of infected sites. (i) Radius of the cluster of infected agents ((beta =1.99)) as a function of time. The red dashed line is a linear fit.

Full size image

The crossover from exponential to power-law growth reflects the drastic change in the structure of clusters of infected agents, as illustrated in Fig. 2b–g, where typical configurations with the same fraction of infected agents are shown ((I/N=0.25, alpha =0.9, gamma =0.025)). As one can see, in the high mobility region ((beta = 1)), infected agents are spread almost everywhere in the system. As (beta) increases, infected sites tend to form a single cluster. This phenomenology is consistent with the literature of mobile agents undergoing SIR dynamics39,40. This structural change is quantitatively documented by the density distribution of infected sites shown in panel (h) of the same figure (see section Methods for details). As one can appreciate, the distribution becomes double-peaked as (beta) increases. The first peak around zero indicates the presence of an extended region of susceptible agents. The peak at high values is due to the growing cluster of infected agents. As highlighted in panel (i), the cluster grows linearly in time and thus the number of infected grows with (t^2).

Another interesting aspect to understand with this model is the trade off between mobility restrictions and and other kind of interventions that have the effect of increasing the removal rate. In particular in Asian countries41, NPIs applied during the COVID-19 waves have relied mostly on contact tracing and/or preventive quarantine, with little mobility reduction, leading to effective and durable control of epidemic spreading, as reviewed by Ref.21. To understand if there is an optimal balance between containing strategies (characterized by (beta)) and efficiency in removing infected agents (denoted by (gamma)), we calculate the fraction of infected population at the epidemic peak (the maximum of I(t)) as a function of the jump parameter (beta) and of the removal rate (gamma). As above, the initial occupation number of each site is, on average, one. The infection rate is (alpha =0.9). The resulting phase diagram is shown in Fig. 3. The color indicates the fraction of infected population: in the violet region, this fraction goes to zero (epidemic is suppressed) while in the yellow region such a value goes to one, indicating an epidemic regime. The phase diagram fully recapitulates the effectiveness of the two strategies used to mitigate the infection spread, a strong lockdown with limited contact tracing, or an efficient contact tracing a moderate reduction of the mobility.

Figure 3

Effect of different containment strategies. The phase diagram is obtained considering as control parameters (beta), that represents mobility restrictions, and (gamma), the efficiency in removing infected agents. The color scale represents the fraction of the initial susceptible population that becomes infected, ranging between 0 (epidemic suppression, violet region) and 1 (fully-blown epidemic, yellow region). Containment is achieved as (beta) increases (corresponding to increasing mobility restrictions) even with low removal rate, or increasing (gamma) (effective removal of infected agents), even with limited mobility restrictions.

Full size image

However, even under the strictest lockdown, several activities could not be stopped (hospitals, food supply chain, …), meaning that a single mobility parameter cannot fully describe this varied situation. To understand what could be the impact of heterogeneous motility patterns on the evolution of the epidemic, we introduce in the model some regions characterized by a high mobility (jump parameter, (beta _2)), while the majority of the the cells have restricted mobility, with a jump parameter (beta _1=1.99) (see Methods for more details). By varying (beta _2) and the density of more mobile cells (parameter (rho)) we are able to draw the phase diagram shown in Fig. 4.

Figure 4

Sites of different mobility affect epidemic spreading. (a) Each cell labelled by (ij) is characterized by its own mobility parameter (beta _{ij}). We consider the special case of a binary mixture ((beta _{ij} = beta _{1,2})) of high and low mobility regions. Changing the density (rho) of (beta _2) sites and the value of (beta _2), we obtain the the phase diagram presented in panel (b), obtained for (beta _1=1.99), (alpha =0.9), and (gamma =0.05), conditions that grant contained epidemic spreading thanks to the low-mobility group. A small amount of sites with small values of (beta _2) can trigger the epidemic spreading.

Full size image

As in the previous case, in the violet area the epidemic spreading is stopped, while in the yellow area the epidemic peak reaches the entire population. Epidemic spreading takes place above a critical curve: for a given value of mobility (beta _2<beta _1), the system can support a maximum fraction of regions with that mobility. Above that fraction, the system falls into an epidemic regime. It turns out that even a small fraction of regions with (beta _2<1.9) triggers the epidemic spreading. Our results confirm that, in the absence of contact tracing able to mitigate the spreading of the infection, only those activities that are strictly necessary should be carried on in order to prevent an enhancement of infections.

The first and second COVID-19 epidemic wave in Italy

As the proposed model is able to describe the transition from an exponential to a power-law epidemic growth, and is capable to provide meaningful predictions on the epidemic trend using only few parameters, it is important to validate it with real data. Therefore, we next tested the model against data from the first and second waves of COVID-19 in Italy. The epidemiological data are very sensitive to the ability of a given country in testing the population, identifying new cases among asymptomatic people, which reflects the capacity of the health system2. During the first wave, testing was restricted by lack of reagents, so the number of positive individuals has been largely underestimated. Daily data for the second wave are much more reliable, thanks to more extensive testing. For this reason, we look at the time evolution of different indexes that are (i) Daily new positive cases (NP), (ii) Positive cases at a given time (P), (iii) Hospitalized cases (H), and (iv) Daily deaths (D). The latter two indicators are independent of testing capacity. The fitting parameters are the mobility parameter (beta), and the removing rate (gamma), while a quantitative measure of the flattening the curve effect, is given by the parameter (delta). The infection rate is kept fixed to the value (alpha = 0.9) for all the simulations. As discussed above, (delta) represents the fraction of agents that do not contract the infection. The mobility parameter (beta) and the removing rate (gamma) can be tuned to reproduce the epidemic spreading. We use official data released daily from the Presidency of the Council of Ministers – Department of Civil Protection42. As shown in Fig. 5 there is a very good agreement between data and model predictions. The model appears to fit particularly well the number of hospitalized cases (H) and the total number of positive cases (P). Both indexes are time-accumulated data, known to lead to an underestimation of uncertainty of fits43. However, daily number of new cases (NP) and deaths (D) are also well fit, and the four indexes considered provide estimates for the fitting parameters (beta), (gamma) and (delta) that are consistent with each other, supporting the validity of the model. Thus, the observed evolution of the indexes is well captured by the model with (beta = 1.95), meaning strong mobility restrictions, and (gamma sim 0.025). The resulting peak reduction ((delta)) is around (50 %).

Figure 5

Model description of the first epidemic wave. (a) Time evolution of the epidemic curves in Italy during the first 150 days for Daily positive cases (NP), Positive cases (P), Hospitalized cases (H) and Deaths (D), as indicated. Black symbols are data, red curves are the fits to the model. Data (from42) are normalized at the peak value. Time is measured in days from the day zero, February 4th, 2020. (b) Mobility parameter (beta), removal rate (gamma), infection rate (alpha) and the corresponding peak reduction (delta) are indicated for the four indexes considered in (a). The shaded area represents the standard deviation, dashed red line is the average.

Full size image

For a further validation of our model, we perform a comparison between the model and the data of the second epidemic wave. As a starting date, we have taken September 20th, and data are updated to December 28th, 2020. For fitting, we have subtracted the baseline and normalized values to the epidemic peak. Results are shown in Fig. 6. As for the first wave, fit to daily (NP and D) and cumulative data (H and P) yield comparable values of the three parameters, (beta), (gamma) and (delta). The second wave is described by a lower value of (beta = 1.55), as a consequence of less severe restrictions, which were less effective, as shown by a smaller reduction of the epidemic peak, described by (delta sim 72.5 %). Note that, also for the fit of the second wave, the infection rate has been kept fixed to (alpha =0.9) while the values of (gamma) found in the two waves are almost the same, indicating that the main difference in the growth of the epidemic can be attributed to a different mobility during the two waves. The values of (gamma) and (alpha) we obtain from the fit are almost the same during the two waves. This fact indicates that (beta) turns out to be the significant fitting parameter. Although at a very coarse-grained level, in the sense that the parameters of the model are not sensitive to the details of the mobility restriction imposed during the two epidemic waves, the model captures the efficiency of the lock-down on March 2020 and correctly returns a small value of (beta) with a substantial peak reduction. The second wave, characterized by mobility restriction heterogeneous in space and time, i. e., different implementations in different regions, are reflected by a higher value of (beta) that suggests the presence of regions of different mobility. Moreover, the strong relationship between peak reduction and mobility restriction provides clear evidence of the crucial role played by NPI in containing the epidemic.

Figure 6

Model description of the second epidemic wave in Italy. (a) Time evolution during the first 75 days of the second wave of Daily positive cases (NP), Positive cases (P), Hospitalized cases (H) and Deaths (D), as indicated. Black symbols are data, red curves are the fits to the model. Data (from42) are normalized at the peak value. Time is measured in days from the day zero, September 20th, 2020. (b) Mobility parameter (beta), removal rate (gamma), infection rate (alpha) and the corresponding peak reduction (delta) are indicated for the four indexes considered in (a). The shaded area represents the standard deviation, dashed red line is the average.

Full size image


Source: Ecology - nature.com

Scientists build new atlas of ocean’s oxygen-starved waters

Dissecting the dominant hot spring microbial populations based on community-wide sampling at single-cell genomic resolution