Epidemic processes on self-propelled particles: continuum and agent-based modelling

Most spreading processes require spatial proximity between agents. The stationary state of spreading dynamics in a population of mobile agents thus depends on the interplay between the time and length scales involved in the epidemic process and their motion in space. We analyze the steady properties resulting from such interplay in a simple model describing epidemic spreading (modeled as a Susceptible-Infected-Susceptible process) on self-propelled particles (performing Run-and-Tumble motion). Focusing our attention on the diffusive long-time regime, we find that the agents' motion changes qualitatively the nature of the epidemic transition characterized by the emergence of a macroscopic fraction of infected agents. Indeed, the transition becomes of the mean-field type for agents diffusing in one, two and three dimensions, while, in the absence of motion, the epidemic outbreak depends on the dimension of the underlying static network determined by the agents' fixed locations. The insights obtained from a continuum description of the system are validated by numerical simulations of an agent-based model. Our work aims at bridging soft active matter physics and theoretical epidemiology, and may be of interest for researchers in both communities.


I. INTRODUCTION
The study of spreading processes on mobile agents is a field attracting growing interest in both communities of epidemiology and active matter physics. On the one hand, human mobility plays a crucial role in the spreading of infectious diseases, as shown by the inclusion of mobility data into epidemic forecasting [1]. Short range mobility -such as individuals walking in a limited spacehas also been taken into account for epidemic modelling [2], in particular by considering a Susceptible-Infected-Recovered (SIR) model in a population of random walk- * matteopaoluzzi@ub.edu † michele.starnini@gmail.com ers [3,4]. Furthermore, the interplay between mobility and spreading dynamics can be used to model behavior change in individuals [5], as well as to show that a feedback mechanism between the epidemic status and the agent's motion can enhance the contagion dynamics, effectively reducing the epidemic threshold [6].
On the other hand, classical spreading process can model well the diffusion of information in a population: individuals aware of the information (infected) transmit it to unaware (susceptible) peers [7]. Such information exchange is mediated locally by social interactions, involving agents with physical proximity, as observed also in the animal reign [8][9][10][11][12]. As a consequence, populations of motile, self-propelled agents self-organise in time and space, with the emergence of coordination and collec-arXiv:2203.12355v1 [cond-mat.stat-mech] 23 Mar 2022 tive behavioral change [13]. Examples range from multicellular organisms to flocks of birds [14], robot swarms [15], or the coherent motion of fish schools avoiding a predator's attack [16]. Also bacteria, which communicate through chemical signals that regulate their motion, show coordinated behavior of the whole population [17,18], a mechanism know as quorum sensing [19].
From a physics standpoint, systems of self-propelled agents are typically modelled as persistent random walkers [20][21][22]. A salient example is the Run-and-Tumble (RnT) walk, which mimics the motion performed by several flagellated bacteria species such as Escherichia Coli [23]. Collectives of such self-propelled particles have been the focus of intense research in the past decades, providing a natural playground to explore living matter from a physics perspective [24][25][26]. These so-called Active systems, made of biomimetic entities, exhibit a remarkable richness of non-equilibrium collective states as a result of different kinds of interactions [27][28][29][30], which can be of very different nature, say, mechanical, 'social', chemical, etc. Understanding how self-propulsion changes the qualitative features of spreading dynamics and how a local information spreads in a collection of moving entities remains a fundamental open problem in the field. A feedback loop that couples motility with local density can be employed in experiments for controlling colloids [31] or photokinetic bacteria [32,33], and it has been found that Susceptible-Infected-Susceptible (SIS) dynamics can be employed for driving pattern formation and collective motion in systems of active particles [6,34,35]. We consider a d-dimensional system where N particles move in a L d box with periodic boundary conditions.
The position of the particles r i (t) evolves according tȯ with e i (t) indicating the unit vector that specifies the swimming direction at time t, which changes randomly at a rate α I,S (see top panel of Figure 1).
The agents' internal states are then subjected to a SIS process [44], where the transition rates between states are defined as usual The S + I β − → 2I reaction takes place with rate β only in regions of space where particles S and I overlap, that is, exclusively for those pairs of particles located within a distance R, chosen to be small with respect to the size of the system R L (see middle panel of Figure 1). Middle panel: A S particle, in blue, located within the interaction radius R of an I particle, in red, gets infected with rate β. Bottom panel: I particles spontaneously recover (become S) with rate µ.
Conversely, infected particles decay spontaneously to the susceptible state with rate µ (see bottom panel of Figure   1).
Here, we focus on the case β > 0 and µ > 0. For µ = 0 and β > 0, the system evolves towards an absorbing state were all particles are infected in the so-called SI dynamics. It has been shown that, even in this casetrivial from the point of view of spreading dynamics -the interplay between the collision rate and the persistence length triggers fractal growth in the diffusive limit [35].
Non-zero values of µ and large β values trigger a mixed phase, where active particles can develop patterns, while for large values of µ (when the recovery rate is faster than the collision rate), the system can eventually evolve into an absorbing configuration of susceptible particles [34].
In the following, we will analyze the system under both a continuum and agent-based perspective. For the latter, we will perform numerical simulations. At each time step t, the interactions between particles are described by a spatial network where nodes represent particles and links represent interactions between them (thus, two particles are linked if they are within a distance R). As particles move, the set of links evolves. This temporal network can thus be described by a sequence of snapshots, each one with mean connectivity k [45].
We run numerical simulations by implementing the Gillespie algorithm [46]. The time required to update the system is in this case a stochastic variable to be sampled. Once the time needed to perform a given update has been computed, the system is updated. The time associated to a given update is generated from an exponential distribution with characteristic time given by the inverse of the sum of all the transition rates involved in the evolution of our system. Then, one of these transitions, or processes (reactions in the spirit of Gillespie), is chosen with a probability proportional to its rate. Such processes are: • Run: We update the particles' positions at a rate λ m = 5000, much higher than any of the other dynamic processes in our model, such that the particle movement can be considered continuous (to be compared with the rates of the other processes involved : α, β, µ 1). The positions of all particles are updated synchronously.
• Tumble: Particles change their direction of motion at a rate α. All the particles' directions are updated synchronously.
• Infection: The rate of infection of a susceptible par- is the number of neighbouring infected particles at time t.
Therefore, the total updating rate of the system is δ(t) = λ m + α + β i∈S I i (t)+µN ρ(t), where the sum runs over susceptible particles and ρ(t) is the prevalence, that is, the fraction of infected particles in the system, at time t.
We fix N = 1024, and the time and length unit by setting µ = 1 and L = 1. We consider that, for each value of β, 1 Monte Carlo (MC) step corresponds to the time it would take to reach the state where all agents are infected under exponential growth ρ(t) = e β k t , leading to 1 MC = log N β k . We start our simulations from a disordered distribution of agents at high values of β. We then let the system relax 60 MC steps and use the final steady configuration as the initial configuration for a new simulation at a lower value of β. We let the system relax 10 MC steps, and then we repeat this procedure by subsequently reducing the value of β to eventually reach the steady configurations at each value of β.
As a reference, we will consider two limit regimes based on a time-scale separation between motion and spreading: (i) the static regime where particles do not move (or move at much longer time scales than the spreading process), and (ii) the homogeneous-mixing regime, where particles move at much shorter time scales than the ones involved in the spreading. In both limits, the only two control parameters are the infection and recovery rates, as the positions are either not updated or updated at random. In the homogeneous-mixing limit the exact steady density of infected particles (or prevalence) is ρ = 1− 1 β k for β k > 1, and ρ = 0 otherwise, resulting in an epidemic threshold β c k = 1 beyond which ρ > 0.

