The influence of active agent motility on SIRS epidemiological dynamics

Abstract

Active Brownian disks moving in two dimensions that exchange information about their internal state stochastically are chosen to model epidemic spread in a self-propelled population of agents under the susceptible-infected-recovered-susceptible (SIRS) framework. The state of infection of an agent, or disk, governs its self-propulsion speed; consequently, the activity of the agents in the system varies in time. Two different protocols (one-to-one and one-to-many) are considered for the transmission of disease from the infected to susceptible populations. The effectiveness of the two protocols are practically identical at high values of the infection transmission rate. The one-to-many protocol, however, outperforms the one-to-one protocol at lower values of the infection transmission rate. Salient features of the macroscopic SIRS model are revisited, and compared to predictions from the agent-based model. Lastly, the motility induced phase separation in a population of such agents with a fluctuating fraction of active disks is found to be well-described by theories governing phase separation in a mixture of active and passive particles with a constant fraction of passive disks.

---

1 Introduction

Epidemiological models have been in use for more than a century to track the spread of infectious diseases.¹ These models assume that the population is compartmentalized into three populations: susceptible (S) individuals, who are prone to the disease, infected (I) individuals who possess the disease and are capable of spreading it, and individuals who have recovered (R) from the disease and may or may not become prone to further infection.² The rates of interconversion between the various populations are specified by epidemiological constants, and the time-evolution of these populations is governed by coupled ordinary differential equations (ODEs).² Several modifications to this susceptible-infected-recovered (SIR) model have been made, to account for policy interventions such as vaccination drives, social distancing protocols, and lockdown mandates.³ SIR models and its variants have been used to describe the spread of several diseases,⁴ such as: the plague,¹ COVID-19,^5,6 varicella⁷ and influenza.⁸ These macroscopic models assume a uniform, well-mixed population, such that the epidemiological constants do not depend on the (heterogeneous) spatial density of the various populations. A fine-grained, discrete treatment of epidemic spread is agent-based modeling^9–12 (ABM), which considers each individual in the population as an autonomous entity, with the disease transmission occurring if a susceptible individual comes within a contagion radius of the infected individual. Another paradigm involves the use of network-based modeling to account for spatial heterogeneities in the population and variabilities in the epidemiological constants.^13–16 More recently, SIR modeling has been coupled^17–19 with advances in the understanding of self-propelled or active matter. Self-propelled entities such as light-activated Janus particles undergo a transition from a homogeneous state to a clustered one, at sufficiently high enough values of packing fraction and self-propulsion speed.^20,21 This experimentally observed phenomenon of motility induced phase separation (MIPS) can be recapitulated in numerical simulations in which the active particles translate with a constant self-propulsion speed, interact sterically, and whose orientations evolve diffusively in time.^22,23 It is therefore possible, using such systems, to examine disease propagation in both homogeneous and spatially heterogeneous configurations, by attaching to the individual active agents an internal state (S, I or R). If the internal state of the particle is coupled to its dynamics, such that the self-propulsion speed of infected particles is different from that of susceptible or recovered particles, then one would expect that MIPS in such a system would be affected. We address this question in the present work, and also aim to draw a connection between the agent-based, or microscopic, and macroscopic descriptions of epidemic spread.

Norambuena et al.¹⁹ simulated a collection of active Brownian particles as self-propelled agents moving in two dimensions in which a susceptible particle gets infected as soon as it comes within a cut-off distance of an infected particle. They considered a low particle density system, meaning that the number of active agents per unit dimensionless area was approximately 0.03, and used a mean-free path approach to determine analytical expressions for the rate of infection in their model. The internal state of the particle does not affect its dynamics, and all particles move with the same self-propulsion speed. Zhao et al.¹⁸ also endow an internal state to active agents and examine how the effectiveness of the disease spread differs from the macroscopic model that assumes a well-mixed population. Here too, a particle's self-propulsion speed is unaltered by its internal state. Useful future directions of research identified in their paper include the study of models in which the internal state is coupled to the motility of the particle, and the simulation of systems at a large enough area fraction so that MIPS can be observed. Forgács et al.¹⁷ study contagion dynamics using an agent-based modeling approach, where the state of the particle determines its motility, i.e., the infected particles move 50% slower than susceptible or recovered particles. Their simulation considers a collection of active Brownian disks moving in two dimensions, the majority of which belong to a large cluster that has undergone MIPS. The evolution of the disease for various values of the epidemiological constants is examined, along with a characterization of the spatial distribution of the susceptible particles around the infected particles. The dynamics of active particles that exchange information about their internal state has been studied in contexts outside of epidemiological modeling as well.^24,25 For example, Paoluzzi et al.²⁴ examine MIPS in a mixture of motile and non-motile particles that can change their identities upon collision. Quantifying phase separation arising in agent-based modeling of epidemiological systems remains an open question, particularly when the internal state is coupled to the motility of the agents, causing the system to behave as a transient mixture of active and passive particles. Also pertinent is establishing the connection (or lack thereof) between the macroscopic and microscopic descriptions of disease spread in populations, in the non-dilute limit. We address these questions in the present work and find that: (a) the analytical theory defining the phase boundary in a system of active and passive particles²⁶ can successfully be adapted to describe phase separation in a collection where the activity of the particles switches transiently, and (b) certain qualitative similarities can be observed in the disease statistics predicted by the macroscopic and microscopic models, although a direct mapping between the two is not found.

The remainder of this manuscript is organized as follows: Section 2 recapitulates salient features of the macroscopic SIRS model, Section 3 describes the numerical implementation of the microscopic model for two protocols of the spread of infection: one in which each infected particle can pass the disease on to exactly one susceptible particle within a contagion radius, and another in which an infected particle can potentially spread the disease to multiple susceptible particles within the contagion radius. A comparison of the contagion dynamics predicted by the one-to-one and one-to-many infection protocols are presented in Section 4. Section 5 presents a comparison of the disease dynamics obtained using the macroscopic and microscopic models. A discussion of phase separation in a mixture of active and transiently passive disks is presented in Section 6, accompanied by a description of the protocol used in the present work for identifying the occurrence of phase separation in such systems, which could be used to analyze MIPS in situations beyond the present study. We conclude in Section 7.

2 Macroscopic model revisited

