Not all interventions are equal for the height of the second peak.

In this paper we conduct a simulation study of the spread of an epidemic like COVID-19 with temporary immunity on finite spatial and non-spatial network models. In particular, we assume that an epidemic spreads stochastically on a scale-free network and that each infected individual in the network gains a temporary immunity after its infectious period is over. After the temporary immunity period is over, the individual becomes susceptible to the virus again. When the underlying contact network is embedded in Euclidean geometry, we model three different intervention strategies that aim to control the spread of the epidemic: social distancing, restrictions on travel, and restrictions on maximal number of social contacts per node. Our first finding is that on a finite network, a long enough average immunity period leads to extinction of the pandemic after the first peak, analogous to the concept of "herd immunity". For each model, there is a critical average immunity duration Lc above which this happens. Our second finding is that all three interventions manage to flatten the first peak (the travel restrictions most efficiently), as well as decrease the critical immunity duration Lc , but elongate the epidemic. However, when the average immunity duration L is shorter than Lc , the price for the flattened first peak is often a high second peak: for limiting the maximal number of contacts, the second peak can be as high as 1/3 of the first peak, and twice as high as it would be without intervention. Thirdly, interventions introduce oscillations into the system and the time to reach equilibrium is, for almost all scenarios, much longer. We conclude that network-based epidemic models can show a variety of behaviors that are not captured by the continuous compartmental models.

Introduction

When recovering from a disease grants temporary immunity against it, it can happen that an epidemic dies out locally, but survives elsewhere, returning at a later point in time. We observe a “second peak”. A second peak can also happen when interventions are effectively applied to slow down the spread of a disease locally, but are then lifted. This phenomenon has a clear geometric component. Standard compartmental models for epidemic curves are inherently a-geometric, because they assume a perfectly mixed population. Agent-based epidemiological models allow for an embedding of the population in a geometric space to capture this effect more clearly.

In this paper we study this geometric effect empirically. We conduct a simulation study about the effect of temporary immunity on spreading processes on different agent-based models. Immunity being only temporary has been observed for various types of coronaviruses, see [28] and references therein. It is currently being debated whether the current Covid-19 epidemic provides lasting or temporal immunity, see e.g. [42].

We also model the effect of three different interventions:

social distancing,

traveling restrictions, and

limiting the number of social contacts.

Our simulation results may serve as a qualitative indication of possible outcomes of epidemic spread and intervention strategies for highly infectious diseases, such as the current COVID-19 pandemic. It also gives qualitative predictions about the ways an epidemic might (or might not) return, depending on the average duration and level of temporary immunity. The focus is twofold: first, on understanding how the underlying space affects the outcome of the simulations, and second, to see and compare the effect of the three intervention strategies on the pandemic. For the first goal, we compare four scenarios with underlying geometries:

Scenario 1: The underlying network is “purely geometric”, ignoring long-range connections. For this scenario we use the square nearest-neighbor torus as the base graph.

Scenario 2: The underlying network is a “mean-field network”, ignoring the spatial component. For this scenario we use the configuration model, which can mimic the local statistical properties of real human contact networks, such as degree distributions.

Scenario 3: The underlying network is a mixture between a geometric and mean-field network. For this we use Geometric Inhomogeneous Random Graphs, which possess geometric features and can match the statistical properties of real human contact networks.

Scenario 4: The epidemic is modeled in a “mean-field continuous space”, ignoring the spatial component and approximating the discrete population by a continuum. For this we use systems of ordinary differential equations (also called: a compartmental model), which are currently the most popular tool for epidemic curve modeling.

Our aim is to compare these models, and see what the effect is of considering more realistic representations of the underlying space and spreading mechanisms. We, however, emphasize that our paper provides qualitative estimates, not quantitative ones. As a result, the study here might underpin or support certain intervention strategies more than others, but we refrain from (and see no justification for) using these simplified models to make numerical predictions with respect to the current COVID-19 outbreak.

Epidemic spread with temporary immunity on a network

We choose the simplest possible model that shows, qualitatively, the behavior that we would like to observe, namely, the effect of temporary immunity. This model can, however, be easily modified to accommodate other compartments, such as an incubation period (exposed state) (as e.g. in [43]), deaths, asymptomatic cases, etc. For simplicity, we describe the dynamics in discrete time. A continuous time version is analogous and shows the same qualitative behavior. The model is similar to the one studied in [18], where the focus is on the critical regime.

The spreading process changes at discrete time steps, each time step corresponds to, say, a day. We fix the network G in advance. We think of nodes in the network as individuals. Each node within the network can be in three possible states: susceptible (S), infected (I) or temporarily immune (T). The neighbors of a node u are nodes with a direct connection (also called link, or edge) to u. A connection may correspond to a friendship or an acquaintance, or simply a contact event. The discrete time dynamics between the three states are described as follows (see also Fig. 2):

Fig. 2

Schematic diagram of node states and their transition probabilities. A susceptible node u may only become infected when it has at least one infected neighbor. Each of the infected neighbors infects u with probability β, independently, in each step.

Infecting:Each infected node, while being infectious, infects each of its neighbors within the network with probability β at every time-step. Infections to different neighbors happen independently.

Healing: When infectious, each node heals with probability γ at every time-step, independently of other nodes, and independently of infecting other nodes. The average infectious period of an infected node is 1/γ time-steps. Upon healing, the node becomes temporarily immune.

Losing immunity: Each temporarily immune node loses its immunity with probability η at every time-step, independently of other nodes. The average immune period of a node is thus 1/η time-steps. After losing immunity, the node becomes susceptible again.

We run these dynamics on three different types of networks, covering the above mentioned scenarios: the purely geometric square grid, the pure network configuration model, and the mixture Geometric Inhomogeneous Random Graph (GIRG) (see Definition 2.1, Definition 2.2, and 2.3) below.

We adapt these dynamics to obtain a continuous S-I-T-S compartmental model (Scenario 4) in Section 2.4. We choose the parameters of the underlying network models in Scenarios 1–3 so that the models have comparable average degree (and when possible, also comparable variance), and then vary the parameters β, γ and η. In the remainder of this section we give a quick summary of our findings, then in Sections 2 and 3 we define the four scenarios and three intervention methods in more detail, and in Section 4 we present our detailed simulation results.

Phases of S-I-T-S on continuous compartmental models.

For Scenario 4, we observe three (commonly known) phases of the epidemic, that we briefly describe here for comparison. For the model description, see Section 2.4.

Subcritical and critical phase: Immediate extinction. Whenever the basic reproduction number, the epidemic dies out without producing a peak. For β/γ < 1, it dies out quickly (logarithmic in the initial number of infected), while for it dies out slowly (polynomial in the initial number of infected). See Fig. 14(b).

Fig. 14

Solutions of the ODE when with 100 initially infected nodes. The infection almost seemingly disappears from the system before the second peak appears when . The infection does not turn into a pandemic when β < γ.

Supercritical phase: Peaks of decreasing magnitude towards a limiting stationary proportion. Whenever β/γ > 1, the proportion of infected population stabilizes at . There are several larger peaks before the equilibrium is reached, see Figs. 14(a)–14(b). The incline and decline of peaks are exponential. See more in Section 2.4.

New supercritical phases in S-I-T-S epidemics on network models.

The presence of the network structure, as well as the underlying space, combined with temporary immunity has a surprising effect on the qualitative behavior of the epidemic. More specifically, new phases arise in the phase diagram compared to those in Section 1.2. The duration of the immune period plays a more profound role than in the continuous compartmental model of Section 2.4, where the single supercritical phase is entirely determined by the ratio β/γ of the infection rate and the healing rate. For finite spatial network models, there is no well-defined definition of R 0: the average node-degree as well as the topology of the network influences the probability of entering each phase. The final size of the epidemic is, at least for the most commonly known S-I-R model, often computable even for network models: See the works [38], [48], [52], [58], [65] for the final size of the epidemic in mean-field type networks, and the phase diagram for S-I-T-S on the lattice [18, Figure 1].

