Using active matter to introduce spatial heterogeneity to the susceptible infected recovered model of epidemic spreading

The widely used susceptible-infected-recovered (S-I-R) epidemic model assumes a uniform, well-mixed population, and incorporation of spatial heterogeneities remains a major challenge. Understanding failures of the mixing assumption is important for designing effective disease mitigation approaches. We combine a run-and-tumble self-propelled active matter system with an S-I-R model to capture the effects of spatial disorder. Working in the motility-induced phase separation regime both with and without quenched disorder, we find two epidemic regimes. For low transmissibility, quenched disorder lowers the frequency of epidemics and increases their average duration. For high transmissibility, the epidemic spreads as a front and the epidemic curves are less sensitive to quenched disorder; however, within this regime it is possible for quenched disorder to enhance the contagion by creating regions of higher particle densities. We discuss how this system could be realized using artificial swimmers with mobile optical traps operated on a feedback loop.


INTRODUCTION
Disease propagation through a heterogeneous environment has become a topic of worldwide interest. Tremendous modeling resources have been applied in efforts to control or at least predict the progress of the global pandemic. The majority of these models have as their basis the conceptually simple yet physically rich compartmentalized susceptible-infected-removed (S-I-R) representation of temporal disease evolution introduced nearly a century ago by Kermack and McKendrick [1]. Under the fundamental simplifying assumption of a mean-field, well-mixed population, in the S-I-R model the population is divided into S (susceptible), I (infected), or R (recovered) individuals, and the dynamic evolution of the epidemic is governed by the transition rates between these categories: a removal rate µ for transitions I → R and an infection rate that relies on on the law of mass action to model transitions from S → I. Individuals in a given bin are indistinguishable, and all spatial details of the system are discarded [2]. Despite their apparent simplicity, S-I-R models and their many variants provide powerful tools for forecasting the general course of an epidemic. Where these models falter is in predicting the specific course of an actual real-world epidemic. This is generally attributed to the lack of homogeneity in individual susceptibility, spatial contacts, and mixing behavior of individuals [3] leading to stochastic effects that can not be averaged away.
Incorporation of heterogeneity has proven to be not at all straightforward, and numerous approaches have been developed over the years. For example, the population can be broken into subpopulations, each with different infection and recovery rates, or the population can be geographically subdivided into regions with diffusive terms to link to the regional S-I-R dynamics [4,5]. Additional heterogeneities in the diffusion can be achieved by incorporating patchiness into the diffusion [6,7]. Much work has been done on connecting individuals via finitedimensional networks rather than through an infinitedimensional mean field [2]; however, the details of the network itself make the problem even more complex since decisions must be made on what is the appropriate degree distribution for the network connectivity as well as whether the network should remain static or should be allowed to evolve either independently or in response to the progress of the disease [8,9]. The impact of heterogeneity in transmission and susceptibility is discussed in [10].
The epidemic model with the ultimate heterogeneity treats each individual as a separate, mobile, interacting unit. Under Agent Based Modeling (ABM), also known as Individual Based Modeling [11], heterogeneity can be included at all levels ranging from varied individual susceptibility and recovery rates, varied contacts between individuals, spatial clustering of individuals in cities or at attractive sites such as bars, and both short and long range transport of individuals such as by bus or airplane [12,13]. The flexibility of these models is also their greatest weakness, since in addition to the computational challenge of tracking potentially millions of individuals on a country-wide scale, there can be a vast number of free parameters that must be painstakingly fitted to real-world data that is not always available at the necessary resolution.
There have been surprisingly limited efforts to address a middle ground of ABM in which many but not all of the details are abstracted away to produce a model that captures spatial heterogeneity in a meaningful way without being swallowed in a proliferation of parameters. This can, in principle, be achieved either by developing more complex analytical models or simpler simulation-based models. One of the earliest approaches for simplifying simulation-based models involved cellular automata, where the mobility of individual agents could be varied up to a level consistent with the mean-field limit [14]. Individuals obeying S-I-R interactions have also been represented as moving particles that are driven and diffusing [15], that never change direction [16], that occasionally make long-range jumps [17], that move at different velocities [18], or that are confined to diffuse only within the region of their 'houses' [19]. To help mitigate the computational expense of such methods, dynamic density functional theory techniques can be applied [20].
The significant progress made during recent years in understanding what are known as active matter models [21,22], where individual particles are self-propelled and interact with each other on a spatial landscape that may or may not include disorder, suggests the natural step of pairing a model of S-I-R type with active particles. The active particles can be of run-and-tumble type [23] or driven diffusive [24]. In a small system of low density, an active matter assembly was able to reproduce the mean field behavior of S-I-R [25]. Generally, however, there has been only limited work on coupling S-I-R modeling with active matter. For example, Paoluzzi et al. considered S-I-R type dynamics to examine information exchange in active clustering transitions [24] but not aspects of the epidemic spreading itself. Recently Zhao et al. studied contagion dynamics in self-propelled flocking models and found that ordered homogeneous states reduce disease spreading while bands and clustering favor the spreading [26].
There are a number of advantages to working with an active matter system. The well-known motility-induced phase separation (MIPS) transition from a low density gas phase to a coexistence between high and low density regions as a function of density and/or mobility of the active particles [27][28][29][30] can provide a natural separation of the particles into clustered communities connected by disordered transport pathways. Contacts between particles can be viewed as an adaptive network that may be tuned to evolve on the same or a different time scale as the progression of the disease. Spatial heterogeneity emerges automatically in the MIPS regime, but can also be inserted using walls, traps, or obstacles. Disease dynamics in such systems can be abstracted by tracking the evolution of the number of S(t), I(t), and R(t) over time, whose temporal behavior will capture the impact of heterogeneities that are averaged out by the mean field approximations of the standard S-I-R model.
In this work, we simulate a large assembly of active matter particles in the MIPS regime where a giant cluster spontaneously emerges. We combine this model with S-I-R interactions in which all particles are initially susceptible (S) but can be infected with probability β when they come into direct contact with an infected (I) particle. Infected particles spontaneously transition to the recovered (R) state at rate µ, and no reinfection is allowed. We study the evolution of epidemics as the ratio of β/µ is varied, and consider the impact on the behavior of adding quenched disorder in the form of immobile obstacles. Increasing the number of immobile obstacles in an active matter system will increase the number of clusters and decrease their sizes. By performing large numbers of realizations, we find that inclusion of quenched disorder increases the number of "failed" outbreaks for small β/µ and increases the average duration of successful epidemics. When β/µ is sufficiently large, the system becomes insensitive to the presence of quenched disorder and approaches the mean field limit, and in this regime the epidemic propagates via spatially well-defined fronts. We also study the average number of susceptible particles surrounding an infective as a function of time, and find that this quantity is modestly altered by the addition of quenched disorder in the mean field limit of high β/µ but becomes strongly affected by quenched disorder as β/µ is reduced.
Our results indicate that for low β/µ, the homogeneous mixing hypothesis breaks down, that is, the infection process departs from mass action and the system becomes much more sensitive to spatial quenched disorder. In the high β/µ regime, the mixing hypothesis is more applicable even though the epidemic is spreading via spatially localized fronts. This implies that localized epidemic mitigation efforts will be more successful at low β/µ but would become ineffective in higher β/µ regimes unless applied to the entire population.
Finally, we discuss how the system we consider could be realized experimentally using feedback control of light activated colloids, where the active behavior of the colloids can be controlled on the individual level. Experiments on this type of system have already demonstrated group formation, responsive states, and predator-prey model realizations [31][32][33]. There are also numerous possible ways to introduce spatial heterogeneities in active matter systems [34][35][36][37][38][39]. Techniques of this type could be used to mimic the S-I-R model with and without spatial disorder. This could permit the creation of table-top epidemic spreading models with active matter.

