key: cord-0720423-ga4equ0z
authors: Ventura, Paulo Cesar; Aleta, Alberto; Rodrigues, Francisco A.; Moreno, Yamir
title: Epidemic spreading in populations of mobile agents with adaptive behavioral response
date: 2021-12-15
journal: Chaos, Solitons and Fractals
DOI: 10.1016/j.chaos.2022.111849
sha: 2ec4c8fc50fc761ece8f566f4af9a641433ee21f
doc_id: 720423
cord_uid: ga4equ0z

Despite the advanced stage of epidemic modeling, there is a major demand for methods to incorporate behavioral responses to the spread of a disease such as social distancing and adoption of prevention methods. Mobility plays an important role in epidemic dynamics and is also affected by behavioral changes, but there are many situations in which real mobility data is incomplete or inaccessible. We present a model for epidemic spreading in temporal networks of mobile agents that incorporates local behavioral responses. Susceptible agents are allowed to move towards the opposite direction of infected agents in their neighborhood. We show that this mechanism considerably decreases the stationary prevalence when the spatial density of agents is low. However, for higher densities, the mechanism causes an abrupt phase transition, where a new bistable phase appears. We develop a semi-analytic approach for the case when the mobility is fast compared to the disease dynamics, and use it to argue that the bistability is caused by the emergence of spatial clusters of susceptible agents. Finally, we characterize the temporal networks formed in the fast mobility regime, showing how the degree distributions and other metrics are affected by the behavioral mechanism. Our work incorporates results previously known from adaptive networks into the population of mobile agents, which can be further developed to be used in mobility-driven models.

Epidemics are of great concern to humankind. While it is possible to construct realistic epidemic models with heavy use of data, the effect of human behavior, from individual to collective level, is crucial to the dynamics but difficult to be incorporated into models. It is known that mobility patterns, as an important part of human behavior, strongly influence the spread of a disease, and mobility data from real world is often incorporated to epidemic forecasting.

Despite the increasing availability of data from human mobility [1, 2] , there are several situations in which such data is not available. Therefore, the use of synthetic mobility models to feed epidemic modeling [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] is a promising, yet understudied topic. Note that in this context, mobility refers to short-range displacement of individuals -such as walking -and not to the transfer of individuals from one place to another. Frasca [4] et al. and Buscarino [5] et al. proposed an SIR model in a population of random walking agents that can also perform long-range jumps. They study the relationship between final outbreak size and some average topological properties of the underlying time-aggregated network of contacts, such as degrees, shortest path lengths and clustering. Interestingly, the long-range jumps cause a small-world effect similar to the Watts-Strogatz model [15] , which seems to explain the increase in the outbreak size with the jumping probability. Their work settled the base for epidemic modeling in populations of random walkers.

One of the major difficulties of working with epidemic spreading with random walking agents is that the spatial component imposes a temporal correlation to the contacts, which strongly influences the epidemic spreading. For very slow mobility, the population can be regarded as a static network, whereas for fast mobility the correlations are broken and the system is essentially homogeneously mixed. For this reason, despite some insights from reaction-diffusion processes [16] , the mid-term between these regimes essentially relies on computational Monte Carlo simulations. More recently, other works have considered spatial heterogeneity and separate communities [7] [8] [9] , heterogeneous interaction radii [11] [12] [13] and different agent's velocities [6, 10] .

Despite the inherent difficulties, mobile agent models for epidemic spreading can be very insightful, as some models can reproduce features of real world human and animal mobility. Starnini and others [17] modified the basic random walk to consider heterogeneous activity time and attractiveness of the individuals, successfully reproducing behaviors of some data sets of the SocioPatterns collaboration [2] . Stehlé and others [18] performed SEIR simulations using interactions between attendees of a conference, and showed that a correct estimation of the epidemic dynamics must account for the duration of contacts, a feature that is in parts reproduced by Starnini's model. Another important example is the Levy walk, a modification of the simple random walk that considers heavy-tailed step size distributions. Despite some recent evidences claim for the use of more complex models [19, 20] , Levy walks have been widely used to describe animal and human mobility patterns [21, 22] .

