key: cord-0916254-cc2ruask
authors: Eastman, Brydon; Meaney, Cameron; Przedborski, Michelle; Kohandel, Mohammad
title: Modeling the impact of public response on the COVID-19 pandemic in Ontario
date: 2021-04-14
journal: PLoS One
DOI: 10.1371/journal.pone.0249456
sha: 50ae2dc616834443392002706d6fcb116b57679a
doc_id: 916254
cord_uid: cc2ruask

The outbreak of SARS-CoV-2 is thought to have originated in Wuhan, China in late 2019 and has since spread quickly around the world. To date, the virus has infected tens of millions of people worldwide, compelling governments to implement strict policies to counteract community spread. Federal, provincial, and municipal governments have employed various public health policies, including social distancing, mandatory mask wearing, and the closure of schools and businesses. However, the implementation of these policies can be difficult and costly, making it imperative that both policy makers and the citizenry understand their potential benefits and the risks of non-compliance. In this work, a mathematical model is developed to study the impact of social behaviour on the course of the pandemic in the province of Ontario. The approach is based upon a standard SEIRD model with a variable transmission rate that depends on the behaviour of the population. The model parameters, which characterize the disease dynamics, are estimated from Ontario COVID-19 epidemiological data using machine learning techniques. A key result of the model, following from the variable transmission rate, is the prediction of the occurrence of a second wave using the most current infection data and disease-specific traits. The qualitative behaviour of different future transmission-reduction strategies is examined, and the time-varying reproduction number is analyzed using existing epidemiological data and future projections. Importantly, the effective reproduction number, and thus the course of the pandemic, is found to be sensitive to the adherence to public health policies, illustrating the need for vigilance as the economy continues to reopen.

On January 11, 2020, China reported the first death due to a novel strain of coronavirus, now officially named SARS-CoV-2, to the World Health Organization (WHO) [1] . To date, the origin of the virus remains uncertain; however, it is thought to have originated in a wet market in Wuhan, China via inter-species transmission to humans. By January 20, 2020, the number of infections in China had sharply risen to 278, including six deaths, and several cases had also been confirmed in Japan, South Korea and Thailand [1] . The Washington, United States on January 21, 2020 [2] , and four days later Canada confirmed its first case in Toronto, Ontario [3] . As of August 2020, epidemiological data curated by Johns Hopkins Center for Systems Science and Engineering indicates that over 19 million people have been infected globally by SARS-CoV-2 and over 700,000 people have died as a result of complications from COVID-19, the disease caused by the virus. The virus is now affecting people in 215 countries and territories around the globe. To attenuate the spread of the virus, several different strategies for reducing transmission have been implemented in countries around the world. These strategies have ranged from quarantining infected individuals and those returning from international travel, to orders to avoid public parks and facilities, to prohibiting gatherings of more than two people (with exceptions for family members), to mandating the use of masks in public spaces. In addition, in some countries, all day cares, educational facilities, and non-essential businesses were ordered temporarily closed by government officials, and citizens were ordered to stay home, except to go shopping for essential supplies, such as groceries and prescriptions. In South Korea, for instance, a recent report indicated that social distancing protocols were responsible for reducing the effective reproduction number in several cities (below 1 in some cases) by estimating transmission-reducing effects via traffic or metro data [4] .

Data obtained by tracking mobile phone locations [5, 6] indicates that a fair proportion of people are abiding by government implemented virus mitigation policies in North America. In particular, data collected for several weeks between March 2020 and mid-April 2020, when government measures were most strict, indicated that mobility toward the workplace decreased between 24% and 52% in different parts of North America while mobility toward retail and recreation decreased between 29% and 70% [5, 6] . These data raise the questions of how these policies have modulated the timeline of SARS-CoV-2 outbreaks, how abiding by guidelines have mitigated the overload to the medical system, and how sensitive these effects are to the proportion of the population abiding by the policies.

To address these and other questions, we developed a mathematical model to study the effects of transmission-reduction guidelines on the spread of SARS-CoV-2, with an emphasis on the policies implemented in Ontario, Canada as a focused example. First, using machine learning techniques, we fit the temporal evolution of infection counts, total recoveries, and total deaths in Ontario, obtained from Ref. [7] , to an SEIRD infection model. Then we use the fitted parameters as the basis for an analysis of the model behaviour, including predicting the occurrence of a second wave, the qualitative behaviour of candidate future strategies, and a prediction of the basic and time-varying reproduction number corresponding to past data and future projections.

The manuscript is organized as follows. In the Materials and methods section, we describe the model used in this work and give a brief discussion of the techniques that were used to estimate the SEIRD model parameters to fit the temporal SARS-CoV-2 data in Ontario. In the Results section, we present our results and give a discussion of the significant model predictions. Specifically, we investigate the conditions that lead to the emergence of a second wave of infections, as well as how the disease behaves under different future transmission-reduction strategies, and how these different strategies affect the reproduction number of the disease. Finally, in the Conclusion we discuss the results and provide some concluding remarks and future research directions.

