Particle velocity controls phase transitions in contagion dynamics.

Interactions often require the proximity between particles. The movement of particles, thus, drives the change of the neighbors which are located in their proximity, leading to a sequence of interactions. In pathogenic contagion, infections occur through proximal interactions, but at the same time, the movement facilitates the co-location of different strains. We analyze how the particle velocity impacts on the phase transitions on the contagion process of both a single infection and two cooperative infections. First, we identify an optimal velocity (close to half of the interaction range normalized by the recovery time) associated with the largest epidemic threshold, such that decreasing the velocity below the optimal value leads to larger outbreaks. Second, in the cooperative case, the system displays a continuous transition for low velocities, which becomes discontinuous for velocities of the order of three times the optimal velocity. Finally, we describe these characteristic regimes and explain the mechanisms driving the dynamics.

Introduction

Spreading processes, such as contagion or rumor transmission, have been modeled using agent-based approaches[1,2]. These dynamics reach a higher level of complexity when, apart from the diffusive nature of the process, the particles move. For example, some oviparous fauna release sperm and eggs in aquatic environments; thus, the fertilization events happen via mobility-driven encounters of these cells[3]. Mobility impacts biological contagion, such as the dispersal of fungi spores in plant fields[4], or the appearance of microcolonies through the attachment to the particles being carried inside the xylem for plants[5] or the veins for animals[6]. Pathogens do not spread only inside a single individual body, rather amongst different individuals too, such that proximity allows the contagion from one organism to another; this has led to a diversity of studies of epidemiology in systems of mobile fauna introducing, for example, the appearance of large outbreaks when a pathogen switches from spreading in one species to another[7], or the particularities of waterborne infections in marine environments, such as fishing or death-based transmission[8,9]. Indeed, aquaculture represents a key industry for nutritional and financial security in developing countries, and can be strongly affected by infections spreading among the confined but mobile fishes[10,11]. Additionally, the coupling between movement and spreading dynamics is also relevant in socio-technological systems, where malware is transmitted through the connection of multiple devices to WiFi and Bluetooth networks, spreading further because as time evolves the devices are connected to different networks due to the mobility of their users[12]. Hence, mobility drives the appearance of temporal connections, which modify the statistical properties of the spreading processes[13,14].

Hosts’ movement facilitated the spread of infections that, throughout the human history, have affected large populations[15-17]. In fact, several infections that emerged in geographically distant locations may converge to the same environment due to their hosts’ mobility, facilitating the interaction between different pathogens[18,19], either cooperating or competing[20-22]. These interactions modify the infection rates, such that the infection rate of a disease depends on the history of other infections of the host. This implies that the primary infection rate, which is the infection rate for a host that has never been infected with any disease, is different from the secondary one, which represents the infection rate for agents that have been previously infected with one disease; and these will be different from the tertiary, quaternary, and higher order infection rates. For example, competition between two different infections can lead to cross-immunity, that is, the particles are immunized against other infections after a primary infection[23,24] (i.e., the secondary infection rates are zero). In contrast, cooperation represents the weakening of the individuals after a primary infection, increasing the secondary infection rates[25-29].

Spreading dynamics coupled with particle movement has been typically modelled with metapopulation approaches, which consider the geographical locations as patches of particles connected through fluxes (which define an adjacency matrix). For example, a recent work showed the presence of a region in parameter space where the mobility has a detrimental effect on the outbreak size[30]. However, metapopulation approaches are limited for representing the temporal aspects of the connecting fluxes and the structured interaction inside each patch. Hence, agent-based approaches, where the particles interact with their proximal neighbors, are useful not only for determining the microscopic mechanisms leading to large outbreaks, but also which movement patterns hinder the spreading process. Our approach is based on a temporal network given by proximal interactions, determined by the distance between mobile particles, in the spirit of temporal contact networks, which consist on empirical measurements of who interacts with whom and when. For example, contact tracing has been used for analyzing the structure of the social structure underlying the spread of sexually transmitted diseases[31]. In fact, correlations and the circadian activity patterns lead to discontinuous transitions in the spreading of two cooperative infections[32]. For single infection dynamics on empirical contact networks, the ratio between the time scales of network dynamics and epidemic dynamics influences the epidemic threshold, such that when the network dynamics is faster, the epidemic threshold is smaller[33].