Modeling and characterization of the S-I-R dynamics
We simulate N = 5000 active particles in a twodimensional system of size L × L where L = 200 and where there are periodic boundary conditions in the x and y directions. The motion of the particles is obtained by integrating the following equation: Here v i = dr i /dt is the velocity and r i is the position of particle i, and the damping constant α d = 1.0. The interaction between two particles is represented with a harmonic repulsive potential F dd i = N i =j k(2r a − |r ij |)Θ(|r ij |−2r a )r ij , where Θ is the Heaviside step function, r ij = r i − r j , andr ij = r ij /|r ij |. We set the spring force to k = 20 and the particle radius to r a = 1.0. Each particle is subjected to a motor force F m i = F Mmi of magnitude F M applied in a randomly chosen direction m during a run time of τ l before instantaneously changing to a new randomly chosen direction during the next run time, producing run-and-tumble dynamics. For each particle, we fix τ l to a value selected randomly from the interval 1.5×10 4 to 3.0×10 4 , and we set F M = 1.5 for susceptible and recovered particles, placing us in the MIPS regime in the absence of quenched disorder [37]. Infected particles have their motor force reduced to F M = 1.0. For some simulations, we include quenched disorder in the form of N obs = 800 obstacles that produce the force F obs . This is taken to be the same as the particle-particle interaction force, with the only difference being that the obstacles are immobile. An image of the system in the presence of obstacles appears in Fig. 1.
Each active particle carries a label marking it as being in one of three states: S, I, or R. If an S particle comes into direct contact with an I particle, for each time step during which the contact persists there is a probability β that the S particle will transition to an I particle. If at a given time step an S particle is in contact with n I particles, the probability of infection is 1−(1−β) n ≈ nβ. Transitions of I particles to state R occur with probability µ at each time step regardless of the state of any particles that may be in contact with the I particle. Thus, the mean time spent in the infected state is 1/µ time steps. The R state is absorbing and R particles experience no further state transitions. In this S-I-R model, the infected I particles are present only as a transient and the system will eventually contain only S and/or R particles. We note that the mean-field rates governing S → I and I → R transitions and determining the basic reproductive number R 0 in classic S-I-R models do not map directly to the values of β and µ that we insert into our model as microscopic parameters. In ABMs, the effective mean-field rates are emergent quantities instead of control parameters.
To initialize the system, we place the particles randomly in the sample and set them all to state S. We allow the system to evolve for 5 × 10 5 simulation time steps until a stable MIPS giant cluster has formed, and define this state to be the t = 0 condition. We then randomly select 5 particles and change their state to I. These particles serve as our index cases, and we choose 5 rather than 1 in order to lower the probability of a failed outbreak. The system continues to evolve under both the motion of the particles and the reactions between states S, I, and R until no I particles remain. We perform 1000 realizations for each parameter set. Since, as is shown in the results, the duration t d = min{t > 0 : I(t) = 0} of an individual epidemic can vary significantly from run to run, we report time in terms of the scaled quantitỹ t = t/t d . As a function of scaled time, we measure the epidemic curves s(t) = S(t)/N , i(t) = I(t)/N , and r(t) = R(t)/N . This rescaling enables us to visually compare features of the progression of the epidemic when we change the ratio β/µ. We also measure the peak infective fraction i max and the final susceptible fraction s ∞ , both of which are commonly used indicators of the severity of an epidemic. To obtain further information on the spatial evolution of the system, we measure the average number of susceptible particles surrounding an infective, where I denotes the indicator function, and the sums over i and j range over the infected and susceptible particles, respectively. For a two-dimensional system of particles with identical radii r a , η cannot exceed the maximum coordination number of z = 6. If the infected individuals are well mixed within the population, the average number of susceptible particles surrounding an infective is η(t) ∝ S(t); more specifically, we would expect η(t) ∝ z c S(t)/N , where z c is the average coordination number of the particles. Departure from this behavior is indicative of a failure of the homogeneous mixing assumption. Low Transmissibility Regime In Fig. 2(a) we show a snapshot of the system at β/µ = 0.5 in the low transmissibility regime in the absence of quenched disorder. The moving particles form a phase separated state of a high density solid and a low density gas. As indicated in the introduction, the relationship between β/µ and the basic reproductive number is an emerging quantity. Since within a cluster the expected number of contacts is z = 6, the expected number of secondary cases from an index case within the cluster will be η = 3, showing that the epidemic will infect a fraction of the cluster. If the index case starts in the gaseous phase, its expected number of contacts is likely z < 1, implying a reproductive number less than one and giving limited cluster-to-cluster transmissions.
The same system in the presence of randomly placed obstacles appears in Fig. 2(b), where the giant dense cluster is now accompanied by numerous smaller persistent clusters that have nucleated around some of the obstacle sites. Since particles within a cluster are locked to one another, the homogeneous mixing assumption fails to hold for the epidemic dynamic within each cluster. Thus by controlling the number and size of the clusters, we explore a range of departures from the mixing assumption, from the most extreme situation when there is only one large cluster, to greater mixing as the number of clusters increases and their sizes decrease.
In Fig. 2(a,b,c) we illustrate the evolution of the S, I, and R particles for the obstacle-free system at times oft = 0.2, 0.3, and 0.4, while in Fig. 2(d,e,f) we show the evolution in the system containing N obs = 800 obstacles. In both cases, when the giant cluster is contacted by an infective, the disease spreads through the cluster, but since the probability of transmission is low, not all of the S particles surrounding a given I become infected, and as a result, a finite number of S remain when the epidemic is complete. See Ref. [40] for a discussion on final epidemic size. When we add quenched disorder to the system, shown as black circles in Fig. 2(d-f), a greater amount of localized clustering occurs in addition to the giant cluster. Since each cluster must be infected separately, this tends to slow the spread of the infection and reduce the peak infective fraction i max , as shown in Fig. 2(e).
Although the dynamics of the spread of the infection is similar with and without quenched disorder, att = 0.4 the number of I particles present is much lower when obstacles have broken the system into smaller clusters, indicating that the epidemic has impacted fewer particles in the system with quenched disorder.
In Fig. 3(a) we plot the epidemic curves showing the fractions of susceptible s, infected i, and recovered r particles as a function of reduced timet for samples with and without quenched disorder. We note that in the pres- ence of obstacles, the duration of the epidemic tends to be longer; however, by plotting the epidemic curves as a function of reduced time it is easier to compare samples with and without quenched disorder. The curves have the shapes expected from the classic S-I-R model. In the system without obstacles, by the end of the epidemic there is still a fraction s ∞ = 0.41 of the population that never became infected, while r ∞ = 1 − s ∞ = 0.59 of the particles have recovered. When obstacles are present, a larger fraction s ∞ = 0.51 of particles have escaped infection. The peak i max in the infected fraction is also considerably reduced in magnitude when obstacles are present. This indicates that the system is sensitive to the presence of spatial heterogeneities introduced by the clustering arising from the presence of fixed obstacles. Within this regime, spatially localized mitigation protocols could be effective, since local quenched disorder can slow the overall mobility of the particles or reduce the effective connectivity among the particles. To further demonstrate the effect of adding obstacles, in Fig. 3(b) we plot η, the average number of S particles surrounding an I particle, versust. Here η is always smaller in the sample containing obstacles.

