Non-analytic behaviour in large-deviations of the susceptible-infected-recovered model under the influence of lockdowns

To allow the sampling of very small probabilities using the large-deviation techniques, a method of manipulating the randomness of the disease outbreaks in a controlled manner is required. For a detailed explanation we refer to a previous publication [35] and just discuss the extension and main idea of the method here.

During a standard simulation the random numbers required to decide whether nodes should become recovered or infected are drawn on demand, usually by calling a pseudo random number generator. Nothing prevents one, however, from drawing these random numbers beforehand and storing them into, here, two random number vectors ξ_λ and ξ_µ , such that for each time step τ and node i there is a corresponding entry in the vectors. Note that we require an estimate of the maximum number of time steps $\tau_{\max}$ that the outbreak is going to last for choosing an appropriate length of these vectors. This will be discussed further in section 5.3.

Additionally we store a vector ξ_η containing all the edges of the investigated graph in randomized order. Let l be the total number of edges in the graph, i.e. the length of the vector ξ_η . Then the pivot point $\theta_\eta = l (1-\eta)$, where η is the fraction of edges affected by the lockdown, can be used to create the locked-down graph by using the first θ_η edges of the vector. A list ξ_P which contains the five initial infected nodes is also maintained.

Thus the outbreak and all measurable quantities are now deterministic functions of the randomness state $\Xi = \left(\xi_\lambda, \xi_\mu, \xi_\eta, \xi_P\right)$.

Our goal is calculating the complete probability distributions P(E), where E is some measurable quantity of our state Ξ, such as C in our case, for a given network and specified parameters. For this purpose we have to employ methods more advanced than typical-event sampling to access states of particularly low probability in numerical simulations. Such large-deviation algorithms have been applied to a variety of models, for example to investigate various graph-[48–51], RNA-[52] and protein-properties [53–55], the Kardar–Parisi–Zhang equation [56–58] and power-grids [59, 60].

Here, we use a large-deviation algorithm, employing a Markov chain Monte Carlo simulation [61] of states given by the random numbers Ξ. The Markov chain $\Xi^{(0)}\to\Xi^{(1)} \to \ldots$ evolves by performing small changes to the given state $\Xi^{(t)}$. The used Markov moves, which we have specifically tailored to our model, are discussed in the appendix. The important part is that the moves allow the system to reach every state from every other state, i.e. ergodicity is fulfiled, and that the moves themselves do not bias the underlying distribution from which the random numbers are drawn. For example if we created a simple sample by drawing all values of Ξ at random and then chose any one of our moves and repeat it k times, the statistical properties our sample would not change.

The probability density functions (pdfs) are then calculated using the $1/t$ Wang–Landau algorithm [62, 63], which is a variant of the original Wang–Landau algorithm [43], that does not suffer from the error-saturation problem of the original algorithm [64, 65].

The basic idea is to initialise a non-normalised probability density estimate $P(C) = 1~\forall C$. Then a Markov-step is performed to generate a trial configuration $\Xi^{\prime}$ from the current configuration $\Xi^{(t)}$. The configurations correspond to the cumulative number of infections $C^{\prime} = C(\Xi^{\prime})$ and $C^{(t)} = C(\Xi^{(t)})$ respectively. They are used in the Metropolis-Hastings probability $\min\{1,P(C^{(t)})/P(C^{\prime})\}$ to decide whether the trial configuration should be accepted, i.e. $\Xi^{(t+1)} = \Xi^{\prime}$ or rejected, i.e. $\Xi^{(t+1)} = \Xi^{(t)}$. A multiplicative factor f > 1 is then used to update the distribution estimate, i.e. $P(C^{(t+1)}) = f P(C^{(t+1)})$, all other entries, i.e. for $C\neq C^{(t+1)}$, are not changed. Note that the Wang-Landau algorithm does not fulfil detailed balance because P(C) changes during the simulation. Still, the $1/t$ variant minimizes this problem [62, 63]. Any remaining bias can be eliminated by subsequent entropic sampling, see below.

During the Wang-Landau simulation, the factor f is iteratively reduced towards 1 following some schedule and the pdf estimate converges to the sought-after density function and just needs to be normalised in the end by demanding $\int_0^1 P(C) \mathop{}\!\mathrm{d}{}C = 1$.

One can split the support of P(C), i.e. interval of allowed values of C, into multiple overlapping smaller intervals and perform an independent simulation for each. This speeds up the simulation [66, 67]. In the end one merges the obtained pdfs, using the fact that the pdfs have to match in the overlapping regions, at least within statistical fluctuations [41, 43].