In spite of the importance of mobility models for disease dynamics, little has been done to include behavioral changes into such models. Human responses crucially influence the spread of a disease, and there is an intense research devoted to account for such effect into epidemic models, which include responses based on risk perception [23] [24] [25] , spreading awareness and individual prevention methods [26] [27] [28] , vaccination games [29, 30] and others [31] .

One particularly popular approach is to consider adaptive contacts in networks, which mimic both the social distancing and isolation of infected individuals. Gross and others [32, 33] proposed a simple model in which susceptible individuals can randomly rewire their links that point to infectious nodes, redirecting them to other susceptibles. This generates a dynamic network that evolves simultaneously to the disease spreading, with susceptibles forming highly connected clusters between themselves to exclude infecteds. This behavior produces rich dynamical phenomena such as bistability, oscillations and hysteresis. The model later received refined analytical treatments [34, 35] and motivated a series of subsequent studies [36] [37] [38] [39] . In particular, Zhou and others [38] proposed a variant of the adaptive rewiring model that considers network growth and isolation avoidance, showing that the combination of these two factors can produce multiple epidemic bursts in an SIS model before the disease is eradicated.

In the present work, we propose to merge the richness of adaptive behavioral responses with the modeling potential of mobile agents. We propose a mechanism through which susceptible agents can avoid contact with infected agents by performing preventive moves. As the model itself considers a very simple adaptive mechanism, in which susceptibles are fully and instantaneously aware of the state of their neighbors, the aim of this work is not to provide results that are directly applicable to the real world. Instead, our goal is to propose the merge between mobility models and adaptive responses, study its basic dynamical properties and motivate future developments.

We also use a simple semi-analytical approach for the case in which the disease dynamics is slow compared to the agents' mobility. While it still relies on computational simulations, it allows to easily extract, among other results, the stationary prevalences and the phase diagrams, including the transcritical and saddle-node bifurcations that our model presents. Ultimately, it also provides an interpolated functional form for the reduction of contacts due to the adaptive mechanism. This can be applied into the local/regional dynamics of more complex models, such as those with metapopulations [40] [41] [42] [43] [44] [45] , to include behavioral responses to epidemic spreading.

The epidemic model employed in this work is the reactive Susceptible-Infected-Susceptible (SIS) with discrete time evolution. Each agent can either be susceptible (S) to the disease or infected (I). At each time step, an infected agent that interacts with a susceptible one can transmit the disease with probability β. In addition, an infected agent can also be healed with probability µ. Infection and healing events are only applied in the next time step, so the order of agent visits does not matter, and reinfection after healing in a single time step is not allowed.

For the baseline population dynamics, we use the simple random walk with hard interaction circles, as in references [4] [5] [6] , yet with no long-range jumps. A population of N agents is initially distributed at random in a square space of length L with periodic boundary conditions. At each time step, each agent can perform a random move (RM) of fixed length v and uniformly random direction θ. Thus, the horizontal (x) and vertical (y) coordinates of the position of the agent at time t + 1 after a random move are given by:

Each agent has an interaction radius r, meaning that if two agents have a spatial distance smaller than r, they interact reciprocally. With this interaction scheme, one can construct a snapshot network at each time step of the model.

To avoid any dependency on the implementation, both epidemic and positional state changes are calculated at each time step, but they are only applied after all changes were calculated.