For the S-I-T-S epidemic spread on all three network models (Scenarios 1, 2, and 3) we observe the following three phases, with 2S below being new compared to Scenario 4: Stochastic transition between phases with a sharp threshold: Due to the stochastic nature of the model, the phases 1, 2S, 2M are, (at least on a finite but large network) intertwined. This means that for a given fixed set of parameter β, γ, η, a single run of the S-I-T-S epidemic may enter any of the three phases, with different probabilities. However, the probability of entering a given phase undergoes a sharp transition in the parameter values: for fixed β, γ, a small change in η results in a shift from almost always seeing phase 2S to almost always seeing phase 2M. This is called a sharp threshold phenomenon, see [26], [27]. The sharp threshold is defined as the value of η(N) where the survival probability crosses 1/2. The same theory suggests that the threshold η(N) converges in the large network limit. The probability of long survival (phase 2M) can be seen in Fig. 3 for eight different models as a function of η, all with the same average degree . We conjecture that the long survival probabilities are monotone in η, see Fig. 20, where the experiment is carried out with more runs for each value of η.

Fig. 3

A comparison of the probability of phase 2M on eight different networks with the same average degree 8 and number of nodes as a function of the η. Probabilities were computed using a 100 runs for all parameter values (η) and models, with . Legend: CM, GIRG: Scenarios 2 and 3, with different parameters and Grid: Scenario 1. For all networks we see a sharp transition at a critical value of ηwhere the system moves from phase 2P (single peak followed by extinction) to 2M (survival with multiple peaks) with overwhelming probability.

Fig. 20

Monotonicity of survival probabilities as a function of η. For the graph a GIRG(3.3, 1.9) (weak no-travel rule), we have tested the monotonicity of the survival probability curve, by increasing the number of runs to 500 for each value of η. The curve appears to be monotonous, with a sharp threshold near η ≈ 0.007.

Subcritical and critical phase: Immediate extinction. The epidemic goes extinct almost immediately (for all runs), and the number of infected decays exponentially in this case. When R 0 ≈ 1, the epidemic may die out more slowly. We emphasize though, that the critical and subcritical phase do not separate so clearly on the stochastic S-I-T-S model on networks: due to the stochastic nature, the epidemic might die out quickly even when R 0 > 1. For the extinction time (logarithmic vs polynomial of the network size), the initial number of infected individuals plays a more important role. This is, however, out of the scope of our current study and we do not pursue this direction further.

Supercritical phases: Possible large outbreak. With positive probability, there is a large outbreak. For network-based S-I-T-S, the supercritical phase separates into two sub-phases:

Extinction after a Single peak: When the duration of the immune period is long, the epidemic has a single peak, after which it immediately goes extinct. Heuristically, the reason for the existence of this new phase is that after the first epidemic peak, earlier infected nodes become immune and maintain their immunity long-enough to provide barriers in the network that the infection cannot pass through. See Fig. 12(a).

Fig. 12

No intervention: The effect of geometry. Each model has average degree 8, and power-law exponent . We see that there is only a match between the mean field network configuration model (CM) when α < 2, hence, there are many long-range connections in the network. The epidemic curve of configuration model with (CM(3.3), orange) matches that of GIRG with (GIRG(3.3,1.3), red). When α is increased to 2.3 (GIRG(3.3,2.3), green). This effects the first peak of epidemic curve to flatten, its magnitude is shrunk by almost 40%, even though the average degree is tuned to remain the same. Top:(500 days immunity) On all three models, a single peak can be observed, and all runs die out within 70 days. Bottom: the first peak has the same height as on to the top picture and is not shown. Without intervention (GIRG(3.3,1.3)), the second and further peaks of the epidemic behave similarly to the mean-field network (CM(3.3)). With intervention, (GIRG(3.3,2.3)), both the oscillation period as well as its amplitude are higher, but the average stationary proportion of infected nodes is lower. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Long-time survival, Many peaks: When the duration of the immunity period is relatively short, the epidemic follows a (qualitatively) similar curve to the one observed in the ODE in Section 2.4. There is a first major outbreak, followed by smaller second, third, etc., peaks, decreasing in magnitude, and eventually settling on a (meta-stable) stationary proportion of infected and immune population. See Fig. 12(b). We mention that this equilibrium is a meta-stable state, since the all-susceptible configuration is an absorbing state. The time to reach that, however, is exponentially long in the network size, see [15].

Height of the first peak: The height of the first peak, in all models, is insensitive to η, just like for Scenario 4, see Fig. 4 A possible explanation for this insensitivity is that during the initial spread, the S-I-T-S process is closely approximated by an S-I-R process, where temporarily immune nodes are simply removed.

Fig. 4

The effect of interventions on the first and second peak, as a function of η (the inverse of the average immunity period). Abbreviations in the legend: ‘hub’: limiting maximal node degree, ‘long’: strong restriction on travel, ‘perc’: social distancing. The first entry is increasing α, a weak restriction on travel. GIRG(3.3,1.3) is the original network. We see that the height of the first peak is insensitive to the duration of immunity, and that cutting long edges is most effective in reducing the first peak. For the second peak, limiting node degree pushes the critical η for appearance of the second peak from 0.002 to 0.006, however, above that the second peak is twice as high as for other interventions. See also Fig. 21, Fig. 22.

Shape and location of the first peak: For Scenario 2 (mean-field network), the shape of the first and later peaks are similar to that of Scenario 4, i.e., exponential incline followed by exponential decline. Scenario 1, the grid, however, has a linear incline followed by linear decline, if there is a first peak, which is not true, when η is large. See Section 4.2. The shape of the curve for Scenario 3 depends on the parameter values: as long as there are many long-range connections, and hubs in the network (nodes with very high degree), the curve resembles that of Scenario 2.

When this is not the case, Scenario 3 starts to gradually resemble the epidemic curve of the grid (Scenario 1): linear incline followed by linear decline, see Fig. 12.

Explanation of the linear peak-shapes: The linear growth was also investigated in a similar setting recently in [64]. There, the authors proposed to model the contact network using a related (one-dimensional) model that is called the Newman-Watts small world [53]. The mathematical explanation for seeing a linear growth curve is as follows: these models are embedded in geometry, and the parameters are such that long edges do not dominate yet, and hence the typical ball-growth is polynomial. This is in accordance with theoretical results (5) related to the average distance in the network being short (poly-logarithmic) when there are many long-connections but long (polynomial) when these are scarce. See [43] and references therein for the role of long-connections.

Topology-dependent fluctuations after the first peak: The fluctuations after the first peak are highly dependent on the network topology. For power-law networks with infinite asymptotic variance, there is no second peak and the system reaches equilibrium very quickly after the first peak. For finite variance networks, several further peaks occur, see e.g. Figs. 23, 25. The transition is most profound when interventions change the topology, see Section 1.4 below.

Fig. 23

No intervention: the epidemic curves for the networks G1, G2 on which we apply the intervention methods. G1 is a GIRG with and average degree 9.6, while G2 is a GIRG with and average degree 8.7. Top: The first peak of the epidemic. Bottom: The second peak and stabilization of the number of infected. The average height and location of the first peak is on day 8 for G1, while it is on day 14 for GIRG with . On G2, we see a second peak at day 79 of average height and the equilibrium number of infected . The equilibrium number of infected is reached without a second peak on G1.

Fig. 25

Later effects of interventions on G1: Top: Social distancing (“GIRG(2.5, 1.3) perc”), as well as the weak no-travel rule (“GIRG(2.5, 2.37)”). No second peak can be observed, just like on the original network. Bottom: Limiting the maximal number of contact (“GIRG(2.5, 1.3) hub”), and the hard no travel rule (“GIRG(2.5, 1.3) long”’): second and further peaks are periodically present, with decreasing magnitude and period roughly the average immunity duration. G1 is a GIRG with and average degree 9.6 on nodes. All networks after intervention have average degree  ≈ 4.9.

We mention that for Scenario 3 the way the second peak happens is entirely different from the first peak: the initial spread is local, emanating from the source, (with some long connections causing non-local new infection centers), while the second peak is spread out, infections occur everywhere in space. This is quite natural, since during the initial wave, the epidemic fills the space, and thus the second peak happens roughly when the immunity arising from the first wave starts to wear off, and at that moment it does not have a well-defined source anymore. This is why we never see the epidemic to die out after a second peak. It either dies out after the first peak, or many peaks occur and the system reaches equilibrium.