In this case, however, some sampling issues were encountered around a non-analytic point ('kink') in the distribution when using multiple sampling intervals. This can be circumvented by making sure the kink is far from the interval boundaries, but we mostly opted to just use one interval for the entire range and let the simulation run longer. Only for N = 6400, where a single interval required too much time for our taste, we used more intervals, six to be precise, and applied sampling using Replica-Exchange-Wang–Landau [43, 68–70], which is similar to the described Wang-Landau algorithm but periodically tries to exchange configurations between overlapping intervals, hence the name Replica-Exchange.

The pdfs were refined using entropic sampling [44]. Note that entropic sampling is known to fulfil detailed balance, so any systematic deviation obtained during the Wang-Landau simulation can be corrected in the final part. For details see [35].

We want to study the behaviour of the model subject to lockdowns, so to obtain non-trivial results one should choose the parameters in such a way that the model would be in the pandemic phase if the lockdowns were absent. Working in discrete time, the parameters relating to disease spread are probabilities rather than rates as in the typical continuous time compartmental models. Following the previous work [35] of two of us, we set the recovery probability µ = 0.14. The actual value is rather arbitrary and basically sets the time scale. What remains then is to choose the transmission probability. For the starting conditions with five initially infected individuals, we measured the epidemic threshold, without lockdowns, in the usual way by finding the value of λ which maximises the variance of C. We consider increasing system sizes N and perform a finite-size scaling analysis. This gives a critical transmission $\lambda_c(\infty) = 0.1186(5)$. Thus, we choose λ = 0.2, which comfortably places the system in the epidemic phase.

Having chosen our parameters pertaining to the spread of the disease itself, we need to choose the parameters governing the lockdown.

When activating the lockdown, a fraction η of edges, is blocked, i.e. removed. The lockdown should have a notable effect, so a natural choice is to thin out the edges to the percolation threshold [71] characteristic for the present Small-World ensemble. Clearly, one could also consider to remove a smaller number of edges, which would change the behaviour only slightly. Given the high numerical effort required, however, we concentrate on the case where the lockdowns have the highest effect. Via finite-size scaling we measured the percolation threshold to be $\eta_c = 0.586(1)$, which is consistent with previous results [72]. Thus, we used this value for the severity η.

As for the lockdown threshold θ_l , which states the fraction of infected individuals above which the lockdown is activated, and the corresponding release threshold θ_r , we have chosen to typically use, unless stated otherwise, a constant ratio $\theta_l = 8 \theta_r$.

Using these choices, we measured the average cumulative fraction $\bar{C}$ of infected as a function of θ_l for increasing system sizes. Each data point is averaged over 100 000 networks. Errors are calculated using bootstrap resampling [73]. An example of such a behaviour for N = 3200 is shown on figure 1. Results for other system sizes look similar.

Zoom In Zoom Out Reset image size

Intuitively, increasing the lockdown threshold will increase the cumulative fraction C of infected. In particular, locking down too late can be seen to have no effect in containing the disease, which also makes sense as the lockdown threshold needs to be reached for the lockdown to have any effect. Still, the behaviour is more complex, as we discuss next.

First, note that the presented result is obtained at the maximum resolution with respect to the possible values of θ_l , i.e. with increment $1/N$. Since the release-threshold is at one-eighth the lockdown-threshold, the release threshold only changes in every eighth data point, due to the integer nature of the measured quantities, giving rise to the apparent discontinuities in the early part of the curve.

Clearly, for small lockdown thresholds the lockdown is able to greatly contain and slow down the disease. Here, increasing the threshold leads to an increasing number of infections, giving rise to a peak in the variance, see inset of figure 1, around $\theta_l = 0.02$ indicating the transition to the epidemic phase.

Interesting behaviour emerges as the lockdown threshold is further increased past $\theta_l = 0.05$. In contrast to the simple second order phase-transition behaviour when increasing λ in the no-lockdown case [35], increasing the lockdown threshold gives rise to rather peculiar behaviour, with notable maximum-minimum pairs in the $\bar{C}(\theta_l)$ curve. The positions of these pairs are of interest, and are determined by fitting a Gaussian or a log-normal function to $\bar{C}$ near these points. Such pairs are seen for $\theta_l = 0.0588$ and 0.0692, as well as 0.0955 and 0.1460. The emergence of these extreme points are most easily rationalised by looking at the infection trajectories. In figure 2 we show 1000 such curves $c(\tau)$, as well as the corresponding $i(\tau)$ curves, for each 0.0955 and 0.1460.

Zoom In Zoom Out Reset image size

