Analytical approximation for invasion and endemic thresholds, and the optimal control of epidemics in spatially explicit individual-based models

Computer simulations of individual-based models are frequently used to compare strategies for the control of epidemics spreading through spatially distributed populations. However, computer simulations can be slow to implement for newly emerging epidemics, delaying rapid exploration of different intervention scenarios, and do not immediately give general insights, for example, to identify the control strategy with a minimal socio-economic cost. Here, we resolve this problem by applying an analytical approximation to a general epidemiological, stochastic, spatially explicit SIR(S) model where the infection is dispersed according to a finite-ranged dispersal kernel. We derive analytical conditions for a pathogen to invade a spatially explicit host population and to become endemic. To derive general insights about the likely impact of optimal control strategies on invasion and persistence: first, we distinguish between ‘spatial' and ‘non-spatial' control measures, based on their impact on the dispersal kernel; second, we quantify the relative impact of control interventions on the epidemic; third, we consider the relative socio-economic cost of control interventions. Overall, our study shows a trade-off between the two types of control interventions and a vaccination strategy. We identify the optimal strategy to control invading and endemic diseases with minimal socio-economic cost across all possible parameter combinations. We also demonstrate the necessary characteristics of exit strategies from control interventions. The modelling framework presented here can be applied to a wide class of diseases in populations of humans, animals and plants.