High Transmissibility Regime
We next consider the case of high transmissibility β/µ = 5.0. In Fig. 4(a,b,c) we plot the spatial evolution of the susceptible, infected and recovered particles in the absence of obstacles. The infection spreads via well defined fronts through the dense region. In Fig. 4(d,e,f) we show the same evolution in the presence of obstacles. There are now multiple dense clusters present, but in each a similar front propagation of the infection appears. In Fig. 5(a) we plot s(t), i(t), and r(t) for the high transmissibility system with β/µ = 5.0 from Fig. 4. Here, s ∞ = 0 and all of the particles become infected regardless of whether obstacles are present. The peak value i max is nearly the same for both cases. An interesting effect appears in which fort < 0.175, adding obstacles depresses i, but fort > 0.185, adding obstacles increases i. This is not merely due to a change in the duration of the epidemic since the curves are plotted in reduced time; instead, it indicates a change in the spatial propagation of the infection, which we will address in Figs. 7 and 8. The crossover in behavior occurs after the initial large infection front has completely swept through either the giant cluster or all of the smaller clusters for the samples with quenched disorder. In Fig. 5(b) we plot the corresponding η versust, which is nearly unchanged by the inclusion of obstacles. These results indicate that under high transmissibility, the system is less sensitive to spatial disorder and the behavior is consistent with the mean field limit. Note that epidemic curves and plots of η(t) for all other values of β/µ can be viewed in the The corresponding η, the average number of S particles surrounding an I particle, versust in the samples without (blue) and with (orange) obstacles. There is a minimal difference in η between the two cases.

