The COVID-19 pandemic constitutes a severe global health crisis. Many countries have implemented stringent non-pharmaceutical control measures that involve social distancing and mobility reduction in their populations. This has led to remarkably successful deceleration and significant "flattening of the curve" of the infection outbreaks, albeit at tremendous economic and financial costs. At this point, societies are in dire need of designing a secure (partial) exit strategy wherein the inevitable recurrence of the infection among the significant non-immune fraction of the population can be thoroughly monitored with sufficient spatial resolution and reliable statistics, provided that dependable, frequent, and widespread virus testing capabilities are accessible and implemented. Until an effective and safe vaccine is widely available, this would ideally allow the localized implementation of rigorous targeted disease control mechanisms that demonstrably protect people's health while the paralyzed branches of the economy are slowly rebooted.
Mathematical analysis and numerical simulations of infection spreading in generic epidemic models are crucial for testing the efficacy of proposed mitigation measures, and the timing and pace of their gradual secure removal. Specifically, the employed mathematical models need to be (i) stochastic in nature in order to adequately account for randomly occurring or enforced disease extinction in small isolated communities, as well as for rare catastrophic infection boosts and (ii) spatially resolved such that they properly capture the significant emerging correlations among the susceptible and immune subpopulations. These distinguishing features are notably complementary to the more detailed and comprehensive computer models utilized by researchers at the University of Washington, Imperial College London, the Virginia Bioinformatics Institute, and others: see, e.g., (1) (2) (3) (4) (5) (6) .
We report a series of detailed individual-based kinetic Monte Carlo computer simulation studies for stochastic variants (7, 8) of the paradigmatic Susceptible-Infectious-Recovered (SIR) model (9, 10) for a community of about 100,000 individuals. To determine the robustness of our results and compare the influence of different contact characteristics, we ran our stochastic model on four distinct spatially structured architectures, namely i) regular two-dimensional square lattices, wherein individuals move slowly and with limited range, i.e., spread diffusively; ii) two-dimensional small-world networks that in addition incorporate substantial long-distance interactions and contaminations; and finally on iii) random as well as iv) scale-free social contact networks. For each setup, we investigated epidemic outbreaks with model parameters informed by the known COVID-19 data (4). To allow for a direct comparison, we extracted the corresponding effective infection and recovery rates by fitting the peak height and the half-peak width of the infection growth curves with the associated classical deterministic SIR rate equations that pertain to a wellmixed setting. We designed appropriate implementations of social distancing and contact reduction measures on each architecture by limiting or removing connections between individuals. This approach allowed us to generically assess the efficacy of non-pharmaceutical control measures.
Although each architecture entails varied implementations of social distancing measures, we find that they all robustly reproduce both the resulting reduced outbreak intensity and growth speed. As anticipated, a dramatic resurgence of the epidemic occurs when mobility and contact restrictions are released too early. Yet if stringent and sufficiently long-lasting social distancing measures are imposed, the disease may go extinct in the majority of isolated small population groups. In our spatially extended lattice systems, disease spreading then becomes confined to the perimeters of a few larger outbreak regions, where it can be effectively localized and specifically targeted. For the small-network architecture, it is however imperative that all longrange connections remain curtailed to a very low percentage for the control measures to remain effective. Intriguingly, we observe that an infection outbreak spreading through a static scale-free network effectively randomizes its connectivity for the remaining susceptible nodes, whence the second wave encounters a very different structure.
In the following sections, we briefly describe the methodology and algorithmic implementations as well as pertinent simulation results for each spatial or network structure; additional details are provided in the Supplementary Materials. We conclude with a comparison of our findings and a summary of their implications.
Our first architecture is a regular two-dimensional square lattice with linear extension ‫ܮ‬ = 448 subject to periodic boundary conditions (i.e., on a torus). Initially, ܰ = ܵሺ0ሻ + ‫ܫ‬ሺ0ሻ + ܴሺ0ሻ = 100,000 individuals with fixed density ߩ = ‫ܮ/ܰ‬ ଶ ≈ 0.5 are randomly placed on the lattice, with at most one individual allowed on each site. Almost the entire population begins in the susceptible state ܵሺ0ሻ; we start with only 0.1 % infected individuals, ‫ܫ‬ሺ0ሻ = 100, and no recovered (immune) ones, ܴሺ0ሻ = 0. All individuals may then move to neighboring empty lattice sites with diffusion rate ݀ (here we set this hopping probability to 1). Upon their encounter, infectious individuals irreversibly change the state of neighboring susceptible ones with set rate ‫:ݎ‬ ܵ + ‫ܫ‬ → ‫ܫ‬ + ‫.ܫ‬ Any infected individual spontaneously recovers to an immune state with fixed rate ܽ: ‫ܫ‬ → ܴ. (Details of the simulation algorithm are presented in the Supplementary Materials.) For the recovery rate, we choose 1/ܽ ≅ 6.667 days (set to 0.15 Monte Carlo steps, MCS) informed by known COVID-19 characteristics (4). To determine the infection rate ‫,ݎ‬ we run simulations for various values, fit the peak height and width of the ensuing epidemic curves with the corresponding SIR rate equations to extract the associated basic reproduction ratio ܴ (as explained in the Supplementary Materials, see Figure S1 ), and finally select that value for ‫ݎ‬ for our individual-based Monte Carlo simulations that reproduces the ܴ ≈ 2.4 for COVID-19 (4). We perform 100 independent simulation runs with these reaction rates, from which we obtain the averaged time tracks for ‫ܫ‬ሺ‫ݐ‬ሻ and ܴሺ‫ݐ‬ሻ, while of course ܵሺ‫ݐ‬ሻ = ܰ − ‫ܫ‬ሺ‫ݐ‬ሻ − ܴሺ‫ݐ‬ሻ and ܴሺ‫ݐ‬ሻ = ܽ ‫‬ ‫ܫ‬ሺ‫′ݐ‬ሻ ‫′ݐ݀‬ ௧ .
The standard classical SIR deterministic rate equations assume a well-mixed population and constitute a mean-field type of approximation wherein stochastic fluctuations and spatial as well as temporal correlations are neglected; see, e.g., (11, 12) . Near the peak of the epidemic outbreak, when many individuals are infected, this description is usually adequate, albeit with coarse-grained `renormalized' rate parameters that effectively incorporate fluctuation effects at short time and small length scales. However, the mean-field rate equations are qualitatively insufficient when the infectious fraction ‫ܫ‬ሺ‫ݐ‬ሻ/ܰ is small, whence both random number fluctuations and the underlying discreteness and associated internal demographic noise become crucial (11) (12) (13) . Already near the epidemic threshold, which constitutes a continuous dynamical phase transition far from thermal equilibrium, c.f. Figure S3 in the Supplementary Materials, the kinetics is dominated by strong critical point fluctuations. These are reflected in characteristic initial power laws rather than simple exponential growth of the ‫ܫ‬ሺ‫ݐ‬ሻ and ܴሺ‫ݐ‬ሻ curves (14) , as demonstrated in Figure S1 (Supplemental Information).
Nor can the deterministic rate equations capture stochastic disease extinction events that may occur at random in regions where the infectious concentration has reached small values locally. The rate equations may be understood to pertain to a static and fully connected network; in contrast, the spreading dynamics on a spatial setting continually rewires any infectious links keeping the epidemic active (6, 15) . Consequently, once the epidemic outbreak threshold is exceeded, the SIR rate equations markedly underestimate the time-integrated outbreak extent reflected in the ultimate saturation level ܴ ஶ = ܴሺ‫ݐ‬ → ∞ ሻ, as is apparent in the comparison Figure  S1 (Supplemental Information).
Once the instantaneous fraction of the population has reached the threshold 10 %, ‫ܫ‬ሺ‫ݐ‬ሻ = 0.1 ܰ, we initiate stringent social distancing that we implement through a strong repulsive interaction between any occupied lattice sites (with ݊ = 1), irrespective of their states ܵ, ‫,ܫ‬ or ܴ; and correspondingly an attractive force between filled and empty (݊ = 0) sites, namely the Ising lattice gas potential energy ܸሺ{݊ }ሻ = ‫ܭ‬ ∑ ሺ2 ݊ − 1ሻ ൫2 ݊ − 1൯ ழ,வ with dimensionless strength ‫ܭ‬ = 1, where the sum extends only over nearest-neighbor pairs on the square lattice. The transfer of any individual from an occupied to an adjacent empty site is subsequently determined through the ensuing energy change ∆ܸ by the Metropolis transition probability ‫ݓ‬ = min{1, exp ሺ−∆ܸሻ} (16, 17) , which replaces the unmitigated hopping rate ݀. As a result, both the mobility as well as any direct contact between individuals on the lattice are quickly and drastically reduced. For sufficiently small total density ߩ = ‫ܮ/ܰ‬ ଶ , most of the individuals eventually become completely isolated from each other. For our ߩ = 0.5, the disease will continue to spread for a short period, until the repulsive potential has induced sufficient spatial anti-correlations between the susceptible individuals. The social-distancing interaction is sustained for a time duration ܶ, and then switched off again. With increasing mitigation duration ܶ, the likelihood for the disease to locally go extinct in isolated population clusters grows markedly. As seen in the bottom row, the prevalence and spreading of the infection thus becomes confined to the perimeters of a mere few remaining centers. Hence we observe drastically improved mitigation effects for extended ܶ: As shown in Figure 2 , the resurgence peak in the ‫ܫ‬ሺ‫ݐ‬ሻ curve assumes markedly lower values and is reached after much longer times. In fact, the time ߬ሺܶሻ for the infection outbreak to reach its second maximum increases exponentially with the social-distancing duration, as evidenced in the inset of Figure 2 (see also Figure 6 below). We emphasize that localized disease extinction and spatial confinement of the prevailing disease clusters represent correlation effects that cannot be captured in the SIR mean-field rate equation description. 
In modern human societies, individuals as well as communities feature long-distance connections that represent `express' routes for infectious disease spreading in addition to short-range links with their immediate neighbors. To represent this situation, we extend our regular lattice with diffusive propagation to a two-dimensional Newman-Watts small-world network (18) , which was previously applied to the study of plant disease proliferation (19) . In contrast to the Watts-Strogatz model (20) , in which the small-world property is generated through rewiring bonds of a onedimensional chain of sites, a Newman-Watts small-world network may be constructed as follows: For each nearest-neighbor bond, a long-distance link (or `short-cut') is added with probability ߮ between randomly chosen pairs of vertices. As illustrated in Figure S2 (Supplemental Information), the resulting network features 2 ߮ ‫ܮ‬ ଶ long-distance links, with mean coordination number < ݇ > = 4 ሺ1 + ߮ሻ.
Again, each vertex may be in either of the states ܵ, ‫,ܫ‬ ܴ, or empty, and each individual can hop to another site along any (nearest-neighbor or long-distance) link with a total diffusivity ݀. A typical snapshot of the SIR model on this small-world architecture is shown in Figure S2 (Supplemental Information). The unmitigated simulation parameters are: ‫ܮ‬ = 1,000, ܰ = 100,000, ‫ܫ‬ሺ0ሻ = 100, ݀ = 1, and ߮ = 0.6. The presence of long-range links increases the mean connectivity, rendering the population more mixed, which in turn significantly facilitates epidemic outbreaks (see Figure S4 in the Supplemental Information). We remark that for the SIR dynamics, the Newman-Watts smallworld network effectively interpolates between a regular two-dimensional lattice and a scale-free network dominated by massively connected hubs; moreover, as the hopping probability ݀ → 0, the small-world network is effectively rendered static. In the two-dimensional small-world network, we may introduce social-distancing measures through two distinct means: i) We can globally diminish mobility by adopting a reduced overall diffusivity ݀ ᇱ < 1; and/or ii) we can drastically reduce the probability of utilizing a long-distance connection to ݀ ఝ ≪ 1. We have found that the latter mitigation strategy of curtailing the infection short-cuts into distant regions has a far superior effect. Therefore, in Figure 3 we display the resulting data for such a scenario where we set ݀ ఝ = 0.05, yet kept the diffusivity unaltered at ݀ = 1; as before, this control was triggered once ‫ܫ‬ሺ‫ݐ‬ሻ = 0.1 ܰ had been reached in the course of the epidemic. The resurgence peak height and growth rate become even more stringently reduced with extended mitigation duration than for (distinct) social distancing measures implemented on the regular lattice.
Finally, we run the stochastic SIR kinetics on two different static structures, namely i) randomly connected and ii) scale-free contact networks. Each network link may be in either the ܵ, ‫,ܫ‬ or ܴ configurations, which are subject to the SIR reaction rules, but we do not allow movement among the network vertices. For the random network, we uniformly distribute 1,000,000 edges among ܰ = 100,000 nodes; this yields a Poisson distribution for the connectivity with preset mean (equal to the variance) < ݇ > = ሺ∆݇ሻ ଶ = 20. For the scale-free network, we employ the Barabasi-Albert graph construction (21) , where each new node is added successively with ݇ = 4 edges, to yield a total of 799,980 edges. The connectivity properties in these quite distinct architectures are vastly different, since the scale-free networks feature prominent `hubs' through which many other nodes are linked. In the epidemic context, these hubs represent super-spreader centers through which a large fraction of the population may become infected (8, 22) .
To implement the stochastic SIR dynamics on either contact network, we employ the efficient rejection-free Gillespie dynamical Monte Carlo algorithm: Each reaction occurs successively, but the corresponding time duration between subsequent events is computed from the associated probability function (23) (for details, see Supplemental Information). The random social network may be considered an emulation of the well-connected mean-field model. Indeed, we obtain excellent agreement for the temporal evolution of the SIR kinetics in these two systems with ܽ = 0.15 MCS (for the scale-free network, a small adjustment to an effective mean-field recovery rate ܽ ≈ 0.18 MCS is required). We implement a `complete lockdown' mitigation strategy: Once the threshold ‫ܫ‬ሺ‫ݐ‬ሻ = 0.1 ܰ has been reached, we immediately cut all links for a subsequent duration ܶ; during that time interval, only spontaneous recovery ‫ܫ‬ → ܴ can occur. In Figure 4 , we discern a markedly stronger impact of this lockdown on the intensity of the epidemic resurgence in both these static contact network architectures, see also Figure 6A below. On the other hand, the mitigation duration influences the second infection wave less strongly, with the time until its peak has been reached growing only linearly with ܶ: ߬ሺܶሻ ~ ܶ, as is visible in Figure 6B . There is however a sharp descent in resurgent peak height beyond an apparent threshold ܶ > 7/ܽ for the random network, and ܶ > 8/ܽ for the scale-free network. For both the two-dimensional regular lattice and small-world structure, a similar sudden drop in the total number of infected individuals ( Figure 6B ) requires a considerably longer mitigation duration: In these dynamical networks, the repopulation of nodes with infective individuals facilitates disease spreading, thereby diminishing control efficacy.
We remark that if a drastically reduced diffusivity ݀ ᇱ ≪ 1 is implemented, the small-world results closely resemble those for a randomly connected contact network ( Figure 6A ). Moreover, we have observed an unexpected and drastic effective structural change in the scalefree network topology as a consequence of the epidemic outbreak infecting its susceptible nodes. Naturally, the highly connected hubs are quickly affected, and through transitioning to the recovered state, become neutralized in further spreading the disease. As shown in Figure 5 , as the infection sweeps through the network (in the absence of any lockdown mitigation), the distribution of the remaining active susceptible-infectious (SI) links remarkably changes from the initial scale-free power law with exponent −1/2 to a more uniform, almost randomized network structure. The disease resurgence wave thus encounters a very different network topology than the original outbreak.
In this study, we implemented social distancing control measures for simple stochastic SIR epidemic models on regular square lattices with diffusive spreading, two-dimensional Newman-Watts small-world networks that include highly infective long-distance connections, and static contact networks, either with random connectivity or scale-free topology. In these distinct architectures, all disease spreading mitigation measures, be that through reduced mobility and/or curtailed connectivity, must of course be implemented at an early outbreak stage, but also maintained for a sufficient duration to be effective. In Figure 6 , we compare salient features of the inevitable epidemic resurgence subsequent to the elimination of social distancing restrictions, namely the asymptotic fraction ܴ ஶ /ܰ of recovered individuals, i.e., the integrated number of infected individuals; and the time ߬ሺܶሻ that elapses between the release and the peak of the second infection wave, both as function of the mitigation duration ܶ. We find that the latter grows exponentially with ܶ on both dynamical lattice architectures, but only linearly on the static networks ( Figure 6B) . Furthermore, as one would expect, the mean-field rate equations pertaining to a fully connected system describe the randomly connected network very well. In stark contrast to the mean-field results (indicated by the purple lines in Figure 6 ), the data for the lattice and network architectures reveal marked correlation effects that emerge at sufficiently long mitigation durations ܶ. For ܶ > 8/ܽ in the static networks, and ܶ > 12/ܽ in the lattice structures, the count of remaining infectious individuals ‫ܫ‬ becomes quite low; importantly, these are also concentrated in the vicinity of a few persisting infection centers. This leads to a steep drop in ܴ ஶ /ܰ , the total fraction of ever infected individuals, by a factor of about 4 in the static network, and 3 in the dynamic lattice architectures. Thus, in these instances, follow-up disease control measures driven by high-fidelity testing and efficient contact tracking should be capable of effectively eradicating the few isolated disease resurgence centers. However, to reach these favorable configurations for the implementation of localized and targeted epidemic control, it is imperative to maintain the original social-distancing restrictions for at least a factor of three (better four) longer than it would have taken the unmitigated outbreak to reach its peak (ܶ ≈ 3/ܽ … 6/ܽ in our simulations) -for COVID-19 that would correspond to about two months. As is evident from our results for two-dimensional small-world networks that perhaps best represent human interactions, it is also absolutely crucial to severely limit all far-ranging links between groups to less than 5 % of the overall connections, during the disease outbreak. The graphs compare the outbreak data obtained without any mitigation (grey) and with social distancing measures implemented for different durations ܶ, as indicated. In all cases, social distancing is turned on once ‫ܫ‬ሺ‫ݐ‬ሻ reaches the set threshold of 10 % of the total population ܰ. The resurgent outbreak is drastically reduced in both its intensity and growth rate as social distancing is maintained for longer time periods ܶ. (The data for each curve were averaged over 100 independent realizations; the shading indicated statistical error estimates.) Inset: time ߬ to reach the second peak following the end of the mitigation; the data indicate an exponential increase of ߬ with ܶ. 
Square lattices with diffusive spreading
On our regular square lattice with ‫ܮ‬ ଶ sites set on a two-dimensional torus, we implement the stochastic Susceptible-Infectious-Recovered (SIR) epidemic model with the following individualbased Monte Carlo algorithm: 1. Randomly distribute ܰ individuals on the lattice, subject to the restriction that each site may only contain at most one individual, and with period boundary conditions. Some small fraction of the individuals will initially be infectious, while the remainder of the population will be susceptible to the infection. 2. Perform random sequential updates ‫ܮ‬ ଶ times in one Monte Carlo step (MCS) by picking a lattice site at random, and then performing the following actions: a. If the selected site contains a susceptible ܵ or a recovered individual ܴ, a hopping direction is picked randomly. If the adjacent lattice site in the hopping direction is empty, then the chosen individual is moved to that neighboring site with hopping probability ݀ that is related to a macroscopic diffusion rate. b. If the chosen lattice site contains an infectious individual ‫,ܫ‬ it will first try to infect each susceptible nearest neighbor ܵ with a prescribed infection probability ‫.ݎ‬ If this attempt is successful, the involved susceptible neighbor ܵ immediately changes its state to infected ‫.ܫ‬ After the originally selected infected individual has repeated its infection attempts with all neighboring susceptibles ܵ, it may reach the immune state ܴ with recovery probability ܽ. Finally, this particular individual, whether still infectious or recovered, tries to hop in a randomly picked direction with probability ݀, provided the chosen adjacent lattice site is empty. 3. Repeat the procedures in item 2 for a preselected total number of Monte Carlo steps. To determine the effective (coarse-grained) basic epidemic reproduction ratio ܴ , we fit the infection curves to straightforward numerical integrations of the deterministic SIR rate equations ݀ܵሺ‫ݐ‬ሻ/݀‫ݐ‬ = ‫ݎ−‬ ܵሺ‫ݐ‬ሻ ‫ܫ‬ሺ‫ݐ‬ሻ/ܰ, ‫ݐ݀/‪ሻ‬ݐ‪ሺ‬ܫ݀‬ = ‫ݎ‬ ܵሺ‫ݐ‬ሻ ‫ܫ‬ሺ‫ݐ‬ሻ/ܰ − ܽ ‫ܫ‬ሺ‫ݐ‬ሻ, ܴ݀ሺ‫ݐ‬ሻ/݀‫ݐ‬ = ܽ ‫ܫ‬ሺ‫ݐ‬ሻ, and adjust the lattice simulation infection probability ‫ݎ‬ ≈ 1.0 and to a lesser extent, the recovery probability ܽ to finally match the targeted COVID-19 value ܴ ≈ 2.4. We note that this slightly `renormalized' value for ܽ is subsequently utilized to set the time axis scale in the figures. On the mean-field level, initially ܴ = ሺ‫ܽ/ݎ‬ሻ ܵሺ0ሻ/ܰ, since all nodes are mutually connected. In spatial settings, ܵሺ0ሻ/ܰ is to be replaced with the mean connectivity (i.e., the coordination number for a regular lattice) to susceptible individuals. The lattice simulation data is fitted with the mean-field result by matching two parameters: the maximum value and the half-peak width of the infectious population curve ‫ܫ‬ሺ‫ݐ‬ሻ, see Figure S1 . The lattice simulation curve digresses from the mean-field curves at low ‫ܫ‬ሺ‫ݐ‬ሻ values, far away from the peak region. In the lattice simulations, the initial rise of the infectious population curve exhibits power-law growths ‫ܫ‬ሺ‫ݐ‬ሻ ~ ‫ݐ‬ ଵ.ସ±.ଵ and ܴሺ‫ݐ‬ሻ ~ ‫ݐ‬ ଶ.ଷ±.ଵ in clear contrast with the simple exponential rise of the mean-field SIR curve as obtained from integrating the mean-field rate equations. We note that these are the standard critical exponents ߠ and 1 + ߠ for the temporal growth of an active seed cluster near a continuous non-equilibrium phase transition to an absorbing extinction state (11) .
Figures S2 about here. Figure S2 shows the dependence of the asymptotic number of recovered individuals ܴ ஶ on the density ߩ for various sets of hopping rates ݀ and initial infectious population values ‫ܫ‬ሺ0ሻ. These data indicate the existence of a well-defined epidemic threshold, i.e., a percolation-like sharp transition from a state when only a tiny fraction of individuals is infected, to the epidemic state wherein the infection spreads over the entire population (11) . As one would expect, this critical point depends only on the ratio ܽ/݀ of the recovery and hopping rates. Varying the lattice simulation parameters just shifts the location of the epidemic threshold. Once the model parameters are set in the epidemic spreading regime, the system's qualitative behavior is thus generic and robust, and only weakly depends on precise parameter settings. For our two-dimensional small-world network, whose construction is schematically depicted in Figure S3 , we employ a similar Monte Carlo algorithm as described above; the essential difference is that individuals may now move to adjacent nearest-neighbor as well as to distant lattice sites along the pre-set `short-cut' links. Figure S4A demonstrates (for fixed diffusivity ݀ = 1) that as function of the fraction ߮ of long-distance links in a two-dimensional small-world network, the epidemic threshold resides quite close to zero: The presence of a mere few `shortcuts' in the lattice already implies a substantial population mixing. The inset, where the ߮ axis is scaled logarithmically, indicates that sizeable outbreaks begin for ߮ ≥ 0.05. Figure S4B similarly shows the outbreak dependence on the diffusion rate ݀ (here for ߮ = 0.6), with the threshold for epidemic spreading observed at ݀ ≈ 0.3. Evidently, prevention of disease outbreaks in this architecture requires that both mobility and the presence of far-ranging connections be stringently curtailed. Figure S4 about here.
For both the randomly connected and scale-free contact networks, we employ the Gillespie or dynamical Monte Carlo algorithm, which allows for efficient numerical simulations of Markovian stochastic processes. It consists of these subsequent steps:
1. Initially, few nodes are assigned to be infected ‫,ܫ‬ while all other nodes are set in the susceptible state ܵ. Each susceptible node ܵ is characterized by a certain number of active links that are connected to infected nodes ‫.ܫ‬ 2. We then determine the rate at which each infected node ‫ܫ‬ will recover, and at which each susceptible node ܵ with a non-zero number of active links becomes infected. From these we infer the total event rate ‫ݎ‬ ௧௧ . 3. Based on this total rate ‫ݎ‬ ௧௧ , we select the waiting time until the next event occurs from an exponential distribution with mean ‫ݎ‬ ௧௧ . 4. We then select any permissible event with a probability proportional to its rate, update the status of each node, and repeat these processes for the desired total number of iterations.