The most common method for the mathematical modelling of infectious disease spread involves dividing the population into distinct compartments and quantifying the flow between the compartments. One example of this is the standard SEIRD model that is illustrated in Fig  1. In this model, it is assumed that all individuals comprising the total population fall into one of five categories: susceptible (S), exposed (E), infected (I), recovered (R), or deceased (D). The "susceptible" compartment refers to individuals who have never been infected by the disease and therefore lack antibodies to defend against it. The "exposed" compartment consists of those who have caught the disease from an infected individual, but have not yet become contagious. The "infected" individuals are those who have the disease and are now contagious. The "recovered" compartment refers to individuals who have gone through the whole life-cycle of the disease and survived. For viral infections, it is often assumed that recovered individuals have a prepared immune response and will not become infected again. Finally, the "deceased" compartment represents the cumulative number of fatalities resulting from the disease. In this model, some individuals may never leave the susceptible compartment during the outbreak of the disease, which represents those who do not contract the disease at all. However, once an individual enters the exposed compartment, they will irreversibly progress to the infected compartment and subsequently to either the recovered or deceased compartment. We assume that anyone who has been exposed to enough pathogen as to become contagious will eventually become infectious and progress through the entire life cycle of the disease. Individuals who are exposed to the virus but not enough to become themselves contagious are still captured by the model, as we assume they have not built up enough antibodies so as to trigger an immune response. As a result, those such individuals remain in the susceptible compartment.

In this work, we use a generalized SEIRD model to quantify the spread of SARS-CoV-2 among the population of Ontario, Canada. The model consists of five compartments whose time evolution is described by Eq (1): S(t) is the population of individuals that are susceptible to the virus, E(t) is the population of individuals who have been exposed to the virus but are not yet contagious, I(t) is the population of individuals that are infected by the virus, R(t) is the population of individuals who have recovered and thus have immunity, and D(t) is the number of fatalities from the virus. We further assume that the S, E, and I compartments are stratified into "distancer" and "mixer" sub-compartments. The "distancer" sub-compartment, denoted with a subscript D, refers to those who are adhering to public health guidelines. The "mixer" sub-compartment, denoted with a subscript M, refers to those who are not adhering to public health policies. Mathematically, this means that S = S D + S M , E = E D + E M , and I = I D + I M . The R and D compartments are not explicitly subdivided into distancer and mixer categories since individuals do not leave these compartments once they enter them. While the stratification devised here can loosely be thought of as merely distinguishing between those who are and are not practicing social distancing guidelines, more strictly speaking, it delineates across all public health transmission-reduction guidelines. These include the wearing of masks in public spaces, adhering to a strict social-bubble, frequent hand washing, and other public health measures. 

To simplify the model, we make the assumption that the proportion of people that are adhering to public health guidelines is given by the time-dependent function θ(t). Mathematically, we impose the following conditions on the stratified S, E, and I sub-compartments:

Making this assumption, Eq (1) is reduced to the following form, which we refer to as the Distancing-SEIRD model: 

where

is a time-dependent transmission rate that depends on the proportion of people who are adhering to public health guidelines. The parameters β DD , β DM , β MD , and β MM have been introduced in Eqs (1)-(4) to quantify the rate of virus transmission from infected to susceptible individuals, while taking into account the transmission-reduction practices. Specifically, β DD describes the rate of transmission from infected distancers to susceptible distancers, β DM from infected distancers to susceptible mixers, β MD from infected mixers to susceptible distancers, and β MM from infected mixers to susceptible mixers. It is expected that individuals who adhere to public health guidelines are less likely to come into contact with others and transmit the virus. Thus, we assume that β DD < β DM = β MD < β MM . These assumptions imply that distancers are least likely to pass the disease to one another and that mixers are most likely to pass the disease to one another. Furthermore, by assigning β DM = β MD , we make the simplifying assumption that transmission between an infected distancer and a susceptible mixer is the same as transmission between an infected mixer and a susceptible distancer. Moreover, we assume that this cross-compartment transmission rate falls between the two inter-compartment transmission rates. Similarly, while the pairs of parameters γ M and γ D , α M and α D , and δ M and δ D are theoretically distinct, we assume that each pair is the same across sub-compartments, i.e. γ M = γ D � γ, α M = α D � α, and δ M = δ D � δ. Effectively, we make the assumption that transmission-reduction strategies primarily impact the transmission rate β, but not the latent period γ −1 or the mean infectious length (α + δ) −1 , which are both instead dictated entirely by the disease itself.

The interpretation of each term in the Distancing-SEIRD model is described in Eq (3). Given the assumptions used in this work, the Distancing-SEIRD model depends upon four distinct kinetic parameters. The time-dependent parameter β controls the rate of virus transmission from infected to susceptible individuals. It is determined by the probability of disease transmission as well as the chance of contact, thus it indirectly incorporates the basic reproduction number, R 0 , and it depends explicitly on the proportion of individuals who are adhering to public health guidelines, as indicated in Eq (4). The mean latent period, given by γ −1 , is the average length of time between exposure to the virus and the point at which an individual becomes contagious. The rate at which infectious individuals are removed from the disease, either via recovery or death, is given by (α + δ). Thus the value of (α + δ) −1 is the mean length of time an infected individual is contagious before they either recover from or die from the disease. The five compartments are subjected to the constraint

where N is the total population. Importantly, the state equations are considered to represent direct counts of the population size in this work, not the proportions of the population in each compartment.

In contrast to other approaches that have been proposed to model the effects of social distancing and self isolation on the spread of SARS-CoV-2, see for example Ref. [8] , the model developed in this work captures two infection pathways that weakly interact with each other. Distinctly, the distancing compartment in the model proposed here does not represent a strict quarantine, since we assume that distancing individuals are still engaging in public life while following public health protocols. Furthermore, the Distancing-SEIRD model developed here is a simplistic, first generation model that does not include age, population density, geographical distribution, environmental conditions, loss of immunity, or multiple strains of the virus. Importantly, recent genomic data indicates that a mutated form of the virus, carrying the Spike protein amino acid change D614G, comprises the majority of cases worldwide and has a higher infectivity than the ancestral strain of SARS-CoV-2 [9] .