We incorporate an adaptive reaction to the local presence of infected individuals. At each time step, each susceptible (S) individual chooses, with probability p a , to avoid its interaction with an infected neighbor by moving away from it, using the following algorithm: (i) if there is exactly one infected neighbor, the susceptible agent moves in its opposite direction with step size v, which we call a preventive move (PM). Here opposite direction means that the displacement vector makes an angle of π with the relative position vector that goes from the S to the I agent (see more infected neighbors, the susceptible chooses one of them at random and promotes a preventive move away from it. (iii) If there are no infected neighbors, it simply promotes a random move (RM) to a uniformly random direction, as already described. Also, with probability 1 − p a , the susceptible promotes a random move regardless of its infected neighborhood. Figure 1 illustrates the PM and RM actions, as well as the SIS epidemic model, while figure 2 shows the agents' algorithm for choosing between PM and RM. Notice that, also as a simplification, the "awareness radius" of the agents is the same as the interaction radius r. Another simplification is that every infected individual is immediately perceived as such, which can be interpreted, for example, as if the disease symptoms are clear and always display immediately after infection.

We perform Monte Carlo simulations of the SIS model under two different motion schemes for the mobile agents: (i) simple random walk (RW), for which no reaction mechanism is considered (i.e., p a = 0) and (ii) random walk with local reactions (LR), for which we use p a = 1, except where explicitly mentioned. We also compare our model with a homogeneously mixed population (HM) with an average number of contacts per unit time given by:

Which is the expected average degree of a population of agents uniformly distributed in space. We call k H the homogeneous degree. We can also define the homogeneous reproduction number (i.e., the epidemic reproduction number if the nodes were homogeneously mixed with k H contacts per time step in average) as:

where λ = β/µ is the infection-to-healing ratio of probabilities. 

We start the analysis by studying the stationary state prevalence ρ * I of the system under a fixed initial infected fraction ρ I (0). We perform simulations under four different regimes with respect to the density of agents and the relative time scale between epidemic and motion dynamics. Henceforth we call slow epidemics to those in which µ v, while for epidemics with similar time scale between the spreading and motion dynamics we have µ ∼ v.

In figure 3 , we show the stationary prevalence as a function of the homogeneous reproduction number R H , for ρ I (0) = 0.30. In each execution, the model is first simulated for a transient period of t tr time steps, and then for t av additional steps during which the data is collected and averaged. Each point is also an average of a number n ex of such independent executions.

First analyzing the low density regime (panels (a) and (b)), we notice that the local reaction mechanism (LR) considerably raises the epidemic threshold, both with respect to homogeneous mixed (HM) and random walking (RW) populations. In the high density regime (panels (c) and (d)), we notice an abrupt phase transition from healthy to endemic stationary states with the LR population, which is not observed within the other models. Such phenomenon was already reported for SIS-like models on adaptive networks [32, 36, 37] , and is actually a fingerprint of another important feature: a bistable phase caused by a saddle-node bifurcation.

Yet from figure 3, we notice that our LR model deviates from homogeneous mixing even in the condition of slow epidemics (panels (a) and (c)), while the simple random walk (RW) can be reasonably described by HM in this regime.

We evaluate the effect of the probability of avoiding infectious contacts p a in figure 4 , which shows the prevalence curves at the low density and similar time scale regime for different values of p a . For intermediate values of the parameter, the curves monotonically fill the range between full avoiding mechanism (LR p a = 1) and simple random walk (RW, equivalent to LR with p a = 0).

To show that our model presents bistability in the high density regime, we plot, in figure 5 , the time series of the LR model for different initial infected fractions ρ I (0), averaged over n ex = 50 executions each, and observe the basins of attraction. The system can converge to two different stationary states, depending on the initial conditions. For this figure, we use N = 1600 individuals in the population to reduce stochastic effects. Yet, the gray-shaded region represents the approximate location of an unstable fixed point where, due to stochastic fluctuations, each execution of the simulation can take different 

As reported in previous works with adaptive networks [32, 37, 38] , susceptible individuals can form highly connected clusters due to the behavioral response mechanism. We report such feature in our model by noticing that, during transient stages of simulations under high density regime, susceptible agents form spatial clusters that are densely connected due to proximity. for a single execution starting with ρ I (0) = 0.13. As the simulation starts, the degree of susceptible agents (blue curve) quickly raises as the S-clusters are formed, whereas the degree of infected agents drops. For this particular execution, the avoidance mechanism was able to reduce and eventually eliminate the disease, but the initial prevalence falls into the shaded region in figure 5 , meaning that the destiny of the system could have been different.

The formation of S-clusters is related to the spatial gaps that are left by the infected agents performing simple random walks. Susceptible agents tend to move into such gaps and form the observed clusters, which move and change through time.

In other works, clusters of susceptibles were observed both when S individuals actively seek connections to other susceptibles [32, 35] and when they simply avoid infectious contacts [37] . Our model is an example of the latter, reinforcing that the S-clusters can occur without explicit preference to susceptibles in the rewiring mechanism.

The S-clusters are transient. If the prevalence at a given time and the disease transmissibility are not high enough, the lack of infective contacts causes the disease to disappear, along with the S-clusters. However, if the prevalence and/or transmissibility are high enough, the infected agents eventually break the S-clusters and the disease takes over the population, which reaches a steady non-clustered regime. The break of the susceptible clusters is similar to that observed by Zhou and others [38] during the epidemic bursts present in their model.

Due to this ambiguity of trajectories from the clustered regime, it is clear that such clusters are closely related to our model's bistability. In the next section, we develop a semi-analytical approach that corroborates this hypothesis.

For slow epidemic evolution, achieved by sufficiently low values of β and µ, we can not only approximate the discrete time dynamics of the disease to continuous, but also assume that any metrics related to the population of agents is a function of the prevalence, as the population quickly responds to the slow changes in epidemic states. With this in mind, the SIS dynamics is given by the following rate equation for the overall prevalence ρ I :

where l SI = L SI /N is the number of links L SI that connect susceptible to infected agents normalized by the population size N , and is a function of the prevalence ρ I . Effectively, this approximation promotes a time scale decoupling of the epidemic and motion dynamics. We know no analytical method to estimate the functional form of l SI (ρ I ) for mobile agents, but we can di-rectly sample it from Monte Carlo simulations of the population with no epidemic dynamics. That is: each agent receives a given state at the beginning, S or I, and holds it during the whole simulation. This way, we can manually set the number of infected agents N I to achieve the desired value of the prevalence ρ I = N I /N . We calculate l SI at each time step in the stationary state, then average its value over time and over independent executions. This process is repeated for different values of ρ I , enough to have a precise shape of the l SI (ρ I ) curve, which is then interpolated to obtain a continuous approximation. For this work, we apply a spline interpolation of third order. Figure 8 shows the data acquired by this method. In the a) panel, we see how the LR model deviates from the mass-action law used for homogeneous mixing [14] . Once obtained, the function l SI (ρ I ) can be used to determine the epidemic dynamics in the slow regime. For instance, the fixed points are the values ρ * I for whicḣ ρ I = 0 in equation 4. This is equivalent to solving the equation:

where λ = β/µ is the infection-to-healing rates ratio. In figure 8 .b), the fixed points can be found as the crossings between λl SI (ρ I ) and the identity line (gray dashed line). Also according to equation 4, the stability of the solution is given by the slope of λl SI (ρ I ) at the fixed point: if it is greater than 1 (i.e., the slope of the identity line), the solution is unstable, and, if it is less than 1, it is stable. Notice from figure 8 .b) that each curve represents a different phase of the LR model. For λ = 0.1, the only solution is the trivial ρ I = 0, i.e., the healthy state. For λ = 0.2, the system presents two more solutions, of which one is stable (the one with greater ρ I ), while the healthy state remains stable. This characterizes the bistable phase, obtained after a saddle node bifurcation. Finally, for λ = 0.4, the healthy solution is unstable (as the initial slope of λl SI (ρ I ) is greater than 1), characterizing the regular epidemic phase. This phase is reached after a transcritical bifurcation.

To compare the semi-analytical approach with Monte Carlo simulations of the epidemic dynamics, we plot the fixed points (both stable and unstable) obtained from the l SI (ρ I ) curves as a function of the disease transmission rate β, rescaled as the homogeneous reproduction number R H = βk H /µ. The results are in figure 9 , where fixed points of the semi-analytical approach are represented as solid (stable) and dashed (unstable) lines, and the stationary prevalences obtained from Monte Carlo simulations are represented as squares (starting from ρ I (0) = 0.30) and circles (starting from ρ I (0) = 0.01). The region where both the healthy and endemic solutions are stable (according to the semi-analytical approach) is shaded in the plot. Figure 9 shows the good agreement between the semianalytical and Monte Carlo formulations. Notice that, in the bistable region, the stationary prevalence of the Monte Carlo simulations depends on the initial conditions, although there is some disagreement at the edges of the bistable region, which may be attributed to stochasticity and population's finite size.