We conclude that the continuous ODE approximations often fail to capture the later behavior of the epidemic, which depends highly on the network topology.

Comparison of intervention methods a) b) c) on scenario 3.

We carry out the intervention methods on Scenario 3, the geometric scale-free network model. We model intervention (A), keeping physical distance by randomly removing connections from the network.

Intervention (B), allowing for less travel, is modeled in two different ways. Once, it is modeled by drawing for each present edge a random exponential threshold with a common mean L, and cutting the edge if it exceeds its threshold value, resulting in a connection-length distribution that decays exponentially above the threshold L. Secondly, it is also modeled by increasing the long-range parameter α to α new.

Intervention (C), limiting the maximal number of contact per person, is modeled by prescribing a maximal node degree M and then for each node u with degree higher than M, randomly chosen connections of u are cut until at most M connections remain, see more details in Section 3.

We choose the parameters M, L, α new so that the average degree is the same after all interventions. See Fig. 1 for a visualization of the intervention methods.

Fig. 1

Visualization of the change in the underlying contact network under interventions for Scenario 3 (see Fig. 26 below). By (w) and (s) we denote the weak and strong travel restrictions, respectively.

Fig. 26

Later effects of interventions on G2: Top: Social distancing (“GIRG(3.3, 1.3) perc”), as well as the weak no-travel rule (“GIRG(3.3, 1.9)”) compared to the original network (“GIRG(3.3, 1.3)”). Second and further peaks are present, equilibrium is reached within 600 days. Bottom: Limiting the maximal number of contact (“GIRG(3.3, 1.3) hub”), and the hard no travel rule (“GIRG(3.3, 1.3) long”’): second and further peaks are more profound, equilibrium is only reached around 1200 days. G2 is a GIRG with and average degree 8.7 on nodes. All networks after intervention have average degree  ≈ 4.7. The first peak can be seen on Fig. 24. The first peak of the red curve (corresponding to truncating long edges) is the first peak of the epidemic. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

We summarize our most important, qualitative findings. For quantitative values see Section 3.

The height of the first peak drops. In all interventions, the height of the first peak is dropping, by at least as much as the shrinkage in average node degree (intervention A) and even more with other interventions. The most effective intervention in this respect is the hard no-travel rule. See left on Fig. 4.

Elongated first peak. The time and duration of the first peak is later/longer than without intervention. For keeping physical distance, the time-shift is the least apparent (a few days), while for intervention B it is the most profound, a factor of 10, see bottom of Fig. 24, as well as Table B.1.

Fig. 24

The first peak under interventions. Top: The effect of interventions on G1: Social distancing (“GIRG(2.5, 1.3) perc”), the hard no travel rule (“GIRG(2.5, 1.3) long”), limiting the maximal number of contact (“GIRG(2.5, 1.3) hub”), and the weak no-travel rule (“GIRG(2.5, 2.37)”). The first peak disappears before day 60 in all cases. The curves for the weak no-travel rule and social distancing are almost identical. Bottom: The effect of interventions on G1: Social distancing (“GIRG(3.3, 1.3) perc”), the hard no travel rule (“GIRG(3.3, 1.3) long”), limiting the maximal number of contact (“GIRG(3.3, 1.3) hub”), and the weak no-travel rule (“GIRG(3.3, 1.9)”). The hard no travel rule is most effective in flattening the curve: both in its height as well as the day of the peak.

Table B1

Time and location of the first and second peak under interventions, for all values of η. O: original graph, A: social distancing, Bw: weak no-travel rule (increasing α, B: hard no-travel rule (cutting long edges), C: limiting maximal node degree. Where the data is missing, the system produces a single peak only. The large numbers for the time of the second peak on G1 are coming from random oscillations of the equilibriated system when there is no clear second peak.

		First peak ( × 103)					Time of first peak					Second peak ( × 103)					Time of second peak
	η	O	A	Bw	B	C	O	A	Bw	B	C	O	A	Bw	B	C	O	A	Bw	B	C
G1	0.001	85	51	52	38	49	8	11	10	25	25
	0.002	84	53	52	39	49	8	10	10	25	25	1.2	1.5	1.7			1103	630	636
	0.003	84	50	51	39	48	8	11	10	24	26	1.7	1.3	1.2			1023	1027	1126
	0.004	84	53	51	39	49	8	11	11	25	26	2.2	1.6	1.6	4.0		609	920	634	193
	0.005	83	51	52	39	49	9	11	10	25	25	2.8	2.0	2.0	5.0		1177	790	697	170
	0.006	85	53	51	39	49	8	11	10	25	26	3.3	2.4	2.3	5.1	9.9	772	680	726	151	184
	0.007	83	53	50	39	49	8	12	10	25	25	3.9	2.7	2.6	5.4	10.1	436	705	708	136	163
	0.008	84	55	50	39	49	8	10	10	25	25	4.4	3.1	3.1	5.7	10.5	1184	1196	460	126	150
	0.009	83	49	49	40	50	8	12	10	24	25	4.9	3.5	3.4	5.5	10.4	308	895	868	117	140
	0.010	86	52	51	39	49	8	11	10	25	25	5.4	3.9	3.8	5.8	10.6	439	321	286	110	130
	0.011	83	51	52	40	48	8	11	10	24	25	6.0	4.2	4.1	5.9	10.5	972	352	244	105	123
	0.012	85	53	52	39	48	8	11	10	24	25	6.5	4.6	4.5	6.0	10.6	1071	225	552	102	116
	0.013	86	53	52	39	50	8	11	10	25	25	7.0	5.0	5.0	6.3	10.9	727	237	246	98	111
G2	0.001	78	46	32	12	37	14	21	30	81	40
	0.002	79	45	33	13	37	14	21	31	80	40
	0.003	80	45	33	13	37	14	21	27	79	42	3.3					154
	0.004	78	43	33	12	38	14	21	26	79	40	4.2	3.1				130	182
	0.005	79	47	33	12	38	14	21	31	78	40	5.1	5.0	4.0			110	171	228
	0.006	80	44	33	13	37	14	22	29	79	41	5.3	4.8	5.4	6.2		101	152	209	332
	0.007	79	46	33	13	37	13	21	28	81	40	5.7	5.2	5.4	6.1	13.5	97	143	183	344	242
	0.008	80	46	34	13	38	14	22	29	82	41	6.1	5.5	6.0	6.7	12.5	89	129	175	314	233
	0.009	78	45	34	13	38	14	21	29	80	39	6.5	5.9	6.4	6.9	13.0	86	121	161	284	208
	0.010	78	47	33	13	38	14	21	29	79	40	7.0	6.0	6.2	7.0	13.8	81	112	147	265	198
	0.011	80	47	34	13	38	13	23	29	80	41	7.4	6.5	6.6	7.0	14.2	78	112	142	236	188
	0.012	77	44	34	13	38	14	21	28	81	40	7.9	6.4	6.5	7.1	14.7	77	107	136	219	174
	0.013	80	45	35	14	39	14	22	26	81	40	8.4	6.8	7.0	7.2	14.7	77	102	131	205	162

Long-survival probabilities drop. Each intervention pushes the probability of the system entering phase (2M) lower for a fixed parameter setting β, γ, η. When increasing η (thus decreasing the immunity duration), at a critical η the system abruptly goes from dominantly phase 2S from dominantly phase 2M, (a sharp threshold). The sharp threshold η of long-survival is increased under all interventions, see Fig. 19. In this respect, interventions B and C perform best.

Fig. 19

The effect of interventions: the long-survival probabilities (phase 2M) are decreasing significantly under each intervention strategy. The four interventions are: (A) social distancing, denoted as ”perc”, (B) limiting maximal degree, denoted as ”hub”, (Ch) a hard no-travel rule, denoted as ”long” on the legend, and (Cw) a weak no-travel rule: increasing α. We see that in all interventions, the long-survival probabilities are below the original curve. The threshold η is lowest for percolation and highest for limiting the maximal degree, while truncating long-edges and increasing α in the network both show an interesting non-monotonous survival probability curve. Top: The initial network is G1, a GIRG with and . Bottom:The initial network is G2, a GIRG with and .

Higher second peak. Intervention B and C result in a higher second peak, or even make a second peak appear where originally it was not present. The worst intervention in this respect is intervention C (limiting the maximum degree), where the second peak can be as high as 1/3 of the first peak. See Fig. 4 and 21. For infinite variance degree GIRG, on Fig. 21) the strong travel restriction also causes a second peak, while social distancing and increasing α does not.