As powerful as mathematical models are at predicting the spread of infectious diseases, all modelling is subject to simplifying assumptions to remain tractable. In addition to the assumptions stated above, the model developed here does not consider how the burden on the healthcare system directly affects the case fatality rate. In reality, it would be expected that if the healthcare system is overburdened by a large number of simultaneous infections, the case fatality rate should increase. This could result from health care practitioners being unable to provide care for all critically-ill patients due to limitations in hospital capacity and staff or medical equipment and treatments. Consequently, even with effective triage procedures, it is anticipated that more fatalities would result from COVID-19 disease in an overburdened healthcare system compared to one with extra resources, even when the total case count is constant over a period of time.

Further, in this work we ignore the role of random events, and assume that the population is well-mixed and that the spread of the virus throughout the province behaves as a single outbreak. In reality, random events can influence the spread of the virus in profound ways, populations are quite heterogeneous and not well mixed, and, especially in the early stages of an outbreak, the spread of the virus is often characterized by multiple smaller outbreaks happening out of sync with one another. As an epidemic outbreak progresses, the validity of these assumptions changes. For instance, in the early days of a viral outbreak in a particular province, we might observe smaller, discrete outbreaks in multiple cities evolving without any influence from one another; however, as those outbreaks grow in size, they eventually overlap and further behave as a single outbreak.

The Distancing-SEIRD model, Eqs (3) and (4), requires calibration with meaningful values of several mathematical parameters. The parameter values can be estimated by fitting to epidemiological data, then the calibrated model can be used to forecast the spread of SARS-CoV-2 in the population. The epidemiological data we use for calibrating the model is collated by the COVID-19 Canada Open Data Working Group [7] . Effectively, this data set has the form I t , R t , D t for t 2 {0 = t 1 , t 2 , . . ., t n−1 , t n = 157}, where I t i represents the active known infections at time t i , R t i the known cumulative recoveries by time t i , and D t i the cumulative known fatalities at time t i . The time points of the data series correspond to each day between January 25, 2020 (t 1 = 0) and June 30, 2020 (t n = 157). To calibrate the model to this data, we introduce an additional non-model parameter x v , which can be interpreted as the visible infective proportion. This value represents the proportion of active infected cases that have a confirmed diagnosis. In particular, 1 − x v represents the proportion of people who are actively infected but who have not received a positive COVID-19 diagnosis. This number could be due to asymptomatic or mildly symptomatic individuals who never approached public health officials for diagnosis or from results that were false negatives due to errors in the testing. In reality, the number of daily tests performed is time-dependent, so we would expect x v to also exhibit time dependence; however, for simplicity, it is fit here to a constant value.

There are many possible functional forms for θ(t), the function which represents the population's adherence to public health guidelines. Here, we explain our rationale for our chosen form. The province of Ontario declared a state of emergency on March 17, 2020. As a result, we assume that θ(t) = 0 for all time points before March 17, 2020; that is, θ(t) = 0 for t � 52 (since the first confirmed infection in Ontario occurred on January 25, 2020 [3] . For the remaining 105 days between March 17, 2020 and June 30, 2020 we expect θ(t) to vary continuously. This represents the varying response of the citizenry to public health guidelines. This variation in public response could be due to new public health guidelines, the re-opening of various aspects of the economy, or even individual variation for personal reasons. To describe θ(t) for these remaining 105 days we subdivide that time into 7 equal intervals of 15 days each and make the simplifying assumption that θ(t) is a linear spline interpolating the points (τ i , θ i ) for i 2 {0, . . ., 7} where τ 0 = 52 and τ 7 = 157. Hence, θ(t) is a continuous, piecewise-linear function for t 1 � t � t n , described by the following form:

To fit the model to the epidemiological data, we define a loss function that is minimized via differential evolution [10] . For a given parameter set ξ = (γ, α, δ, θ 1 , . . ., θ 7 , β DD , β DM , β MM , x v ), the loss function χ(ξ) is defined as the weighted average of the following quantities: the mean squared error (MSE) between visible infectious cases in the model and in the data, the MSE between recovered cases in the model and the data, and the MSE between fatalities in the model and the data. Explicitly, χ(ξ) has the form

where I(t), R(t), D(t) are trajectories of Eq (3) parameterized by ξ. The weights W I , W R , and W D were chosen to be W I = 10, W R = 1, W D = 1, giving primary importance to deviations from the reported infection counts in the fitting process.

In Fig 2, we present a plot of the epidemiological data along with the fit trajectories of the Distancing-SEIRD model, Eqs (3) and (4), acquired by this fitting procedure. The figure depicts a noticeable difference between visible and total infections, where the visible infections is simply equal to the number of total infections multiplied by the proportion of visible cases, x v . Importantly, in the fitting process, the epidemiological time-series I t is considered to be strictly the visible infections. The full set of parameter values, ξ � , resulting from this fitting procedure is presented in Table 1 . The corresponding predicted public-response function is illustrated in Fig 3 and depicts two prominent peaks. A local sensitivity analysis of the model, conducted around the fit parameter values, is presented in Appendix 1.