It can be seen that the lockdown is triggered once in both cases. After the lockdown is lifted, the infection curves both show a notable second wave of infections. Interestingly C reaches a higher final value for $\theta_l = 0.0955$ even though $i(\tau)$ peaks at a higher value for $\theta_l = 0.146$. This is due to the fact that an earlier lockdown ensures a larger proportion of the population is susceptible in the second wave, and hence the disease can spread to more nodes. This can be seen from the $c(\tau)$ curves, where the earlier lockdown leads to a larger cumulative number of infected nodes in the long-run.

Also note the second peak in the variance around 0.17. This is the transition from the lockdown having some effect at containing the disease to having no effect at all.

We want to address these parameter-space regions within the large-deviation simulations later on. To obtain precise limiting values of θ_l for the first and second peak of the variance, we also performed finite-size scaling. For this purpose, we defined the finite-size transition points in the usual way as the peak locations of the variance $\sigma^2(\theta_l)$. We measured their positions by fitting Gaussian-shaped functions around the maxima. Figure 3 shows the position of these maxima in σ² as a function of system size. The main plot corresponds to the thresholds of the second peak in σ², while the inset corresponds to the first peak. To actually perform the finite-size scaling, the function

$$\begin{equation} \theta_l^c\left(N\right) = \theta_l^c\left(\infty\right) + a N^{-b} \end{equation} \tag{ 3 }$$

is fitted to the positions of the second maxima, as shown in figure 3, giving a value $\theta^c_l(\infty) = 0.1683(5)$ for the critical threshold. The other parameters are found to be $a = 0.47(15)$ and $b = 0.52(5)$.

Zoom In Zoom Out Reset image size

The positions of the variance maxima corresponding to the initial transition are shown on the inset of figure 3 as a function of system size. The data is nicely linear on a log-log scale, indicating that here the model follows equation (3) with $\theta^c_l(\infty) = 0$. The other parameters are found to be $a = 12.0(5)$ and $b = 0.778(5)$. Thus, this measurement predicts that the threshold capable of completely containing the outbreak is zero for an infinitely-sized system. We rationalise this by considering the long-range links of the Small-World network. With the system nicely in the pandemic phase with $\lambda = 0.2 \gt 0.1186$, the long range links of the Small-World network allow the disease to spread rather quickly from the beginning to all regions of the network. For this reason, the spread could only be prevented in the thermodynamic limit $N\to\infty$ if the system was already in lockdown when the initial five nodes are infected, i.e. $\theta_l = 0$.

The employed large-deviation algorithm requires a good estimate of the length of the disease outbreak, because this determines the amount of random numbers that need to be controlled by the Markov chain. If the duration was chosen to short it would result in too many unfinished outbreak dynamics which in turn would lead to underestimated values of C. This would lead to a skewed and incorrectly measured density of states. For this reason, before we set up the large-deviation simulations, we investigated the life-span of the disease for various system sizes using standard Monte Carlo sampling.

For each considered parameter set $(N,\theta_l,\theta_r)$, we measured the time τ it took until $i(\tau) = 0$ was reached. We did this for 100 000 networks, respectively, measuring one disease outbreak dynamics each time. From this raw data, we extracted the time at which $98\%$ of the outbreaks are completed, which we denote by $\tau_{98}(N,\theta_l,\theta_r)$. This value is plotted as a function of $\theta_l = 8\theta_r$ for four different system sizes in figure 4.

Zoom In Zoom Out Reset image size

This data is used to decide the duration of the large-deviation simulations as follows. For a given parameter set $(N,\theta_l,\theta_r)$, we use a maximum time of $\tau_{\text{max}} = g \tau_{98}(N,\theta_l,\theta_r)$, where $g \in [2.7,3.0]$ is a factor that we chose on a case by case basis.

Our results presented below show that we are able to sample the distribution P(C) up to a value of C = 1, i.e. including the cases where all nodes are infected. This implies that the chosen time $\tau_{\text{max}}$ is actually large enough. Note that we kept track of the duration of the outbreaks encountered during the large-deviation simulations and counted how often the simulations were not finished when we ran out of numbers. From this we could calculate the frequency $f_{\neq}$ of observing unfinished outbreaks. For the vast majority of simulations this frequency was $f_{\neq} = 0$, i.e. the simulations were long enough. The worst simulation exhibited a frequency of $f_{\neq} = 10^{-5}$, which we still deemed small enough.

Firstly, we observe the pdfs for increasing system size for a lockdown threshold $\theta_l = 8\theta_r = 0.1683$. For system sizes below (and including) N = 1600, we sample the histograms at the highest possible resolution, i.e. with a bin-size of $1/N$. For larger system sizes, we increase the bin-size to $2/N$ for computational efficiency. The probability densities P(C) for a few system sizes are shown in figure 5.