Fig. 21

The effect of interventions on the first and second peak on G1, as a function of η (the inverse of the average immunity period). Abbreviations in the legend: ‘hub’: limiting maximal node degree, ‘long’: cutting long edges, ‘perc’: social distancing. The first entry is increasing α. GIRG(2.3,1.3) is the original network. We see that the height of the first peak is insensitive to the immunity duration, and that cutting long edges is most effective in reducing the first peak. For the second peak, limiting node degree pushes the critical η for appearance of the second peak from basically 0 to 0.005, however, above that the second peak is twice as high as for other interventions. Cutting long-edges also has a similar effect, although to a lesser extent. See also Figure 4. Star indicates that there is no second peak, and the height we see is the height of the equilibrium, see the top figure on Fig. 25.

A possible explanation for a higher second peak is the following. Under intervention C (limiting node degree) and Intervention B (cutting long-edges), there are no more hubs in the network: nodes with very high degree. Once infected, hubs infect a large proportion of their neighbors roughly at the same time. Moreover, hubs are much closer to each other (in terms of graph distance) than the average distance in the graph, thus the infection can travel quickly between hubs. To summarize, hubs synchronize the system. Once the hubs are removed, this synchronization is also removed from the system and oscillations similar to the grid appear. Limiting node degree without travel restrictions has the strongest effect on the removal of hubs. Limiting travel reduces the degree of hubs, but do not remove them completely, since nodes that were hubs before the intervention, have typically many local connections as well, see Fig. 16.

Fig. 16

Visualization of the change in the underlying contact network under interventions. The network on the top left is similar to G2: a GIRG with average degree 4.8 on nodes. All interventions result in an average degree  ≈ 2.6, just like on Fig. 15, even though the pictures here look much more sparse.

More oscillations to reach equilibrium. Interventions C and B, while they are most effective regarding the first peak, introduce high oscillations in the system, and the time to reach equilibrium can be much longer than in the original model (from a few days to more than a thousand days). Again, these oscillations are explained by the lack of synchronization. See Figs. 25 and 26 for these later effects.

To summarize, we conclude that the effect of interventions vary to a high extent and depend on the precise way the intervention changes the network topology. Again, ODE-approximations fail to capture these various effects.

Network models

Scenario 1: Lattice network models

Arguably the simplest geometric network model is the nearest neighbor lattice. To avoid boundary effects, and make the network completely homogeneous, we shall study epidemics on the torus instead of a square that we define here:

Nearest neighbor tori

Arrange nodes on a square grid, and label them with two coordinates for . To obtain connect to when either and or and . To obtain connect to when and . The networks have nodes, see Fig. 6.

Fig. 6

Left: The two dimensional torus on nodes. Each node has four neighbors. Right: The modified two dimensional torus on nodes. Each node has 8 neighbors. The nodes in the bottom row are also connected to the nodes in the top row, and the nodes in the very left column are also connected to the nodes in the very right column.

Lattice models are homogeneous and the average distance within the network is large, it grows polynomially with the number of nodes. Writing for the number of connections on the shortest path between nodes u and (also called: graph distance), the average distance becomes

Scenario 2: The configuration model and other a-geometric network models

A-geometric random network models are often used as null-models to compare network data to purely random networks. Here we describe a commonly used model, the configuration model, and mention similar alternatives. The configuration model is a well-studied object in mathematical network science, dating back to Bollobás [10] and Molloy and Reed [46], [47]. The main advantage of this simple model is that it can mimic the degree distribution of real-life networks.

Configuration model

Fix N ≥ 1 the number of nodes. Prescribe to each node its node-degree so that the total degree is even. To form a graph, to each node u we assign number of half-edges and the half-edges are then paired uniformly at random to form edges. The resulting random multi-graph is the configuration model. The erased configuration model can be obtained by erasing all self-loops and multiple edges between nodes.

See Fig. 7 for an illustration of the algorithm, and Fig. 8 for two realizations. Observe that we pair half-edges randomly, hence, even though the prescribed degrees might be deterministic, the configuration model is a random (multi)graph. It is a multi-graph since multiple connections between nodes and even connections to the same node may occur, however, these become proportionally insignificant when the size of the graph is sufficiently large. An easy way to obtain a simple graph from the configuration model is to erase second, third, etc edges between nodes, as well as self-connections. It is shown that the erasure does not affect the empirical distribution of the model in any significant way [34]. For a mathematical overview on the properties of the configuration model and erased configuration model, see [34].

Fig. 7

The algorithm producing the configuration model. Left: First, for each node in the network we prescribe its degree and draw it as half-edges. Middle: Then, we match half-edges randomly to form edges (connections). This may be done in a sequential way, by always choosing a uniform pair from the remaining half-edges. Here, an intermediate stage is shown when there are five edges formed. Right: All half-edges are matched. The output is a random graph called the configuration model.

Fig. 8

Two examples of the configuration model. The nodes have average degree 4, with power-law degree distribution, with exponent (left) and (right).

In real-life networks, an extreme variability of node degree is often observed, see [1], [25], [51]. Extreme node degree variability results in the presence of a few individuals with extreme influence on spreading processes, the hubs or superspreaders [51]. Mathematically, this extreme degree variability can be expressed using the empirical distribution of node degrees, that follows a power-law:with exponent τ > 2. In our simulations we study the case when the empirical degree distribution follows such a power law, i.e., the node-degrees satisfy (2), for some parameter τ > 2. In this case, the node-degrees in the erased configuration model version also follow (2), with the same exponent τ.

Other popular a-geometric models include the Chung-Lu and the similar Norros-Reittu model [16], [17], [54], and preferential attachment models (also known as the Barabási-Albert model) [3]. In the Chung-Lu and the Norros-Reittu model, only the expected degrees of nodes are prescribed, rather than their exact degree. It is shown that these models, with similar parameter settings, behave qualitatively similarly to the configuration model, see [34] for references. The preferential attachment model is a growing network model, that is best suitable to study the evolution of networks on a longer time-scale. For a relatively short time-frame typical to epidemic spread, a time snapshot of the preferential model could be used. Mathematical results suggest that snapshots of this model behave, again, qualitative similarly to the configuration model with similar degree structure, although numerical values might differ [5], [21], [23], [25]. As a result, we choose the configuration model for baseline comparison.

An important feature of the configuration model is the small-world property. Heuristically, this means that two arbitrary nodes in the network can be connected via very short paths, using only a few connections. Mathematically, the small world phenomenon means that the average graph distance in the network is small, at most logarithmic in the network size. For the configuration model, when τ > 3, the network is proven to be a small world [35], [66], while, when τ ∈ (2, 3), the average graph distance is ultra small, meaning a double-logarithm of the network size:In the former case, when τ > 3, the shape of the epidemic curve is shown to converge in the large network limit, see [4], [6], [39].

While the configuration model easily accommodates power-law node degrees, and has the small-world property, it neither contains communities nor clustering [25], [63]. Empirically, clustering quantifies the effect commonly known as “a friend of a friend is also likely to be my friend”. Mathematically, clustering means the presence of triangles in the network. Communities are parts of the network that have significantly more connections than the same number of randomly selected nodes. One could argue that clustering and communities in human contact networks often arise due to spatial effects: people living nearby tend to know each other with higher probability. To accommodate clustering and communities, as well as keep node degree variability high, we introduce our last model of investigation.

Scenario 3: Geometric inhomogeneous random graphs

The last network class we run the S-I-T-S epidemic on is a mixture of pure geometric and purely random network models. For this we use a general Geometric Inhomogeneous Random Graph (GIRG) model, that is a state-of-the-art model for real-world social and technological networks, embedded in geometric space. The presence of underlying geometry underpins the model: individuals are embedded in space, just like in real life, allowing for local community structures to be present in the contact network, leading to strong clustering [13], [63]. The GIRG incorporates edges bridging spatial distance on all scales, ranging from short to long-range edges, as well as a high variability of node degree: in fact, the model is scale-free in two respects: both in spatial distance that edges cover, and in node-degree variability [13].