As mentioned in the Mathematical Model section, some of the mathematical parameters have a direct biological interpretation. For example, the reciprocal of the parameter γ is an estimate of the latency period of the disease. Based on our fitted estimate for γ, this corresponds to roughly four days, in agreement with other mathematical models [11] and slightly below a report that puts the 95% confidence interval for the latency period at 4.5-5.8 days [12] . Similarly, the procedure predicted values for α and δ that correspond to an infectious period (the mean time from infection to either recovery or death) of roughly 18 days. In Ontario, a nonhospitalized case is marked as recovered exactly 14 days past symptom onset. As a result, the infectious period reflected in the epidemiological data is expected to be greater than 14 days, especially since more serious hospitalized cases can take longer than 14 days to resolve in some individuals. For instance, the World Health Organization reports that the time from onset to recovery is around two weeks for mild cases but is 3-6 weeks for patients with severe or critical disease [13] . Furthermore, the fitting procedure also determined a value for θ 1 on the order of 10 −5 suggesting that there was not a large reduction of transmission during the first 15 days after the province declared a state of emergency. There are a variety of explanations for this beyond just accounting for the latent period of the disease. For example, while March 17, 2020 was the date Ontario entered state of emergency, non-essential businesses were closed in waves on March 24, 2020 and April 4, 2020. The procedure predicted an x v value of roughly 0.73, corresponding to 27% of cases not being reported. Finally, the predicted values of β DD and β DM differ by roughly, 5%, whereas there is a predicted nearly nine-fold reduction in transmission between β DM and β MM . As such, the calibrated model appears to suggest that as long as an individual is mindful of transmission-reduction guidelines, they are protected against transmission even from those who eschew such public health mandates.

One of the key goals in planning future disease-control measures is to prevent the occurrence of a second wave of infections. Moreover, once appropriate public health measures are established, prevention of a second wave requires public adherence to these policies. Thus far, the public response function θ(t), presented in Eq (5), is only defined for t 2 [t 1 , t n ]. However, to mathematically predict a second wave, it is necessary to consider θ(t) for t > t n .

More generally, given a prescribed function of θ(t) for t � t c , we wish to predict the effects of θ(t) for t > t c on the control of an additional wave of infection. While the form of θ(t) for t > t c (hereafter denoted θ(t > t c )) is unlikely to be constant (as human behaviour is unlikely to be constant), to begin, we assume that future θ(t) takes a constant value, denoted θ(t > t c ) = θ f . Then, to impose the absence of a second wave, we search for a condition on θ f that ensures that the number of active cases do not increase for t > t c . That is, we wish to show to find a condition on θ f such that I 0 (t > t c )�0. In Appendix 2 we demonstrate that such a condition is only possible if I 0 (t c )�0. Moreover, we show that the condition separates the interval [0, 1] into Table 1 . Further details of this process and the resulting algorithm for determining the conditions under which such a θ critical exists, and how to arbitrarily approximate this θ critical value, can be found in Appendix 2.

For the fitted parameter values, see Table 1 , that correspond to Ontario's epidemiological data as of t c = τ 6 (June 15, 2020), the critical value of θ is given as θ critical � 0.5326. Importantly, this analysis is not unique to the case of Ontario-a similar number could be obtained for any population, given estimates of the model parameters for that population and recent infection data.

Here we have chosen to take the infection as of June 15, 2020 and project forward as, at that time, the infection count in the model was trending downward. Importantly, as can be seen in Appendix 2, if I 0 (t c )>0, then another wave of infection is inevitable. (3) and (4), and the publicresponse function, Eq (5), obtained by fitting to COVID-19 epidemiological data from Ontario between January 25, 2020 and June 30, 2020 inclusive. To demonstrate the sharp sensitivity of this condition, we plot visible infection trajectories of the Distancing-SERID model, Eqs (3) and (4) Table 1 and the form of θ(t) from Eq (5) with the addition that θ(t) = θ f for t > t c . However, in the first trajectory, we consider PLOS ONE θ f = 0.975 θ critical , in the second trajectory we consider θ f = θ critical , and in the third trajectory we consider θ f = 0.1.025 θ critical .

The scenarios presented in Fig 4 demonstrate the sharp, and non-symmetric, response to perturbations in θ f about θ critical . That is, the number of additional active cases observed after a 2.5% perturbation below θ critical far outweigh the number of cases avoided after a 2.5% perturbation above θ critical . This fragility in the development of a second wave of infections demonstrates the need to err on the side of caution while in social and public settings. The long timescale of the simulations in Fig 4 is a result of how close the value of θ f is to the value of θ critical . For instance, if θ f was instead taken to be 20% larger than θ critical (resulting in θ f � 0.64) the downward slope of the infection trajectory is much sharper, and the model predicts less than 10 active cases in the province by June 2021.

An effective strategy for combating the spread of the virus is one that causes the current infection number to decrease over time and the cumulative deaths and recoveries to plateau. By solving the Distancing-SEIRD model, Eqs (3) and (4), for future θ(t) schedules, where t > t n = 157, we can make qualitative, short-term predictions about the nature of the pandemic beyond the current epidemiological data. In Fig 5, we present the temporal evolution of the infections and deaths for six possible future θ(t) schedules. The schedules were chosen to represent prospective plans that have been considered for the coming months. As discussed above, θ critical can be found at any given point in time to determine whether the infection numbers will increase or decrease beyond that point. For t = t n = 157 days and the parameters in Table 1 , we find θ f � 0.5326, and the importance of this cutoff is reflected in Fig 5. Particularly, as illustrated in Fig 5A and 5B , when the public response is above θ critical , infections die out over time with no second wave, and the decay rate positively correlates with the value of θ(t). Conversely, when there is no adherence to public policy, i.e. θ(t > t n ) = 0, the infections dramatically increase, leading to a large second wave of infections, see Fig 5C. Similarly, we see that when θ(t > t n ) decays exponentially over time, as in Fig 5D, infection counts initially decrease, but begin to increase once θ(t > t n ) has dropped below a new θ critical threshold. An interesting case that has been considered by governments is to "pulse" disease control measures at regular intervals between two fixed values. From Fig 5E and 5F , we see that the effect of pulsing depends on the specific details of the pulsing strategy. For example, when pulsing between θ = 0 and θ = 0.8, overall, there is a net increase in the infections over the span of three months. However, pulsing between θ = 0.4 and θ = 0.8 at two week intervals leads to the infection numbers decreasing over time.