We consider the SIRS model in this work, in which the immunity gained by the recovered particles can be lost, causing them to become susceptible to the infection again. The terms “macroscopic” model and “well-mixed” model are used interchangeably, to refer to the set of ODEs given by eqn (1). This model is typically used to capture the contagion dynamics of diseases like influenza^27,28 and the Omicron variant of SARS-CoV-2,²⁹ where the immunity gained is temporary. The population of each compartment is denoted by N_S, N_I and N_R, and evolve according to the following coupled ODEs:where β* indicates the rate of infection or transmission of the disease from an infected to a susceptible individual, γ* represents the rate of recovery of an infected individual, and α* represents the rate of relapse or the loss of immunity of a recovered individual, resulting in its conversion to a susceptible agent. We also define the relative transmission rate ω* ≡ β*/γ* and relative relapse rate μ* ≡ α*/γ* for convenience. The asterisks on the macroscopic (population-level) rate constants are used to distinguish them from their counterparts in the microscopic model, which is discussed later in the paper.

Depending on the context, the labels {S,I,R} could refer to the type of population being discussed, or the normalized value of that population, e.g., S = N_S/N. The governing equations [eqn (1)] are subject to the initial conditions S₀ ≡ S(t = 0), I₀, R₀, and the normalization condition

S + I + R = 1. (2)

Two possible steady-state solutions to the system of eqn (1) emerge: the first corresponds to epidemic extinction, with S_∞ = 1, I_∞ = 0, R_∞ = 0, and the second isValues of ω* < 1 correspond to epidemic extinction, while ω* > 1 results in the endemic state given by eqn (3). Interestingly, the steady-state susceptible population depends only on the rates of infection and recovery, while the populations of the other-two species (I and R) depend on all the three rate constants.

Fig. 1 illustrates the effect of the rate constants on the time-evolution of the infection, for different initial populations of the various populations. The rate constants and the initial conditions uniquely determine the dynamics of population evolution. The steady-state value of the various populations (S_∞, for instance), however, are solely determined by the ratios ω* and μ*, and are independent of the initial conditions. In the language of dynamical systems theory,³⁰eqn (3) is an attracting state for ω* > 1.

Fig. 1 Population evolution as a function of time for two parameter sets with ω* = 2, μ* = 1.5 and different values of the initial populations of susceptible and infected individuals, i.e., (a) S₀ = 0.7, I₀ = 0.3 and (b) S₀ = 0.9, I₀ = 0.1.

Fig. 2 illustrates a plot of the steady-state populations as a function of the rate constants. Fig. 2(a) illustrates the case of varying the infection rate fixed values of the recovery and relapse rate. Low values of the relative transmission rate, i.e., ω* < 1 result in epidemic extinction, in which there are no infected particles in the long-time limit. As ω* ≥ 1, however, the steady-state population of susceptibles declines with an increase in the infection rate. The transition of the steady-state numbers from the low ω* branch to the high ω* branch appears to follow a transcritical bifurcation.³⁰ Holding the infection and recovery rates constant, while varying the relapse rate, as shown in Fig. 2(b), has no effect on the population of susceptibles. There are more recovered than infected individuals for values of μ* < 1, while the balance is reversed as the value of μ* crosses unity. Lastly, Fig. 2(c) examines the consequences of varying the recovery rate as the infection and relapse rates are held constant. At both small and large values of the recovery rate, the recovered population is nearly zero. This is because for ω* < 1 the population is solely made of susceptible individuals, while the infected individuals dominate the population for ω* ≫ 1.

Fig. 2 Steady-state population as a function of (a) relative transmission rate at constant values of recovery and relapse rates, (b) relative relapse rate at constant values of infection and recovery rates, and (c) relative transmission rate at constant values of infection and relapse rates. Data in Fig. 1 and 2 have been generated by solving the coupled ODEs given by eqn (1), using a code shared on Mathworks File Exchange.³¹

3 Numerical simulations of microscopic model for SIRS dynamics in active Brownian particles

We simulate a system of N active Brownian particles (ABPs) of unit radius each (d = 2r = 2) moving in a periodic square box of side L. The position of the particles evolves in time according to

ṙ_i = v^HM_i + U₀e_i (4)

where U₀ is the self-propulsion speed of the disks. The particles move in the direction e_i ≡ [cos θ_i,sin θ_i], where θ_i denotes the particle's orientation measured with respect to the positive x-axis. The particle orientation evolves in time according to a rotational diffusive process, such that

〈[small theta, Greek, dot above]_i(t)[small theta, Greek, dot above]_i(t′)〉 = 2D_rδ(t − t′) (5)

where D_r denotes the rotational diffusion constant. The positions and orientations of the particles are updated using a forward Euler algorithm with a timestep Δt. The v^HM_i term on the RHS of eqn (4) represents a harmonic interaction that operates only when the centre-to-centre separation of disk i and k is smaller than their diameter d. The steric interaction is proportional to the extent of overlap of the disks, and acts to alter the position of each disk in an overlapping pair such that they are just in contact. The functional form of v^HM_i is given bywhere r̂_ik = r_ik/r_ik is the unit vector along the line joining the particle centres, Θ denotes the Heaviside function, the stiffness K = 0.5, and n_o denotes the number of overlapping particles in the neighborhood of the ith particle.

We next describe the update rules to simulate the spreading of diseases in this microscopic, agent-based model, which largely follows the algorithm outlined by Forgács et al.¹⁷ The simulations are initialized in a homogeneous configuration, with all the agents arranged in regularly spaced intervals on a square lattice. Each particle has an associated internal state, i.e., it could be a susceptible (S), infected (I), or recovered (R) agent. The contagion dynamics (change in the internal state of the agent, i.e., S, I, R) as well as the physical movement of the agents happen simultaneously. This approach differs from that adopted by Forgács et al.,¹⁷ in which the contagion dynamics starts from a phase separated state.

A census of the number of particles in each sub-category is taken at the beginning of each timestep. The allowed transitions are S → I, I → R and R → S, with the rate constants associated with the transitions given by β, γ and α, respectively. The protocol governing the spread of the infection is as follows: a loop is run over all the N_I infected particles in the box at a given time instant, and the number of susceptible and infected particles within a cut-off radius r_c of an infected particle I⁽ⁿ⁾ is recorded as N⁽ⁿ⁾_S and N⁽ⁿ⁾_I, respectively. Note that the index n runs over all the infected disks in the system at a given time instant. By this definition, the smallest allowable value for N⁽ⁿ⁾_S is zero, while that for N⁽ⁿ⁾_I is unity. The probability of infection is calculated as p_inf = βN⁽ⁿ⁾_IΔt. The probability of recovery is given by p_rec = γΔt. The values of the three probabilities {p_inf,p_rec,p_rel} are capped at unity. If any of these probabilities exceed unity during a timestep, it is reset to unity.

