Integrated vaccination and non-pharmaceutical interventions based strategies in Ontario, Canada, as a case study: a mathematical modelling study

Recently, two coronavirus disease 2019 (COVID-19) vaccine products have been authorized in Canada. It is of crucial importance to model an integrated/combined package of non-pharmaceutical (physical/social distancing) and pharmaceutical (immunization) public health control measures. A modified epidemiological, compartmental SIR model was used and fit to the cumulative COVID-19 case data for the province of Ontario, Canada, from 8 September 2020 to 8 December 2020. Different vaccine roll-out strategies were simulated until 75% of the population was vaccinated, including a no-vaccination scenario. We compete these vaccination strategies with relaxation of non-pharmaceutical interventions. Non-pharmaceutical interventions were supposed to remain enforced and began to be relaxed on 31 January, 31 March or 1 May 2021. Based on projections from the data and long-term extrapolation of scenarios, relaxing the public health measures implemented by re-opening too early would cause any benefits of vaccination to be lost by increasing case numbers, increasing the effective reproduction number above 1 and thus increasing the risk of localized outbreaks. If relaxation is, instead, delayed and 75% of the Ontarian population gets vaccinated by the end of the year, re-opening can occur with very little risk. Relaxing non-pharmaceutical interventions by re-opening and vaccine deployment is a careful balancing act. Our combination of model projections from data and simulation of different strategies and scenarios, can equip local public health decision- and policy-makers with projections concerning the COVID-19 epidemiological trend, helping them in the decision-making process.


Introduction
The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a novel, emerging coronavirus, belonging to the class of enveloped, single-stranded, positive-sense RNA viruses [1]. SARS-CoV-2 is the infectious agent responsible for the still ongoing coronavirus disease 2019 (COVID-19) pandemic. This respiratory disease is generally mild and even asymptomatic, but in a proportion of cases can worsen, becoming particularly severe and even life-threatening [2].
Owing to the lack of effective vaccines and drugs that could prevent the infection and properly treat patients, as well as given the highly contagious nature of the virus and its quick spread, public health agencies have decided to implement non-pharmaceutical interventions (NPIs) [3]. NPIs have become necessary also due to the unprecedented strain on the health systems

Non-pharmaceutical interventions in Ontario
Canada belongs to those countries that have decided to implement a mitigation rather than an eradication strategy. The first confirmed COVID-19 case was reported on 25 January 2020 and, since then, as of 2 January 2020, Canada has experienced a burden of 590 280 cases, with 15 715 deaths. The Federal government and the province of Ontario have enforced an integrated package of public health interventions. In particular, the Ontario province has closed schools and universities on 14 March, declared emergency on 17 March, restricting public events and recreational venues, and closed non-essential workplaces on 24 March. These control measures have been escalated and gradually become more and more stringent, significantly contributing to mitigate against COVID-19 and saving lives.

Vaccine deployment in Ontario
Currently, two mRNA-based COVID-19 vaccines have been authorized by the Canadian regulatory body 'Health Canada'.
On 14 December 2020, Canada started administering the first vaccine doses, according to the federal vaccine plan shown in figure 1. The Ontario Government plan consists of three phases: during the phase I, vaccine doses are limited and healthcare workers working in hospitals, long-term care and retirement homes, congregate care settings and other health facilities as well as members of remote, under-served, either urban or rural, Indigenous communities (such as First Nation communities, Métis and Inuit adults) are being prioritized. During phase II, which is likely to start in Winter 2021, the list of people who can get vaccinated will be further expanded. Finally, during phase III, any person willing to be immunized will receive their vaccine doses. Specifically concerning the Ontario province, more than 20 hospitals will be involved in the immunization campaign.