We handpicked several future θ(t) schedules in the previous section. Here, we use the functional form of θ(t), obtained by fitting to Ontario's epidemiological data for t 2 [t 1 , t n ], to predict potential functional forms for θ(t > t n ). To begin, we recognise that Fig 3 illustrates θ(t) exhibits oscillatory behaviour for t > τ 2 . Hence, we first consider the case of maintaining this oscillatory behaviour for future times t > t n . Specifically, we assume that θ(t) = θ(t−T) 8 t > t n where the period T = t n −τ 2 , as illustrated in Fig 6A. This functional form could represent general adherence to public health policies, periodically and temporarily undone by spurious social mixing behaviour. The corresponding infection and fatality trajectories, depicted in Fig  6B, illustrate an oscillatory decay of active cases toward zero in the province over time.

While this type of periodic schedule is intuitive, it is also interesting to consider the case where the oscillations in θ(t), observed in Fig 3, are instead driven by oscillations in the active Table 1 . The blue lines represent the active number of infections, the solid red lines represent the cumulative number of (daily) deaths, and the dashed red lines represent the monthly fatalities. The vertical dashed line is placed at June 30th to differentiate between past and future data. 

infection numbers. This scenario could correspond to the case where citizens hear news about a dip in active infections, which in turn causes the public to eschew public health transmission reduction policies in favour of social mixing behaviour, thus resulting in an increase in active infections. This new increase in active infections is in turn reported on by the news media resulting in increased vigilance in the populous and stricter adherence to public health policies, eliciting a decrease in active infections, thus causing the cycle to continue. In this hypothetical scenario, we expect the function θ(t) to effectively depend on the infections, such that θ(t) = θ(I(t)). Moreover, given the delay between exposure and symptom onset and the delay between symptom onset and confirmed diagnosis, we would anticipate a temporal delay in the dependence of θ on I. That is, θ(t) = θ(I(t−L)) for some time lag L.

To determine the value of the time lag, we calculated the Pearson correlation coefficient between θ(t) and I(t−L) for several values of L, using the epidemiological data for t < t n . The results are presented in Fig 7 and show that the maximum correlation is obtained when L = 6, which results in an R 2 value of approximately 0.753. We next assumed a linear dependence for θ, with the form θ(t) = min(max(0, mI(t − 6)+ b), 1). Using linear regression, we fit this functional form to θ(t) for t � t n , which gave approximately m = 1.9949 × 10 −4 and b = −3.9818 × 10 −1 , with a correlation coefficient of approximately 0.868. The resulting response function is depicted in Fig 8A, Importantly, we observe that while this response function is adequate for flattening the active infection curve, it is not enough to fully eradicate the virus in the population on its own.

The basic reproduction number, R 0 , is the number of expected secondary cases caused by a single infectious individual placed in an otherwise wholly susceptible population. Broadly, the effective reproduction number, R t , represents the expected number of secondary cases seeded by an initial case at time t during the pandemic. The effective reproduction number depends upon properties of the virus and disease itself as well as upon the behaviour of the citizenry in response to public health policy. This number is often argued to be the most important number to track in order to manage a pandemic outbreak. If R t < 1 is sustained, then transmission is slowing and the virus is likely to die out, however if R t > 1 then transmission is increasing within the population. As the pandemic has progressed, local public health officials have enacted various policies to control this number, including social distancing procedures, mandatory face-covering bylaws, and business closures.

Here we used a Bayesian method that was developed by Bettencourt and Riberio [14] to estimate the effective reproduction number R t from Ontario's epidemiological data. Following this approach, we assumed that new infection cases arise according to a Poisson distribution with mean value equal to the inverse of the serial interval. We further assumed that the distribution of R t is a Gaussian centred around R t−1 with standard deviation σ = 0.1. Additionally, we made the assumption that the serial interval, for the purpose of calibrating the Poisson distribution, is 4 days [15] . The estimated R t , obtained using this method, is presented in Fig 9 for the SARS-CoV-2 outbreak in Ontario between February 25th, 2020 and June 30th, 2020. While Ontario had its first confirmed case of COVID-19 on January 25th, 2020 the province did not see a sustained increase in case count until late February. The method for estimating R t 

that we employed requires more than 1 new case a day in a 7-day moving average, as a result we cannot provide estimates for R t for dates earlier than February 25th, 2020 in Ontario.

We also used the Bayesian method, along with the Distancing-SEIRD model, to predict the trends in R t for the future θ(t) functions presented in Fig 5. The results are depicted in Fig 10 and show that the effective reproduction number is quite sensitive to the public response. Specifically, in Fig 10A, which corresponds to maintaining θ(t) = 0.6 for t > t n , R t ends up slightly below 1 as time progresses, meaning that active infections decay to zero. This result is expected since θ critical = 0.6. In comparison, in Fig 10B, the public response schedule is θ(t) = 0 for t > t n , corresponding to limited public health containment measures. We see that in this case, R t spends the entire months of July and August 2020 at around R t = 1.5. In September, R t starts decreasing until crossing the R t = 1 line around October 2020. In Fig 5B, there is a rapid increase in active cases before a peak is reached around October 2020. Thus, the qualitative change in R t around October 2020 in Fig 10B does not correspond to a change in public behaviour, since the public response function remains constant at θ(t) = 0, but rather corresponds to a sufficient (but not total) depletion of the susceptible pool resulting in a decrease in transmission below the R t = 1 bifurcation. Finally, in Fig 10C, corresponding to the public response schedule of θ(t) = 0.8 for t > t n , the behaviour of R t is similar to the case of θ(t) = 0.6 in Fig 10A. Notably, with higher public response, long-term R t values The cases where θ(t) was taken to be time-dependent for t > t n are presented in Fig 10D-10F . We see from Fig 10D that an exponential decay in θ(t) with a half-life of two months leads to a steady increase in R(t) above R(t) = 1 after mid-July 2020. However, the predicted R(t) values in this case remain smaller than those corresponding to θ(t) = 0, Fig 10B, until approximately September 2020. Fig 10E and 10F both correspond to oscillatory schedules for θ(t) for t > t n , and this periodic behaviour is reflected in the predicted trend in R t . However, while the tails of both R t plots are oscillatory, only the public response schedule in Fig 10F is capable of reducing the transmission sufficiently to cause a long-term decay in cases, see Fig 5. In contrast, the oscillations in the tail of R t in Fig 10E result in enough time intervals for which R t > 1 and not enough time intervals for which R t < 1, culminating in an overall long-term increase in cases.