We consider two protocols by which the infection might spread:

Protocol A: one-to-one

For each value of the index n during the loop over N_I, if the outcome of a binomial trial with probability p_inf is non-zero, then one susceptible particle from N⁽ⁿ⁾_S is picked at random for conversion to an infected particle. Similarly, if the outcome of a binomial trial with probability p_rec is non-zero, and there are multiple infected particles within a cut-off radius r_c of an infected particle I⁽ⁿ⁾, then the state of one infected particle changes to recovered. An infected particle in this protocol can spread the disease to only one other particle in a timestep. Having multiple infected particles in the vicinity of a susceptible particle only increases the probability of infection.

Protocol B: one-to-many

For each value of the index n during the loop over N_I, the number of susceptible particles that get infected is decided by drawing a sample from a binomial distribution of probability p_inf, N⁽ⁿ⁾_S number of times, and calculating S⁺, the total number of successful outcomes. A total of S⁺ particles from amongst N⁽ⁿ⁾_S are then randomly selected to be infected. Once the loop over all N_I is completed, the number of infected particles that recover is decided by drawing a sample from a binomial distribution of probability p_rec, N_I number of times, and calculating I⁺, the total number of successful outcomes. A total of I⁺ particles from amongst N_I are then randomly selected to undergo recovery. In this protocol, therefore, an infected particle could potentially transmit the disease to multiple susceptible particles in its neighborhood.

For both the protocols discussed above, the relapse of recovered particles into the susceptible category is governed by a binomial process of probability p_rel = αΔt, and is invoked once the loop over all N_I is completed. Furthermore, the infected particles become immobile (U₀ = 0), while both the susceptible and recovered particles retain their activity. The infected particles regain their mobility upon recovery. The instantaneous fraction of active disks in the system is therefore given by x_A ≡ (N_S + N_R)/N.

The population and system size in all cases are chosen to be N = 1600 and L = 100, corresponding to a number density of ρ = N/L² = 0.16, and an area fraction of ϕ₀ = ρπd²/4 ≈ 0.5. The self-propulsion speed is fixed at U₀ = 0.1. A discrete timestep width of Δt = 0.1 is used in all the simulations, such that the displacement of a disk over a single timestep is smaller than its diameter. The contagion radius is chosen to be r_c = 3 in all the simulations, unless specified otherwise. The motility of the system is quantified using the Péclet number, Pe ≡ 3U₀/dD_r. The persistence length [small script l] ≡ U₀/D_r of the active particles must be smaller than the box dimensions, to minimize finite-size effects.^32,33 We adjust the rotational diffusion constant at each Pe so that this condition is met. The effects of thermal noise on the translational motion are ignored in the present work, but may be included by adding a white-noise process to the RHS of eqn (4). In Section 4 that compares the contagion dynamics predicted by the two protocols (A and B), and in Section 5 that explores the connection between the microscopic and macroscopic models, the numerical results are obtained from simulations that are n = 2 × 10⁴ steps long. In Section 6 that charts the phase behavior of a collection of active and transiently passive disks, simulations of at least [scr O, script letter O](10⁵) steps are used. The total length of the simulation is simply the product of the number of simulation steps and the discrete timestep, and is denoted by t_sim = nΔt.

4 Comparison between protocols A and B

In this section, we compare the dynamics of infection spread, and the steady-state statistics for numerical simulations performed using the two protocols (A and B) discussed above. Given the broad span of the parameter space, we restrict ourselves to studying two values of the relative transmission rate ω ≡ β/γ, and perform a scan across a range of relative relapse rates, μ ≡ α/γ, by holding the recovery rate constant at γ = 0.1, and varying α and β appropriately. We divide the parameter space into four quadrants as shown in Fig. 3 and examine the simulation results accordingly.

Fig. 3 Schematic of parameter space explored. Parameters in quadrant I correspond to the case of high infectiousness (or transmission), and high loss of immunity (relapse), while quadrant III corresponds to the case of low infectiousness and low relapse. Other quadrants could interpreted in a similar manner.

Fig. 4(a) explores the first and second quadrants of the parameter space and illustrates that, at a higher value of the relative relapse rate, μ = 2, we notice two different behaviors, based on the value of the relative transmission rate. In quadrant I, both the relative transmission and the relapse rates are high, implying that the disease spread is more probable, as is the replenishment of the numbers of the susceptible population. In this quadrant, we observe that the higher transmission rate plays a more dominant role, leading to a faster spread of the disease, and a rapid decrease in the number of susceptibles. Furthermore, it is immaterial if the infection spreads via the one-to-one or the one-to-many route, as they both result in a nearly identical prediction for the steady-state susceptible population. In quadrant II, the relapse rate is higher than the infection rate. The steady-state in such a case is decided by the protocol for disease spread. Stipulating a one-to-one spreading protocol causes the infection to die out faster, so that the population is entirely composed of susceptibles in the long time limit. Allowing for a one-to-many spreading protocol for the infection compensates for the low value of the transmissibility in this quadrant, resulting in a faster spread of the disease, and a lower value of the steady-state susceptible population as compared to the one-to-one protocol. Fig. 4(b) explores the third and fourth quadrants of the parameter space, for the smallest value of the relative relapse rate considered in this section, μ = 0.1. Over the time window considered in the figure, the numbers predicted by the one-to-many protocol are comparable to or lower than that predicted by the one-to-one protocol, indicating, unsurprisingly, that the former is more efficient in spreading the disease. This effectiveness of the one-to-many protocol, however, is manifest only transiently, as the steady-state values of the susceptible population are independent of the relative transmissibility, and the route for disease spread.

Fig. 4 Time-evolution of the susceptible population predicted by the two protocols at (a) high and (b) low values of the relative relapse rate μ ≡ α/γ, for relative transmission rates of ω = 0.5 and ω = 2.

From the above analysis, it is clear that the sharpest contrast between the steady-state outcomes predicted by the two protocols occurs in quadrant II. Given the fluctuations in the values of the susceptible population, the average steady-state value, S_∞, is estimated by computing the mean of the last 20% of the time series. Fig. 5 illustrates the steady-state susceptible population as a function of the relative relapse rate, for two different values of the relative transmission rate. The effectiveness of the one-to-many protocol in governing the spread of infection is most evident from the low transmission regime as identified in Fig. 5(a). For μ > 0.2, the S_∞ resulting from the one-to-many protocol is significantly lower than that predicted by the one-to-one protocol, indicating a spread of the disease amongst a larger fraction of the population. In the high transmission regime, however, the protocol for disease propagation has a less pronounced effect on the steady-state statistics, as evinced by Fig. 5(b). Having compared the outcomes from the two protocols, we now present results obtained with protocol A only for the rest of the paper.