Contact- or activity networks of humans have been found to show similar behavior, including heavy-tailed degree distributions, strong clustering, and community structures [1], [8], [49], [50], [51], as well as heavy-tailed distance-distribution for edges, see references in [36].

Recently, several spatial random graph models were developed to incorporate these features: hyperbolic random graphs [9], [33], [56], scale-free percolation [19], and GIRGs [11], [12], [13]. These models can be unified into a general model containing all these three models as special cases. The qualitative behavior of the three models is the same. The following definition is general, and in its last sentence we specify it for GIRG. The underlying space can be the earth’s surface, or the two-dimensional Euclidean space. We denote by x∧y the minimum of two numbers x, y.

Geometric Inhomogeneous Random Graph (GIRG)

Fix N ≥ 1 the number of nodes. Assign to each node a fitness and a spatial location Φ(u). Fix α > 0. For any pair of nodes with fixed connect them by an edge with probabilitywhere dist is a distance function1 , i.e., it measures the distance between the locations . When Φ(u) is chosen uniformly in a box of volume N, we obtain the GIRG.

Any or all of the vertex set, fitnesses, and spatial locations can be random (see Fig. 9). GIRGs have a natural interpretation: the fitnesses express the ability of nodes to have many connections, Φ embeds them in space, and α is the long-range parameter: the smaller α is, the more the model favors longer connections. The parameter space of GIRG is rich enough to model many desired features observed in real networks:

Fig. 9

Two examples of Geometric Inhomogeneous Random Graphs (GIRGs). The nodes are placed randomly into a square of area N. Each node draws a random fitness from a power law distribution with exponent (left) and (right). We used the same location for nodes and the same underlying uniform variables to simulate fitnesses in both cases: for a uniform variable we set the fitness of node to on the left, while on the right. Each pair of nodes with positions x1, x2 and weights respectively, is connected with probability where . Connections are again generated in a coupled way, using the same set of uniform variables for the two pictures, thresholded at p(2.95)and p(3.3), respectively.

extreme variability of the number of neighbors (degrees),

connections present on all length-scales,

small and ultra-small distances,

strong clustering,

local communities.

In real-life networks, (a) corresponds to the presence of a few individuals with extreme influence on spreading processes, the hubs or superspreaders. Mathematically, (a) corresponds to the empirical distribution of node degrees following a power-law, as in (2) above, for some τ > 2. Setting a power-law fitness distribution for in GIRG yields that the degrees satisfy (2), while setting results in a (fairly) homogeneous network with some longer edges, that is known in the literature as (continuum) long-range percolation [7]. The abundance of long-range connections is tuned by α in (4): the smaller α, the more likely are long-range connections (b). The power-law exponent τ and the long-range parameter α tune the average graph distance see [7], [11], [19], [20]:Comparing this to the average distance in the configuration model in (3) and to the lattice models in (1), one sees that geometry and long-range connections play a role when τ > 3, and the model interpolates between the small-world configuration model and the lattice.

Empirically, clustering quantifies the effect commonly known as “a friend of a friend is also likely to be my friend”. Mathematically, clustering means the presence of triangles in the network. Communities and clustering are naturally present in SCLM, because the model favors connections between nodes that are close to each other in space (see Fig. 9). Thus, GIRG incorporates all five desired features.

For the spread of information or infections in such networks, it is known that large-degree nodes and many triangles have opposing effects. On the one hand, nodes of large degree (also called hubs, super-spreaders, or influencers) contribute to fast dissemination, and foster explosive propagation of information or infections [24], [32], [57], [59]. On the other hand, clustering and community structures provide natural barriers that slow down the process, while long-range edges accelerate the spread [2], [37], [40], [41], [45].

Scenario 4: S-I-T-S compartmental model

For baseline comparison to continuous epidemic models, we use a continuous S-I-T-S model, where there are only three possible states, susceptible (S), infected (I) or temporarily immune (T), in a total population of N individuals. Susceptible individuals may become infected via a contact to infectious individuals, while infectious individuals heal and thus become temporarily immune. The assumption is that individuals make contact with a random other individual at rate β. Immune individuals loose their immunity at a certain rate and become susceptible again.

The deterministic ordinary differential equation (ODE) model is governed by the three main parameters:

β: rate of infecting a susceptible individual when being infected,

γ: rate of healing, and becoming temporarily immune,

η: rate of losing temporal immunity and becoming susceptible again.

See Fig. 5 for a schematic diagram. We assume the population size is N, denote by S(t), I(t), T(t) the number of susceptible, infected and temporarily immune individuals at time t, then the corresponding ordinary differential equation becomes:For the infection to survive, β > γ is necessary, otherwise . When β > γholds, this ODE has a stationary state, which is explicitly computable:Observe that the equilibrium proportion of infected is roughly linear in η, as long as η ≪ γ. One can also compute the basic reproduction number, and observe that . Qualitatively, the ODE exhibits either (sub)critical behavior, implying extinction, or supercritical behavior, as described in Section 1.2and depicted in Figs. 14(a)–14(b). On these figures, for the supercritical case one can observe that

Fig. 5

Schematic diagram of states and their transition rates in the continuous compartmental model. Infected individuals make contact at rate β to randomly chosen individuals, as a result, the number of susceptibles increases at rate β · S · I/N. This is a mean-field approximation of the graph model in Section 1.1, see Fig. 2.

the height of the first peak does not depend on the duration of the immune period; only on β, γ.

for fixed β, γ, increasing the duration of the immune period results in deepening the valley between the first and the second peak. The minimum between the first and second peak is always strictly positive, however, it approaches 0 as 1/η tends to infinity. For somewhat realistic values (the minimum gets near 0 around 1/η ≈ 210). Since this is a continuous system, the second peak is always present, explaining the absence of phase 2S in this compartmental model.

Fig. 22

The effect of interventions on the first and second peak on G2, as a function of η (the inverse of the average immunity period). Abbreviations in the legend: ‘hub’: limiting maximal node degree, ‘long’: cutting long edges, ‘perc’: social distancing. The first entry is increasing α. GIRG(3.3,1.3) is the original network. We see that the height of the first peak is insensitive to the immunity duration, and that cutting long edges is most effective in reducing the first peak. For the second peak, limiting node degree pushes the critical η for appearance of the second peak from 0.002 to 0.006, however, above that the second peak is twice as high as for other interventions. See also Figure 21.

Modeling intervention methods on GIRGs

A) Social distancing measures: percolation

In the current Covid-19 outbreak, many governments issued a collection of social distancing measures: keeping 1.5 meter distance from each other, wearing face masks, hand-gloves, washing hands frequently, and so on. Such measures result in reducing the number of contact moments between individuals and thus resulting the chance for the virus to spread. When the social distancing measures are not combined with travel restrictions, this results in the underlying contact network to lose its connections roughly randomly. This has been observed also in [22].

Hence, we model the effect of social distancing by randomly removing a certain proportion of edges in the network. This method is called percolation in the mathematical literature [31]. The same method is used to model social distancing in e.g. [55], where percolation is applied to the configuration model. We denote the resulting network by G perc.

B) Travel restrictions

In the current Covid-19 outbreak, many governments issued everybody to stay at home and only travel when necessary. Outside country travel is not allowed for some countries, e.g. Spain, Italy, Germany, Hungary. This results in the underlying contact network to lose many of its long connections.

We model the removal of long-connections in the network in two different ways: a (seemingly) milder and a stronger intervention measure.

B(w): Weak travel restriction. One way to reduce the number of long-range connection in the network is by increasing the parameter α in GIRG. This intervention does not cause a ‘hard cutoff’ for the length of connections: some long-range connections might remain, but the intervention is strong enough to push the distances in the network from poly-logarithmic to polynomial, see (5).

Increasing α, when τ > 3, decreases the probabilities of long edges: indeed, in (4), for two nodes with locations the ratiois raised to the power α, so as long increasing α reduces the chance of a long connection being present.