In a scenario where particles move, we study the contagion dynamics of either a single or two cooperative infections arguing about the role of particle velocity as a control parameter. First, we introduce the model describing the contagion process and the particle movement. Second, we analyze the behavior of a single infection spreading among a population of mobile particles. Finally, we extend our analysis to the dynamics of two cooperative infections.

Model

As initial condition, we distribute N particles randomly in a two-dimensional space of size L × L with periodic boundary conditions. Each particle i is assigned a direction of movement randomly from a uniform distribution. This direction of movement does not change with time. The movement is uniform rectilinear with velocity v, such that the time evolution of the coordinates of particle i, x(t) and y(t), is given bywhere the motion direction is represented by the angles and v is the traveled distance per unit of time. The Susceptible-Infected-Recovered (SIR) model describes the dynamics[34] between three states: an infected particle transmits the infection with probability p to each particle located at a distance smaller than d, and recovers after one time step. The time scale of the dynamics is set by the time it takes one infected particle to recover, which we define as a time unit. For two cooperative infections, and , the dynamics is similar to the SIR model[25], with primary infections from an infected neighbor happening with probability p, while secondary infections occur with probability q.

For a single infection, the initially infected population does not transmit the infection for probabilities below a critical value, . However, for , the infection affects a finite fraction of the system[35]. For two cooperative infections, the case represents two independent contagion processes, while for both infections have a cooperative interaction. Particularly, focusing on an infected particle EGO, each neighbor particle i updates its state according to the following rules (Fig. 1): a) if i is in state S, it will be primarily infected with probability p; b) if i is infected or recovered with/from the same infection as EGO, nothing will happen; c) if i is infected/recovered with/from the other infection, it will be secondarily infected with probability q. In our simulations, we initially assign the state AB (I for single infection) to a randomly chosen particle and the state S to the rest. At each time step, first, the state of all the particles is updated (Fig. 1); and then, their positions are updated synchronously according to Eq. 1.

Figure 1

Scheme of the model. (a) Contagion of infection from EGO to its neighbors (particles at a distance smaller than d). Labels on the links indicate the infection probabilities. (b) Dynamics of primary and secondary infections, happening with probabilities p and q, respectively.

At each time step, the set of interactions is described as a random geometric graph[36], because the initial positions and the movement directions are uniformly distributed. In that topology, the expected number of neighbors is

We identify two referential spatial scales in random geometric graphs. First, random geometric graphs with periodic boundary conditions display a percolation transition[37] in 2D at . Thus, fixing and (same particle density as in ref.[38]), we obtain . Second, the system is fully connected for .

The temporal network describing the contacts in our system displays temporal correlations. In fact, temporal correlations are related to the probability that, if there is a contact between particles at time t, this contact remains at time . Considering the particle velocity as a control parameter, the temporal correlations have two limits. One limit is the static case, where the particles do not move, such that the set of contacts does not change in time (see Supplementary Fig. S1), and there will be finite outbreaks only if (necessary but not sufficient condition). In the other limit, particles move so fast such that the set of contacts between two consecutive updates are uncorrelated, and can be represented by generating a random location for each particle at every update. Thus, the probability to observe the same interaction or new interactions in two consecutive time steps can be used to link the velocity (one of our control parameters) to empirical temporal networks, such as those specifying human contacts measured through radio frequency identification devices[39,40].

Recapitulating, our model has three control parameters (d, v, p). We will keep one of these parameters constant for studying the behavior of the order parameter, which will be the density of (doubly) recovered particles () in the final absorbing configuration for single (two cooperative) infection(s).

Results

Single infection

We start studying the case of maximum infection probability () in the space. For a given particle velocity v, there is a threshold interaction range d* (v) such that the infection spreads to a finite fraction of the population for (Fig. 2a). For very high particle velocities, this threshold approaches an asymptotic value d (circle in Fig. 2a) while, as the velocity is decreased, d* increases, reaching a maximum for and then decreasing to (triangle in Fig. 2a). Thus, we identify four different regimes when, for a fixed d, the particle velocity is varied: i) : the infection does not spread, independently of the velocity; ii) : there is a threshold velocity v, such that for the infection affects a finite fraction of the population; iii) : the infection spreads in the static case (), but the low velocities are detrimental for the contagion process, implying that increasing v leads to smaller . In this regime, for a given d, there are two critical velocities, and , such that there are no macroscopic outbreaks for ; iv) : the outbreak reaches all the system, independently of v.