COVID-19 dynamics model
We use the model and fitting techniques developed in [6]. All data and code are available on a GitHub repository [7]. The modified SIR model is developed to forecast epidemic trajectories using minimal cumulative case data of an outbreak. The model in [6] is used to estimate and quantify the rate of underreporting, the efficacy of preventative measures in healthcare settings, the basic reproduction number and the efficacy of NPIs. The model uses re-scaled time; time is measured in infection incubation periods,t to further reduce the required parameters.
The model is purposefully kept simple to facilitate datafitting. The main purpose of this model, as outlined in [6], is to provide confident estimates of parameter values from which we can forecast forward. In order to do this effectively with modest datasets during an outbreak, it is necessary to maintain control of the number of parameters.
The compartments we are interested in are the number of mild cases, I m , the number of severe cases, I s , the cumulative known cases C K and cumulative incidence C I . Mild cases are those which present no or mild symptoms; severe cases are those which present obvious symptoms or require hospitalization. We make the following assumptions in the model: -The total population is constant.
-Acquired immunity lasts longer than the outbreak.
-There is no co-infection or super-infection.
-The testing/reporting rate in the population is relatively constant. -The probability of a case being severe versus mild is constant.
where N is the total susceptible population, R p is the population reproduction number which differs from the basic reproduction number as we calculate this value in the middle of an outbreak for a specific population; for a discussion on the basic reproduction number and its estimate see [6], p s is the probability of severe infection, r is the average testing/reporting rate, M(t) is a mitigation function, andt is time scaled by the average infectious lifetime 1/μ. Note that Model (2.1) accounts for, but does not assume, that severe infections are less infectious than mild or asymptomatic infections by a factor of p. We include p since in novel or dangerous epidemics severe cases are often managed and quarantined once known. Figure 2 shows a compartmental diagram of Model (2.1).
Model (2.1) is coupled with four initial conditions, derived from three parameters: where R is the number of cases that have already recovered or died by the time we start fitting data.
The function M(t) is a predetermined function that captures the effects of non-pharmaceutical interventions on rates of infection. When fitting, we use a functional form of This allows us to quantify both socially driven relaxation and adherence to non-pharmaceutical interventions.

Vaccine deployment versus social distancing
The primary aim of vaccination is to reduce the number of susceptible individuals in the community. This is accounted for with the change N → N(t). N(t) is defined as a function which starts at N(ŝ) ¼ N, where N is the total population andŝ is the day of vaccine deployment, and approaches the expected vaccinated population N(s*) = vN for 0 ≤ v ≤ 1.
As an example, we employ where L is the proportion of the population we expect to be vaccinated,ŝ and d are shape parameters that determines how quickly the vaccine is deployed.

Data fitting
We fit our model to the cumulative case data for the province of Ontario from 8 September 2020 to 8 December 2020 for the 2020 COVID-19 pandemic; 8 September 2020 was chosen because of the observed significant increase in the number of reported infected cases. Fits were done until 8 December 2020 as at the time of the analyses this was the most current data available. Of course, projections can be adjusted by taking more recent data as it is made available.
To reduce the number of scenarios considered, we fit system (2.1) to cumulative case data prior to vaccine deployment. Using the fitted parameters, we can then explore different functional forms of N(t) to determine the effects of different vaccination strategies on the epidemic curve. We fit the vector of parameters: ( 2 :6) using least squares on the cumulative reported cases per day and the new daily reported cases Figure 2. A compartmental diagram of Model (2.1). From a population N, we assume infection presents as either mild or severe with probability p s . Infection spreads at rate R p M(t). Infection clears at rate 1. All infections are added to two compartments: C I , the cumulative case load of the outbreak and C K the known case load. All severe infections are assumed to be known and reported, and a fraction r of mild infections are known and reported. The compartment C K is highlighted as it is the compartment we use to fit to data.
royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20210009 As described in [6], we bootstrap the fitting by running 1000 fits starting from randomized initial guesses and over random subsets of the data. We then average over all simulations. This process has two benefits: it allows us to perform sensitivity analysis while simultaneously avoiding problems caused by complex surfaces with many local minima.