Fig. 5 Steady-state susceptible population as a function of the relative relapse rate μ for (a) low and (b) high values of the relative transmission rate ω. Error bars represent standard deviation of data obtained from the time-averaging procedure used to estimate the mean steady-state population. Where invisible, error bars are smaller than symbol size.

5 Connection between microscopic and macroscopic models

A major distinction between the microscopic agent-based model (ABM) and the macroscopic population-based model for disease spread is the specification of a contagion radius r_c for the former, which makes the spread of the infection depend not only on the number of infected and susceptible individuals at a given time, but also on their locations. We illustrate below a few salient features of the steady-state numbers predicted by the microscopic model, and how they compare to the macroscopic model predictions.

Fig. 6(a) illustrates that the steady-state population for the microscopic model is independent of the initial fraction of the infected population I₀, at a fixed value of the contagion radius (r_c = 3), for various values of the relative transmission and relapse rates. The independence from initial conditions, over the range examined in this figure, is a trait shared by the microscopic and macroscopic models.

Fig. 6 Steady-state susceptible population as a function of (a) initial fraction of infected individuals, (b) relative relapse rate at fixed values of infection rate, (c) relative transmissibility at fixed values of relapse rate, and (d) contagion radius at fixed values of the epidemiological constants. All results obtained using protocol A (one-to-one) for disease spread.

Fig. 6(b) shows that the steady-state susceptible population decreases as a function of the ratio of the relative relapse rate (μ), for fixed values of ω and the contagion radius (r_c = 3). This marks a crucial departure from the macroscopic model (see Fig. 2(b)) in which S_∞ is solely a function of the relative transmission rate.

The variation of the steady-state population as a function of the relative transmission rate ω is shown in Fig. 6(c), for fixed values of μ and contagion radius (r_c = 3). The susceptible population is independent of the relative transmission rate for small values of the latter. Beyond a threshold value of the relative transmission rate, however, the susceptible population decreases with ω. The crossover value depends on the ratio α/γ, unlike in the macroscopic model where the transition occurs at ω* = 1 and is independent of the relapse rate. While the existence of a bifurcation in the macroscopic model predictions (a system of coupled ODEs) is unsurprising,³⁰ it is remarkable that an evidence of bifurcation is also seen in the agent-based model. This also indicates that the stochastic update rule for the various compartments is faithful to the contagion dynamics as predicted by the ordinary differential equations of the macroscopic model.

Lastly, Fig. 6(d) illustrates the dependence of S_∞ on the contagion radius, for various values of the epidemiological rate constants. For all the cases examined in the figure, the steady-state susceptible population decreases as a function of r_c.

We briefly revisit the comparison between protocols A and B before concluding this section. The results reported in Section 4 all come from numerical simulations with a contagion radius of r_c = 3. Fig. 7 illustrates the effect of using a larger value of the contagion radius, r_c = 10, on the steady-state numbers obtained using protocol A, over a range of relative relapse rates and a fixed relative transmission rate of ω = 2. Keeping the contagion radius fixed at r_c = 10, using a lower value of ω, or a different protocol results in a nearly identical plot. We note that at a large value of the contagion radius, there is essentially no distinction between the steady-state predictions of the two protocols. The long-time population of susceptibles is practically zero and the particles are split between the infected and recovered categories. This indicates that when an infected particle has an abundance of susceptible neighbors to transmit the disease, it is immaterial if the disease spreads via the one-to-one or one-to-many route, and there are negligible susceptible individuals remaining in the long time limit.

Fig. 7 Steady state population as a function of the relative relapse rate, as obtained using protocol A.

In the limit of a large contagion radius, each particle can “see” all the other particles in the box, and one could therefore expect that the effect of spatial heterogeneity is reduced, bringing the microscopic model predictions closer to that obtained from the macroscopic model. Probing this line of thought, we note a qualitative similarity between Fig. 7 and the macroscopic model results given by Fig. 2(b), in that the steady state numbers are independent of the relative relapse rate. A distinction between the microscopic and macroscopic model predictions is that while the former predicts a vanishing of the susceptible population across the range of the relative relapse rates considered, the latter predicts a finite non-zero value for the steady-state susceptible population.

In Fig. 8, the macroscopic and microscopic model results are plotted simultaneously, with the caveat that even though the same numerical values have been used for the epidemiological constants (e.g. α = α* = 0.05), there is no direct mapping between the two models. Keeping the relative relapse rate μ fixed, increasing the relative transmission rate ω drives the spread of infection from a state of extinction (S_∞ = 1) to one in which the fraction of susceptible individuals has reduced considerably. As noted in the discussion of Fig. 6(c), the relative transmissibility at which the transition away from the epidemic extinction state occurs depends on the value of μ in the microscopic model. For the macroscopic model, however, the location of this transition is fixed at the analytically determinable value of ω* = 1, and is independent of the relative relapse rate μ. Fig. 8 illustrates the effect of the contagion radius on the location of this bifurcation: higher values of r_c push the transition to lower values of ω. This makes intuitive sense: a smaller infection rate is required when the infected particles can see a larger number of the susceptible population, resulting in a more effective spread of the disease.

Fig. 8 Steady-state epidemic statistics as a function of the relative infection rate, at a contagion radius of (a) r_c = 3 and (b) r_c = 10. The same legend scheme is followed in both the subfigures, with filled symbols denoting the microscopic model results, obtained using protocol A, and hollow symbols representing the well-mixed model results. The numerical values of the parameters used in both the models are the same, although there is no direct mapping between the two.

We have examined additional factors which could determine the bifurcation point in the microscopic model. The contagion dynamics for agents moving with a reduced self-propulsion speed U₀ = 0.05, for two cases is analyzed. In the first case, the rotational diffusivity is set to D_r = 10⁻³ so that the Péclet number remains at Pe = 75, at which MIPS is observed for an area fraction of ϕ₀ = 0.5.³⁴ In the second case, the rotational diffusivity is left unchanged at D_r = 2 × 10⁻³, so that the Péclet number falls below the threshold required for observing MIPS at ϕ₀ = 0.5. From Fig. 12, we see that the crossover point in all the cases is not only a function of the relapse rate and the contagion radius, but also depends on the mobility of the agent (U₀ and D_r) which indirectly determines the local packing fraction in the box.