Mathematical models of pandemic progression can serve as useful tools for policy-makers who are crafting a societal pandemic response. If presented in a widely accessible nature, they can also be used to persuade the general public on best practices for controlling the spread of the disease. One of the main tactics the public can employ to combat the SARS-CoV-2 virus is social distancing: limiting the amount of social interaction in order to lower the effective transmission rate of the virus. Indeed, we saw in this work that the effective reproduction number of the virus is highly sensitive to the public response. Unfortunately, the widespread use of social distancing and other mitigation techniques can cause many economic, social, and personal challenges, leading some to question the value of social distancing in achieving pandemic control. Naturally, questions arise regarding how many people must social distance, and for how long, in order for these tactics to be effective. Effective reproduction number, R t , in Ontario. Plotted for t � t n = 157 days, estimated directly from epidemiological data using the Bayesian method described in Ref. [14] . Grey area represents the 95% confidence interval about the mean values in black.

https://doi.org/10.1371/journal.pone.0249456.g009

In this work, we presented a mathematical model that is capable of describing the infections, recoveries, and deaths caused by SARS-CoV-2, and can therefore shed light on a number of these questions. The model developed here is a generalized SEIRD model that includes segmentation of the total population into two groups, where one group is adhering to public health measures, while the other is not. The proportion of the population that is following 

public health guidelines can be controlled so that time-varying strategies can be examined to make short-term predictions on the spread of the pandemic. This model can be extended to incorporate several different sub-populations, each adhering to a different degree of public health measures. In addition, stochastic effects can be incorporated into the model to account for random events and fluctuations due to individual differences in transmissibility and predisposition toward severe disease conditions.

A key result of this work is the development of a method to predict the appearance of a second wave of infections, given the current infection numbers and the disease-specific parameters. Importantly, the appearance of a second wave is found to depend sensitively on the public response. Indeed, we showed in this work that a slight decrease from the critical value of θ(t) leads to a significant second wave of infections that continues to peak for many months into the future. Importantly, the method for calculating θ critical developed here can be used as a benchmark to calculate the approximate level of adherence to public health policies that is necessary in a population to effectively control the viral spread and lead to a prolonged decrease in active infections.

Importantly, the significance of θ critical in distinguishing the nature of the future pandemic progression was reflected in our projections of the pandemic. We saw that maintaining supercritical public adherence, i.e. θ(t) > θ critical , leads to a decrease in active infections, while reducing public adherence below θ critical can lead to an exponential increase in infections and a surge in case fatalities over a period of several months. Increasing θ(t) well beyond θ critical leads quickly to a decline in active cases and a saturation in case fatalities. This means that a near perfect adherence to public health advice could efficiently dampen the pandemic over a period of months. It would be interesting to consider how these results are impacted by spurious injections of infected individuals into the population. This extension of the model could correspond to, for example, infections that are imported into the population due to national or international travel.

In this endeavour, we also saw that time-varying public health measures can be effective at reducing the active infections and flattening the case fatality curve. One particularly effective strategy supported by this work is the pulsing of public response measures, at two-week intervals, between two different values of public adherence. We showed that while this strategy leads to an oscillatory response in the effective reproduction number of SARS-CoV-2 in Ontario, by a suitable combination of the two values, R t can on average remain below 1. This situation leads to a decay in active case numbers over a period of a few months, while still permitting a degree of socialization. These results are significant because, as we saw in the calculated effective reproduction number that was determined directly from Ontario's epidemiological data, R t has been hovering around 1 since approximately April 2020. Public adherence to health guidelines will therefore be important, more than ever before, as schools reopen and child care centers move to full capacity, in the midst of the flu season this Fall. If an increase in active infections should ensue, our work provides a method to project the progression of the pandemic under different public responses. These results can be used to devise an efficient plan that can dampen infections, while still permitting the economy to operate at some capacity.

Given the impact of the public response on the spread of the pandemic, we also investigated potential methods for predicting future population behaviour from previous behaviour. We considered the case in which the oscillations that were observed in the epidemiological data were spurred by actions that would remain constant even as active cases decrease. For instance, the oscillations in θ(t) from Fig 3 could be explained as a natural waxing and waning of responsible pandemic response from the populace. Notably, during the time period over which these oscillations occur, the official public health policy in Ontario was constant. That is, the alternating peaks and valleys could represent general adherence to public health policies periodically undone by spurious social mixing behaviour as the public's resolve weakens instead of being driven by direct changes to public health policy. This lack of adherence is only temporary however, as the public eventually rebounds back to following public health policies. Interestingly, in this case, maintaining the oscillatory trend was found to cause a slow oscillatory decay in active cases and a saturation in case fatalities in our pandemic projections for Ontario.

