EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS Vol. 13, No. 3, 2020, 710-729 ISSN 1307-5543 – www.ejpam.com Published by New York Business Global Scenario tree and adaptive decision making on optimal type and timing for intervention and social-economic activity changes to manage the COVID-19 pandemic K. Nah1, S. Chen2, Y. Xiao3, B. Tang1, N. L. Bragazzi1, J. M. Heffernan1,4, A. Asgary5, N. H. Ogden6, J. Wu1,∗ 1 Laboratory for Industrial and Applied Mathematics, Department of Mathematics and Statistics, York University, Toronto, Ontario, Canada 2 Department of Mathematics and Statistics, York University, Toronto, Ontario, Canada 3 Department of Mathematical Sciences, University of Cincinnati, Cincinnati, OH, USA 4 Modelling Infection and Immunity Lab, Centre for Disease Modelling, Department of Mathematics and Statistics, York University, Toronto, Ontario, Canada 5 Disaster & Emergency Management, School of Administrative Studies & Advanced Disaster & Emergency Rapid-response Simulation (ADERSIM), York University, Toronto, Ontario, Canada 6 Public Health Risk Sciences Division, National Microbiology Laboratory , St-Hyacinthe/Guelph Public Health Agency of Canada, Guelph, Canada Abstract. We introduce a novel approach to inform the re-opening plan followed by a post- pandemic lockdown by integrating a stochastic optimization technique with a disease transmission model. We assess Ontarios re-opening plans as a case-study. Taking into account the uncertainties in contact rates during different re-opening phases, we find the optimal timing for the upcoming re-opening phase that maximizes the relaxation of social contacts under uncertainties, while not overwhelming the health system capacity before the arrival of effective therapeutics or vaccines. 2020 Mathematics Subject Classifications: 92D30, 90C15, 37N40 Key Words and Phrases: Scenario tree, COVID-19 social distancing, lockdown exit strategy, re-opening, transmission dynamics model, stochastic optimization 1. Introduction and Background Planning re-opening after several phases of social distancing escalation in order to manage the Coronavirus disease 2019 (COVID-19) pandemic has been challenging for policy- and decision-makers globally. Re-opening too soon could cause an immediate local resurgence; re-opening too late would lead to huge unnecessary economic losses and other ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v13i3.3792 Email address: wujh@yorku.ca (J. Wu) https://www.ejpam.com 710 c© 2020 EJPAM All rights reserved. K. Nah, S. Chen et al. / Eur. J. Pure Appl. Math, 13 (3) (2020), 710-729 711 adverse effects on health. Most governments around the world have adopted a stage-wise re-opening strategy, in which different social-economic activities are permitted in different stages while the epidemic is closely monitored to ensure the needs for intensive care units (ICU) beds and other medical resources do not exceed the health system capacity. COVID-19 transmission dynamic models [18, 20, 22–24] have been used to estimate the peak demand of ICU beds. These models, however, have ignored the inherent stochasticity in nature and human behaviour (i.e., contact rates) that will affect future social distancing relaxation stages, which will ultimately affect the magnitude and timing of the epidemic peak in future outbreak scenarios, including ICU bed needs. This uncertainty poses great challenges to public health decision- and policy-makers in the design, i.e., timing, type of change, that future phases of re-opening of the economy should include. To meet these challenges, here we develop a quantitative stochastic programming model to support a multi-stage decision-making process under substantial uncertainty. Our key idea is to make the best possible decision in re-opening given currently avail- able information that may be weak. Such available weak information is provided by a disease transmission dynamical system that incorporates intrinsic uncertainty in individ- ual behavior, resulting in randomness of the daily average contact rate. Indeed, the exact contact rate for different future stages of re-opening are not known. However, they can be informed by, for example, mathematical modelling studies of historical measures of COVID-19 public health mitigation strategies during different closure stages – bounds of contact rates in the future may be reasonably established in terms of the economic sec- tor and social, economic and health considerations, from which ranges in the timing and magnitude of peak demand of ICU units can be inferred. Furthermore, we can consider different scenarios of the contact rate distributions. Consequently, the distribution of peak demands can be established as well. This provides a scientific basis to calculate the ex- pected loss of contacts, rather than considering only the worst- and/or best-case scenario. A natural approach is to calculate the peak demand based on the average of the future contact rate, however the dynamics of the disease transmission renders dependence of the peak demand on the contact rate highly nonlinear, and it is known that the expectation of a nonlinear function of a random variable is in general different from the nonlinear function of the expectation of the random variable. A critical consideration in re-opening is the recourse option. It is possible that, once re-opened, or when re-opening is initiated, outcomes may be worse than expected, and corrective actions must be taken. A complete re-implementation of lockdown would be the safest solution, but in the long-term it is not acceptable, feasible or sustainable. We, therefore, must consider re-opening plans that allow for some variability in outcomes so that severe recourse actions will not be needed. In this type of multistage decision problem, the first stage decision shall therefore leave enough room for future recourse in unfavourable cases. An aggressive first stage decision might lead to a better result in a favorable case, but could cause failure in an unfavourable one. On the other hand, a pessimistic first stage decision incurs very high social-economic costs. An optimal first stage decision shall take into account the recourse options and the distribution of future uncertainties, and therefore, will yield a proactive and reactive re-opening policy. The K. Nah, S. Chen et al. / Eur. J. Pure Appl. Math, 13 (3) (2020), 710-729 712 policy is not the best in the favorable case, and is not the best in the worst case either, but it is overall optimal in a holistic sense. Here we develop a stochastic programming model based on an adopted disease trans- mission dynamical system [18]. The model formulation and its numerical simulations take full consideration of: (1) the adaptive decision process; (2) the contact rate uncertainties; (3) recourse actions and consequences; and (4) the underlying dynamics of COVID-19 which is described by an ODE system. Stochastic programming is a well established model and is widely used in finance, transportation, agriculture, etc., see [23]. To our best knowledge, this is the first attempt to integrate stochastic optimization and disease trans- mission models for informing re-opening optimal strategies. It is also innovative in that it is the synthesis of an underlying dynamical process described by a disease transmission model involving other public health interventions. In our numerical simulations, we take the province of Ontario, Canada, as a case study. The optimal re-opening plan is achieved by maximizing the total activities from the beginning of the re-opening process to an assumed date when an effective therapeutics or vaccine against COVID-19 is available, subject to the capacity of the health system and the uncertainties of contact rates during the re-opening phases. Using the tools of stochastic programming and the disease transmission model, we find optimal timings for each re-opening phase, that maximizes the relaxation of social contacts under uncertainties while not overwhelming the health system, awaiting the arrival of an effective SARS-CoV-2 therapeutics or vaccine. 2. Materials and Methods 2.1. The Transmission Dynamics Model We use a compartmental model describing COVID-19 transmission dynamics that in- corporates feasible public health measures, including detection of cases by testing and their isolation, and tracing of contacts with cases and their quarantine. In the model formula- tion, the population is divided into susceptible (S), exposed (E), asymptomatic infectious (A), infectious with symptoms (I), and recovered (R) compartments according to the clinical-epidemiological status of individuals. We also include diagnosed and isolated (D), isolated susceptible (Sq), and isolated exposed (Eq) compartments based on control inter- ventions. Within the modelling framework, we take contact tracing into account, where a proportion, q, of individuals exposed to the virus are traced and isolated (referred to as the quarantine proportion). Quarantined individuals can either move to the compartment Eq or Sq, depending on whether transmission has occurred (with probability q), while the other proportion, 1 − q, consists of individuals exposed to the virus who are missed from contact tracing and, therefore, move to the exposed compartment E if infected, or stay in K. Nah, S. Chen et al. / Eur. J. Pure Appl. Math, 13 (3) (2020), 710-729 713 the compartment S otherwise. The transmission dynamics model is S′ = −βc+ cq(1 − β)S(I + θA)/N + λSq, E′ = βc(1 − q)S(I + θA)/N − σE, I ′ = σρE − (δI + α+ γI)I, A′ = σ(1 − ρ)E − γAA, S′q = (1 − β)cqS(I + θA)/N − λSq, E′q = βcqS(I + θA)/N − δqEq, D′ = δII + δqEq − (α+ γD)D, R′ = γII + γAA+ γDD, (1) where N denotes the total population. The model was developed in a series of studies, and parameterized by fitting to the incidence data in a number of settings including the first reported epicenter Wuhan [17, 19, 21], in other provinces in China [20], in South Korea [20], in Italy [14] and in the Province of Ontario, Canada [18, 24]. In what follows, we will use the model fitted to the incidence data in Ontario, however, the model structure and optimization techniques we develop here can be applied to other settings and regions. In order to develop optimal strategies to gradually relax the social distancing subject to key medical resources (hospital wards, ICU beds, for example), we divide compartment D(t) (those diagnosed but not-yet-recovered individuals) as follows: D = Dmild +DICU +Dward. In other words, we distinguish the confirmed cases into those who show only mild symp- toms and do not require hospitalization (Dmild), those who are hospitalized in the ICU with severe symptoms (DICU ), and those who are hospitalized in non-ICU units (Dward). The disease progression and medical treatment process is described as follows: D′mild = (1 − h)(δII + δqEq) − (α+ γmild)Dmild, D′ICU = h(1 − w)(δII + qEq) + bwardDward − (α+ γICU )DICU , D′ward = hw(δII + δqEq) − (α+ γward + bward)Dward. (2) Here, h is the proportion of hospitalized cases among newly confirmed cases, among which a proportion w are hospitalized in non-ICU units directly. In our simulations and analyses, we use values of h and w estimated in [2], namely, h = 0.16 and w = 0.77. The recovery rates for individuals in each state are denoted by mild, ICU and ward. Here, we use values from [10, 15, 16], namely, γmild = 1/5, γICU = 1/14 and γward = 1/12. The value of bward we will use is estimated to be 0.26 × 1 3 which is the product of two quantities: the proportion of patients who are moving to ICU beds among those who are originally admitted to non-ICU units, 0.26 [22], and the average duration of the pre-ICU period of 3 days [22]. Specifically, after March 2, the model was used to determine a time varying contact rate c(t) and quarantine proportion q(t), that were assumed to be exponentially decreasing and increasing with exponential rates r1 and r2, respectively. That is, by introducing T0, T1 and Ts corresponding to March 14 (School closure began), March 18 (Emergency declaration, K. Nah, S. Chen et al. / Eur. J. Pure Appl. Math, 13 (3) (2020), 710-729 714 Figure 1: The flowchart of the transmission dynamics model, where the population is stratified by susceptible (S), exposed (E), asymptomatic infectious (A), symptomatic infectious (I), and recovered status (R), and by quarantine and isolation status (i.e. Sq and Eq). The compartment of diagnosed but not yet resolved cases is further stratified by the need for the use of hospital wards (Dward) and/or ICU beds (DICU ), and mild (Dmild). with closure of public events and recreational venues; Escalation), and March 24 (Closure of all non-essential workplaces), respectively, a full paramterization of the model for each phase of the escalation process in Ontario, Canada was determined (see [18]). Table 1 contains the values of parameters used for the simulation. Parameters estimated from [18] are presented with standard deviations. 2.2. Stochastic Optimization of De-escalation Plans We investigate a social distancing de-escalation plan in which the escalation process experienced in the Province Ontario is reversed through a staged approach. Therefore, four main de-escalation phases will be considered: De-escalation-Phase 0 and 1: Opening of workplaces; De-escalation-Phase 2: Resumption of public events and activities; De-escalation-Phase 3: School opening. We set up the initiation time of de-escalation phase 0 to be May 19, the first day of the (partial) workplace re-opening. We look for the optimal initiation time of de- escalation phase 1, 2 and 3 while keeping the number of diagnosed cases below a certain threshold. This threshold is set by the health system capacity, especially the number of ICU beds in the province. We aim to minimize the reduction in contacts compared to the contacts of the pre-pandemic level until an expected arrival time of a COVID-19 therapeutics or vaccine, assumed, herewithin, to be January 2021. In the final discussion section, we discuss how our approach can be modified to obtain results subject to different assumptions on public health measures and the arrival time of vaccines. We aim to find the optimal de-escalation strategy (the optimal time to switch between contact rates) within the health system capacity, taking into account the various scenarios of social-economic activities at each de-escalation phase. We emphasize that in our study, we assume that the contact rates at each de-escalation phase (1, 2, and 3) are random variables. K. Nah, S. Chen et al. / Eur. J. Pure Appl. Math, 13 (3) (2020), 710-729 715 Table 1: Parameters for COVID-19 transmission dynamics in Ontario, Canada. Parameter Definition Mean (STD) c0 Contact rate before March 14 11.5801 (0.3456) c1 Contact rate between March 14 to March 18 10.1202 (0.9185) c2 Contact rate between March 18 to March 24 8.0495 (0.2787) r1 Exponential decrease of contact rate 0.0466 (0.0152) cb Minimum contact rate after March 24 2.1987 (0.2400) β Probability of transmission per contact 0.1469 (0.0023) q0 Fraction of quarantined exposed individuals before March 24 0.1145 (0.0114) r2 Exponential increase of quarantine fraction 0.1230 (0.0123) qb The maximum quarantine fraction 0.3721 (0.0371) σ Transition rate of exposed individuals to the infected class 1/5 λ Rate at which the quarantined uninfected contacts were 1/14 released into the wider community ρ Probability of having symptoms among infected individuals 0.7036 (0.0261) δI Transition rate of symptomatic infected individuals 0.1344 (0.0134) to the quarantined infected class δq Transition rate of quarantined exposed individuals 0.1237 (0.0086) to the quarantined infected class γI Recovery rate of symptomatic infected individuals 0.1957 (0.0111) γA Recovery rate of asymptomatic infected individuals 0.139 γmild Recovery rate of diagnosed individuals 0.2 with mild symptoms γICU Recovery rate of individuals in ICU units 1/14 γward Recovery rate of individuals in non-ICU units 1/12 α Disease-induced death rate 0.008 θ Modification factor of asymptomatic infectiousness 0.0275 (0.0128) h Proportion of hospitalized cases among newly confirmed cases 0.16 w Proportion of cases admitted to non-ICU units 0.77 among newly hospitalized cases bward Rate of transfer from non-ICU units to ICU units 0.26 × 1 3 Initial values Definition Mean (STD) S(0) Initial susceptible population 1.471107 E(0) Initial exposed population 8.9743 (0.6558) I(0) Initial symptomatic infected population 5.3887 (0.9442) A(0) Initial asymptomatic infected population 19.4186 (3.9406) Sq(0) Initial quarantined susceptible population 0 Eq(0) Initial quarantined exposed population 0 Dmild(0) Initial diagnosed population with mild symptom 4 DICU (0) Initial diagnosed population in ICU units 0 Dward(0) Initial diagnosed population in non-ICU units 1 R(0) Initial recovered population 0 We consider four time-points, t0, t1, t2 and t3, each of which corresponds to the first re-opening date of some workplaces (May 19, 2020), de-escalation-phase 1, de-escalation- phase 2 and de-escalation-phase 3. Each de-escalation phase is characterized by a contact rate, with cr,0, cr,1, cr,2 and cr,3 denoting the contact rates in each de-escalation phase K. Nah, S. Chen et al. / Eur. J. Pure Appl. Math, 13 (3) (2020), 710-729 716 Table 2: Parameters for COVID-19 transmission dynamics post-re-opening in Ontario, Canada. Parameter Definition Value cp Contact rate between April 22 and May 19 4.5191 cr,0 Contact rate during de-escalation phase 0 (since May 19) 5.3736 [cr,1, cr,1] Range of the contact rate during de-escalation phase 1 [6, 9.394] [cr,2, cr,2] Range of the contact rate during de-escalation phase 2 [9.394, 10.5] [cr,3, cr,3] Range of the contact rate during de-escalation phase 3 [10.5, 11.58] qr Quarantine fraction after re-opening (since May 19) 0.3741 δI,r Transition rate of symptomatic infected individuals 0.1586 to the quarantined infected class after re-opening (since May 19) βr probability of transmission per contact (since May 19) 0.13 (0, 1, 2 and 3). Symmetric triangular distributions are defined on the ranges presented in Table 2. We set up the maximum level of contact rate to be the one at the pre- pandemic level and assume that cr,3 = 11.58. The ranges of each triangular distribution are chosen from [25]. They are calculated by modifying weights of four different contact matrices at households, workplaces, schools and communities and others according to strategies of relaxing social distance in different de-escalation phases in Ontario, Canada. In [25], we used age- and location-specific contact matrices to calculate the average total contact rate (per day) during different phases of COVID-19 mitigation in Ontario, Canada. Using the transformations between different populations developed in [3], we calculated the contact matrix for Ontario, Canada based on the POLYMOD 2015s survey data ([11, 13]) for Canada. We then constructed contact matrices for four different location settings: households, workplaces, community-and-others, and schools among 18 age groups. These age groups have different daily schedules during weekdays and weekend, for example, young kids may go to daycares, students go to primary schools, high schools and colleges, working adults and seniors go to work during the weekdays. Assuming the overall contact mixing is a function of the four baseline matrices, base on possible age groups-specific schedules, we calculated different weights for these four baseline matrices from May 19, 2020 for the three de-escalation phases. In order to obtain the contact rate cr,0, the quarantine rate qr and case detection rate δI,r post re-opening (i.e. after May 19), we have first fixed all the other parameter values (all the parameters are constants) to those estimated in the previous study [18]. Note that, we also introduced a new parameter cp describing the contact rate during the pre- re-opening stage (i.e. April 22 to May 19) as the estimation in the previous study [18] only used the data till April 21. Using the least square method, we fitted the model to the data of cumulative reported cases between April 22 to June 9 in Ontario, where the initial conditions (i.e. April 22) are obtained by solving the model in the previous study [18] based on the estimated parameter values. The estimated values for cp, cr,0, qr, δI,r are displayed in Table 2. Further, the transmission probability after re-opening is assumed to be lower than that of the pre-re-opening stage (βr > β) due to increased use of PPE (masks) by the public and personal distancing. The estimate of βr is taken from the result of the similar modeling study in China [20]. K. Nah, S. Chen et al. / Eur. J. Pure Appl. Math, 13 (3) (2020), 710-729 717 We consider the acceptable de-escalation strategies which meet the constraints DICU (t) ≤ DICU for t0 < t < t0 + T, where, t0+T is the time corresponding to the end of cost-evaluation and DICU is a capacity of the ICU beds available for COVID-19 patients. The total number of existing ICU and acute beds in Ontario, Canada as of April 16, 2020 is 3504 and 20,354, respectively [1]. We assume that only a portion of these resources is available for COVID-19 patients. A de-escalation strategy refers to the vectors (ε0, ε1, ε2) determined by εi = ti+1 − ti, for i = 1, 2, 3. Here, ε0 represents the length of de-escalation phase 0 starting from t0 (May 19), 1, 2 and 3 are the duration of the de-escalation phases 1, 2 and 3, respectively. The initiation time of de-escalation phase 3 is at the time t3 = t2 + ε2. Note that ε3 is determined by the equation ε0 + ε1 + ε2 + ε3 = T and we denote a de-escalation strategy corresponding to the scenario (cjr,1, c j r,2, c j r,3) to be (εj0, ε j 1, ε j 2). We consider the situation that re-opening phases are separated by at least two weeks so that the minimal 2 weeks phase switching time is assumed εi ≥ 14, for i = 0, 1, 2. The days of cost-evaluation are set so that the end of the evaluation time corresponds to January 2021. We solve the following scenario based stochastic programming model [4] to minimize the intensity of reduced contacts (cost) during de-escalation phases 0, 1, 2 and 3,∑ j wj ( 3∑ i=0 (cr,3 − cjr,i)ε j i + ug(εj0, ε j 1, ε j 2) ) , where g(εj0, ε j 1, ε j 2) = −log(DICU − max t0