The connection between the well-mixed model and the agent-based model has been examined in detail by Paoluzzi et al.³⁵ They consider mobile agents on a two-dimensional lattice (in a periodic box of size L) that undergo SIR dynamics. The length of the steps is governed by the Lévy exponent (called β in their work, but we will use the symbol λ, to avoid confusion with the rate of infection), and the step direction is chosen from a uniform random distribution. In the limit of large mobility coefficient λ → 2, the motion of the agents is akin to Brownian motion, while the λ → 1 corresponds to a Lévy flight where the agents can take steps whose lengths are picked at random from the interval [0,L/4] using Mantegna's algorithm.³⁶ The lower values of λ are seen to agree with the analytical results for the well-mixed SIR model. Paoluzzi et al.³⁵ also studied the effect of a mixture (high and low λ) of the mobility coefficients on the contagion dynamics. They find that even a small number of sites with a small value of λ (meaning higher mobility) can trigger epidemic waves. Note that the agents have no finite-size and hence no steric repulsion exists between them. The step sizes are entirely user-defined, and drawn from a known distribution. In our work, although the step size is uniform by design [with a value of U₀Δt per unit time, as seen from eqn (4), where U₀ is the self-propulsion speed of the ABP], the actual sizes of the steps vary due to steric repulsions between the agents. In any case, the maximum size of the step taken in any timestep is smaller than the particle diameter (d = 2), and our simulation box dimensions (L = 100) are such that d ≪ L. These step sizes are far smaller than the ones encountered by Paoluzzi et al.³⁵ The absence of such long steps is perhaps the reason why the microscopic model in our case does not completely converge to the well-mixed model results, although qualitative similarities are observed when the contagion radius is increased.

A common paradigm to study the spread of epidemic is to use a network-based approach,^13–16 in which the members of the population are represented by nodes, and their connectivity denoted by the edges that join them. The infectiousness of the disease-causing vector is allowed to be different for each node in the network. This approach allows the decoupling of the connectivity of the agents from the probability of disease transmission. In the active particle-based model considered in our paper, although the input parameters for the agent mobility and the epidemiological constants are picked independently, the coupling of their effects is an emergent phenomenon, due to the protocol of disease spread which depends on the spatial positioning of the various disks. Großmann et al.¹⁴ consider SIR-dynamics on a static network where the infectiousness of the nodes can take on a distribution of values. They find that a large variation in the infectiousness leads to a smaller final size of the epidemic, stemming from an increased probability of epidemic extinction, and therefore a lowering of the herd immunity threshold.³⁷ In our paper, increasing the contagion radius appears to result in a similar outcome, as evinced in Fig. 8, as larger r_c results in smaller values of the transmission rate for making the infection endemic, i.e., S_∞ < 1.

We conclude this section by noting that the precise mapping (if indeed one exists) between the rate constants used in the ABM and those appearing in the ordinary differential equations at the population-level remains unknown, we predict that the contagion radius r_c and the transmission protocol would crucially affect this relationship.

6 Phase separation in motility-modified SIRS model

A collection of self-propelled (or active) particles interacting sterically undergo a motility induced phase separation (MIPS) at large enough values of the Péclet number and the area fraction ϕ₀ of the particles. In this transition, the particles go from being in a gas-like, single-phase to a phase-separated state consisting of a dense large cluster that dynamically exchanges particles with the surrounding dilute phase [see ref. 22 and 23 for an extensive review of the topic]. The boundary separating the homogeneous state from the phase-separated one in the Pe–ϕ₀ plane has been determined through direct numerical simulations^34,38,39 and analytical theory.^38,40–43

The effect of the presence of passive particles – ones that move translationally under the effect of Brownian noise or not at all – on the phase separation behavior of an active–passive mixture has also received interest.^26,44–47 Stenhammar et al.²⁶ studied such a mixture with a total area fraction of ϕ₀, of which a number fraction x_A is active. They derive the following analytical expression for the phase boundary in the Pe–x_A plane:

using the kinetic model introduced by Redner et al.³⁸ Here κ is an empirical fitting parameter³⁸ that represents the average total number of particles that are lost from a phase-separated cluster in an escape event. The boundary obtained using κ = 4.05 is seen to accurately demarcate the homogeneous and demixed states in the phase diagram generated from numerical simulations at ϕ₀ = 0.6. A value of κ = 4.5 accurately predicts the phase boundary in the Pe–ϕ₀ plane for a system composed solely of active disks (x_A = 1). Takatori and Brady⁴⁴ derive an expression for the phase boundary in an active–passive mixture using an alternative approach that relies on the concept of active swim pressure in a collection of self-propelled swimmers. They obtain an agreement with the predictions of Stenhammar et al.²⁶ without any fitting parameters.

We summarize the common metrics used in the literature to quantify motility-induced phase separation, before describing the algorithm introduced in the present work to identify systems that have undergone phase separation. The fraction of particles in the largest cluster N_LC/N is a popular metric^32,48,49 to track the approach to MIPS. When the majority of the particles in the system belong to a single cluster, it is taken to be an indication of MIPS. Large fluctuations in the number of particles within a subregion are taken to be a sign of inhomogeneity and the onset of MIPS.^50,51 Another signature for the occurrence of MIPS is the appearance of a bimodality in the probability distribution of the local area fraction.^52–54 As a system evolves from a homogeneous state to a phase separated one, the size of the domains (calculated from either the static structure factor or the density correlation function) grows with time as a power law^34,38,55,56 [script L](t) ∼ t^1/3. These metrics have also been used to analyze systems containing a mixture of active and passive particles,²⁶ in which the relative populations of the two species remain constant in time.