Supplemental Material. Duration of Epidemic
Our simulations reveal a strong stochasticity of the behavior. Depending on the particular randomly chosen locations of the index cases, the duration t d of the epidemic can vary widely. In particular, for some realizations the outbreak fails to take hold and is extinguished without affecting a significant fraction of the particles. To illustrate this, in Fig. 6 we plot the distribution P (t d ) of the epidemics measured in simulation time steps with and without obstacles for 1000 realizations. In Fig. 6(a-e), we show the low transmissibility regime with β/µ = 0.4, 0.45, 0.5, 0.6, and 1.0. Here the distribution is bimodal and there is a clear division between small t d , where we find failed outbreaks that do not infect a significant fraction of the particles, and larger t d , where successful epidemics occur that involve a substantial fraction of the particles. This behavior is similar to what has been observed in other studies [41]. Addition of quenched disorder in this regime increases the probability that the outbreak will fail, but also increases the average duration of successful epidemics. In contrast, for high transmissibility, as shown in Fig. 6(f,g) at β/µ = 2.0 and 3.0, P (t d ) is unimodal since all outbreaks produce successful epidemics. Additionally, there is no longer a significant difference in the distribution for systems with and without quenched disorder.

Ability of I to Contact S
We can also distinguish the two regimes of behavior us-t d 0.00 0.05 ing features in η by comparing the value of η in samples with and without quenched disorder. In Fig. 7(a) we plot η versust in samples with and without obstacles. The time scalet reaches a value oft = 1.0 when the number of recovered has increased to 95% of its maximum value, r(t = 1.0) = 0.95r ∞ = 0.95(1 − s ∞ ). Use of this time scale allows us to exclude the stochastic late time behavior when the last few straggling infectives are recovering. Att = 0, η is always high since the initial seed I particles are surrounded only by susceptible particles. As the epidemic spreads, the average number of S particles around I particles decreases. When β/µ ≤ 1.5, the curves monotonically decrease to a saturation value between η = 2. (a) η vst for varied β/µ in samples without obstacles (thin lines) and with obstacles (thick lines). When β/µ is large, a local minimum in η appears neart = 0.1 due to the formation of a propagating front. (b) The difference ∆η between the value of η in the sample with disorder and the value in the sample without disorder as a function oft . For β/µ ≤ 1.5, there is no front propagation and the addition of quenched disorder always reduces the value of η. For β/µ > 1.5, a front appears, and once the front has passed, ∆η drops below zero, indicating an enhancement of the infection rate when quenched disorder is present.
form of a front, which is visible as the appearance of a local dip in η centered neart = 0.2. As the front moves rapidly through the largest cluster, most of the infected particles are surrounded by I particles behind the expanding front, while only I particles at the edge of the front are adjacent to S particles. This depresses the value of η. Once the front has passed through the cluster, the mobility of the particles bring more S from the gas phase into contact with the remaining I, and η recovers somewhat before saturating to a low value between η = 1.5 and 2.0. In Fig. 7(b) we plot the difference ∆η = η obs=800 − η obs=0 between η for the samples with and without quenched disorder from Fig. 7(a). When β/µ ≤ 1.5, ∆η reaches a constant value and is always negative, ∆η ≈ −0.25, indicating that the quenched disorder is always reducing the effectiveness of the spread of the epidemic. Once the system enters the front propagation regime for β/µ > 1.5, ∆η becomes nonmonotonic and shows local peaks and dips. Fort < 0.4, there is a dip when the front is passing through the largest clusters. In this regime, ∆η is negative, indicating that the quenched disorder slows the front to some extent. For times above the minimum of the dip, ∆η increases and becomes positive, indicating that the addition of quenched disorder is actually increasing the effectiveness of the epidemic spread. This can also be seen in Fig. 5(a), where i is reduced in the presence of quenched disorder fort < 0.175 but is slightly increased fort > 0.175, indicating that the disorder can accelerate the infection at later times. This enhancement of the epidemic arises after the largest cluster has become fully infected and some of the infected particles break away from the cluster and enter the gas phase. Within the gas phase, the quenched disorder induces the formation of smaller localized clusters, as shown in Fig. 1. These smaller clusters, once contacted by an infective, undergo the same rapid front propagation as the initial infection wave. When quenched disorder is not present, there are no smaller clusters and the infection must propagate through the gas phase and infect the remaining S particles one by one, an inefficient process. Epidemic Phase Diagram Based on the features in Fig. 7, we can construct a phase diagram of the behavior of the system as a function of β/µ versust , as shown in Fig. 8. When β/µ > 1.5, s ∞ = 0 and the entire system becomes infected, while the initial invasion of the infection occurs by front propagation. For β/µ ≤ 1.5, the low transmissibility regime marked LT, the infection spreads much more homogeneously, as illustrated in Fig. 2, and s ∞ > 0 so that not all of the particles have been infected by the end of the epidemic. Within this regime, the addition of quenched disorder always reduces i max and increases s ∞ . In the high transmissibility regime with β/µ > 1.5, at smallt the infection propagates as a front through the largest cluster, defined as the front propagation regime FP. Here the addition of quenched disorder can slow the front propagation but does not stop it. After the front has crossed the entire largest cluster, we find the CP regime in which the secondary clusters start to show front propagation. Here the quenched disorder can increase the effectiveness of the spread of the infection by increasing the number of secondary clusters that are present. At larger values oft , all of the clusters have been infected and the epidemic is making its way through the gas phase. In this regime, which we call Diff for diffusive, there is little difference between the systems with and without quenched disorder. The locations of the phase boundaries should depend on the amount of quenched disorder and the activity level of the particles.