Zoom In Zoom Out Reset image size

We are able to calculate probabilities as low as 10⁻⁸⁵ in the case of N = 6400. For this system size, we needed to calculate C roughly 10⁹ times during the large-deviation simulation. Standard typical-event Monte Carlo sampling, which directly generates random numbers when needed instead of employing a Markov chain, deals with common events. Using a comparable computational effort as the large-deviation method, the typical-event sampling is only able to estimate probabilities on the order of 10⁻⁹, as illustrated in the plot where we have also included results from typical sampling with the same sample size. In the range accessible by typical-event sampling we see a good agreement with the large-deviation data.

The distributions displayed in figure 5 exhibit three peaks. One around C ≈ 0, where the disease quickly dies out, and two for high C. This pair of peaks appears, as we are right at the transition (variance peak on figure 1) controlled by lockdown effectiveness η in parameter space. The system exhibits the highest probabilities for a somehow reduced epidemic with C ≈ 0.7 and for an unaffected spread with a peak around C ≈ 0.9.

Notable are also the seemingly non-differentiable points around C ≈ 0.35 in the pdfs. Investigation of the disease trajectories revealed this to be the point where at least one lockdown takes place almost for sure. This will be discussed in the next section where the number of lockdowns as a function of C is also presented.

Further, we connect our model to large-deviation theory by calculating the empirical rate functions [45, 47], defined as

$$\begin{equation} \Phi\left(C\right) = -\frac{\ln{P\left(C\right)}}{N} + \Phi_0, \end{equation} \tag{ 4 }$$

where Φ₀ is a shift that ensures each of the $\Phi(C)$ have their minimum at $\Phi = 0$.

The obtained rate functions are displayed in figure 6. A convergence with increasing system size is visible, indicating that the large-deviation-principle indeed holds for this model.

Zoom In Zoom Out Reset image size

The previous work without lockdowns [35] found this to also be the case, and clearly the lockdowns retain this behaviour. Therefore, the pdfs are given to first order by $P(C) \propto \exp\left\{-N\Phi(C) + o(N)\right\}$ with o(N) some sub-linear term. Hence, one could potentially make analytical process in regards to P(C) through application of the Gärtner–Ellis theorem [45–47, 74], at least in the region of the convex envelope and where the rate function is differentiable.

Secondly, we now investigate the effect of varying the lockdown and release thresholds. For different values of θ_l and θ_r , we calculated and compared the probability density functions. The particular values of the parameters correspond to the regions visible in figure 1, of which we choose the positions of the minima and maxima of $\bar{C}$. Furthermore, we investigated relatively early lockdowns by choosing points where $\bar{C}(\theta_l)$ is on the first rise. For comparison we also considered the critical value of $\theta_l = 0.1683$ as well as the case with disabled lockdowns. The investigated values are presented in table 1. For actually presenting the pdfs, we have selected a subset of this set for clarity.

With these values in mind, we ideally should only vary one parameter at a time. For this reason, we first fixed the lockdown threshold at the critical value of 0.1683, and varied the release threshold θ_r according to the values presented in table 1, where we also include a simulation where releasing is disabled completely, i.e. non-lifted lockdowns. This is the subject of section 6.2.1. Secondly, we fix the release-threshold at 0.0210 and vary the lockdown threshold according to the (admissible) values from table 1. Finally, we 'compile' these results together by investigating the pdfs of the constant ratio $\theta_l = 8\theta_r$ in section 6.2.3.

6.2.1. Varying the release threshold

Here we study the effect on the pdfs of varying the release threshold. The pdfs are shown on figure 7. It can be seen that the distributions actually align with the no-lockdown case before deviating away at a particular C value specific to each parameter set. These are the C values where the lockdowns become relevant for that particular parameter. The first to deviate is the $\theta_r = \text{disabled} \equiv 0$ case, with the point of deviation, also marked by a non-analyticity, visibly increasing with θ_r .

Zoom In Zoom Out Reset image size

This can be explained as follows. If the lockdown is released the disease will likely be able to propagate through the system better and therefore infect more people compared to the case where the lockdown is not lifted. For rather larger values of C this results in a higher probability, and, due to normalization, for medium values of C in a lower probability. For very small values of C, the lockdown is never triggered.