In the present work, however, the fraction of active and passive disks fluctuate in time due to the disease spreading, and we seek to identify an appropriate metric for the identification of MIPS in such systems. To that end, Fig. 9 illustrates the time evolution of N_LC/N in a system consisting entirely of susceptible particles, with the epidemiological constants set to zero. There is no spread of infection in such a system, and the fraction of active disks is therefore unity at all times. For the values of Pe and the area fraction considered in Fig. 9, the system undergoes MIPS, as evinced by the snapshots recorded at the various time instances [Fig. 10]. Our goal is to define metrics for the identification of MIPS based on this reference time series, for application to other systems in our work in which the fraction of active disks fluctuate in time. Firstly, we note that the normalized size of the largest cluster reaches a steady state value of N_LC/N ≈ 0.83 following an initial transient. Secondly, after t ≈ 5520, there are no significant dips in the value of N_LC/N, and the fluctuations in this quantity are minimal. We use these two observations to devise a methodology (Fig. 11 and 13) for ascertaining if a system has undergone MIPS or not, given the time-series of the largest cluster. In case a system has phase separated, this algorithm also estimates the time needed for MIPS. The various parameters needed by the algorithm are {M_cut,M_min,τ,σ}, and a brief explanation is as follows. If the fractional size of the largest cluster μ ≡ N_LC/N does not exceed M_cut at any point in its time series, then we consider that MIPS has not occurred. The algorithm searches for a chunk of data in the time series in the interval [t_s,t_s + τ], such that each data point in the interval exceeds M_min. If such a chunk is not found in the input data series, the algorithm concludes that MIPS has not occurred. Provided such a data chunk is found, we then test if the standard deviation of the data series (normalized by the total number of particles), is smaller than σ. If this requirement is met, the algorithm concludes that MIPS has occurred and returns t_s, the starting point of the data chunk, as the time at which MIPS has occurred. If the standard deviation of the data series (normalized by the total number of particles) exceeds σ, then we conclude that MIPS has not occurred. The parameter values used in this algorithm are listed in Table 1. The time to MIPS as estimated by the algorithm for the timeseries indicated in Fig. 9 and 11 are 3730 and 5527, respectively.

Fig. 9 Normalized size of the largest cluster as a function of time, for the case of N = 1600 active disks in which all the epidemiological constants have been set to zero and there is no contagion dynamics. Snapshots of the system at the three indicated time-instances are given in Fig. 10.

Fig. 10 Snapshots of the system in Fig. 9, recorded at (a) t = 1000, (b) t = 5527, and (c) t = 8000, respectively.

Fig. 11 Cluster size evolution as a function of time for two representative systems. Parameters needed for the algorithm that detects motility induced phase separation are identified here, and their numerical values provided in Table 1.

Fig. 12 Steady-state value of susceptibles as a function of the relative recovery rate, for agents moving at a slower self-propulsion speed U₀ = 0.05, with (a) Pe = 75, at which MIPS is observed, and (b) Pe = 37.5, at which no MIPS is observed for ϕ₀ = 0.5.³⁴ The contagion radius used in both (a) and (b) is r_cut = 3. Each data point in the figure was obtained from simulations of n = 5 × 10⁴ steps.

Fig. 13 Flowchart illustrating the algorithm for ascertaining if a system has undergone MIPS, and to evaluate the time needed for phase separation in case it has.

Table 1 Parameters used in the algorithm (Fig. 13) to determine if phase separation has occurred or not

Parameter	Value used in present work
M cut	0.8
M min	0.6
τ	0.1tsim
σ	0.1

---

Fig. 14 explores the phase behavior of the system in the Pe–x_A plane, obtained from simulations using a wide range of the epidemiological constants (see Table 2 for values). At the steady-state value of the fraction of active disks for any given value of the Péclet number, we denote if MIPS has occurred or not using the algorithm described above. The boundary between the homogeneous states and the MIPS states is well described by eqn (7) when a value of κ = 1.875 is used. Snapshots of the system at a fixed value of Pe and varying fractions of the active disks are given in Fig. 15A–D. We notice that there is no preference for disks with identical internal states to cluster together. A more quantitative analysis would involve the calculation of the pair correlation function for the various populations.

Fig. 14 Phase separation in a mixture of active and passive disks. Open circles indicate a homogeneous phase, while closed diamonds represent a phase-separated system. Dash-dotted redline represents eqn (7) with κ = 1.875. Snapshots of the system at the locations A, B, C, D are given in Fig. 15.

Table 2 Data points in the Pe–x_A phase plane (Fig. 14) and the epidemiological constants used in the simulations for obtaining them

Pe	x A	β	γ	α
50	1.0	0	0	0
	0.84	0.1	0.1	0.05
	0.62	0.05	0.1	0.2
	0.53	0.2	0.15	0.2
	0.38	0.2	0.1	0.2
	0.22	0.2	0.05	0.2
	0.004	0.2	0.001	0.2
60	1.0	0.05	0.1	0.05
	0.82	0.1	0.1	0.05
	0.63	0.2	0.2	0.2
	0.55	0.2	0.15	0.2
	0.42	0.2	0.1	0.2
	0.23	0.2	0.05	0.2
	0.03	0.2	0.005	0.2
75	1.0	0	0	0
	0.78	0.1	0.1	0.05
	0.59	0.2	0.2	0.2
	0.51	0.05	0.1	0.2
	0.37	0.2	0.1	0.2
	0.22	0.2	0.05	0.2
	0.006	0.2	0.001	0.2
90	1.0	0.05	0.1	0.05
	0.76	0.1	0.1	0.05
	0.56	0.2	0.15	0.2
	0.41	0.2	0.1	0.2
	0.21	0.2	0.05	0.2
	0.03	0.2	0.005	0.2
100	1.0	0	0	0
	0.78	0.1	0.1	0.05
	0.59	0.2	0.2	0.2
	0.53	0.05	0.1	0.2
	0.39	0.2	0.1	0.2
	0.23	0.2	0.05	0.2
	0.008	0.2	0.001	0.2

---

Fig. 15 Representative snapshots of the system for Pe = 75 at the active-disk fraction (A–D) indicated in Fig. 14 and recorded at the halfway point of the total simulation. The black disks denote susceptible agents, while blue and orange disks represent infected and recovered individuals, respectively.

In Fig. 16, the time taken for motility induced phase separation, as identified using the algorithm described in Fig. 13, is plotted as a function of the fraction of active disks in the system, for a range of Péclet numbers. We observed no definitive trend, implying that we cannot conclude if the presence of transiently immobile disks helps to aid or suppress phase separation in an active–passive mixture. Forgács et al.¹⁷ observe that the presence of quenched disorder, or immobile obstacles, in an active matter system causes the formation of numerous small clusters in addition to the large cluster that is characteristic of MIPS. Additionally, MIPS is a re-entrant phenomenon,^41,48 meaning that increasing the Péclet number in an already phase separated system can cause the system to go back to being in a homogeneous phase.

Fig. 16 Time taken for motility induced phase separation, as a function of the active disk fraction.