III. CONTINUUM MODEL
In order to gain insight into the large-scale behavior of the system described so far, we adopt a continuum approach by considering the SIS dynamics on top of the runand-tumble master equation [43]. We introduceS(r, e, t) andĨ(r, e, t) as the probability density function of, respectively, susceptible and infected particles at location r with orientation e at the time t. We are interested in the time-evolution of S(r, t) = deS(r, e, t) and I(r, t) = deĨ(r, e, t). We can associate to S(r, t) and I(r, t) the currents J I,S (r, t) that are defined as J S (r, t) = v S de eS(r, e, t) and J I (r, t) = v I de eĨ(r, e, t). Focusing our attention on the case where the motility parameters α I,S and v I,S might dependent on the internal state of the agents but homogeneous in space, we obtain that the dynamics of the concentration fieldsS andĨ is governed by the following equations where we have introduced the projector operators Q ≡ 1 − P. The operator P acts on a generic function, integrating over all possible directions e of motion [43]. We denote it as follows where f = f (e, r, t). Within this notation, the average prevalence (the concentration of infected particles) in the system at time t is given by ρ(t) = I(r, t)dr.
In this way, we can write the following set of equations for the densities and their currentṡ We now focus our attention on equations for the currents. Without loss of generality, let us discuss the equation for the current of S. As a first choice, we can consider that the active dynamics define the relevant time scale through the tumbling rate α S and thus write Active particles will reach stationarity on time scales t α −1 S . We can thus consider safely a diffusive limit obtained by considering α −1 [47]. Following the same trend of ideas for the current of I, in this limit of vanishing currents, i. e.,J S,I = 0, we get the following constitutive relations In this limit, we are assuming active particles reach a stationary state before the spreading process. In this picture, we obtain that the dynamics of the system is captured by a two-component reaction-diffusion process conserving the total mass ρ(t) = dr [S(r, t) + I(r, t)] [48,49] and is described by the following equations where suitable boundary conditions have to be taken into account.
Another limiting case can be obtained considering that the spreading process is much faster than the RnT motion. In this limit, the spreading dynamics reaches a stationary state before active particles are able to reach the diffusive limit. For studying this limiting situations it turns out convenient to write the equations for the currents in the following way In this case, once we defineD S,I ≡ v 2 S,I /βd, we obtain the following equations As one can see, the equation for I has the form of a backward diffusion equation that tends to make the I profile less smooth as time increases. In the limit β → ∞ or   where the number of active links reaches its maximum.
In the following, we address each case in detail.