Using the sampled and interpolated approximation for l SI (ρ I ), we can numerically calculate the critical points. In figure 10 , we show two phase diagrams: one for the agents' step size v (upper panel) and another for the local reaction parameter p a , both as a function of the homo-geneous reproduction number R H . As the velocity v increases, the bistable region shrinks and disappears, as the fast motion prevents the agents from forming S-clusters. For greater velocities, the critical value of R H approaches that of a homogeneously mixed population (i.e., R H = 1). From the p a phase diagram (lower panel), we infer that the size of the bistable region grows with the parameter p a that controls the intensity of the local reaction, as expected.

We can further study the structure of the networks that are formed by this model in the clustered regime by looking at the degree distributions. To simplify the execution, we also remove the epidemic dynamics for this measurement, so the S and I agent states are static. We consider the connectivity between different classes of agents, so for example k SS is the number of links of a susceptible agent that point to other susceptibles, k SI is the number of links of a susceptible that point to infected agents, and so on. Also k S and k I are the total degrees of susceptibles and infecteds, respectively. Figure 11 shows each of these degree distributions, grouped by the state of the link targets. As a reference, we also plot a Poisson distribution (gray x symbols) f (k; s) = e −s s k /k! with s = k H (1 − ρ I ) (a), s = k H ρ I (b) and s = k H (c). These are, in each case, the expected degree distributions if the agents were homogeneously distributed at random in the space. Each vertical line also shows the average of each distribution of the corresponding color.

From figure 11 .a), we can see how susceptibles are highly connected to each other, but weakly connected to infected agents, as the k SS distribution spans over higher values than the k IS . In panel b), we see that k II is well described by the equivalent Poisson distribution because infected agents perform simple random walks, whereas k SI is slightly reduced due to the local reaction mechanism. Finally, panel c) shows how the overall degree distributions are distorted and broadened from the basic Poisson curve, similar to the effect observed by Gross and others [32] in adaptive networks. It also shows that, on average, susceptible agents are more connected than infected ones.

Yet under the same framework that considers static disease states, we can analyze how the average degrees, as well as some other network metrics, vary with the prevalence, in order to determine where the S-clustering regime occurs. Figure 12a shows the average degrees of susceptible ( k S ), infected ( k I ) and ( k ) agents of the snapshot network, normalized by the homogeneous degree k H , as a function of the prevalence. The peak on the k S curve is a consequence of the S-clusters, and it is also visible in the total average degree k . The average degree of infected agents k I , on the other hand, is always smaller than k H , as a consequence of the local reaction mechanism itself, which reduces the contacts with infected agents.

In figure 12 .b), we show the average clustering coefficient C and the degree assortativity r deg of the snapshot networks also as a function of the prevalence. The value a) b) c) FIG. 11 . Degree distributions for snapshots of the dynamic network, averaged over time steps and executions, using simulations of the local reaction (LR) model with static disease states. We group the degrees according to the state of target nodes: susceptibles (a), infecteds (b) and both (c). Each dashed line is the average of the distribution with the same color. As a null model, we show a Poisson distribution ("x" symbols) that represents homogeneously random interactions in each situation, as explained in text. Other parameters are: N = 400, r = 1, v = 0.3, pa = 1, kH = 12.57 and a fixed prevalence ρI = 0.2.