Interestingly, for all cases the non-analytic point where the deviation between the pdfs for lockdown and not lockdown appears is considerably larger, at least $C\unicode{x2A7E} 0.3$, than the lockdown threshold $\theta_l = 0.1683$ here! This means, the location of this points is not only determined by the value of θ_l but also by the complex network topology.

Furthermore, it can be seen that increasing the release threshold from zero deforms the pdf by increasing the height of the first peak at intermediate C, with alignment around the second peak, i.e. for C > 0.85, indicating the behaviour at high C is independent of the release threshold. We then deduce the behaviour of the pdf is, for large C, likely dictated by the lockdown threshold.

This is emphasised by investigating the average number $\bar{L}$ of lockdowns, see figure 8. Interestingly, for high values of C, lockdowns were not triggered. This means they exhibit values for the numbers of infected below the lockdown threshold, i.e. extremely severe outbreaks are slow-spreading. Also, for the present case with the lockdown threshold fixed at $\theta_l = 0.1683$, the system was found to exhibit typically at most a single lockdown. Note the appearance of a transition between zero and one lockdown. This transition corresponds very well to the nonanalytic point seen in figure 7 and the transition point increases monotonically with θ_r .

Zoom In Zoom Out Reset image size

This can be explained as follows. To reach the lockdown threshold, the system needs a substantial amount of simultaneously infected nodes. If the lockdown is then quickly lifted due to high value of θ_r the resulting value of C will be larger than when the lockdown is lifted late $\theta_r \to 0$. Thus, the $\bar{L}(C)$ curves are shifted to the right. Note that C is not an independent parameter, i.e. cannot be directly controlled, although it is plotted on the x-axis here! For lower values of θ_r the system spends more time in lockdown where the disease spread is greatly limited. Thus, the only way for the system to achieve higher value of C is to have never locked down in the first place, which can only occur inspite of large values of C if the spreading is slow.

6.2.2. Varying the lockdown threshold

The pdfs P(C) for varying the lockdown threshold θ_l , are shown in figure 9. Note that the release threshold is fixed at $\theta_r = 0.0210$. As before, normalisation ensures the pdfs align with the no-lockdown case for small values of C, because as long as $i(\tau) \unicode{x2A7D} C$ holds no lockdown can be triggered. The curves start deviating at a value of C that here depends on θ_l , again accompanied by a non-analyticity of P(C). These points exhibit values of C that are clearly larger than the respective values of θ_l . It is worth mentioning that for low values of θ_l more than one single non-analytic point appears.

Zoom In Zoom Out Reset image size

Also, increasing the lockdown threshold mainly moves the dominant peak to the right. However, with increasing lockdown threshold, we approach the case where the lockdowns fail to affect the system as the threshold is never reached. This gives rise to a small second peak at first, which is visible in the pdf for $\theta_l = 0.146$ in the logarithmic scale. At the critical threshold $\theta_l = 0.1683$ the second peak is of notable magnitude and also visible in the linear scale. Finally, with disabled lockdown, the once dominant peak has vanished completely and only this second peak remains, apart from the peak near C = 0 of course.

We also notice that in contrast to the previous case of varying the release threshold, varying the lockdown threshold is still having a strong impact on the probabilities even for C > 0.85. In particular the behaviour gets richer because low lockdown thresholds result in triggering several lockdown-release pairs. This can be seen in figure 10, where we display the average number $\bar{L}$ of triggered lockdowns as function of C. Again, the points where $\bar{L}$ changes strongly correspond to the non-analytic points of P(C). Indeed a system with a low lockdown threshold must undergo multiple lockdowns in order to attain high C values. By contrast, higher lockdown thresholds allow the disease to spread a bit more rapidly, although to obtain large C values it must still spread slow enough to stay below the lockdown threshold. It is also seen that intermediate thresholds can still suffer a severe outbreak with a single lockdown on average.

Zoom In Zoom Out Reset image size

6.2.3. Fixed-ratio variation

In figure 11 we investigate how the probability distribution functions change with the lockdown and release thresholds simultaneously varied, but with their ratio kept at a constant factor of eight.

Zoom In Zoom Out Reset image size

The left part of the curves aligns with the no-lockdown curve, for the same reasons as discussed before. For low values of θ_l , i.e. quickly triggered lockdowns, the system heavily favours low value of C, with probabilities as small as 10⁻⁸⁰ for high values C. Increasing the lockdown threshold of course allows the system to exhibit an increasingly higher fraction C of infected.

With respect to the pdf's shape, the result for the lowest threshold, i.e. $\theta_l = 0.021$, exhibits multiple peaks of similar heights. This is due to multiple lockdowns and is discussed below. Note that this value of θ_l is close to the peak location of the variance, which corresponds to the epidemic threshold, see figure 1.