A. The static limit
It is worth noting that, in the case of RnT walkers, the static limit can be obtained in two different ways, i. e,.
α → ∞ at fixed v, or v → 0 at fixed α. For D → 0, Eqs. (10) and (11) become ∂ t S(r, t) = −βI(r, t)S(r, t) + µI(r, t) and ∂ t I(r, t) = βI(r, t)S(r, t)−µI(r, t), that have to be solved with the initial conditions S 0 = S(r, 0) and i.e. particles are initially picked from a uniform distribution, the interaction network is a random geometric graph [50]. The latter displays a percolation transition as the interaction radius R is increased, with many microscopic connected components for low R and the emergence of a macroscopic connected component at R = R c . The interaction radius R is related to the average connectivity k (average number of links per node) by where A d (R) is the area of the d-dimensional sphere of radius R, being A 2 (R) = πR 2 , A 3 (R) = 4 3 πR 3 . The connectivity of the graph plays a key role on the emergence of an endemic phase: if the average degree of the graph k is above the percolation threshold k c (d) (which crucially depends on the dimension d of the system), then the graph displays a giant connected component that allows the emergence of macroscopic outbreaks. For k k c (d), the epidemic threshold will be the same as in the homogeneous-mixing regime, β c /µ = k −1 [39,51].
The choice of the interaction radius R is also relevant in the study of dynamical processes on motile agents as it defines the instantaneous underlying network structure.
We set the interaction radius R for all numerical simulations such that, in the static limit, the average connectivity of the network is slightly above the percolation threshold, k k c (d). Hence, using eq. (17) and considering the critical connectivities ( k c = 4.52 for d = 2 this case corresponds to an underlying contact network evolving much faster than the spreading process on its top, known as fast-switching or annealed network limit.
In this limit, the underlying structure can be effectively approximated by a fully-connected graph, in which at each time step all particles may interact with every one else with a probability proportional to the average connectivity k . In this mean-field regime the SIS dynamics can be solved exactly [52], and the epidemic threshold is C. Static to mean-field limit crossover To quantify this phenomenon, we consider the continuum model in the diffusive limit in a one-dimensional space, whose dynamics is governed by the following equa- with x ∈ [0, L] and where we are assuming that the two species diffuse with the same diffusion constant D.
We thus solved Eqs. (18) numerically using Euler explicit scheme using periodic boundary conditions. We discretized the equations on a grid of N g = 10 3 points, with ∆x = 1 (L = N g ), and ∆t chosen in way such that ∆t ≤ 1/2D. As control parameters, we move in D ∈ [0, 10 3 ] and β ∈ [10 −1 , 10 2 ]. For the discrete Laplace operator, we adopted a standard finite difference method.
Moreover, we added to both initial condition I(x, 0) and S(x, 0) a small amount of noise and we averaged over N s = 100 independent noise realizations.
The results from the numerical solution of the continuum model are shown in Fig. 3 (here the prevalence is ρ = dx S(x, ∞)). As one can see, the functional form of ρ as a function of β depends on D in a non-trivial way.
We obtain that larger values of ρ are reached sooner for large D values. We can thus define β c (D) as the value of β such that ρ > 5 × 10 −2 . Once we identify the static limit β c (0), we obtain that β c (D) = β c (0) for D/β < 1, while β c (D) starts to decrease (shown in Fig. 3   In Section IV C we showed that, in the one-dimensional continuum model, the static to mean-field crossover is exclusively governed by the diffusion coefficient D (Fig. 3).
This picture is confirmed from the agent-based numerical simulations, see Fig. 4 (b), which shows that the curves ρ(β, v, α) with the same diffusivity D = v 2 2α but different values of α collapse to a single curve ρ(β, D) (green, yellow, and red dots in Fig. 4 (b)). Notably, this is true only in the diffusive regime, while for very small values of α (blue dots) the prevalence is different.
Indeed, when β, the typical scale of the infection process, is higher than α, infections occur more frequently  (Fig. 4 (b)). Com-paring the blue symbols in Fig. 4

(b) for ballistic and
RnT agents moving at same velocity, we confirm that for α ≤ 1, the SIS process spreads at the time scales of the ballistic component of RnT motion. On the contrary, comparing the red symbols in Fig. 4 where brackets in the right hand side denote averages over independent steady-states. the epidemic threshold β c (D), diffusivity appears to control the crossover between different universality classes: the static one, corresponding to the SIS spreading on a finite dimensional network, and the homogeneous-mixing, mean-field one, corresponding to the SIS spreading in an infinite dimensional structure. Interestingly, beyond D ≈ R 2 , the emergence of a finite fraction of infected agents is controlled by mean-field behavior and is thus largely independent of the underlying dimension of the space in which agents move. Figure 5 shows that for large enough D, the curves ρ vs. β (once rescaled by β c ) are indistinguishable from the mean-field one within our numerical accuracy, giving support to the fact that the spreading mechanism in populations of fast enough diffusive agents is homogeneous, generically well captured by the homogeneous-mixing approximation.
In this work, we have studied the impact of motility on spreading dynamics. We established a general framework to tackle this problem, based on two paradigmatic models of both motility and spreading, namely, Run-and-  [29,34,53]. We focused in the diffusive limit, As a limiting situation, one has the so-called static limit, that is reached for vanishing values of the diffusion constant D → 0. This static limit is quite intuitive: in absence of diffusion, regions with overlap of different species are the only ones where spreading can occur. This situation is well-represented by an epidemic process on random geometric graphs [52], where the outbreaks will only occur in the connected components with at least one initially infected particle. On the other hand, as the diffusion becomes important, the model undergoes a crossover towards mean-field behavior, and thus the system reaches the so-called homogeneous-mixing limit, in which all particles can, on average, interact with all the others. We tested and documented the presence of this crossover by solving numerically the continuum model in one spatial dimension. The mean-field regime is obtained because the initial density profiles of infected and susceptible agents relax towards the homogeneous profile on a typical time scale smaller than the one associated to the SIS process. We observed that the critical value β c is diffusion-dependent, i. e., β c = β c (D), and continuously switches between two regimes: (i) for small D values, β c is D-independent, basically controlled by the behavior of the SIS model on a static short-ranged network, while (ii) it decreases to smaller β c as D is increased.
We tested the predictions of the continuum model against numerical simulations of SIS epidemic process on top of run-and-tumble walkers in two and three dimensions, showing the presence of the crossover from a mean-field regime, that is reached at high diffusion constant, to the static limit for small diffusion values. The nature of such crossover is revealed by the behavior of the local prevalence that is characterized by localized spots at zero diffusion that become more and more extended for increasing value of the diffusion constant. Above a threshold D c ≈ R 2 , the transition towards a state with non-zero prevalence appears to be of the mean-field kind. While the specific location of the epidemic threshold does depend on D, at large enough diffusivities, the transition is controlled by the mean-field, homogeneousmixing, regime.
As a future direction, it might be interesting to explore different motility regimes, e. g., when the persis-tence time is of the same order as the SIS time scales, allowing to investigate the role of persistent motion on the spreading dynamics [54]. It might be also interesting to study pattern formation driven by feedback between the SIS state variable and motility parameters [34,35].