When τ > 3, the node fitnesses are, typically, not high enough2 so that, typically, hence the intervention of increasing α is very effective. The average distance increases from poly-logarithmic to polynomial in the network size. On the other hand, when τ ∈ (2, 3), node pairs with ratio happen frequently enough (even on global scale, nodes with fitness are present), hence, increasing α turns out to be ineffective in reducing the average distance in the network, since many long-range edges remain. This is in good correspondence with (5), where the average distance is doubly-logarithmic regardless of the value of α when τ ∈ (2, 3). We denote the resulting network by G.

B(s): Strong travel restriction. Another way to reduce the number of long-range connection in the network is by removing long-edges in the network. Mathematically, we set a threshold distance T, and for each edge present in the original network we draw an independent Exponential variable with mean T; and remove the edge e whenThis intervention does cause a ‘hard cutoff’ for the length of connections at the threshold T. The intervention is also strong enough to push the distances in the network from poly-logarithmic to polynomial. We denote the resulting network by G long.

C) Limiting the maximal number of contacts

Another intervention measure that one can imagine is to limit the maximal number of contact a person can have. On the contact level, this rule barely affects individuals with low number of contacts, but it aims to decrease the number of contact of ‘superspreaders’ in the network: nodes that have a relatively high degree.

We model this effect by choosing a node-degree threshold M, and for each node u that has, in the original network, more than M connections, we remove many randomly chosen connections from node u. The removal of these connections, in turn, also reduces the number of connections of other vertices. Thus, the resulting network has maximal node degree at most M. Since we removed connections randomly, some long connections still remain in the network if they were originally present. We denote the resulting network by G hub.

Simulation results

In this section we present our results for the network based S-I-T-S epidemic model. We first comment on the general description of the types of phases that may occur that are mentioned in Section 1.2, then we compare our findings for each Scenario (the purely geometric lattice models, the mean field network configuration model, and their mixture geometric inhomogeneous random graphs) with respect to the various intervention strategies.

Throughout this section, these underlying networks for the different scenarios (and interventions) are sampled once per parameter setting and kept constant over the simulation runs. On every graph, for every described parameter setting of the S-I-T-S model, the results are based on 100 simulations of the epidemic spreading starting from a single node. This node is sampled uniformly at random at the beginning of every simulation, and differs per run.

Recall from Section 1.3 that the epidemic on networks can have two supercritical phases, one with a single peak, and one with many peaks, where the system survives for longer. The transition between these phases is stochastic, but undergoing a sharp threshold.

Long-time survival probabilities. The probability of long survival (phase 2M) can be seen in Fig. 10 for eight different models as a function of η, all with the same average degree : two configuration models with and respectively, three geometric inhomogeneous random graphs (GIRG) with and (3.3, 2.3) respectively, and the modified lattice (called ‘grid’ on the figure). All networks have the same number of nodes. We see that the configuration model with the smallest τ has the smallest threshold η, followed closely by GIRG with . The grid and long range percolation (LRP) with have the highest η ( ≈ 0.085 and  ≈ 0.01), these are networks where the average distance is linear. Theory [26], [27] suggest that the threshold η(N) converges in the large network limit.

Fig. 10

A comparison of survival probability of S-I-T-S on eight different networks with the same average degree, as a function of the η, the rate of losing immunisation. Survival probabilities were computed using a 100 runs for all parameter values (η) and models. The infection and healing rates for all cases. CM: the configuration model with GIRG: GIRG power-law fitnesses with and Grid: the modified lattice . All networks have the same number of nodes, and the same average degree . For all networks we see a sharp transition at a critical value of ηwhere the system moves from phase 2P (single peak followed by extinction) to 2M (survival with multiple peaks) with overwhelming probability.

The effect of geometry: long-range edges matter Comparing the S-I-T-S epidemic curves on the configuration model and on GIRGs, with matching parameters (degree distribution and average degree) can be seen in Fig. 11. When we run the epidemic with the same parameters β, γ, η on configuration models and on GIRGs, with the same degree power-law exponent and and same mean degree, we find that the epidemic curves are in good match, as long as the long range parameter α of GIRG is in the interval (1,2) when τ > 3. For τ ∈ (2, 3), hubs dominate the network (see left picture on Fig. 9), and the role of the parameter α is insignificant. This is in accordance with theoretical results (5), where distances are doubly-logarithmic regardless of the value α as long as τ ∈ (2, 3).

Fig. 11

No intervention. Mean field network vs spatial scale-free networks with many long-range edges. A comparison of the number of infected on the configuration model versus GIRG with matching parameters. The continuous line shows the median of infected individuals, the shaded area is covers 95% of all runs for which the number of infected nodes was positive. Each model has average degree 8. Observe that the epidemic curve of configuration model with (CM(2.5), dark blue) matches that of GIRG with (GIRG(2.5,2.3), light blue), while the epidemic curve of the configuration model with (CM(3.3), yellow) matches that of that of GIRG with (GIRG(3.3,1.3), red) Top: (500 days immunity), a single peak can be observed. Bottom: (111 days immunity), many peaks can be observed. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

In Fig. 12, we see that there is only a good match between the a-geometric configuration model and GIRG when τ < 3 or when τ > 3 and α ∈ (1, 2). When τ > 3, the effect of the long range parameter α plays a crucial role. When α > 2, even when keeping the average degree and the degree distribution the same, the curve flattens by roughly 40%. The explanation for this is that in this parameter range, even though hubs are present in the network, they mostly connect to nearby nodes, and thus average distance is significantly larger (see (5)), thus the infection needs more time to spread, see right on Fig. 9.

We explain more about this effect under Intervention B(w), the weak travel restrictions in Section 4.1 below. For the first two scenarios, the height of the first peak is 80 · 103 while in the second scenario, the first peak only has magnitude 50 · 103 nodes, an almost 40% decrease.

The paper [44] studies the effect of the exponent α in more detail on an S-I-R epidemic (corresponding to the first peak here) on a related model, when (4) is replaced byThis is called the max-kernel in the literature, while (4) corresponds to the product-kernel. We mention, however, that the topology of networks generated by the max-kernel is significantly different from that of the product kernel, and the model in [44] falls under a different universality class, see [29, Theorem 1.1].

Intervention strategies.

In this section we study the effect of four different intervention strategies: social distancing, a weak and a strong travel restriction, and limiting the maximal number of contacts per node. We apply each intervention to two initial graphs on nodes: G 1, an GIRG with and average degree 9.7, and G 2, an GIRG with and average degree 8.7. To be able to compare the effect of the different intervention strategies, we choose the parameters of each intervention so the intervention removes about 40% of all connections, resulting in four sub-networks of G 1 and G 2, respectively.

After obtaining the sub networks, we run S-I-T-S on each of them with parameters and varying the values of η from 0.001 to 0.013 at step-size 0.001, corresponding to average immunity period 1/η decreased from 1000 days down to 77 days.

We visualize and study below the effect of each intervention on the S-I-T-S epidemic in five aspects:

the change in the underlying contact network,

a visualization of the spread of the initial infection,

the probability of long-survival as a function of the immunity parameter η,

the height and location of the first and second peak,

the full epidemic curve.

Without intervention

For a baseline comparison, we carry out the above listed four points on the initial networks:

On the top left picture of Figs. 15 and 16, we see similar networks than G 1 and G 2: GIRGs with and (3.3,1.3), respectively. The size is smaller, nodes. Both graphs have approximately the same average degree 4.8.

Fig. 15

Visualization of the change in the underlying contact network under interventions. The network on the top left is similar to G1: a GIRG with average degree 4.8 on nodes. All interventions result in an average degree  ≈ 2.6.

On the top left picture of Figs. 17 and Figs. 18, the visualization of the first 19 days of the epidemic are shown: the epidemic spreads very quickly on G 1, and also quite fast on G 2: (almost) the whole network is infected in both cases. This is explained by the next point, the time of the first peak is about 8 and 14 days, respectively.

Fig. 17

Visualization of the initial 19 days on G1 under interventions. The network on the top left is G1, a GIRG with and average degree 9.6 on nodes. All interventions result in an average degree 5.7. The darker the color, the earlier the infection. We chose typical runs to illustrate the effect of interventions.

Fig. 18