We also considered the scenario in which the observed oscillations in the Ontario public response up to June 30, 2020 were driven by oscillations in the infection numbers. As discussed above, this scenario could correspond to the case where people hear news about a dip in active infections, which in turn causes them to engage in social mixing behaviour. This causes an increase in active infections, which is in turn reported on by the news media. The result is a temporary increase in vigilance in the populous and stricter adherence to public health policies, which further elicits a decrease in active infections, thus causing the cycle to continue. We saw that this scenario quickly leads to an increase in active infections, ultimately saturating at a value that is higher than Ontario's peak value prior to June 30, 2020. The pandemic is sustained in this case and the predicted fatalities continue to increase over time. These results reflect the importance of remaining vigilant in following public health guidelines, even when the infection numbers are down and give the false impression that it is safe to engage in social mixing behaviour.

We hope that this work can assist public health officials in crafting policies to lessen the severity of the pandemic, to inform the public on the potential outcomes of the pandemic, and to serve as an impetus for those who are skeptical of following evidence-based public health guidelines. Importantly, we stress that the results presented in this work are estimates that are based upon a simplistic mathematical model that relies on imperfect and incomplete data. Our results should be interpreted qualitatively as an approximation of the possible outcomes of the pandemic given different levels of public adherence to public health measures. A more comprehensive mathematical analysis of the pandemic, which could potentially lead to more accurate quantitative predictions, would require a significantly more complex model and a larger set of data to determine key disease parameters. Indeed, above we have suggested several potential extensions to the model. In addition, the model could be further extended to incorporate the age and sex of people in the population, as well as to account for underlying medical conditions. It would also be interesting to extend the model to include influenza infections, since the flu season typically begins in the Fall and will likely lead to unforeseen challenges during the SARS-CoV-2 pandemic response.

where subscripts denote nominal values and R x,p depends on the time through the population compartments.

In Fig 11, we present the calculated relative sensitivities for each of the fitted parameters and non-zero initial conditions in the model. The sensitivities are calculated at the end of the 52 day simulation in Fig 11(A) , which corresponds to the date that state of emergency was declared in Ontario. We see here that the susceptible population S(t) is least sensitive of all population compartments to small changes in the parameter values and non-zero initial (3) and (4), to perturbations in the parameter values and initial conditions. Each parameter was perturbed by + 1% and initial conditions were perturbed by + 1 from their nominal values in Table 1 . In all plots, the public response function θ(t), for t � t n = 157 days, is fit from epidemiological data, see 

conditions. Moreover, the other populations exhibit a fairly strong overlap in their sensitivity values, with the smallest sensitivity to changes in the initial values of E(t) and I(t). The deceased population D(t) is most sensitive to the parameter δ at this time point. In contrast, we see from Fig 11(B) that further into the pandemic, the sensitivity of D(t) to changes in δ is lessened and approximately equal for the E, I, R and D populations. Apart from this parameter, α and β MM have the highest sensitivities among all population compartments. Moreover, due to its large size, S(t) is still the least sensitive population compartment by orders of magnitude. Interestingly, the E and I compartments show at least approximately one order of magnitude increase in sensitivity to θ 6 and θ 7 at this time point, compared to the R and D populations.

We consider sensitivity projections, beyond the epidemiological data, in Fig 11(C) and 11(D), where the future response function is taken to have a constant value. We see that in the case where future θ(t) values are taken to be zero, corresponding to no adherence to public health advice, the susceptible population is thoroughly depleted and becomes more sensitive to perturbations in parameter values and initial conditions than the deceased population. In addition, the exposed population is the most sensitive of all population compartments for all parameters. In comparison, a high adherence to public health policies, such as future θ(t) = 0.8 as in Fig 11(D) , leads to a sensitivity pattern that is nearly identical to Fig 11(B) , where θ(t) was fit from the epidemiological data. This pattern is marked by a highly insensitive susceptible population, shown closest to the center of the plot. Unlike Fig 11(B) , the E, I, R and D populations nearly overlap for all parameters, including θ 6 and θ 7 .

Importantly, the results in Fig 11 suggest that the relative sensitivities exhibit a qualitatively different trend when the pandemic is winding down compared to when the infection and death counts are high. Particularly, when the pandemic is under control, the susceptible population is nearly insensitive to changes in the parameter values and sits at the center of the radial plot, distinct from the E, I, R and D populations. However, when the infections are high, the depletion of the susceptibility pool drives up its sensitivity values, resulting in a near overlap in sensitivity trends between all population compartments. This trend is also observed when the future response functions are taken to be time-dependent, as in Fig 12. Specifically, inspection of the results in Fig 12 reveals that only plots (B) and (D) exhibit the qualitative trend in which the susceptibility pool is robust to changes in the parameter values. Indeed, only the future public response functions chosen in plots (B) and (D) lead to containment of the pandemic, see Fig 5. On the other hand, the other two public response functions lead to an explosion of infections, see Fig 5 , and this behaviour is reflected in the trends in the relative sensitivity. Notably, the susceptible population exhibits a much higher sensitivity to all parameter values in plots (A) and (C). We thus see that the qualitative trends in the relative sensitivity of the susceptible population provides a way to determine the state of the pandemic.