Figure 2

Single infection: fraction of recovered particles in the final absorbing configuration for: (a) maximum infection probability () in the (d, v) space; (b) , in the (p, v) space. Symbols depict the cases a, dmin (circle) and d (triangle); b, (square).

In fact, there are two cases that help us understand the limit cases of d*:

In the static limit (), particles do not move. Hence, the infection spreads along a random geometric graph and, considering , it affects a finite fraction of the population only for interaction ranges in the supercritical percolation regime, (triangle in Fig. 2a).

The uncorrelated contact sequence limit () represents a sequence of spatial configurations in which the positions of the particles are randomly reassigned after each time step. This implies that, for a high enough number of particles and , there are no common contacts in two consecutive time steps. For an infection probability p, the dynamics in this limit, known as annealed network, is described by the mean-field approximation for SIR dynamics[35], where follows the implicit equation . This equation shows the dependence of on , and has a non-trivial solution for . Hence, the epidemic threshold is . In fact, considering this limit for different interaction ranges, has a universal behavior for all d when p is rescaled by (see Supplementary Fig. S2). From the mean-field solution, and considering that p is a probability (), we deduce that there is a unique solution (the trivial solution, ) for . Thus, we estimate dmin substituting in Eq. 2, where we obtain (circle in Fig. 2a).

Finite particle velocity: non-monotonic behavior

Previously, we have described a set of interaction ranges () for which there is a non-monotonic behavior of the outbreak size with the velocity for (Fig. 2a). Now, we want to study the behavior with velocity for and for a fixed interaction range . Considering the uncorrelated contact sequence limit as a point of reference, where our order parameter depended only on , we normalise the infection probability with the average degree  , such that the case in Fig. 2a () corresponds to the case in Fig. 2b. We observe that, for a given velocity, there is an epidemic threshold p(v), that is, there are macroscopic outbreaks for , and the epidemic threshold has also a non-monotonic behavior with v, approaching the limit as (Fig. 2b).

For simplicity, we explain the observed non-monotonic behavior for the case (Fig. 2a). For , the infection spreads following an expanding front (see Supplementary Fig. S3); however, as the velocity is increased, the infected particles are further from those that infected them, such that for high velocities there is not a well-defined expanding front and the dynamics is described by a mean-field time evolution (global mixing). We quantify the presence or absence of an expanding front measuring the average connected component size of recovered particles at each time step (i.e., considering only the recovered particles and the contacts between them) as a function of the number of recovered particles R. If the dynamics is driven by an expanding front, , while if the spreading is driven by global mixing for low R, and for higher values of R there is a percolation transition at (Fig. 3a). Hence, as , the contagion among the population is driven by an expanding front in the static case. Moreover, the spatial configurations for and are ordered, as SIR dynamics belongs to the universality class of directed percolation, implying that there are not any contacts between recovered and susceptible particles, and then the system is composed by coherent regions of recovered (the core of the expanding front), infected (the border of the front) and susceptible particles (those out of the front). Hence, the contagion process in the static case represents an expanding front with ordered contacts.

Figure 3

Single infection: mechanism for a non-monotonic outbreak size with particle velocity. (a) Average connected component size of the network composed by the recovered particles and the contacts between them, as a function of the number of recovered particles R. The black vertical line represents the critical value R for percolation in a random geometric graph. (b) Time evolution of the ratio between the number of IS links and the sum of IS and IR links. Inset: the three possible scenarios for a mobile system in 1D, with susceptible and infected particles depicted, respectively, in white and red colors. Velocity is changed from (blue) to (cyan) with , and from to (red) with , for and .

However, for low velocities, the contagion among the population is still described by an expanding front (Fig. 3a), but the movement of the particles leads to local mixing at the edge of the front. For example, there are contacts between recovered and susceptible particles, which did not occur in the static case. As the infection is transmitted through the contacts between infected and susceptible particles (IS links), we compare its abundance with the number of contacts between infected and recovered particles (IR links). After an initial transient, the number of IS and the number of IR links are balanced in the static case. However, when velocity is slightly increased, the IS links are less frequent (Fig. 3b). Considering that the average degree is fixed by Eq. 2, this change in the relative abundance implies that each infected particle will have a smaller amount of contacts with susceptible particles. Indeed, as , this decrease represents a dynamical behavior similar to a lower (i.e., ), such that the observed discontinuity may represent a dynamical percolation transition.