Introduction
For emerging diseases, such as COVID-19 [1][2][3][4][5][6][7], and re-occurring diseases typified by annual influenza [8][9][10][11][12][13], there is usually a large variety of control measures of different effectiveness and different socio-economic cost. Given limited resources, it is crucially important to have a good understanding of what combination of control measures would constitute an optimal control strategy [14][15][16][17][18][19][20][21][22][23][24][25], and to be able to identify such an optimal strategy quickly. This requires two key elements: (i) an accurate description of epidemics that allows a general understanding and rapid exploration of different scenarios and (ii) a reliable model of control measures to quantify their relative cost and impact on epidemics.
Rapid and reliable identification of an optimal control strategy for an emerging epidemic is still a major challenge mainly due to theoretical and computational difficulties in describing the spatial stochastic spread and persistence of diseases. Simulations of individual-based models [15,[25][26][27][28] and network models [3,[29][30][31] often are highly detailed and may provide reliable predictions, but can take a long time to provide results. Complex simulation models do not allow quick exploration of different epidemic and control scenarios. Analytical models, in contrast, can provide insights about epidemic principles, e.g. regarding the invasion threshold for a pathogen [32][33][34][35][36][37][38][39], or the level of vaccination needed to stop the spread of the pathogen in a population [40,41]. However, most analytical insights are based on classical, non-spatial epidemiological models [32] that do not account for spatial dynamics. Analytical approximations in the form of 'moment closure' [42][43][44] for individual-based models and 'pair-approximation' [29] in lattice and networks models currently have two important limitations. One limitation is that the approximation schemes are uncontrolled and not guaranteed to provide an exact result in any particular limit [45,46]. An alternative more reliable approximation scheme introduced in Ovaskainen & Cornell [47] and Ovaskainen et al. [48] for spatial point processes and individual-based models provides asymptotically exact results when interactions between individuals are sufficiently long-ranged. The approximation can be applied to a wide class of individual-based models [49], but until now has rarely been applied to epidemiological problems (but see [50] for applications to spatially explicit metapopulations).
A second limitation of moment closure and pair approximation techniques is revealed when estimating the invasion threshold for infectious diseases [43,44]. Bolker, in particular, has shown [43] that spatial moment equations cannot be used to compute invasion eigenvalues for a dynamic system such as an epidemic. Instead, as stated in Brown & Bolker [44], 'to compute the threshold, one needs to compute the spatial structure of the initial phase of the (potential) epidemic'. Accordingly, it is necessary to estimate or assume the local spatial structure [29,44,50] and use that to determine whether a global invasion can proceed. The invasion threshold is therefore estimated not at the start of the epidemic, but at a certain later time when the so-called 'pseudoequilibrium' at the local scale is achieved [44]. Examples of such estimates can be found in individual-based models [44], in network-based models [29] as well as in spatially explicit metapopulation models [50].
In this paper, we adapt the method of Cornell et al. [49] to present a reliable analytical description of epidemics in spatially explicit epidemiological models. Our solution permits rapid exploration of the relative impact of different control scenarios with a view to the identification of the optimal control strategy. Our approach also overcomes the limitations associated with uncontrolled approximation schemes. By exploiting the fact that the invasion is a localized phenomenon when interactions are localized, we estimate the invasion threshold at the start of an epidemic. In the case when re-infection is possible, we also quantify the fraction of infected individuals at the endemic equilibrium. We show that our analytical model allows a classification of a wide class of control measures into spatial and non-spatial measures based on their influence on spatially explicit transmission of infection. This classification leads to a quantified estimate of the impact of control measures on disease incidence, and, therefore, to the identification of the optimal control strategy subject to the known relative cost of different control measures. This paper was motivated by the current COVID-19 pandemic [51,52], and, therefore, for numerical results, we used parameters of the COVID-19 epidemic in the UK as the most relevant example. However, our approach is applicable to a wide class of diseases in populations of humans, animals and plants.

Model
We use the modelling framework introduced in Cornell et al. [49]. The modelling framework [49] formulates individual-based models by stochastic spatio-temporal point processes, derives an exact expression for the moment equations to all orders and, using a perturbation scheme [47,48], provides equations that reliably approximate the effects of space and stochasticity.
We consider the simplest example of a spatially explicit host population where all individuals are stationary and distributed randomly (Poisson process). Within this modelling framework, individuals are distributed with a constant spatial density, n, in infinite two-dimensional Euclidean space. However, one can think of this as representing a population of total size N distributed over a large area A, so that n = N/A. Individuals can be susceptible, infected or recovered, described by their expected densities Q S , Q I and Q R , respectively. Correspondingly, for a real-life population of total size N, these densities are calculated as Q S = S/A, Q I = I/A and Q R = R/A, where S, I and R denote the total number of individuals in the susceptible, infected or recovered compartments.
Individuals change type by one of the following transitions (figure 1a): one infected individual interacting with N susceptible individuals creates βN newly infected individuals per unit of time, where β is the infection rate per contact; each infected individual becomes a recovered and immune individual with rate μ; the immunity to a pathogen disappears with rate γ, i.e. the average duration between the end of infectiousness and the loss of immunity is 1/γ. The SIRS model reduces to a SIR model, when γ = 0. For diseases with confirmed short-lasting immunity, one should use γ > 0, and the model corresponds to an SIRS epidemiological model.
The infection transmission in the model is spatially explicit: the infection rate β is distributed in space according to a dispersal kernel, b(r). This means that an infection of a susceptible individual located at the point x S by an infected individual located at the point x I will occur with the rate βb(r), where r is the distance between individuals, r = |x S − x I | (figure 1b). The dispersal kernel b(r) is defined as a normalized ( Ð 1 0 2prb(r) dr ¼ 1), nonincreasing, non-negative function of the distance r between infected and susceptible individuals. Moreover, the kernel, b(r), is a finite-ranged kernel characterized by the length scale L (e.g. a radius of the Tophat kernel, or a standard deviation of the Gaussian kernel). Thus, within the spatially explicit model with finite-ranged kernel b(r), one infected individual has direct contact only with a finite number of susceptible individuals and does not contact individuals located at distances much larger than L. By contrast, in non-spatial models because of homogeneous mixing each infected individual can infect each susceptible individual in a population directly with the same constant infection rate. In our model, which considers individuals in infinite space, homogeneous mixing is obtained in the limit L → ∞. In such a case, when L = ∞, our model provides the same results as a corresponding non-spatial epidemiological model; see electronic supplementary material, Supplementary Note S1.
Dynamical equations for the expected densities Q S , Q I and Q R are constructed from equations presented in Cornell et al. [49] for the spatial point process 'Infection' (to model the dispersal of infection, i.e. when infected individuals infect susceptible individuals) and for the spatial point process 'Change in Type' (to model the recovery of individuals, and the loss of immunity). The complete closed system of equations describing the dynamics of Q S , Q I and Q R are presented in the royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200966 electronic supplementary material, Supplementary Note S2. According to the definition of the modelling framework introduced in [49], the expected densities of individuals are given by an infinite perturbation series in the approximation of a long-ranged dispersal kernel, i.e. in the approximation of large L. We use only the leading and sub-leading contributions: the leading contribution is identical to a non-spatial analogue of this model, i.e. where the dispersal kernel is infinitely long-ranged (L = ∞); the sub-leading contribution becomes relevant when infection occurs over local spatial scales, i.e. the sub-leading contribution accounts for spatial and stochastic effects in a population of discrete individuals interacting over finite spatial scales.

Basic reproduction number and the invasion threshold
The basic reproduction number R 0 is a powerful but also intuitively clear concept [33][34][35][36]38,39]. It defines the invasion threshold for a pathogen and hence for the invasion of disease: when R 0 > 1 an initial introduction of a pathogen grows, and the pathogen invades the population causing an epidemic; when R 0 < 1 an initial introduction of a pathogen decays. Also, in non-spatial models, R 0 determines the so-called 'herd immunity threshold' [40,41]: the pathogen cannot invade the population of N individuals where N(1 − 1/R 0 ) individuals are immune to the infection. Herd immunity may be achieved by the accumulation of natural immunity or by vaccination. However, it has been shown that the estimates of R 0 based on non-spatial models may break down when stochasticity and discreteness of individuals are taken into account [53]. R 0 is defined as 'the expected number of secondary cases produced by a typical infected individual during its entire period of infectiousness in a completely susceptible population' [33]. Thus, R 0 can be calculated as 'the product of the infection rate and the mean duration of the infection' [35]. Assuming the invasion of a pathogen occurs at time t = 0, the infection rate can be expressed as Q À1 where the term in square brackets represents the rate of emergence of new infected individuals, which can be computed as the net rate of change of infected individuals (infections minus recoveries) plus the rate of recoveries; the pre-factor Q À1 I (0) guarantees that the infection rate is calculated for a unit initial introduction (in non-spatial systems an initial introduction corresponds to a single infected individual per whole finite population). Multiplying the infection rate by the infectious period 1/μ yields the basic reproduction number R 0 : Expressions for Q I (0) and (dQ I /dt)| t=0 are non-trivial in spatially explicit models. In contrast with non-spatial models, a single infected individual within a spatially explicit model interacts not with the whole infinite-sized population, but only with individuals that can be contacted directly via a finite-ranged dispersal kernel, b(x). Even if mathematically the dispersal kernel b(x) is positive everywhere in infinite space, the infection transmissions over sufficiently large distances have an extremely low probability and therefore can be omitted from consideration. Thus, practically, the invasion is a localized phenomenon. Therefore, it is reasonable to consider a local neighbourhood area around an introduced infectious individual, such that all individuals within this neighbourhood area have a sufficiently high probability of having direct contact with the initially introduced infected individual during a single period of infectiousness. For clarity and in order to provide a mathematically transparent and intuitively clear estimate for R 0 , we define a neighbourhood area around an introduced infectious individual as a disk that includes a newly infected individual from that initial infected individual with probability P = 0.95. The radius L * of the neighbourhood area is determined by the equation  royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200966 is also large. Next, considering the whole space as such that is covered by non-overlapping neighbourhood areas, the introduction of a single infected individual to each neighbourhood results in the initial density Q I (0) ¼ 1=(p L 2 Ã ). We assume that the total size of the population as well as the density, n, of the population are conserved; therefore, an infection can be introduced to each neighbourhood area by converting a single susceptible individual in each neighbourhood area into an infectious individual. Consequently, the initial values of other densities are as follows: Q S (0) = n − Q I (0), Q R (0) = 0, i.e. the total density equals n. The initial rate of change (dQ I /dt)| t=0 consists of leading and sub-leading contributions and is expressed straightforwardly from dynamical equations for these contributions, see electronic supplementary material, Supplementary Note S2. Substituting the formal expressions into the definition of R 0 one obtains: whereb(k) andg SI (k,0) represent a Fourier transform of the dispersal kernel, b(x), and the second-order spatial cumulant g SI (x, t) between susceptible and infected individuals in the initial configuration in the system at time t = 0, respectively. We estimateg SI (k,0) approximately by constructing an auxiliary dynamical system where an equilibrium state has the same spatial structure as the initial condition in the main system, i.e. there is a single infected individual per neighbourhood area, on average. Details of this estimation and resulting analytical expression forg SI (k,0) are presented in electronic supplementary material, Supplementary Note S3. By exploring the dependence of expression (2.3) on the range L of dispersal kernel, we find that, in the limit of infinitely long range L → ∞, the only non-vanishing contribution is equal to the basic reproduction number in non-spatial models, R 0 = βN/μ. Notice that equation (2.3) does not depend on the rate γ of the immunity loss; therefore, equation (2.3) is applicable for both SIR and SIRS models.

Endemic equilibrium and endemic threshold
For diseases with short-lived immunity, i.e. when γ > 0, an important characteristic of disease incidence is provided by the density of infected individuals at the endemic equilibrium, which we denote as Q endemic I . The disease is endemic when Q endemic I . 0, and we refer to the condition Q endemic I ¼ 0 as 'the endemic threshold'. Within the analytical framework used in the current analysis, the analytical expression for Q endemic I is obtained as a stable fixed point of the corresponding dynamical equations for the densities of individuals. The resulting expression is shown below: where R (0) 0 ¼ bN=m is the basic reproduction number in nonspatial models.

Control measures
It is natural to classify many control interventions into spatial and non-spatial interventions based on their influence on spatially explicit pathogen transmission. Spatial control interventions reduce the distance over which individuals mix (e.g. the characteristic length scale L of the dispersal kernel b(x)), i.e. they reduce the number of contacts. However, we assume that this does not change the transmission rate at small distances, i.e. βb(0) = constant; figure 1c. Spatial interventions include many forms of restricted movement of individuals including quarantine. Within the context of the current pandemic of SARS-CoV-2, spatial interventions include lockdown [7,9,15] and social distancing acting on a large spatial scale, i.e. restriction of long-distance travelling [54,55]. Within agricultural systems, involving the spread of crop disease, spatial interventions include restriction of movement of inoculum, for example, on farm machinery [56], as well as quarantine in the movement of contaminated produce [57]. The 'strength' of the spatial control is denoted by the parameter s L , where 0 ≤ s L ≤ 1. Thus, we model spatial control measures as such that transform model parameters in the following way: where the scaling of β follows from the requirement βb(0) = constant and the normalization of the kernel b(x). Non-spatial control measures reduce the infection rate, β; they do not change the characteristic length scale, L, of the dispersal kernel b(x); figure 1c. Within the context of an epidemic of SARS-CoV-2, and similar highly transmissible human pathogens, non-spatial control measures include using facemasks [4,5] and other personal protective equipment (PPE), good sanitation, hand washing, surface disinfecting [52,58], distancing between individuals by 1-2-3 m in public referred to as social distancing [31]. Contact tracing and subsequent isolation of infected people [2,3] can be considered as non-spatial control measures since this intervention does not reduce the distance over which undetected infectious individuals mix in the population. Within agricultural systems, non-spatial control measures include the use of pesticides and roguing (removal of symptomatically infected hosts). The 'strength' of the non-spatial control is denoted by the parameter s β , 0 ≤ s β ≤ 1. Non-spatial control transforms the parameter β in the following way: ð2:6Þ The impact of control measures on disease incidence is calculated by their impact on R 0 and Q endemic I according to definitions (2.5) and (2.6), i.e. R 0 → R 0 (s L ,s β ) and Q endemic We incorporate vaccination into the model by removing a fraction of the population from the susceptible compartment. For example, when 70% of the population is vaccinated, 0.7n individuals per unit area are no longer susceptible and play no further role in the epidemic, so the initial spatial density of susceptible individuals is equal to 0.3n. We explore how the impact of control measures depends on the range of the dispersal kernel and the shape of the dispersal kernel (e.g. Gaussian or Tophat, or intermediate case; see electronic supplementary material, Supplementary Note S4 for further details).

Socio-economic cost of control measures
The socio-economic cost of applying control measures during a fixed time interval is assumed to be a function of the control strengths s β and s L . To demonstrate that an optimal control strategy depends on the relative cost of control measures, we consider three hypothetical scenarios with different costs for spatial and non-spatial control measures. Scenario (a): spatial control is moderately cheaper than non-spatial control. Scenario (b): spatial control is moderately more expensive than non-spatial control. Scenario (c): spatial control is much more expensive than royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200966 non-spatial control. For illustration purposes, we choose the following cost functions: In each scenario (a)-(c), we considered all possible combinations of control measures with a constant total cost (equal to the sum of the cost of s L and the cost of s β ) and identified a combination under which an invasion threshold is achieved with the minimal total cost. Details of the mathematical formulation of the problem and its solutions are shown in electronic supplementary material, Supplementary Note S5.

Parameters for numerical calculations
For numerical results, we used parameters of the COVID-19 epidemics caused by SARS-CoV-2 in the UK as the most relevant example: the population size N = 60 million corresponds to the approximate mainland GB population size; the basic reproduction number for SARS-CoV-2 R (0) 0 ¼ 4 [1], where the definition based on non-spatial models was used, R (0) 0 ¼ bN=m; the average duration between onset of asymptomatic infectiousness and the end of infectiousness μ −1 = 4.4 days [5,59]; the infection rate per individual β = 0.15 × 10 −7 , inferred from the equation R (0) 0 ¼ bN=m using default values of N, R (0) 0 and μ; the average duration between the end of infectiousness and the loss of immunity is assumed to be infinite for SARS-CoV-2, γ −1 = ∞, subject to the current absence of evidence for wide-spread loss of immunity to SARS-CoV-2. (We also considered γ −1 = 1 year which is consistent with data for other coronaviruses in human populations [60].) To study spatial effects, we assumed that individuals are distributed in space and interact with each other as follows: the average number of individuals per unit area n = 1; individuals are distributed randomly (Poisson process); individuals interact with each other via a Gaussian dispersal kernel, b(x), depending upon the distance, x, between interacting individuals, bðxÞ ¼ 1=ð2pL 2 Þ expðÀx 2 =ð2L 2 ÞÞ. The main results are shown using the length scale L = 3 (here L is measured in units where the population density is equal to 1); however, we also considered L = 1, 2, and L = ∞.

Results
First, to derive an accurate analytical description of epidemics in spatially explicit individual-based models, we used the unified framework for the analysis of individual-based models introduced in Cornell et al. [49]. We obtained the analytical expression for the basic reproduction number, R 0 , shown in equations (2.1)-(2.3); for diseases with short-lived immunity (γ > 0), we obtained the analytical expression for the density of infected individuals Q endemic I at the endemic equilibrium, equation (2.4).
Second, to quantify the relative impact of control interventions on epidemics, we classified control interventions as spatial and non-spatial based on their impact on dispersal kernel (figure 1c, equations (2.5) and (2.6)), and calculated R 0 and Q endemic I for a range of possible combinations of spatial (with strength s L ) and non-spatial (with strength s β ) control measures. We used numerical values of parameters that correspond to ongoing COVID-19 epidemics in the UK, see results of calculations shown in figure 1c,d.
To identify the optimal strategy for controlling a disease, we first defined a disease as being controlled when the invasion threshold R 0 (s L ,s β ) = 1 is achieved; an alternative definition for the case of endemic diseases is when the endemic threshold Q endemic I (s L ,s b ) ¼ 0 is achieved. Next, we explored properties of these thresholds. We find that the thresholds for R 0 (s L ,s β ) = 1 and Q endemic I (s L ,s b ) ¼ 0 differ in spatially explicit models (figure 1e). Thus, in spatially explicit models the pathogen (and hence disease) can invade a population, but then can fail to become endemic even when re-infection is possible, contrary to predictions for non-spatial models. This difference between thresholds occurs due to the lower density of susceptible individuals when the disease is endemic, making it less probable to find a susceptible individual near an infected individual. However, when the range of the dispersal kernel, L, increases, it is expected that the difference between invasion and endemic thresholds decreases and eventually vanishes. This was confirmed by our model: figure 2a where both thresholds coincide when L = ∞. Earlier studies [44,61] showed that in spatially explicit models where the host density is uniform, the invasion threshold is unchanged whether the interaction is short-or long-ranged. In line with those findings, we find only a weak dependence of the invasion threshold on L (figure 2a(i)). Such a dependence is weak because changing only L increases the local transmission rate due to the normalization of the dispersal kernel, thus the reduction in the number of contacts is compensated by the increase in the intensity of transmission. In contrast with the invasion threshold, the dependence of the endemic equilibrium on L is much stronger ( figure 2a(ii)). This is again explained by the lower density of susceptible individuals at the endemic equilibrium. In principle, a disease can also be controlled by increasing the fraction of the population that has immunity to a pathogen, e.g. this can be achieved by vaccination. We have considered the trade-off of two types of control measures and the level of vaccination required to achieve invasion and endemic thresholds (figure 2b). Vaccination has a strong effect on both thresholds. Unsurprisingly, if a larger fraction of the population is vaccinated, then weaker control measures are required to control the disease. However, the inverse is practically useful when available vaccines or resources are limited: if stronger control measures are in place, then a lower fraction of the population needs to be vaccinated to control the disease. The shape of the dispersal kernel (e.g. Gaussian or Tophat, or intermediate case) has a weak effect on the invasion threshold, but a stronger effect on the endemic threshold, see electronic supplementary material, Supplementary Note S4 for further details.
Finally, the identification of the optimal combination of spatial and non-spatial control interventions depends on the accurate estimation of the relative socio-economic cost of control interventions. We demonstrated this by using three hypothetical scenarios with different relative costs of control interventions (see equations (2.7)). The static optimal control strategy for each scenario is shown in figure 3a-c. It is highly likely that the optimal control strategy would change during the course of an epidemic, due to inevitable changes in the cost of control measures. To illustrate this, we assumed that the socio-economic cost of spatial control measures increases in time (e.g. due to the strong economic impact of lockdown on the economy). Correspondingly, we assume that the socio-economic cost of non-spatial control would slightly reduce in time (e.g. due to the development of new production facilities, or the increase in compliance with non-spatial control regulations). Under these circumstances, the time-dependent  Figure 3. Optimal control strategies. (a-c) Optimal control in scenarios with different socio-economic cost of control measures defined in equation (2.7) and shown here as inset plots. The green curves are the curves of minimum total cost that intersect the R 0 = 1 contour, the point of intersection (located in the centre of the red circle) provides the optimal combination of control interventions, i.e. optimal values for s L and s β . (d ) Assuming socio-economic cost changes in time from scenario a to b to c, the optimal control strategy is expected to change in time as shown by blue arrows; (e) when reducing spatial control (such as exiting a lockdown strategy) the optimal strategy (green arrow) requires simultaneous increase of non-spatial control measures, while a non-optimal control strategy (red) keeps the strength of non-spatial control constant.
royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200966 optimal combination of control interventions shifts towards using only non-spatial control measures (figure 3d), details are discussed in electronic supplementary material, Supplementary Note S5. If, eventually, the cost of spatial control measures, such as lockdown, become unacceptably high, lockdown may be relaxed or lifted entirely (i.e. decreasing the strength, s L , of spatial control measure). Our results show that the reduction of spatial control without a simultaneous increase of non-spatial control inevitably increases disease incidence by increasing the reproduction number; see the 'red' strategy in figure 3e. Alternatively, it is possible to keep the reproduction number constant or even reduce it if non-spatial control interventions increase in strength while spatial control measures are being lifted; see the 'green' strategy in figure 3e.

Discussion
We applied the unified framework of Cornell et al. [49] for the analysis of individual-based models to assess the optimal combination of spatial and non-spatial control interventions taking into account the socio-economic cost of control and the epidemiological impact on epidemics. In the current paper, we considered a simple example of a spatially distributed population of hosts. However, the methods of our paper could be used to model epidemics in populations in which individuals move or aggregate in clusters. The approach can also be applied to a wide class of diseases in populations of humans, animals and plants. For example, the methods we describe here can be readily applied to allow consideration of a large number of different types of individuals and interactions, therefore permitting analysis of, for example, stochastic multitype epidemics [37] of multiple infective strains [62] in heterogeneous spatially explicit populations with age structure [63] and different levels of mixing [63,64]. The work also contributes towards a new generation of analytical epidemic models that combine the following three characteristics that are necessary for realistic and reliable predictions: (i) stochasticity, (ii) spatially explicit dynamics and (iii) reliable mathematically controlled approximations. The strong advantage of analytical models is that they provide insight and a general understanding of what combination of control measures would constitute an optimal control strategy taking account of the relative socio-economic cost of the measure. In addition, analytical models allow a rapid exploration of different intervention scenarios. This has a practical value: for example, the trade-off derived here between the two kinds of control measure and vaccination can suggest an optimal strategy for resource allocation. In particular, we showed that the application of stronger control measures reduces the level of vaccination required to achieve the herd immunity threshold [41,65]. Therefore, when the required number of doses of vaccine is not available or unacceptably costly, the application of stronger control measures can help to achieve the herd immunity threshold. This is valuable in the case of the current COVID-19 pandemic with the availability of vaccines [66]. The results of the methods explored here are also important in light of preparation for future epidemics especially in identifying ways to minimize the cost of future optimal control strategies [67][68][69].