Increasing the threshold beyond the critical threshold shifts the behaviour of the system such that it is in a strong epidemic phase. Most of the peaks of the pdf at small values of C disappear somewhere between $\theta_l = 0.021$ and $\theta_l = 0.0588$.

Still, most of the pdfs also exhibit a peak around C = 0.7, which corresponds to a rather large outbreak, that is affected by the lockdowns. A notable exception for this is $\theta_l = 0.0955$, where this peak is slightly shifted a bit to higher values of C, i.e. here this typical outbreak is larger. Starting with $\theta_l = 0.146$ we see a peak around C = 0.9 emerging, which corresponds to the peak without lockdowns. Thus, the lockdown threshold for these systems is so high that the high values of C are reached without triggering the lockdown anymore.

To understand the shapes of the pdfs better, we look at the average number $\bar{L}$ of lockdowns as shown in figure 12. Indeed, not more than one lockdown is triggered for $\theta_l\unicode{x2A7E} 0.146$.

Zoom In Zoom Out Reset image size

Consistently with the two previous cases, low C values exhibit relatively few lockdowns because the disease goes extinct rather early. Still, the case of the relatively low lockdown threshold $\theta_l = 0.021$ creates multiple lockdowns, i.e. several infection waves.

As was the case previously, the maximum number of lockdowns typically occurs around intermediate values of C. For high values of C the slowly spreading infections dominate the dynamics, which leads to relatively few lockdowns again. For high threshold values, even no lockdowns are triggered there.

To understand what is going on for the lowest considered threshold value θ_l in more detail we now look at the conditional probability $P(L|C)$ that L lockdowns are triggered given a value of C. The normalization is such that $\sum_L P(L|C) = 1$ for all values of C. This is displayed in figure 13.

Zoom In Zoom Out Reset image size

Interestingly we can see that for about $C\unicode{x2A7D} 0.45$ the dynamic is dominated by a specific number L. In contrast, for values around C ≈ 0.6, the distributions exhibit has a notable probability for 6, 7 or 8 lockdowns. This region is also characterized by large fluctuations. For even larger values of C, the typical number of lockdowns decreases again.

To analyse the outbreak dynamics further [35], we store a number of the outbreak trajectories during the entropic sampling. We elected to store 200 000 such curves for each pdf. The trajectories are binned according to their corresponding value of C. Let T denote a time series, T = i or T = c, with

$$\begin{equation} T = \left(T\left(0\right), \dots T\left(\tau_{\text{max}}-1\right)\right). \end{equation} \tag{ 5 }$$

Three examples of such curves are shown on figure 14 for $\theta_l = 8\theta_r = 0.021$ and for low, medium and high value of C, respectively. They present the general behaviour in each of the three regimes. For low values of C, typically very short-lived trajectories with one or two lockdowns appear, as seen as the previous section. For intermediate values of C one may experience many lockdowns. For high C typically only a few or no lockdowns occur, depending of course on the parameter set.

Zoom In Zoom Out Reset image size

In the following subsection we use these time series to construct heatmaps to investigate the similarity of pairs of time series. Secondly, heatmaps pertaining to other properties of the outbreaks are presented and discussed.

6.3.1. Disparity maps

To compare the similarity of two outbreaks, i.e. of two time series T₁ and T₂, we first normalise by their respective maxima. The length of the time series is denoted by l₁ and l₂ respectively. Let now $l_{\max} = \max\{l_1,l_2\}$. A distance d between two time series is defined [35] as

$$\begin{equation} d\left(T_1,T_2\right) = \frac{1}{l_{\text{max}}} \sum_{\tau = 0}^{l_{\max}-1} \left|T_1\left(\tau\right)-T_2\left(\tau\right)\right|\,. \end{equation} \tag{ 6 }$$

The disparity $V_T(C_1,C_2)$ is the averaged distance for all pairs of time series T with bin values C₁, C₂, respectively.

For brevity, we only present the disparity V_i of the infection curves $i(\tau)$ here. Figure 15 shows V_i as a heatmap for $\theta_l = 8\theta_r = 0.0210$. Note that using the large-deviation approach has allowed us to create such a heatmap over the entire allowed range of C; particularly also in the range which is inaccessible by standard Monte Carlo sampling.

Zoom In Zoom Out Reset image size

