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.

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.
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][16][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][21][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][26][27][28][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 www.nature.com/scientificreports www.nature.com/scientificreports/ 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 ξ π ∈ [0, 2 ) i 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 i (t) and y i (t), is given by where the motion direction is represented by the angles ξ i 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, < p p c . However, for > p p c , the infection affects a finite fraction of the system 35 . For two cooperative infections, the case = q p represents two independent contagion processes, while for > q p 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.
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 www.nature.com/scientificreports www.nature.com/scientificreports/ We identify two referential spatial scales in random geometric graphs. First, random geometric graphs with periodic boundary conditions display a percolation transition 37 Thus, fixing = N 2 12 and = L 1280 (same particle density as in ref. 38 . Second, the system is fully con- 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 + t 1. 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 > d d c (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 ρ (ρ ab ) in the final absorbing configuration for single (two cooperative) infection(s).

Results single infection.
We start studying the case of maximum infection probability ( = p 1) in the d v ( , ) 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 > ⁎ d d (Fig. 2a). For very high particle velocities, this threshold approaches an asymptotic value d min (circle in Fig. 2a) while, as the velocity is decreased, d* increases, reaching a maximum for ≈ . v d 0 5 and then decreasing to Fig. 2a). Thus, we identify four different regimes when, for a fixed d, the particle velocity is varied: i) < d d min : the infection does not spread, independently of the velocity; ii) : there is a threshold velocity v c , such that for > v v c the infection affects a finite fraction of the population; iii)  d d c : the infection spreads in the static case ( = v 0), 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, − v c and + v c , such that there are no macroscopic outbreaks for the outbreak reaches all the system, independently of v.
In fact, there are two cases that help us understand the limit cases of d*: • In the static limit ( = v 0), particles do not move. Hence, the infection spreads along a random geometric graph and, considering = p 1, it affects a finite fraction of the population only for interaction ranges in the supercritical percolation regime, > d d c (triangle in Fig. 2a). • The uncorrelated contact sequence limit ( → ∞ v ) 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 〈 〉  k N, there are no common contacts in two consecutive time steps. For an infec- www.nature.com/scientificreports www.nature.com/scientificreports/ tion probability p, the dynamics in this limit, known as annealed network, is described by the mean-field approximation for SIR dynamics 35 p k . This equation shows the dependence of ρ on 〈 〉 p k , and has a non-trivial solution for 〈 〉 > p k 1. Hence, the epidemic threshold is = 〈 〉 p c k 1 . In fact, considering this limit for different interaction ranges, ρ p ( ) has a universal behavior for all d when p is rescaled by 〈 〉 k 1 (see Supplementary Fig. S2). From the mean-field solution, and considering that p is a probability ( ≤ p 1), we deduce that there is a unique solution (the trivial solution, ρ = 0) for 〈 〉 < 〈 〉 = k k 1 min . (circle in Fig. 2a).
Finite particle velocity: non-monotonic behavior. Previously, we have described a set of interaction ranges (  d d c ) for which there is a non-monotonic behavior of the outbreak size with the velocity for = p 1 (Fig. 2a) Fig. 2b. We observe that, for a given velocity, there is an epidemic threshold p c (v), that is, there are macroscopic outbreaks for > p p c , and the epidemic threshold has also a non-monotonic behavior with v, approaching the limit 〈 〉 = p k 1 c as → ∞ v (Fig. 2b). For simplicity, we explain the observed non-monotonic behavior for the case = p 1 (Fig. 2a). For = v 0, 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 〈 〉 C R 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,  (Fig. 3a). Hence, as , the contagion among the population is driven by an expanding front in the static case. Moreover, the spatial configurations for = v 0 and = p 1 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.
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  www.nature.com/scientificreports www.nature.com/scientificreports/ 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  d d c , this decrease represents a dynamical behavior similar to a lower 〈 〉 k (i.e., < d d c ), 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 d 0 , at the next time step their distance will be = + (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 = p 1, most of the particles in the system are infected during this transient, leading to ρ ≈ 1 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 ( = q 1). Focusing on the static particles case ( = v 0), there are no finite outbreaks for < d d c . For connected systems with short-range interactions (  d d c ), the epidemic threshold p c is non-linear with 〈 〉 − k 1 (see Supplementary  Fig. S5). While for  d d c 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 〈 〉 = p k 1 c to 〈 〉 = . p k 0 5 c (see Supplementary Fig. S5). However, in the case  d d c , p c has a linear growth with the inverse average degree 42 〈 〉 − k 1 , and the transition is continuous (see Supplementary Fig. S7), due to the immediate secondary infection of the particles after a primary infection (  q p c ); hence, the dynamics is driven by primary infections, which lead to continuous transitions. However, if we decrease q for high  d d c , p c grows, shifting from 0.5 to 1, and having a discontinuous transition (see Supplementary Figs S8 and S9).
In the uncorrelated contact sequence limit ( → ∞ v ), we find a universal shape for all d, when p and q are normalized with 〈 〉 − k 1 (Fig. 4). For non-interacting infections ( = q p), there is a continuous transition at 〈 〉 = p k 1. However, as 〈 〉 q k 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 c for low 〈 〉 q k , but leading to a change in the nature of the transition when q is varied. In fact, for = q 1, the behavior displayed by these annealed networks is understood taking into account the results on mean-field (see Supplementary Fig. S11) 25 . Specifically, for < < d d d c min , the system displays a  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  = d d 30 c with = q 1, which has a discontinuous transition for → ∞ v (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 c following a non-monotonic behavior with v ( Fig. 5a,b). However, in the region associated with global mixing effects, the transition point p c does not vary too much, but the lines of constant ρ ab 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).
The appearance of a discontinuous transition is explained taking into account previously reported mechanisms for this dynamics 26 . The gap in ρ ab at p c 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 = q 1, 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 = p p c , if the two infections meet, there will be a high fraction of doubly recovered particles ρ ab 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  d d c , 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  d d c . 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][45][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.

Data Availability
The codes used for generating the results are publicly available in https://github.com/jorgeprodriguezg/Particle-velocity-controls-phase-transitions-in-contagion-dynamics.