Another common metric to track motility induced phase separation is the local area fraction ϕ_loc. This is evaluated by dividing the periodic box into multiple smaller boxes and measuring the area fraction occupied by disks in each sub-box, to get a distribution of ϕ_loc values.⁵⁴ In Fig. 17 and 18, the probability and cumulative distribution functions of the local packing fraction (respectively), for the data series represented in Fig. 9 of the manuscript, are plotted at various timepoints. The onset of MIPS is indicated by the appearance of bimodality in the probability distribution function. The CDF peaks sharply around the average packing fraction (ϕ₀ = 0.5) in the absence of MIPS, and is seen to broaden as MIPS progresses (Fig. 18).

Fig. 17 (a) Probability and (b) cumulative distribution functions of the local area fraction, for the motility induced phase separation process represented in Fig. 9.

Fig. 18 (a) Probability and (b) cumulative distribution functions of the local area fraction, for the snapshots given in Fig. 15.

We have also calculated the average number of mobile (susceptible and recovered) agents that are nearest neighbors to a an immobile (infected) agent at a given time (t), similar to the methodology adopted by Forgács et al.¹⁷ This metric is calculated as follows:

where the indicator function returns 1 (0) if its argument is true (false). The numerical implementation of eqn (8) allows for a 0.5% tolerance, and uses a value of 1.005d in place of d. In Fig. 19, the average number of neighbors is evaluated for the Pe = 75 case, for steady-state active fractions (x_A) on either side of the phase boundary. We note that point A in the phase diagram is obtained for a system in which all the particles are of the susceptible type, with the epidemiological constants set to zero. There are no infected particles in this case, and η(t) is therefore not defined. It is clear from the figure that the average number of neighbors drops with the decrease in the fraction of active disks in the system. We performed a similar analysis for all the other Péclet numbers we have examined in the current study, and focused on those points at which MIPS has not occurred. It is clear that the infected particles precipitate the formation of microclusters, even in the absence of a global motility induced phase separation (Fig. 20).

Fig. 19 Average nearest neighbors η as a function of time for Pe = 75 at various values of the steady-state active fractions. The alphabets in parentheses in the legend entry correspond to the points indicated in Fig. 14.

Fig. 20 Average nearest neighbors η as a function of time for (a) Pe = 50, (b) Pe = 60, (c) Pe = 90, and (d) Pe = 100, at various number fractions of the active disks (x_A).

7 Conclusions

We performed an agent-based modeling of disease spread according to the SIRS model using a collection of active Brownian particles moving in two dimensions whose internal state encodes their state of infection. Two protocols for infection were considered, and their efficacies for the spread of the disease were analyzed for various combinations of the epidemiological constants. The coupling of the particle's internal state to its motility causes the population to behave as a collection of particles in which the fraction of active disks is time-dependent. We developed an algorithm to determine the occurrence of motility induced phase separation in such systems with transient activity, and find that it is well-described by the theories for phase separation in active–passive mixtures where the fraction of active disks remains constant in time. Although a direct mapping between the agent-based (microscopic) and macroscopic model is not found, several common features between the contagion dynamics predicted by the two models are noted. We see evidence for a transcritical bifurcation in the microscopic model where the agents are modeled as active Brownian particles. The use of active Brownian disks permits a tractable method to tune the density distribution of the system by changing the Péclet number. Humans in general, however, do not move at a constant self-propulsion speed with randomly varying orientations. Simulating the dynamics of individuals in a crowd has typically relied on the use of social forces that describe the interaction between the individual members.^57,58 The use of such pedestrian models to describe the motion of individual agents could permit the extension of the present work to model epidemic spread in human populations. Another interesting exercise for future work could be the effect of the type of motility on the nature and location of the bifurcation point, i.e., would a system of agents modeled using social forces exhibit a different kind of bifurcation when the steady-state numbers are plotted as a function of the relative rate of transmission. The concept of over-dispersion has been observed in the case of the COVID-19 pandemic,¹⁴ in which a few members of the population transmit the infection to many, while most individuals infect only a few or none at all. This aspect could be included in our framework by prescribing that certain infected agents in the system to follow protocol A (one-to-one), while a few others follow protocol B (one-to-many).

Author contributions

R. K. and A. K. designed research, performed research, contributed new analytic tools, analyzed data, and wrote the paper.

Data availability

MATLAB codes for the solution of the macroscopic epidemiological model and HOOMD-blue⁵⁹ scripts for agent-based modeling are available freely on GitHub.⁶⁰ Visualization of simulation data was performed using OVITO.⁶¹ Simulation data are available upon request from the authors.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We gratefully acknowledge the support of the Charles E. Kaufmann Foundation of the Pittsburgh Foundation (Grant No. 1031373-438639). R. K. thanks Ligesh Theeyancheri for valuable insights regarding the calculation of the local area fraction.

References