Outcome of different vaccine roll-out strategies
Using the estimated parameters, we explore different scenarios that could result with vaccine roll-out. We start by defining a logistic roll-out strategy We choose a logistic form of N(t) as it is the simplest functional form that allows for (a) a ramping up of vaccine procurement and deployment and (b) saturation. Later on, we also consider a linear vaccination strategy where daily deployment is constant until 75% of the population is vaccinated.
Using public data about the vaccine roll-out [5,8], we set parameters of N(t). We use L = 0.75 based on a STATCAN survey [8]; meaning we are limited to 75% of the population of Ontario getting vaccinated. Using the projections from [5], we set d and T such that approximately 10% of the population will be vaccinated by March 2021 and all voluntary vaccination will be completed by the end of 2021. This gives us d ≈ 0.433 and T = 24.73, keeping in mind our re-scaled time from the fitting. We first study the outcome of this vaccination strategy assuming that the non-pharmaceutical interventions are to remain in place at the current levels of acceptance and adherence. These parameters are summarized in table 1 and the vaccination curve is shown in figure 3.
Of course, at some point non-pharmaceutical interventions must be lifted. We can compete these two factors by modifying M(t). If we allow for t > t* where t* is when we begin relaxation of nonpharmaceutical interventions. The parameter θ is a shaping parameter that controls the speed of relaxation. Relaxation from the current value of R p to the accepted average basic reproduction number of 2.5 is shown in figure 4 for different values of θ. A Holling Type II function is the simplest function that allows for relaxation over time as well as saturation.
We set τ to be the date relaxation begins and set θ accordingly so that by the end of 2021, non-pharmaceutical interventions are completely relaxed. We maintain that M(t) is a function of time alone so that we are can remain general about how and when these relaxations may manifest themselves. For instance, nonpharmaceutical interventions may relax through policy, fatigue, or apathy as interventions remain in place for long periods of time.       figure 5 are the total cases ( purple), known cases (green), active mild cases (red) and active severe cases (blue). We see that all active cases are resolved by September 2021, given the vaccination strategy used in the simulation. As shown in figure 6, when considering the current infection rate, with the assumption that will approach 75% vaccination of the population by September 2021, we observe that the effect of the vaccine will only be seen after the peak has occurred where the province will be able to ensure that the outbreak ends around July 2020 ( pink line) rather than January 2022, as seen in the case of no vaccine roll-out (blue line).

Relaxation
For the sake of comparison, we present a scenario whereby NPIs relaxation starts on 8 December 2020 without any vaccination plan. This scenario is presented in figure 7 as a worst case scenario stemming from the model. As expected, without any mitigation strategies we would observe a greater number of infected cases, compared to the case where a 75% vaccination rate is assumed with NPIs remaining in place at the current levels of acceptance and adherence.
The figure shows that by removing NPIs without any vaccination, we will see an approximately 80% attack rate in Ontario. This is in line with estimates provided by other models [11,12].  Figure 9 shows that phasing in relaxation starting on 31 March 2021 allows the province of Ontario to maintain control of vaccination roll-out and reopening. In this case, we vary θ and set T = 20.4. We seen in figure 9a,d that if we phase-in reopening slowly then we can maintain control of the outbreak and allow for the realization of vaccination to be seen. Figure 9b,e,f,c shows that increasing the speed of reopening can be detrimental to the control of the outbreak and can outpace the vaccine roll-out. Figure 10 shows the effects of removing all restrictions, immediately (θ = 0) on 1 May 2021. We see that by May 2021, we are able to reopen without risk of significant increase of new cases. In figure 10a, we see that the cumulative caseload does not change significantly, but figure 10b shows a non-negligible increase in new cases as reopening occurs.

A linear vaccination strategy
We can also look at a more hypothetical, but equally valid vaccination strategy employed by the equation which suggests a constant vaccination rate across the population. Figure 11a we see that active cases with such a strategy are not resolved until December 2021, and in figure 11b we see that we continue to observe new cases until August 2021.   (2.8). Shaded areas represent the 95% confidence interval obtained from the fits. The current vaccination roll-out plan will not take hold in the population until after we reach the peak of the outbreak. Known data points are given as black dots.