The disparity appears to form regions in the heatmap. Looking at the diagonal, which represents comparing the curves in one bin to one another, for $0\lt C\lt0.2$ the curves are evidently rather similar. These low-outbreak curves are characterized by early lockdowns stopping the spread of the disease. Following the diagonal further the disparity increases and one can no longer clearly distinguish regions from one another. This also makes sense as an increasing number of lockdowns makes it increasingly unlikely that the lockdowns of two infection curves occur at the same time step. Since the lockdowns trigger rapid changes in the time evolution of the disease this increases the disparity between two such time series.

Comparing the time series of any fixed C value, e.g. C = 0.2 with the other time series we can see that they quickly become dissimilar to one another when the two corresponding values of C differ. Looking at figure 13 this makes sense, as the number of lockdowns triggered varies a lot, but has a high correlation with the C value. Also the visual 'steps' in similarity in the range of small values of C correspond to changes in the typical number L of lockdowns.

We also show the disparity V_i for a higher value $\theta_l = 8\theta_r = 0.1683$, around the critical threshold, in figure 16. Here, three regions can be distinguished. Looking at figure 12 we see that the first region $0\unicode{x2A7D} C \unicode{x2A7D} 0.34$ corresponds to no lockdowns, whereas $0.34 \lt C \unicode{x2A7D} 0.81$ corresponds to 1 lockdown and the last region $C\unicode{x2A7E} 0.81$ exhibits no lockdowns again. Note that the border of region one and two corresponds to the seemingly non-analytic behaviour seen on the pdf in figure 11.

Zoom In Zoom Out Reset image size

The time series in the region of $0.34\lt C\lt0.81$ exhibit extremely low disparities to one another, which shows that there is low variability in the time series. This is also true, although to a slightly less degree, for the third region, i.e. $C\unicode{x2A7E} 0.81$. Interestingly, also off-diagonal parts representing the disparity between time series from different regions exhibit an internal structure to some degree. This indicates that the three regions are somehow subdivided further. We do not go into details here.

6.3.2. Conditional density

Other properties of the outbreaks can be studied by considering the conditional densities $\rho(Q|C) = P(Q,C)/P(C)$, with Q some measurable quantity. The time series are binned according to the cumulative fraction C of infections, and then a histogram of Q given C is constructed. These are presented as heatmaps. In this paper, we consider the following quantities for Q:

• $Q = \tau_{\max}$, that is the number of time steps for the infection trajectories to reach their global maximum of $i(\tau)$. If a trajectory reaches the same maximum multiple times we take the time it took to reach the first one,
• Q = M, the relative height of the global maximum, i.e. $i(\tau_{\max})$.

Across the various data sets, some similar patterns emerge in the heatmaps. To illustrate the main points, we only present those for $\theta_l = 8\theta_r = 0.0421$ and 0.1683 respectively.

The conditional density $\rho(\tau_{\text{max}}|C)$ for $\theta_l = 0.0421$ is shown on figure 17. As we saw previously, with these values for the parameters, the system has a tendency to experience multiple lockdowns. The lockdowns will lead to a steep decline in $i(\tau)$. Therefore, if lockdowns occur, the global maximum will be very close to one of the times where the lockdown was triggered, which can be at different times, but typically not at all times. This is revealed in this heatmap by the multiple 'bands' of high probability which are visible for large values of C. Note that with the propagation of the disease, with increasing time, less and less susceptible nodes remain. Therefore, the earlier lockdowns have a higher chance to lead to the global maximum, which is apparent by the colour of the earlier 'bands'.

Zoom In Zoom Out Reset image size

There are multiple discontinuities in $P(\tau_{\max}|C)$. They, like above, correspond to the C values where the dominant number of lockdowns changes.

Figure 18 shows the conditional density $\rho(M|C) = \rho(M,C)/P(C)$ for $\theta_l = 8\theta_r = 0.0421$. Note the phase-transition like discontinuity for C ≈ 0.05, where the maximum M is somehow constrained by the lockdown threshold. For low values of C clearly no lockdowns are triggered. When C exceeds θ_l significantly, lockdowns are triggered and the number of infections is drastically reduced. This leads to a significant increase of P(C) for values of C just above θ_l , as we saw before and to a significant increase of $P(M,C\unicode{x2A7E} \theta_l)$ for values of M also just above of θ_l , i.e. a peak near M = 0.045. Note that in $P(C,M)$ also the small peak for M near 0.01 continues to exist for C larger than θ_l , but due to the normalization by a much larger value of P(C), as compared to $C\lt\Theta_l$, the small peak is not visible any more in figure 18. Thus, the other peak is dominant for C larger than θ_l and appears as a discontinuity.

Zoom In Zoom Out Reset image size

By contrast, to attain higher C the infection curves must spread with more vigour, triggering the lockdowns and having their global maximum restricted by the lockdown threshold.