1. W. O. Kermack and A. G. McKendrick, Proc. R. Soc. London, Ser. A, 1927, 115, 700–721.
2. R. M. Anderson and R. M. May, Infectious Diseases of Humans, Oxford Science Publications, Oxford, 1991.
3. M. Cabrera, F. Córdova-Lepe, J. P. Gutiérrez-Jara and K. Vogt-Geisse, Sci. Rep., 2021, 11, 10170.
4. O. N. Bjørnstad, K. Shea, M. Krzywinski and N. Altman, Nat. Methods, 2020, 17, 557–558.
5. Y. Wu, Y. Sun and M. Lin, Comput. Electr. Eng., 2022, 102, 108230.
6. T. Liu, J. Huang, Z. He, Y. Zhang, N. Yan, C. J. Zhang and W. K. Ming, Sci. Rep., 2023, 13, 9164.
7. J. O. Giraldo and D. H. Palacio, Epidemiol. Infect., 2008, 136, 679–687.
8. D. Osthus, K. S. Hickmann, P. C. Caragea, D. Higdon and S. Y. Del Valle, Ann. Appl. Stat., 2017, 11, 202–224.
9. F. Peruani and G. J. Sibona, Phys. Rev. Lett., 2008, 100, 168103.
10. F. Peruani and G. J. Sibona, Soft Matter, 2019, 15, 497–503.
11. M. C. González and H. J. Herrmann, Phys. A, 2004, 340, 741–748.
12. J. P. Rodríguez, F. Ghanbarnejad and V. M. Eguíluz, Sci. Rep., 2019, 9, 6463.
13. I. Z. Kiss, J. C. Miller and P. L. Simon, Mathematics of Epidemics on Networks, Springer, 1st edn, 2017.
14. G. Großmann, M. Backenköhler and V. Wolf, PLoS One, 2021, 16, e0250050.
15. C. Nowzari, V. M. Preciado and G. J. Pappas, IEEE Control Syst., 2016, 36, 26–46.
16. R. Pastor-Satorras, C. Castellano, P. Van Mieghem and A. Vespignani, Rev. Mod. Phys., 2015, 87, 925–979.
17. P. Forgács, A. Libál, C. Reichhardt, N. Hengartner and C. J. Reichhardt, Sci. Rep., 2022, 12, 11229.
18. Y. Zhao, C. Huepe and P. Romanczuk, Sci. Rep., 2022, 12, 2588.
19. A. Norambuena, F. J. Valencia and F. Guzmán-Lastra, Sci. Rep., 2020, 10, 20845.
20. I. Buttinoni, J. Bialké, F. Kümmel, H. Löwen, C. Bechinger and T. Speck, Phys. Rev. Lett., 2013, 110, 238301.
21. J. Bialké, T. Speck and H. Löwen, J. Non-Cryst. Solids, 2015, 407, 367–375.
22. M. E. Cates and J. Tailleur, Annu. Rev. Condens. Matter Phys., 2015, 6, 219–244.
23. A. Zöttl and H. Stark, Annu. Rev. Condens. Matter Phys., 2023, 14, 109–127.
24. M. Paoluzzi, M. Leoni and M. C. Marchetti, Soft Matter, 2020, 16, 6317–6327.
25. W. Zhong, Y. Deng and D. Xiong, Soft Matter, 2023, 19, 2962–2969.
26. J. Stenhammar, R. Wittkowski, D. Marenduzzo and M. E. Cates, Phys. Rev. Lett., 2015, 114, 018301.
27. S. Edlund, J. Kaufman, J. Lessler, J. Douglas, M. Bromberg, Z. Kaufman, R. Bassal, G. Chodick, R. Marom, V. Shalev, Y. Mesika, R. Ram and A. Leventhal, Epidemics, 2011, 3, 135–142.
28. G. Bucyibaruta, C. B. Dean and M. Torabi, Infect. Dis. Model., 2023, 8, 471–483.
29. F. Nill, Chaos, Solitons Fractals, 2023, 173, 113678.
30. S. Strogatz, Nonlinear Dynamics and Chaos, CRC Press, 2nd edn, 2015.
31. G. Valentini, SIR Epidemic Spread Model, 2020, https://www.mathworks.com/matlabcentral/fileexchange/75100-sir-epidemic-spread-model.
32. A. Patch, D. Yllanes and M. C. Marchetti, Phys. Rev. E, 2017, 95, 012601.
33. A. Patch, D. M. Sussman, D. Yllanes and M. C. Marchetti, Soft Matter, 2018, 14, 7435–7445.
34. J. Stenhammar, D. Marenduzzo, R. J. Allen and M. E. Cates, Soft Matter, 2014, 10, 1489–1499.
35. M. Paoluzzi, N. Gnan, F. Grassi, M. Salvetti, N. Vanacore and A. Crisanti, Sci. Rep., 2021, 11, 24467.
36. R. N. Mantegna and H. E. Stanley, Phys. Rev. Lett., 1994, 73, 2946.
37. M. G. M. Gomes, M. U. Ferreira, R. M. Corder, J. G. King, C. Souto-Maior, C. Penha-Gonçalves, G. Gonçalves, M. Chikina, W. Pegden and R. Aguas, J. Theor. Biol., 2022, 540, 111063.
38. G. S. Redner, M. F. Hagan and A. Baskaran, Phys. Rev. Lett., 2013, 110, 055701.
39. D. Levis, J. Codina and I. Pagonabarraga, Soft Matter, 2017, 13, 8113–8119.
40. J. Stenhammar, A. Tiribocchi, R. J. Allen, D. Marenduzzo and M. E. Cates, Phys. Rev. Lett., 2013, 111, 145702.
41. J. Bialké, H. Löwen and T. Speck, Europhys. Lett., 2013, 103, 30008.
42. S. C. Takatori and J. F. Brady, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91, 032117.
43. V. M. Worlitzer, G. Ariel, A. Be’er, H. Stark, M. Bär and S. Heidenreich, New J. Phys., 2021, 23, 033012.
44. S. C. Takatori and J. F. Brady, Soft Matter, 2015, 11, 7920–7931.
45. D. Rogel Rodriguez, F. Alarcon, R. Martinez, J. Ramírez and C. Valeriani, Soft Matter, 2020, 16, 1162–1169.
46. N. K. Agrawal and P. S. Mahapatra, Phys. Rev. E, 2021, 104, 044610.
47. J. Zhang, T. Huang, G. Xu and Y. Chen, Commun. Theor. Phys., 2022, 74, 075601.
48. J. Su, M. Feng, Y. Du, H. Jiang and Z. Hou, Commun. Phys., 2023, 6, 58.
49. R. Kailasham and A. S. Khair, Soft Matter, 2023, 19, 7764–7774.
50. Y. Fily and M. C. Marchetti, Phys. Rev. Lett., 2012, 108, 235702.
51. C. Desgranges, M. E. Tuckerman, M. Ferrari, P. M. Chaikin, S. Sacanna and J. Delhommelle, Soft Matter, 2023, 19, 7334–7342.
52. G. Gonnella, D. Marenduzzo, A. Suma and A. Tiribocchi, C. R. Phys., 2015, 16, 316–331.
53. M. Knezevic and H. Stark, Sci. Rep., 2022, 12, 19437.
54. L. Theeyancheri, S. Chaki, T. Bhattacharjee and R. Chakrabarti, Phys. Rev. Res., 2024, 6, L012038.
55. C. B. Caporusso, L. F. Cugliandolo, P. Digregorio, G. Gonnella, D. Levis and A. Suma, Phys. Rev. Lett., 2022, 131, 068201.
56. C. B. Caporusso, P. Digregorio, D. Levis, L. F. Cugliandolo and G. Gonnella, Phys. Rev. Lett., 2020, 125, 178004.
57. K. Teknomo, Transp. Res. F: Traffic Psychol. Behav., 2006, 9, 15–27.
58. J. C. Miller, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 010101(R).
59. J. A. Anderson, J. Glaser and S. C. Glotzer, Comput. Mater. Sci., 2020, 173, 109363.
60. R. Kailasham and A. S. Khair, Contagion dynamics in active Brownian particle systems, 2024, https://github.com/khair-group/abp_sirs_contagion_dynamics.
61. A. Stukowski, Modell. Simul. Mater. Sci. Eng., 2010, 18, 015012.

---

Footnote

† Present address: Department of Chemical Engineering, Indian Institute of Technology Indore, Khandwa Road, Simrol, Madhya Pradesh 453552, India.

---