Visualization of the initial 19 days on G2 under several intervention measures. The network on the left is G2, a GIRG with and average degree 8.7 on nodes. All interventions result in an average degree 4.9. The darker the color, the earlier the infection. Color scaling is the same as on Fig. 17.

The probability of long-survival (phase 2M) as a function of the immunity parameter η can be seen in Fig. 19: the top curves on the top/bottom picture show the curves for G 1 and G 2, respectively. Probabilities are calculated based on a hundred runs for each model and value of η. Observe that the survival curves drop for all intervention types, and stay below 0.7 even after the sharp increase.

Table B.1 shows the height and location of the first and second peak for G 1, G 2 (column O, original graph), while Fig. 21, Fig. 22 plots these heights as a function of η.

We see that the location and height of the first peak is insensitive to the change in η, it is  ≈ 84*103 on day 8 for G 1, while  ≈ 79*103 on day 14 on G 2.

There is no second peak on G 1, the ‘star’ shaped dots on the curve indicate the height of the equilibrium. Observe that this grows linearly with η, in accordance with the ODE-solution in (6).

On G 2 there is a second peak, and its height is increasing with η, and its location is decreasing with η. However, there is no clear relation between the location and 1/η, the average immunity duration, see Table B.1.

On Fig. 23, the full epidemic curve () can be seen, on top, the first 60 days while on the bottom picture the first 500 days. We indeed see that for G 1, when τ < 3, the system reaches equilibrium without a well-defined second peak.

Before we start with details, a short summary: social distancing seems to be the least effective intervention in terms of reducing peaks, but at least it does not influence the presence or the height of the second peak.

On the top right picture of Figs. 15 and 16, we see the change in the network after percolation: even though the edge density and the average degree drops, long-range edges still remain present, and the graph looks qualitatively the same as the ones without intervention (top left). The long edges will carry the majority of the infections.

On the top right picture of Figs. 17 and Figs. 18, the visualization of the first 19 days of the epidemic are shown under social distancing: the epidemic spreads somewhat slower, but it still manages to spread to fill the space and infect a large part of the network. Observe that both top-right pictures are qualitatively similar to the top left pictures, where the same spread is shown without intervention.

On Fig. 24, Fig. 25, Fig. 26 and 27 we see the comparison of the epidemic curves: the first peak appears later on the percolated graphs, and their height is roughly 40% less than the height of the original peak. On there is no second peak, (just like on the original G 1), while the second peak on G 2 only drops by 10%, though it appears later. The equilibrium number of infected also drops, roughly by 33%.

Fig. 27

Later effects of interventions on G1, G2, for : Top: G1and all interventions. Bottom G2 and all interventions. The curves are qualitatively similar to the period between second and further peaks differ from those on Figs. 25 and 25. The first peak is not affected by the change of η (below 2% change in height, that could be also due to statistical error margins, so we omit to plot it. The first peak of the red curve (corresponding to truncating long edges) is the first peak of the epidemic. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

On Fig. 19, the long-survival probabilities of epidemics are shown as a function of η. Social distancing (called “perc” in the legend) has its sharp threshold for G 1 and 0.007 for G 2, in this respect it is one of the less effective measures even with respect to survival probability.

Table B.1 shows the height and location of the first and second peak under social distancing (column A), while Fig. 21, Fig. 22 plots these heights as a function of η.

We see again that the location and height of the first peak is insensitive to the change in η, it is  ≈ 52*103 on day 11 ± 1 for G 1, while  ≈ 45*103 on day 22 ± 1 on G 2. Thus the first peak appears later on the percolated graphs, and their height is roughly 40% less than the height of the original peak.

There is no second peak under social distancing on G 1 (the large data are due to fluctuations in the equilibrium system), while on G 2 its height is increasing with η, and roughly of the original graph’s second peak.

On Fig. 24, Fig. 25, Fig. 26 and 27 we see the comparison of the epidemic curves, just like on the original graph, there is no second peak on either, while the second peak on G 2 only drops by 10%, though it appears later. The equilibrium number of infected also drops, roughly by 33%.

Before we start with details, a short summary: strong travel restriction is the most effective interventions in terms of reducing the first peak, but it elongates the epidemic enormously and causes later oscillations and often a higher second peak. The effectiveness of increasing α depends on the node-degree power-law exponent: when τ > 3, its effect is comparable to the strong travel restriction, but when τ ∈ (2, 3), it is not effective or rather, it does not model travel restrictions well.

B(w): Weak travel restriction.

On the middle left picture of Figs. 15 and 16, we see the change in the network after increasing α to a value that results in the average degree to drop to 2.6: for G 1 (when Fig. 15), this intervention only effect node-pairs with and most long-range edges remain present, while many short and medium-scale edges are deleted. For τ < 3, increasing α is not a good model for the no-travel rule. However, for τ > 3 (on Fig. 16), most node-pairs have and thus increasing α does have a similar effect than cutting long-range edges (compare to middle right picture). Observe that both pictures have the same average degree 2.6, yet the middle left picture on Fig. 16) looks more sparse: this is because edges are much more localized in space.

On the middle left pictures of Figs. 17 and Figs. 18, the visualization of the first 19 days of the epidemic are shown under the weak no-travel rule: for G 1 (Figs. 17), the remaining edges are mostly long, and thus the picture looks qualitatively similar to the top left, without intervention. However, for G 2 (Fig. 18), the intervention is very effective, and the epidemic spreads much slower. Within the first 19 days, it only reaches a small fraction of the network localized in various patches.

On Fig. 19, the long-survival probabilities of epidemics are shown as a function of η. The weak no-travel rule, increasing α (last row in the legend) has its sharp threshold for G 1 and 0.007 for G 2, in this respect it is also among the less effective measures for both values of τ.

Table B.1 shows the height and location of the first and second peak under the weak no-travel rule (column Bw), while Fig. 21, Fig. 22 plots these heights as a function of η.

We see again that the location and height of the first peak is insensitive to the change in η, it is  ≈ 50*103 on day 10 ± 1 for G 1, while  ≈ 33*103 on day 29 ± 2 on G 2. We conclude that increasing α only model travel restrictions well when τ > 3. In this case, its effective, the only more effective measure is the strong travel restriction. On G 2, the height of the first peak is 44% of the original peak, while on G 1, its 60%.

We conclude that this intervention is very effective in terms of reducing the first peak on G 2.

On there is no second peak, (just like on the original G 1), while the second peak on is about 80% of the height of the second peak on G 2, though it appears later.

On Fig. 24, Fig. 25, Fig. 26 and 27, we see the comparison of the epidemic curves, On there is no second peak, just like on the original G 1.

We mention that increasing α already has an effect on the epidemic curve when the average degree is kept constant, see Fig. 12 that shows the effect of increasing α on the epidemic curve: for the configuration model () and for GIRG with and herd immunity is reached first when 80000 nodes are infected, while in the second scenario, the first peak only has magnitude 50000 nodes, an almost 40% decrease.

For short average immunity duration, (large η), later peaks can be seen at the bottom picture of Fig. 12, where we observe that the intervention is effective in the sense that the stationary proportion of infected nodes is lower (roughly, 15% decrease: 4.1k vs 4.8k nodes). However, the amplitude of the oscillations do increase when (vs ) and thus the second peak is higher by  ≈ 15% (7.3kvs 6.3k). This is in accordance with intuition: flattening the curve results in an elongated first peak of the epidemic and a higher second peak.

B(s): Strong travel restriction.

On the middle right picture of Figs. 15 and 16, we see the change in the network after applying the no-travel rule with a parameter L so that the average degree drops by 40%: for G 1 (when Fig. 15), already resulted in a 40% drop, while for G 2, was needed. (when Fig. 16). Again, observe that both pictures have the same average degree 2.6, yet the middle left picture on Fig. 16) looks more sparse: this is because edges are much more localized in space.

On the middle right pictures of Figs. 17 and Figs. 18, the visualization of the first 19 days of the epidemic are shown under the hard no-travel rule: for both networks, the effect of cutting long-range edges is drastic: during the course of the first 19 days, the epidemic barely manages to leave a small area near its source (for G 2) and only infects a smaller part of the network (for G 1).