On this note, with higher lockdown threshold a third phase in the $\rho(M|C)$ surface begins to appear. This is shown on figure 19 for $\theta_l = 8\theta_r = 0.1683$.

Zoom In Zoom Out Reset image size

It can be seen that the behaviour is qualitatively similar to that of figure 18, with the emergence of a new phase around C > 0.8. This is because the relatively high lockdown threshold allows for slowly spreading but long-lasting trajectories with relatively high value of M yet still lower than the lockdown threshold. This is seen from the non-zero probability for M lower than the lockdown threshold in this regime.

Here we detail how we perform the Markov moves. Note that the probabilities with which the moves are performed are not important for the correctness of the algorithm although they do affect the efficiency. During Wang-Landau a Metropolis-Hastings criterium is used to accept the steps. It is a standard rule of thumb to choose the move probabilities such that about 50% of Markov steps are accepted.

With a probability of $92\%$, a standard move is performed, i.e. changes are made to the values of the elements of ξ_λ and ξ_µ . Typically, with in 99% of all standard-move cases, this is done as follows: One of the vectors ξ_λ and ξ_µ is chosen, and a random index k and a random number $\chi \in [0,1]$ is drawn uniformly to set $\xi[k] = \chi$. This choice of the vector and corresponding re-drawing of the random number is repeated B times. The choice of B does not determine the correctness of the algorithm, but rather the efficiency. The convention is to choose B such that roughly ${\sim}50\%$ of trial configurations are accepted.

For the remaining 1% cases of standard moves, the random numbers corresponding to the infection and recovery of the initial infected as well as their immediate neighbours at τ = 0 are all re-drawn uniformly from $[0,1]$. We observed that this strong special move improved the convergence of the Wang-Landau simulation. Since all of the random numbers are uniformly drawn, these moves do not skew any of the underlying statistics.

With a probability of $1\%$, a lockdown move is performed. First we decide how many edges we want to change by drawing a random integer between 1 and 15. For each edge we want to change we then randomly and uniformly choose an index $i\in [0,..,\theta_\eta)$ and $j\in [\theta_\eta,.., l-1]$ and swap the respective edges $\xi_\eta[i] \leftrightarrow \xi_\eta[j]$.

With probability $1\%$ a rotation is performed. The elements of ξ_λ and ξ_µ are shifted by N elements to the right (50%) or to the left (50%) with periodic boundary conditions. This approximately reflects a shift of the trajectory by one time-step in either direction. Note that this can be done very efficiently by not actually shifting the vectors in RAM but storing the offset instead.

With the remaining probability of $7\%$ a patient move is performed. There are two types of patient moves. With a probability of $3/7$ a random patient move is performed, where one of the entries of ξ_P is redrawn by uniformly drawing a new node as initial patient. Note that duplicates are not allowed in ξ_P . Otherwise, i.e. with probability $4/7$ a neighbour patient move is performed. One of the initial infected nodes is chosen randomly. Next one neighbour is chosen, each with probability $1/D$, where D is the maximum degree of the (not-locked-down) network. It is worth mentioning that for nodes that have less than D neighbours not choosing any neighbour is also possible. This ensure detailed balance, and thus, in the long-time limit, all nodes are selected as being infected at the begin of the outbreak with the same probability.

For our simulations we decided to draw the lockdown-graph once per outbreak. This made it easier for the large-deviation simulations, because otherwise we would need to know the maximal number of lockdowns that we can encounter in the rare-event simulation beforehand.

Alternatively one could create a lockdown-graph once per lockdown, i.e. successive lockdowns have independent lockdown graphs. This was implemented for our typical-event sampling simulations, allowing us to compare with the single locked-down network strategy. Naturally, this change only affects outbreak simulations where at least two lockdowns occur.

In figure 10 we see that intermediate values of C exhibit many lockdowns for $\theta_l = 0.0421$, $\theta_r = 0.0210$ and N = 3200 (remaining parameters at default values). Note that this interval is accessible via typical-event sampling, allowing us to compare with the complete pdf presented previously in the region where many lockdowns are experienced. Then, using a sample size of 10⁸ we used this typical-event sampling to measure the pdf of the other lockdown strategy for these parameters. The comparison with our large deviation data is displayed in figure 20.

Zoom In Zoom Out Reset image size

It can be seen that the two pdfs agree very well in the regime accessible via typical-event sampling, showing that any effect of redrawing the locked-down network averages out when calculating the pdf P(C). This is of great advantage—as drawing the locked-down graph only once per outbreak is considerably less expensive than redrawing it for every lockdown.