Suppose that I 0 (t c )>0, then there exists an open interval O containing t c such that I 0 (t)>0 for all t 2 O. As a result, I(t)>I(t c ) for some t > t c . Hence, if I 0 (t c )>0, then there will always be another wave of infection regardless of the choice of θ(t) for t > t c . For the remainder of this section we suppose that θ(t > t c ) = θ f for some θ f and that I 0 (t c )�0. We wish to show that if another wave of infection can be prevented for t > t c , then it can be prevented entirely by (3) and (4), to perturbations in the parameter values and initial conditions. Each parameter was perturbed by + 1% and non-zero initial conditions were perturbed by + 1 from their nominal values in Table 1 . Simulation time in all plots is 300 days. Moreover, in all plots, the public response function θ(t), for t � t n = 157 days, is fit from epidemiological data, see 

Novel Coronavirus (2019-nCoV) SITUATION REPORT-1; 2020

Novel Coronavirus (2019-nCoV) SITUATION REPORT-3; 2020

Novel Coronavirus (2019-nCoV) SITUATION REPORT-7; 2020

Behaviour of infections for θ f 2 {0.55, 0.7, 1.0} when I 0 (t c )>0. This plot demonstrates that no choice of θ f is sufficient at preventing an additional wave of infection as even maximum adherence, θ f = 1, is not enough to prevent I(t) from increasing. Moreover

Potential roles of social distancing in mitigating the spread of coronavirus disease 2019 (COVID-19) in South Korea

COVID-19 Community Mobility Report, Canada

United States

Open access epidemiologic data and an interactive dashboard to monitor the COVID-19 outbreak in Canada

Epidemic analysis of COVID-19 in China by dynamical modeling

Tracking changes in SARS-CoV-2 Spike: evidence that D614G increases infectivity of the COVID-19 virus

Differential Evolution-A Simple and Efficient Heuristic for Global Optimization over Continuous Spaces

Novel coronavirus 2019-nCoV: early estimation of epidemiological parameters and epidemic predictions

The incubation period of coronavirus disease 2019 (COVID-19) from publicly reported confirmed cases: estimation and application. Annals of internal medicine

Report of the WHO-China Joint Mission on Coronavirus Disease

Real time bayesian estimation of the epidemic potential of emerging infectious diseases

Serial interval of novel coronavirus (COVID-19) infections. International journal of infectious diseases

The authors want to thank the editor and the reviewers for their many helpful comments.

Conceptualization: Brydon Eastman, Cameron Meaney, Michelle Przedborski. 

The authors have declared that no competing interests exist.

Local sensitivity analysis. Here we investigate the sensitivity of the outputs for the Distancing-SEIRD model, Eqs (3) and (4) , to perturbations in the nominal parameter set presented in Table 1 . To compute the sensitivities, we change the parameter values and the initial values of E and I one-at-a-time by a small amount, Δp. Here we take Δp to be + 1% of the nominal parameter value p 0 . Since the quantities E 0 and I 0 represent discrete counts of the population, we take the change Δp to be + 1 for these quantities. Then the relative sensitivity of each model population x = {S, E, I, R, D} for the parameter or initial condition p is calculated as follows:

Details for the prediction of another wave of infection and an algorithm for determining θ critical . If the function θ(t) is discontinuous, then the functions S(t), E(t), and I(t) remain continuous. However, as can be seen in Eq 3, the functions S 0 (t) and E 0 (t) would be discontinuous for discontinuous θ(t). In the case when γ D = γ M , α D = α M , and δ D = δ M then by Eq (3) the function I 0 (t) = γE − (α + δ)I is entirely independent of θ(t). Hence, for discontinuous θ(t) we see that S(t) and E(t) are continuous functions and I(t) is a continuous differentiable function. This is relevant because we may introduce a jump discontinuity in θ(t) at the point t c in order to prevent another wave of infection, however I will remain C 1 even in the presence of this discontinuity.choosing θ f � θ critical for some unique value of θ critical 2 [0, 1] and, conversely, if we choose θ < θ critical , then another wave of infection would be inevitable.ConsiderAs a result, for fixed t we have @ θ E 0 (t) = β 0 (θ)SI < 0 and @ E I 0 (t) = γ > 0, hence @ θ I 0 (t)<0. So if a θ critical exists, then it is unique. To see if a θ critical can exist for an arbitrary parameter set, we first take θ = 1, solve the system in Eq (3), and use the result with the equation for I 0 (t) from 3, to see if I 0 (t) � 0 for all t > t c . If not, then we have no hope of preventing a second wave as @ θ I 0 (t)<0. Similarly, we should check if θ = 0 will result in I 0 (t)�0 for all t > t c . If so, then we have the opposite case where any choice of θ f will end the infection. In this situation, no optimization on θ f is needed. Instead let us suppose that θ = 0 results in an additional wave of infection and θ = 1 results in no additional wave of infection, then, since @ θ I 0 (t) < 0, we are guaranteed the existence of an optimal θ critical 2 (0, 1) that acts as a separatrix between these two behavioural sets.To find such a θ critical we can use a variant of the bijection method as presented in Algorithm 1. While there are many variant methods that could be used that may be more efficient, the bijection method variant presented in Algorithm 1 is just one example that has the benefit of being simple and is guaranteed to terminate. Moreover, while faster methods can certainly be implemented, Algorithm 1 will approximate θ critical within machine epsilon accuracy in 53 iterations (as after 53 iterations, the interval [θ 0 , θ 1 ] will have width 2 −53 � 1.11 × 10 −16 hence the midpoint is at most distance 2 −54 � 5.55 × 10 −17 from θ critical ).

Input: Parameters t c , γ, α, δ, β DD , β DM , β MM , and function θ(t) defined on [0, t c ] for solving the ODE system (3) and a tolerance �. Output: An approximationŷ 2 ½y critical À ∊; y critical þ ∊� of the value θ critical such that θ > θ critical avoids a second wave of infection and θ < θ critical does not or an error indicating that such a θ critical does not exist.