On Fig. 19, the long-survival probabilities of epidemics are shown as a function of η. The hard no-travel rule, (denoted by “long” in the legend) has its sharp threshold for G 1 and 0.01 for G 2, in this respect it is one of the most effective measures in containing the epidemic: even with average immunity duration 111 days, we only observe a single peak for half of the runs.

Table B.1 shows the height and location of the first and second peak under the weak no-travel rule (column B), while Fig. 21, Fig. 22 plots these heights as a function of η. We see again that the location and height of the first peak is insensitive to the change in η, it is  ≈ 39*103 on day 25 ± 1 for while  ≈ 13*103 on day 80 ± 2 on and so the height of the first peak is only 45% of the original peak, while on its only 15%.

However, the first peak appears much later, by a factor 3 later on G 1, and a factor  ≈ 5 later on G 2.

We conclude that this intervention is the most effective in terms of reducing the first peak on G 2, but it elongates the epidemic.

This intervention does cause a second peak to appear also G 1. The height of second peak increases only slightly with η, and thus, for low values of η, the second peak is higher on than on G 2, while for high values of η it is lower. On the height of the second peak is above the equilibrium proportion on G 1 for lower values of η, but not for higher values.

So, the intervention causes the second peak to be rather insensitive to the value of η. It would be interesting to investigate this further.

On Fig. 24, Fig. 25, Fig. 26 and 27, we see the comparison of the epidemic curves: this intervention is the most effective on both G 1 and G 2, the peak is much lower and appears later.

Thus, the price for a (much) lower peak in case of the strong travel restriction is a very elongated epidemic curve.

Before we start with details, a short summary: even though this intervention has comparable effects to social distancing in reducing the first peak, it is causing very high second peak. We gave a plausible explanation for this: removing hubs removes the synchronization from the system, causing high oscillations and a long time to reach equilibrium.

On the bottom picture of Figs. 15 and 16, we see the change in the network after applying the rule of limiting the maximal node-degree. The maximal degree M is chosen so that the average degree drops by 40%: on the bottom picture of Figs. 15, while on the bottom picture of Figs. 15. Observe that limiting the maximal degree, but removing edges randomly to achieve that, results in a network where there are a still many long-range edges, but no hubs remain. Due to the presence of the remaining long edges, the network looks more densely connected compared to the hard and weak no-travel rules.

On the bottom pictures of Figs. 17 and Figs. 18, the visualization of the first 19 days of the epidemic are shown when the maximal node-degree is limited: for both networks, the effect is apparent, but more drastic on G 2, when τ > 3 During the course of the first 19 days, the epidemic grows in a scattered way, mostly using the remaining long edges. The effect is somewhere between the weak and the hard no-travel rule in terms of its capability to restrain the spread.

On Fig. 19, the long-survival probabilities of epidemics are shown as a function of η. Limiting the maximal degree, meaning no hubs, (denoted by “hub” in the legend) has its sharp threshold for G 1 and 0.01 for G 2, so in terms of long survival it is the most effective measure in containing the epidemic (together with the strong travel restriction for G 2).

Table B.1 shows the height and location of the first and second peak under the weak no-travel rule (column B), while Fig. 21, Fig. 4 plots these heights as a function of η. We see again that the location and height of the first peak is insensitive to the change in η, it is  ≈ 49*103 on day 25 ± 1 for while  ≈ 37*103 on day 41 ± 2 on . Interestingly, this intervention has roughly the same effect as social distancing in terms of reducing the first peak, especially on G 1, the difference is a few percentages only: On the height of the first peak is roughly 60% of the original peak, while on its only 50%.

This intervention, unlike social distancing, causes a second peak to appear on that again increases very slowly with η and is roughly 1/5 of the height of the first peak on . The height of second peak on is surprisingly high: roughly a factor 2 higher than that on G 2. Its height, compared to the first peak on is also significant: as big as 1/3 of the first peak.

So, the intervention causes the second peak to be rather insensitive to the value of η. It would be interesting to investigate this asynchronization, the surprisingly high second peak, and its insensitivity to η further.

On Fig. 24, Fig. 25, Fig. 26 and 27, we see the comparison of the epidemic curves: also further peak are quite high, and the system takes a very long time to equilibriate.

Thus, the price for a lower first peak is a very high second peak, followed by quite high later peaks and a long time for the system to equilibriate.

Monotonicity of survival probabilities

On Figs. 10 and 19, some curves do not appear to be monotone increasing, but rather, after the sharp increase they seemingly develop a small peak of survival probabilities. We believe that this is only due to the increased fluctuations around the critical value η. For one model, namely, G 2 under weak no-travel rule (which is an instance of a GIRG with ), we have tested this on a larger scale simulation. We have run 500 runs for each value of η, and plotted the survival probabilities as a function of ηon Fig. 20.

On this larger scale, the curve appears to be monotone, and we believe that this is indeed the case.

S-I-T-S on lattice network models.

Epidemic models on lattice networks have been under investigation for decades, see e.g. the book [61], or the survey [62]. There are close relations to percolation cluster growth models see e.g. [14], [30], and stochastic lattice gas models [18]. In the more applied sciences, similar epidemic models have been developed under the name agent based modeling, see e.g. [60] and references therein, or the survey [62]. The paper [18] conducts a simulation study of a similar model3 to the S-I-T-S model studied here, but focusing on the critical regime (and critical exponents) rather than the supercritical regime investigated here.

Arguable, due to their homogeneous nature, lattice models do not represent realistic epidemics well. Nevertheless, they provide an interesting comparison, especially, when we apply homogenization methods to GIRG, such as the strong no travel rule or limiting the maximal node degree. The S-I-T-S on lattice networks has again subcritical, critical and supercritical phases. Here we briefly summarize our findings for supercritical S-I-T-S on lattice networks. Within the supercritical phase, we distinguish the two sub-phases 2S and 2M as before, see Fig. 13. Between these phases, there is a sharp transition between the two phases.

Fig. 13

A comparison of the number of infected on the torus with in total nodes, for two values of η, see in subfigure captions. The continuous curve shows the median, the shaded area 95% of all runs of all runs for which the number of infected nodes was positive.

However, the shape of the peak is different from the mean-field configuration model and that of GIRGs: the number of infected grows linearly, not exponentially, followed by a linear decline and extinction. The linear growth corresponds to the linear circumference of graph distance balls in two dimensions. The epidemic proceeds like a wave. The number of nodes ever infected grows quadratically at the beginning, corresponding to quadratic volume growth in the model.

Conclusion and outlook

We have studied the first and second peak in an epidemic spreading model with temporary immunity on three different underlying networks, and compared them to the classical compartmental model that is based on ordinary differential equations (ODE). On network models, a new supercritical phase arises: the epidemic dies out after a single peak when the average immunity duration is long enough. When this duration is short enough, the model survives, and reaches a meta-stable equilibrium, similar to the ODE solution. However, as the variability of the node degrees in the underlying network increases, the amplitude of the oscillations after the first peak become smaller, or even vanish, showing that the ODE approximations might not capture some aspects of epidemic spread on networks well enough.

We have used scale-free networks embedded in a geometric space to study interventions methods. The geometry allows for intuitive modeling of several intervention methods: social distancing, travel restrictions, and meeting a limited number of people. We compare the effect of these interventions by ‘permanently’ modifying the underlying contact network, while maintaining the same average node-degree across interventions. We find that the strong travel restrictions are most effective in elongating and diminishing the first peak: the shape changes from exponential to linear. However, travel restrictions and meeting a limited number of people result in a higher second peak, where the latter restriction yields the highest second peak. Our main conclusion is that the behavior of the epidemic highly depends on the topology of the network.

Several future adjustments can be made to make the model more realistic. Other compartments could be added, such as an exposed, quarantined, and death state. Moreover, one could study a scenario where interventions are only applied once certain thresholds in the number of infected nodes are exceeded, or are being lifted. Our framework could also be used to investigate the effect of geometry for other intervention methods, such as a mobile app.

CRediT authorship contribution statement

Joost Jorritsma: Methodology, Software, Validation, Visualization, Data curation, Writing - original draft. Tim Hulshof: Conceptualization, Software, Writing - review & editing. Júlia Komjáthy: Conceptualization, Methodology, Writing - original draft, Writing - review & editing, Supervision, Funding acquisition.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.