DISCUSSION
Our model suggests that active matter can be used to capture a variety of epidemic behaviors. There are a number of active matter systems, such as active colloids, in which the activity of the particles can be controlled on an individual basis using optical rastering methods. Experiments of this type have been developed in order to use active colloids to mimic group formation, to introduce an effective visual perception mobility, and to produce other kinds of collective behaviors such as quorum sensing [31,32]. In order to implement an S-I-R model, individual active colloids could be traced and tagged according to their infective state, and when they interact with other colloids, there can be a probability that the infection will pass to a susceptible colloid. This can be done in a motility induced phase separated regime or in a diffusive regime for varied β/µ. The experiments could then be repeated multiple times to obtain the average behavior. Within a given sample, certain colloids could remain inactive and be counted as passive or obstacle particles, or actual obstacles could be put in place on the substrate. Additionally, a wealth of rules could be introduced, such as the inclusion of hyperactive particles that could serve as superspreaders, as well as possible mitigation effects. This could potentially position active matter as a table-top experimental system for modeling epidemics. Our work indicates that active matter can be used as a simulation tool to study epidemics in a system that can be tuned readily between states that are sensitive to spatial disorder and states that are insensitive to disorder.
In conclusion, we have shown how an active matter system of self-propelled particles can be used to model spatial heterogeneity in an S-I-R epidemic spreading model. The active particles naturally form spatial clusters in the motility induced phase separated regime. For low transmissibility, the epidemic spread is percolative and the system is sensitive to the addition of quenched disorder, which both increases the probability of a failed outbreak and increases the average duration of successful epidemics. In this regime, the mixing hypothesis of classical S-I-R models breaks down. For high transmissibility, all of the particles are eventually infected and the epidemic spreads in well defined fronts. In this case, the addition of quenched disorder can slow the spread of the epidemic at early times by slowing the propagation of the initial front. At later times, however, since the quenched disorder introduces a larger number of small clusters in the gas phase, the epidemic can spread more efficiently compared to a system without quenched disorder. Our results indicate that spatial disorder can impact epidemic spreading in both the high and low transmissibility regimes. Our system could be realized experimentally using light activated colloidal particles with specified feedback rules to mimic the S-I-R model, and our results suggest that active matter systems could provide a new way to create table-top epidemic experiments.
Supplemental material for "Using Active Matter to Introduce Spatial Heterogeneity to the Susceptible-Infected-Recovered Model of Epidemic Spreading" Here we present the full epidemic curves of susceptible s(t), infected i(t), and recovered r(t) along with the corresponding η(t), the average number of S particles surrounding an I particle, for additional values of β/µ that were not included in the main text. Each curve is averaged over 1000 realizations. Att = 1.0, the epidemic is over and i = 0.