The decrease in the relative abundance of links between infected and susceptible particles (see Supplementary Fig. S4) is illustrated with a similar set-up in a 1D configuration. In this case, at the edge of an expanding front, there are three possibilities (Fig. 3b inset): i) the susceptible and the infected particles are approaching: the susceptible particle will get infected, but at next time step it will be surrounded by many recovered particles; ii) they are getting away: if initially they are separated by d0, at the next time step their distance will be such that, if , (optimal velocity); iii) they move in the same direction, leading to the same microscopic behavior as in the case of static particles. In conclusion, this non-monotonic behavior is explained through the cases i) and ii) that, for an expanding front, imply a detrimental effect of low velocities on the contagion.

Finally, for high velocities, the system approaches asymptotically the uncorrelated contact sequence limit. This means that an infected particle at time t is not in contact with the particle that infected it (at time t − 1), such that the dynamics cannot be described by an expanding front. As a consequence, the initial transient lasts longer (high abundance of IS links, Fig. 3b) and, as , most of the particles in the system are infected during this transient, leading to in the final absorbing configuration.

Two cooperative infections

In this section, we consider the contagion of two cooperative infections which, unless otherwise stated, have a maximum cooperative interaction ().

Focusing on the static particles case (), there are no finite outbreaks for . For connected systems with short-range interactions (), the epidemic threshold p is non-linear with (see Supplementary Fig. S5). While for the transition is continuous (see Supplementary Fig. S6), long-range interactions appear as d is increased, facilitating the presence of discontinuous transitions[41], with the critical point shifting from to (see Supplementary Fig. S5). However, in the case , p has a linear growth with the inverse average degree[42] , and the transition is continuous (see Supplementary Fig. S7), due to the immediate secondary infection of the particles after a primary infection (); hence, the dynamics is driven by primary infections, which lead to continuous transitions. However, if we decrease q for high , p grows, shifting from 0.5 to 1, and having a discontinuous transition (see Supplementary Figs S8 and S9).

In the uncorrelated contact sequence limit (), we find a universal shape for all d, when p and q are normalized with (Fig. 4). For non-interacting infections (), there is a continuous transition at . However, as is increased, a discontinuous transition appears, with a characteristic gap that grows for higher values of q (see Supplementary Figs S10 and S11). This behavior is qualitatively similar to that already reported for Erdös-Rényi networks[26], with the same p for low , but leading to a change in the nature of the transition when q is varied. In fact, for , the behavior displayed by these annealed networks is understood taking into account the results on mean-field (see Supplementary Fig. S11)[25]. Specifically, for , the system displays a continuous transition at ; in fact, is not higher enough than p for observing a discontinuous transition. As d is increased (), a discontinuous transition appears (see Supplementary Fig. S12) and shifts towards lower values. Finally, for , , exhibiting a continuous transition, in agreement with the results found for the static case with , as the random redistribution of the particles, with , does not imply large changes on the topology.

Figure 4

Cooperative contagion with interactions defined by uncorrelated contact sequences. The system has a phase transition at , which becomes discontinuous as the secondary infection probability q is increased. Symbols denote different values of q: (squares), (diamonds) and (circles). Colors represent different interaction radii d, changing gradually from (blue) to (red) with .

Finite particle velocity

We vary v to describe how the differences between the limits of static particles and the uncorrelated contact sequence arise. We focus on the interaction range with , which has a discontinuous transition for (see Supplementary Fig. S12), while the static case leads to a continuous transition (see Supplementary Fig. S6), and corresponds to the regime associated with the non-monotonic behavior of with v for a single infection (Fig. 2b).

There is a continuous transition for low velocities, with p following a non-monotonic behavior with v (Fig. 5a,b). However, in the region associated with global mixing effects, the transition point p does not vary too much, but the lines of constant are shifting to lower values of p as v is increased, leading to the appearance of a discontinuous transition with a gap which grows with the velocity (Fig. 5a,c).