of C, which is naturally high on dense random geometric networks, is enhanced at the range in which the S-clusters occur, having a good correlation with the k S curve. The assortativity r deg is also enhanced by the avoidance mechanism, displaying a broader peak, which means that the effect of the mechanism into the assortativity prevails even when the S-clusters are not very expressive. As also reported in [32] , the increase in degree correlations may be an effect of the segregation between S and I agents, as reported in the degree distributions in figure 11 .

We finally notice that the range of unstable prevalences in figure 8 .b) is compatible with the region at which k S and C have their peaks in the plots of figure 12 meaning that, as observed, the S-clusters are unstable when the epidemic dynamics is active, causing the disease to either spread globally by breaking the S-clusters or be eradicated, as explained in section III. This reinforces the relationship between the S-clusters and the model's bistability.

We propose a simple mechanism to include a form of behavioral response to epidemics in mobile agents, based on the avoidance of contacts with infective individuals. We show that, with such mechanism, we can merge the rich dynamical features of adaptive contacts, initially studied on networks, with the overlooked potential of mobile agents for epidemic modeling. Although this work is focused on the high density and slow epidemics regime, we show that different outcomes can be obtained in each regime.

For the low density regime, which is often used to reproduce empirical data [17] , the local reaction mechanism considerably suppresses the infectious contacts and thus the stationary prevalence, besides increasing the epidemic threshold. This is because, when the agents are spatially sparse, the susceptibles easily find the direction to avoid infected agents. These results, however, might be sensibly affected by considering a "watch radius" different from the disease transmission radius, which is proposed as a future work.

In the high density regime, the stationary prevalence is not strongly reduced from the homogeneous mixing scenario, but new dynamical features are introduced: the bistability and the clusters of susceptibles. In the bistable regime, the transient S-clusters can either succeed to eradicate the disease or permit it to spread, depending on the initial prevalence, disease infectiousness and stochastic factors. The bistability is inherited from the adaptive contact changes, but the collective motion of susceptibles between the gaps left by infecteds is an interesting feature that is exclusive to our spatial model.

We also apply a simple semi-analytic approach to describe the dynamic features of the LR model in the slow epidemics regime, based on simulating the population with static epidemic conditions. Due to spatiotemporal correlations of the random walk network, this approach is less powerful than those proposed on adaptive networks [33] [34] [35] , being unable to capture higher order phenomena such as hysteresis and oscillations that are possibly present in our model too. A more powerful analytical approach can be built by using simulation bursts [33, 46] to capture higher order moments. Nevertheless, our simple framework can be generalized to a variety of other dynamic population models, provided that the epidemic dynamics is slow, and a functional form of the behavioral reduction of contacts can be extracted from it.

Finally, we characterize the networks obtained by snapshots of the population with static epidemics in the steady state. With the degree distributions and averages, we show how the local reaction mechanism deviates the behavior from the simple random walk, pointing how the infective contacts are avoided while susceptibles are joint into highly connected clusters. From a practical point of view, this represents a situation in which space is limited and infected individuals do not change their behavior, and is not realistic due to several aspects. However, and as our main goal, we seed the idea of merging adaptive reactions with mobile agent models, calling for further works in the topic. M. also acknowledges partial support from the Government of Aragón, Spain through grant E36-20R, and by MCIN/AEI and FEDER funds (grant PID2020-115800GB-I00). The funders had no role in study design, data collection, and analysis, decision to publish, or preparation of the manuscript.

Communications in Nonlinear Science and Numerical Simulation

Small worlds: the dynamics of networks between order and randomness