for different vaccination percentages and non-pharmaceutical interventions relaxation measures
The effective reproduction number is given by It shows the average number of cases one active case will create. Figure 12 shows the effects of the different relaxation/ vaccination strategies on the effective reproduction number. The different coloured lines represent the different scenarios described above. The dashed lines show when we have an R eff < 1 for each strategy. The figure suggests that reopening in May 2021 presents the least risk for a new outbreak to occur. By reopening earlier, we risk R eff increasing above 1 for a longer period of time, creating increased risk for large outbreaks.

Interpretation
We explored the effects of an integrated/combined package of both pharmaceutical and non-pharmaceutical public health control measures. We demonstrated that these two interventions need to be fully harmonized in order to be able to properly control the spread of the outbreak. In the literature, there are some models modelling the impact of COVID-19 vaccination in terms of the epidemiological trend of the pandemic. For example, Libotte and coauthors [13] designed a mathematical compartmental SIR model and explored two major scenarios, according to the optimal control theory, and Jentsch et al. [14] devised a coupled social-epidemiological model, simulating schools and unessential workplaces closures, using an evolutionary game theory-based approach and exploiting mobility data. Some vaccination strategies were age-based, whereas others targeted social contacts, or particularly frail and vulnerable groups. Immunization could curb deaths by 22-77%. Vaccine availability impacted the proportion of preventable deaths: in case of vaccines already available since January 2021, targeting seniors first would be the most optimal strategy, whereas, in case of availability from July 2021 on, targeting social contacts first would result in better control of the outbreak. Similar results were obtained by Grundel et al. [15], who combined an optimization-based control approach with an age-dependent epidemiological model. To be able to relax physical/social distancing measures, a vaccination strategy     targeting high-risk subjects was more effective in the short-term, whereas an approach based on the reduction in social contact rates was the optimal choice in the long-term. In terms of vaccine doses and effectiveness, authors found that greater amounts of a less effective vaccine resulted in more relaxed NPIs with respect to less amounts of a more effective vaccine.
We are able to leverage the current outbreak trajectory through data-fitting to refine scenario-based approaches that compete reopening with vaccination. Current vaccination strategies employed throughout the country are focused on high-risk rather than high-transmission populations. As reopening is, generally, not age-dependent (except school re-openings), and the focus of vaccination is on reducing risk, not necessarily transmission, the questions addressed herein would not benefit from an age-structured model. Moreover, an age-structured model would necessarily increase the number of parameters fitted to unreasonable numbers, given the data available.
The goal is to give potential dates to consider reopening with lowest risk, and timelines of potential herd immunity (in most cases, our results estimate that we will reach minimal cases per day by the end of second quarter 2021). While an age-structured model would be beneficial for developing policies in greater detail, the age-structured case data are sparse at best. In [16], age-related data are given at a weekly resolution and are updated infrequently. While the model and methods introduced in [6] and used here are able to fit agestructured data (even those sparsely provided) it is difficult to extend scenarios from such data at said resolution.
Moreover, our study focused on the provincial level (Ontario), as a case study. At the global level or regional level, few predictive models exist. For instance, by incorporating expert opinions and mathematical models, McDonnell et al. [17] have estimated that due to the time of manufacture and distribution of the vaccine, sufficient doses globally are likely not to be available before September 2023. While our model equips local decision-and policy-makers with projections about the future epidemiological trend of the COVID-19 outbreak. It would be extremely useful to provide public health agencies across the world with a tool that is flexible enough to explore the effect of different mitigation strategies.

Conclusion
Reopening and vaccine deployment is a careful balancing act. Using our mathematical model, we extensively explore the effect of adopting various vaccination and relaxation strategies on the COVID-19 epidemiological long-term projections. Our model shows that while we need not wait for herd immunity to begin relaxation of NPIs, careful monitoring of the situation is required. The model shows that relaxation in Ontario, Canada could be safely implemented as early as May 2021, given a strong vaccination effort. Our linear strategies for vaccination show that reopening would require a significant delay. Our findings are able to provide public health bodies with important insights on the effect  Figure 12. A reduction of R eff with target dates for different scenarios described above. We see that opening too early will cause R eff to increase above 1 which will increase the risk of localized outbreaks. If relaxation is delayed until May 2021 and phased-in, reopening can occur with very little risk.