Figure 5

Cooperative contagion on mobile particles with velocity v for and . (a) The fraction of doubly recovered particles in the final absorbing configuration has a first order phase transition for high velocities v. The dashed red line indicates the discontinuous transition point, while the square symbol represents . The nature of the transition is confirmed by the time evolution of the number of infected particles, which is (b) continuous for and (c) discontinuous for . The plotted curves correspond to (brown), 1.6 (dark violet), with .

The appearance of a discontinuous transition is explained taking into account previously reported mechanisms for this dynamics[26]. The gap in at p represents the occurrence or not of avalanches of secondary infections in our system. The avalanches occur when the two infections and meet after following different relatively long paths in the network. As , when this meeting event occurs, an avalanche of secondary contagion events spreads over the paths that the infections followed independently previously. The requirements for this phenomenon are that the contact sequence describing the interactions has loops, such that this will not happen in static trees, and that the infections do not meet after short paths, implying that network has a low level of local clustering. Specifically, in static networks, this dynamics leads to continuous transitions in low-dimensional lattices (1D and 2D), and to discontinuous transitions for higher dimensions (4D lattices and Erdös-Rényi networks)[26,41].

In our system, the mobility allows that the infections spread to a higher (and located further) fraction of particles and also decreases the effect of the local clustering, as the movement hinders the meeting events after a short path. Then, the two infections follow different and long paths before meeting: when , if the two infections meet, there will be a high fraction of doubly recovered particles in the final absorbing configuration; otherwise, the fraction of singly recovered particles will belong to a continuously growing branch, but the fraction of doubly recovered particles will tend to zero. Hence, there will be a discontinuous transition with two branches, with the probability of being in the upper branch representing the probability of the two infections meeting after a long path of independent contagion events.

Discussion

The dynamics of a single infection illustrated how the epidemic threshold can be controlled with the particle velocity, leading to a non-monotonic behavior with the particle velocity. For low velocities and , the dynamics lead to higher epidemic thresholds than in the static case. Specifically, we have described the mechanism leading to a dynamical fragmentation, which has already been reported in the context of disease spreading in temporal networks[32], but also in other dynamics such as the coevolving voter model[43]; in contrast to the later, where the evolution of the topology depends on the dynamical configuration, our approach leads to a dynamical fragmentation even when the particles’ movement is independent of their states. In the infinite velocity limit, we found a scaling relationship between the epidemic threshold and the topological parameters (for , ). These results suggest that blood or sap velocity inside different species may have evolved in order to minimize the risk of the microcolonizations by pathogens in their environment.

For the case of two cooperative infections, not only the epidemic threshold but also the order of the phase transition varied with the particle velocity. Previous work studying this dynamics reported that the nature of the phase transitions was influenced both by the topology and the level of cooperation[26,41]. Specifically, intermediate levels of cooperation in mean-field approximations were leading to abrupt transitions[25]. However, although a qualitatively similar behavior was observed for two-dimensional lattices with long-range interactions, the phase transitions were continuous with short-range interactions. We analyzed the interplay between this dynamics and the mobility focusing on the case . While in the case of low velocities there was a continuous transition, the global mixing effects dominated for higher velocities, making the system evolve towards discontinuous transitions. This behavior may arise from the dynamical fragmentation of the short-term loops due to the high particle velocity: in contrast to the static case, where the short loops are relatively abundant, the dynamical short-term loops do not appear in the high-velocity regime. This leads to the nucleation, which happens only when the population of singly recovered (a, b) has grown enough such that the proximal interaction between one infected A (B) and one recovered b (a) is more likely, leading to a cascade of secondary infection events.

We anticipate that the mechanism leading to a non-monotonic behavior with the particle velocity, arising due to the differences between front dynamics and global mixing, may appear for different dynamics on systems of mobile particles, such as synchronization or evolutionary game theory[38,44-46]. For example, the non-monotonic behavior of the synchronization with the velocity of the particles is shown for low connectivity, but it disappears when the number of interacting neighbors is increased[47,48]. Additionally, future research will explore other movement approaches, such as the Vicsek model[49], which leads to collective motion through the coupling of the movement direction between proximal particles, heterogeneity in the particle velocities, or the coupling between the particle velocity and its dynamical state.