key: cord-0855322-mzi8504f
authors: Lazebnik, T.; Bunimovich-Mendrazitsky, S.
title: The signature features of COVID-19 pandemic in a hybrid mathematical model - implications for optimal work-school lockdown policy
date: 2020-11-04
journal: nan
DOI: 10.1101/2020.11.02.20224584
sha: be708619c5f6de75c0f6ade949462522ce904de3
doc_id: 855322
cord_uid: mzi8504f

Background: The coronavirus disease 2019 (COVID-19) first identified in China, spreads rapidly across the globe and is considered the fastest moving pandemic in history. The new disease has challenged policymakers and scientists on key issues such as the magnitude of the first-time problem, the susceptibility of the population, the severity of the disease, and its symptoms. Most countries have adopted lockdown policies to reduce the spatial spread of COVID-19, but they have damaged the economic and moral fabric of society. Timely action to prevent the spread of the virus is critical, and mathematical modeling in non-pharmaceutical intervention (NPI) policy management has proven to be a major weapon in this fight due to the lack of an effective COVID-19 vaccine. Methods: We present a new hybrid model for COVID-19 dynamics using both an age-structured mathematical model and spatio-temporal model in silico, analyzing the data of COVID-19 in Israel. The age-structured mathematical model is based on SIRD two age-class model. The spatial model examines a circle of day and night (with one-hour resolution) and three main locations (work / school or home) for every individual. Results: We determine mathematically the basic reproduction number ( R_0 ) via the next-generation matrix based on Markov chain theory. Then, we analyze the stability of the equilibria and the effects of the significant differences in infection rates between children and adults. Using the hybrid model, we have introduced a method for estimating the reproduction number of an epidemic in real time from the data of daily notification of cases. The results of the proposed model are confirmed by the Israeli Lockdown experience with a mean square error of 0.205 over two weeks. The model was able to predict changes in ( R_0 ) by opening schools on September 1, 2020, resulting in ( R_0 ) = 2.2, which entailed a month quarantine of all areas of life. According to the model, by extending the school day to 9 hours, and assuming that children and adults go to school and work every day (except weekends), we get a significant reduction in ( R_0 ) of 1.45. Finally, model-based analytical-numerical results are obtained and displayed in graphical profiles. Conclusions: The use of mathematical models promises to reduce the uncertainty in the choice of Lockdown policies. Our unique use of contact details from 2 classes (children and adults), the interaction of populations depending on the time of day (the cycle of day and night), and several physical locations, allowed a new look at the differential dynamics of the spread and control of infection. Using knowledge about how the length of the work and school day affects the dynamics of the spread of the disease can be useful for improving control programs, mitigation, and policy.

Introduction and related work 10 At the beginning of 2020, the novel severe acute respiratory syndrome coronavirus 11 2 (SARS-CoV-2), also known as COVID-19 reached Europe and the western world 12 from China [1] . Similar to other diseases from the coronavirus family, COVID-19 is 13 transmitted human-to-human, but it turned out that COVID-19 is more infectious and 14 transmissible than previous coronavirus [2] . 15 The World Health Organization (WHO) has declared COVID-2019 a public health 16 emergency of international concern [3, 4] . Currently, due to a lack of an efficient vaccine or 17 clinical treatment to COVID-19, policy-makers are forced to rely on non-pharmaceutical 18 intervention (NPI) policies to reduce the infection rate and control the epidemic. A few 19 examples of NPI policies are masks, social distancing, work capsules, and partial to a full 20 lockdown of central locations (restaurants, malls, offices, etc.). These policies reduce the 21 infection rate but have an influence on the economy, healthcare system, and social life. 22 Multiple studies have been conducted to study epidemics in general and the COVID-19 23 epidemic in particular from both biological and epidemiological perspective providing the 24 information about the epidemics behavior with relevant data for later modeling [1, [5] [6] [7] [8] [9] . 25 In addition, these studies examine the clinical characteristics of COVID-19, enabling a 26 deeper understanding of the disease spreads in the population. 27 Mathematical models are shown to be a useful tool for policy-makers to make data-28 driven decisions based on investigation of different scenarios and their outcomes in a 29 controlled manner [5, 10, 11] . These mathematical models can be divided into two main 30 groups. 31 First, models aim to predict the different parameters such as the total COVID-19 32 related death and peak in hospitalized individuals, given the historical data up to some 33 point. For example, Nesteruk [12] used the data from January 16 to February 9 (2020), 34 from mainland China with the continuous SIR (S-susceptible, I-infected, R-recovered) 35 model. Nesteruk [12] fitted the SIR model using the least mean square method, which 36 resulted in poor fitting to later officially confirmed infected cases [4] . This can be 37 explained by the fact that the SIR model with no modification is too simplistic to 38 correctly represent the dynamics of the COVID-19 epidemic. In [12] , a simple calculation 39 method was proposed that allows quick results which can predict the COVID-19 spread 40 with a fine accuracy if calculated using a very accurate data of infections and recoveries. 41 Model proposed by Tuite et al. [11] which used data from January 25 to March 1 42 (2020) Ontario, Canada, is SEIR (E-exposed) model, the extension of the SIR model. 43 Authors of [11] took into consideration four levels of infection severity, social isolation in 44 the exposed and infected states, hospitalization dynamics, and death state. [11] estimated 45 that 56% of the Ontario population would be infected over the course of the epidemic 46 with a peak of 55500 cases in intensive care units (ICU) and 107000 total cases. At the 47 same time, in all of Canada there are only half of this number of infections [4] . This 48 error in their prediction is associated with incomplete information due to lower frequency 49 of tests and inaccurate ICU records. Nevertheless, Tuite et al's [11] model presents a 50 more detailed dynamic between the infected population and the healthcare response 51 which policy-makers can take into consideration. 52 Ivorra et al. [13] proposed nine states model divides individuals into isolated, hospi-53 talized, or dead. In addition, authors of [13] take into consideration the fact that the 54 recorded confirmed data is under-sampled. 55 Furthermore, a machine learning-based model was proposed by Allam et al. [14] using 56 several machine learning algorithms (Knn, Random Forest, Support Vector Machine) on 57 data of 53 clinical cases. This model was used to predict infection severity and spread. 58 However, this approach suffers from an imbalanced sample of the distribution of severity 59 in the population, as simple or asymptomatic cases are not reported, leading to poor 60 representation of the real dynamic. 61 Second, models aim to analyze and optimize NPI policies. A model proposed by 62 Zhao et al. [15] is an extension to the SEIR model where the susceptible population 63 is separated into two groups: individuals not taking infection-prevention actions and 64 people taking infection-prevention actions as an NPI policy. In addition, authors of [15] 65 included a probability of willingness to take infection-prevention actions with changes 66 over time. They predicted 148.5 thousand infections by the end of May 2020 in Wuhan 67 alone while in all of China there were 84.5 thousand infections at the same time [4] . The 68 authors introduces a stochastic element to the SEIR model making it more robust for 69 social changes happened during the epidemic. 70 Di Domenico et al. [16] used data from March 17 to May 11 (2020)Île-de-France with 71 a stochastic age-structured transmission extension of the SEIR model integrating data 72 on age profile and social contacts of four age-based classes. In this model, hospitalization 73 dynamics with ICU cases are taken into consideration in the model. Model shows that 74 during full lockdown the reproductive number is estimated to be 0.68, due to an 81% 75 reduction of the average number of contacts. These results show that dividing the 76 population into several age-based classes better represent the spread of COVID-19 from 77 an epidemiological perspective [17, 18] .

In this paper, we provide and study a more accurate spatio-temporal model for 79 the COVID-19 transmission by using individual two age classes SIRD named hybrid 80 model (D-death) model. We study two important factors concerning the diffusion of 81 COVID-19: schooling/working hours and physical location of infected population has 82 provided new insights into the epidemic dynamics. Based on the different impact of 83 COVID-19 to the immune response, severity of infection and transmission disease in 84 different age groups (mainly children and adults) [17, 19] , we proposed a two classes 85 age-structured SIRD epidemic model dividing the population into children and adults. 86 Moreover, we developed a numerical, stochastic simulator based on this hybrid model 87 (https://teddylazebnik.info/coronavirus-sir-simulation/index.html) for COVID-19 popu-88 lation spread in addition to the analytical examination of the epidemic dynamics.

This paper is organized as follows: First, we introduce our mathematical hybrid 90 model based on the SIRD model with dual age-structured. Second, we present the 91 model's equilibria, stability analysis, and asymptotic form. Third, we present the spatial 92 model extending the SIRD model by introducing a day-night circle and three locations of 93 disease transmission (home, work, school). Fourth, an analysis of NPI policies including 94 optimal lockdown, optimal work-school duration is presented, and a comparison to 95 the Israeli historical data from August and September. Finally, we discuss the main 96 epidemiological results arising from the model.

Hybrid model 98 As shown in [20] and [21] , the spatial model plays an important role in describing the 99 spreading of communicable diseases, because individuals move around inside a zone or 100 habitat in time.

In this section, we present a hybrid model (Fig 1) that is based on a SIRD model 102 for two age classes using eight populations (dynamics are shown in Fig 2) and a spatial 103 model where these populations are distributed in space and time between work, school 104 and home (Fig 5) . We will define these sub-models in the following sections. Two class age-structured epidemic model 106

The SIR model is proven to be a meaningful mathematical tool for epidemic analysis [22] . 107 This model with the needed modifications has already been shown to predict the epidemics 108 such as COVID-19 [23] , influenza [24] , Ebola [10] , and others.

dynamics. An extension of the SIR model to two age-classes has been investigated for 114 explanation of Polio outbreak by Bunimovich-Mendrazitsky and Stone [25] .

Model definition 116 The model considers a constant population with a fixed number of individuals N . Each 117 individual belongs to one of the three groups: susceptible (S), infected (I), and recovered 118 (R) such that N = |S| + |I| + |R|. When an individual in the susceptible group (S) 119 exposed to the infection, it is transferred to the infected group (I). The individual stays 120 in this group on average d I→R days, after which it is transferred to the recovered group 121 (R).

In addition to these groups, we define a death group (D) which is associated with 123 individuals that are not able to fully recover from the disease or succumb to death. 124 Individuals from the infected group (I) recover from the disease in some chance and 125 move to (R), while the others do not and move to (D). Therefore, in each time unit, 126 some rate of infected individuals recover while others die or remain seriously ill. 127 We divide the population into two classes based on their age: children and adults 128 because these groups experience the disease in varying degrees of severity and have 129 different infection rates. In addition, adults and children are present in various discrete 130 locations throughout many hours of the day which affects the spread dynamics. Individ-131 uals below age A are associated with the "children" age-class while individuals in the 132 complementary group are associated with the "adult" age-class. The specific threshold 133 age (A) may differ in different locations but the main goal is to divide the population 134 into two representative age-classes. Since it takes A years from birth to move from a 135 child to an adult age group, the conversion rate is set as α := 1 A .

By expanding the designation to two age-classes, we let S c , I c , R c , D c and S a , I a , R a , D a 137 represent susceptible, infected, recovered, and death groups for children and adults, 138 respectively such that

The model does not take into consideration death during the epidemic unrelated 140 to the disease itself because in the United States in 2018 the birth rate was 11.6 for 141 every 1000 individuals [26] while the mortality rate in 2017 was 8.6 for every 1000 142 individuals [27] , resulting in around 0.3% increment of the population size which is 143 assumed to be small enough to be neglected. In addition, we introduce two death states 144 for children and adults, respectively, that died from the epidemic.

Eq (1-10) describes the epidemic's dynamics.

In Eq (1), dSc(t) dt is the dynamical amount of susceptible individual children over time. 147 It is affected by the following three terms. First, with rate β cc each infected child infects 148 susceptible children. Second, with rate β ca each infected adult infects the susceptible 149 children. Finally, children grow and pass from the children's age-class to the adult's 150 age-class with transition rate α, reduced from the children's age-class. N c is the size 151 of the children population and used to take all variables as fixed proportions of the 152 population N .

In Eq (2), dSa(t) dt is the dynamical amount of susceptible adult individuals over time. children's age-class to the adult's age-class with transition rate α, added to the adult 156 age-class. Second, with rate β ac each infected child infects the susceptible adult. Finally, 157 with rate β aa each infected adult infects a susceptible adult. N a is the size of the adult 158 population and used to take all variables as fixed proportions of the population N .

is the dynamical amount of infected individual children over time. 160 It is affected by the following four terms. First, with rate β ca each infected child infects 161 the susceptible adult. Second, with rate β cc each infected child infects a susceptible child. 162 Third, individuals recover or die from the disease after a period γ c . Finally, children 163 grow and pass from the children's age-class to the adult's age-class with transition rate 164 α, reduced from the adult's age-class. While the last process has a minor impact relative 165 to the first three processes, we do not neglect it in order to count edge-cases.

November 2, 2020 6/28 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted November 4, 2020. ; https://doi.org/10.1101/2020.11.02.20224584 doi: medRxiv preprint

In Eq (4), dIa(t) dt is the dynamical amount of infected adult individuals over time. It 167 is affected by the following four terms. First, with rate β ac each infected child infects the 168 susceptible adult. Second, with rate β aa each infected adult infects a susceptible adult. 169 Third, individuals recover or die from the disease after period γ a . Finally, children grow 170 and pass from the children's age-class to the adult's age-class with transition rate α, 171 added to the adult's age-class. While the last process has a minor impact relative to the 172 first three processes, we do not neglect it to count edge-cases.

In Eq (5), dRc(t) dt is the dynamical amount of recovered individual children over 174 time. It is affected by the following two terms. First, in each point, a portion of the 175 infected children recover after period γ c which is multiplied by the rate of children 176 that do recover from the disease ρ c . Second, children grow from birth and pass from 177 the children's age-class to the adult age-class with transition rate α, reduced from the 178 children's age-class.

In Eq (6), dRa(t) dt is the dynamical amount of recovered adult individuals over time. 180 It is affected by the following two terms. First, in each point, a portion of the infected 181 adults recover after period γ a which is multiplied by the rate of adults that do recover 182 from the disease ρ a . Second, children grow from birth and pass from the children's 183 age-class to the adult age-class with transition rate α, added to the adult age-class.

In Eq (7), dDc(t) dt is the dynamical amount of dead individual children over time. It 185 is affected by the portion of the infected children that do not recover after period γ c 186 which is multiplied by the rate of children that do not recover from the disease ψ c .

is the dynamical amount of dead adult individuals over time. It is 188 affected by a portion of the infected adult that do not recover after period γ a which is 189 multiplied by the rate of adults that do not recover from the disease ψ a .

It should be noted that both Eq (7) and (8) can be presented as one equation

as D a and D c are accumulative populations which do not infect the dynamics separately. 192 Nevertheless, to provide a better representation of the age-based death in the epidemic, 193 we chose to divide the dead individuals into the same age-groups which allows later 194 analysis of the death of children and adults separately. In summary, the interactions 195 between disease stages presented in Fig 1 are modelled by the following system of eight 196 coupled ordinary differential equations (ODE) Eq (9) with initial conditions at t = 0, 197 Eq (10) as :

November 2, 2020 7/28 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted November 4, 2020. ;

.

The parameters used in the calculation of the model throughout the paper are 199 presented in Table. 1. The values are cited from the sources themselves except A and β ca 200 calculated from the data of the correlated source [6, 17] . The threshold of the children's 201 age to become adults in parameter A is set to 13 year as the mean value of the group 202 of ages in which the percentage of critical cases relative to all cases is the highest as 203 reported by Dong et al. [6] ( Table 2 ). The susceptible contacts in children who become 204 infected due to direct disease transmission from an adult β ca are calculated based on 205 data of 10 children. Eight children out of the 10 have been exposed to 30 adults in total 206 and later found to be infected, as reported by Cai et al. [17] . Therefore, the infected 207 rate is set to 8 30 = 0.266.

Numerical solution

To obtain a better understanding of how different parameter values influence the system 210 dynamics, we illustrate the behavior of the system using numerical analysis.

and D a (t)). We processed eight-order Runge-Kutta integration on 214 the differential system to enable numerical simulations. Runge-Kutta integration of 215 the equations was implemented by Octave programming (version 5.2.0), using standard 216 program lsode, for a set of initial conditions described in Eq (10).

The solutions of the system (9-10) are shown in Fig 3. In addition, the population 218 size, adults, and children are set to be N = 8000000, N c = 2240000, N a = 5760000 to 219 present the distribution of the Israeli population in 2017 as published by the Israeli 220 November 2, 2020 8/28

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted November 4, 2020. ; 

where the x-axis in all eight graphs represents the 222 time (in days) that has passed from the beginning of the dynamics and the y-axis is the 223 size of each population, respectively. A maximum in the percent of infected children and 224 adults (39%, 89%) is reached in the 8th and 11th days as shown in Fig 3. Therefore, the 225 maximum infected population was 89% of the whole population. Besides, all children 226 infected and recovered after 17 days while all adults recovered after 82 days.

In epidemiology, it is essential to quantify the severity of outbreaks of infectious diseases. 229 The standard parameter indicates the severity called the basic reproduction number R 0 . 230 In the standard SIR model [22] , R 0 is defined to be the ratio between individual infection 231 rate and recovery rate, where almost everyone is susceptible (namely, S c + S a ∼ N ). For 232 any start condition and parameters, at t → ∞ the model's state is an instance of Eq 233 (10) according to Lemma 1.

Definition 1. The model's state at time t is defined by

where g is a function g :

November 2, 2020 9/28 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted November 4, 2020. ; Table. 1. The graphs show the evolution in time (days) of S c (t), S a (t), I c (t), I a (t), R c (t), R a (t), D c (t), and D a (t). Adult and children graphs are presented with a dotted and solid lines, respectively. Susceptible, infected, recovered, and dead are shown in green, red, blue, and black, respectively. November 2, 2020 10/28

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted November 4, 2020. ; https://doi.org/10.1101/2020.11.02.20224584 doi: medRxiv preprint Proof. Define a non-homogeneous continuous-time Markov chain M with state space 240 with eight values in respect to S c , S a , I c , I a , R c , R a , D c , and D a . Table 2 represents 241 the transition matrix of M where

In order to model the dynamics as a Markov chain one needs to show that

where {X} t is a stochastic process with values in the state space for all t ≥ 0 and all 244 states i 0 , ...i t , j and h is an arbitrary small step in time [31] . M satisfies Eq (12) 

such that (= * ) stands for I c (t+h) = I c (t)+dI c (t), I a (t+h) = I a (t)+dI a (t). From Eq (13) 247 and (14), v 1 (t + h) and v 2 (t + h) depends only on the previous state X t = i t , respectively. 248 Let π ∈ N 8 be the initial condition for the system. Therefore, the model's state at some 249 time t is defined by A t π. For any initial condition π, the vector lim t→∞ (A t π) takes the 250 form 

and for lim t→∞ (1 − α) t = 0 meaning that asymptotically the population in S c decreases 258 to zero. Therefore, the asymptotic state of the model takes the form:

, for any initial condition π and parameters such that α, γ a , γ c , β aa , β cc > 0.

In order to find g (from Lemma 1) for given initial conditions and parameters, it is 262 possible to use the next generation matrix approach [32] . Let V + i be the rate of transfer 263 of individuals by all means except of appearance of new infections in compartment i, 264 and V − i be the rate of transfer of individuals out of compartment i. Let F i be the rate 265 of appearance of new infections in compartment i. Eq (1 -8) takes the form:

November 2, 2020 11/28

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted November 4, 2020. ; https://doi.org/10.1101/2020.11.02.20224584 doi: medRxiv preprint Table 2 . Markov chain transformation matrix representation of Eq (1) (2) (3) (4) (5) (6) (7) (8) .

In equilibrium, I c and I a are both equal to zero so the derivatives at equilibrium, 270 focusing on I c and I a from Eqs. (3) and (4), are mapped to the third and forth elements 271 in vectors F and V giving the matrices F and V.

The next generation matrix is defined as FV −1 [32] . 

After one time unit, the amount of infected individuals can be 275 calculated using v = Gv. Therefore, for any start condition and model parameters one 276 can find g by calculating

where C ∈ N is the first index that satisfies I c (C) + I a (C) = 0. 

Both γ c and γ a are biological properties of the disease; on the other hand, β aa , β ac , β ca , 286 and β cc can be managed using social distance, quarantine, masks, and other meth-287 ods. Figs 4a -4e presents five projections of stability from the parameter space 288 November 2, 2020 12/28

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted November 4, 2020. ; https://doi.org/10.1101/2020.11.02.20224584 doi: medRxiv preprint {β aa , β ac , β ca , β cc } calculated using Eq (16) . Fig 4a shows the β aa − β ac projection 289 where β cc = 0.308, β ca = 0.266. There are no values such that R 0 < 1 and from the 290 color gradient, it is possible to see that β ac has slightly more influence on R 0 relatively 291 to β aa . Fig 4b shows the β aa − β ca projection where β cc = 0.308, β ac = 0. There are no 292 values such that R 0 < 1 and from the color gradient, it is possible to see that β aa has 293 a minor to no influence on R 0 relatively to β ca . Fig 4c shows the β cc − β ca projection 294 where β aa = 0.308, β ca = 0.266. R 0 < 1 for any combination of β cc , β ac such that 295 β ac < 0.08 − 1.6β cc . Fig 4d shows the β cc − β ac projection where β aa = 0.308, β ac = 0. 296 R 0 < 1 for any combination of β cc , β ac such that β cc ≤ 0.3. Fig 4e shows the β aa − β cc 297 projection where β ca = 0.266, β ac = 0. R 0 < 1 for any combination of β aa , β cc such that 298 β cc < 0.9 − 0.155β aa .

In order to obtain the equilibria points of the model and their stability properties, 300 the Jacobian (J) is obtained by linearizing Eqs (3) and (4) in the system (Eqs (1 -8) ). 301 The calculations show that 

The system (Eqs (9-10)) has four non-trivial equilibria which may be found by setting 303 all rates in Eq (9) to zero; marked by EQ 1 , EQ 2 , EQ 3 , and EQ 4 . Thus equilibria are 304 provided as the state vector of the eight population states in Table 3 , where

There are other equilibria for trivial cases such as S a (0) = 0 or S c (0) = 0 which are 306 unrealistic for real-world dynamics. Below we investigate the stability of four nontrivial 307 equilibria. 308 EQ 1 : is the case for I c (t) = I a (t) = 0 for every t ∈ N (see Table 3 ). This is a trivial 309 equilibrium where there is no epidemic at all. Because the model does not take into 310 consideration the birth of new children, after A days all children become adults and EQ 1 311 is derived. By settings J(EQ 1 ), all eigenvalues are negative for all parameter values 312 except λ = −γ a + β aa N a+N c N a . Therefore, this equilibrium is stable if β aa N a+N c N a < γ a .

EQ 2 : is the case for I a (t) = 0 for every t ∈ N (see Table 3 ). This equilibrium can be 314 achieved if β ac = 0 because in any other case (assuming S a (0) > 0) exists time t such 315 that dIa(t) dt > 0 and therefore I a (t + 1) = 0. 316 EQ 3 : is the case for I c (t) = 0 for every t ∈ N (see Table 3 ). This equilibrium can be 317 achieved if β ca = 0 because in any other case (assuming S c (0) > 0) exists time t such 318 that dIc(t) dt > 0 and therefore I c (t + 1) = 0. 319 EQ 4 : is the case for ∃t such that I c (t) = 0 and I a (t) = 0 (see Table 3 ). This 320 equilibrium is the generic case for the model.

According to Lemma 1, EQ 2 , EQ 3 , and EQ 4 are not possible asymptotic forms of 322 the system and are therefore not stable. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted November 4, 2020. ; https://doi.org/10.1101/2020.11.02.20224584 doi: medRxiv preprint . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted November 4, 2020. ; Table 3 . The four equilibria solutions for Eqs (9) (10) . 5. t d a hours of the day that adults are at home and t n a = 24 − t d a hours of the day that 334 adults are at work. 335 We assume the transition from home to either work or school and back is immediate and 336 that everybody is following the same clock. Each simulation step simulates one hour. 337 The population size is constant during the simulation and initialized in the beginning of 338 each iteration by setting children population size N c and adult population size N a . In 339 each simulation step, the following three actions take place: 340 

group there is a change of β aa , β ac , β ca , or β cc according to the age-class of the 342 two members that the first will be infected. 343 2. Each infected child or adult that was infected for 1 γc , 1 γa simulation steps becomes 344 either recovered or deceased, respectively. 345 3. According to the hour of the day, the adults transition to home or work and the 346 children to home or school.

The spatial model adds day-night circle and three main locations to the dynamics of the 348 hybrid model.

November 2, 2020 15/28 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted November 4, 2020. ; Table. 1. Children and adult graphs are presented with a dotted and solid lines, respectively. Susceptible, infected, recovered, and dead are shown in green, red, blue, and black, respectively. 

A maximum in the percent of infected children and adults (38.5%, 100%) is reached 355 on the 9 th and 12 th day, as shown in Fig 6. In addition, all children were infected and 356 recovered after 12 days while all adults recovered after 28 days. The simulation shows 357 similar results to the results derived by solving the model (Eqs (9-10)) as presented in 358 Fig 3. The main difference between the two is that the simulation predicts 4.43 times 359 shorter duration from the time of maximum infected adults to the time that all adults 360 are either dead or recovered, as the infected adult population equals zero on the 82 nd 361 and 28 th days. In addition, the maximum of infected adults is reached on the 11 th and 362 12 th days resulting in 71 and 16 days for full recovery. One of the major hopes of politicians is for a NPI policy in which the epidemic does 368 not reach an outbreak at any time after they initialized a given decision. We define this 369 November 2, 2020 17/28 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted November 4, 2020. ; https://doi.org/10.1101/2020.11.02.20224584 doi: medRxiv preprint condition mathematically as follows: 370 Definition 2. For given initial condition and model's parameter (P). A solution of 371 Eq (9-10) is defined as outbreak dynamics if ∃t ∈ N : R 0 (t) > 1. 372 We will examine two policies, based on this condition, to determine if each one is 373 possible to fulfill the condition. If so, the optimal NPI policy is based on the parameter-374 space criteria: 375 1. The influence of the duration of the work/school day. 376 2. Lockdown in homes with partial to full separation between individuals.

Duration of working and schooling day 378 We will assume that children meet only children in school and adults meet only adults at 379 work and at home adults and children meet each other. This is a good approximation of 380 the real dynamics as a relatively small percent of adults work with children during the 381 day which keeps the interactions between children and adults at this time relatively small 382 to the extent of interaction when both adults and children are at home and therefore 383 can be neglected.

Based on the proposed spatial model, we change the number of hours children and 385 adults spend in school and work each day, respectively. Fig 7a shows the average R 0 as 386 a function of the duration in hours of the working (d W ) and schooling (d S ) day. The 387 dots are the calculated values from the simulator and the surface is a fitting function 388 R 0 (d W , d S ). The fitting function is calculated using the least mean square (LMS) 389 method [33] . In order to use the LMS method, one needs to define the family function 390 approximating a function. The family function 391 f (a, c) = p 1 + p 2 a + p 3 c + p 4 ac + p 5 a 2 + p 6 c 2 has been chosen to balance between the accuracy of the sampled data on the one hand 392 and simplicity of usage on the other [34] ,

and was obtained with a coefficient of determination 

obtained with a coefficient of determination R 2 = 0.748.

The duration of either working or schooling day has a local minimum at (d W , d S ) = 399 (18, 14) with R 0 = 0.808, as shown in Fig 7a. The optimal point from Eq (17) is 400 (d W , d S ) = (9, 17.5) with R 0 = 0.967 obtained by setting the gradient of R 0 (d W , d S ) 401 to zero. The inconsistency in the values is resulted by the error in the fitting function 402 but both show the same behavior in which longer working hours reduce significantly 403 the infection rate but too many hours increase the infection rate back as the lack of 404 circulation of contagious adults and children during the day helps to reduce the infection 405 rate. Furthermore , Fig 7b presents a binary classification where red cells are cases 406 with an outbreak and green cells are cases without outbreak during the simulation as a 407 function of the working (d W ) and schooling (d S ) duration in hours each day.

November 2, 2020 18/28 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted November 4, 2020. ; https://doi.org/10.1101/2020.11.02.20224584 doi: medRxiv preprint (a) Average R0. (18) . There is a sharp decline between a shorter 412 working-schooling duration of 10 hours or less and longer than that. This is associated 413 with the fact that with more than 12 hours (half-day) the dynamics that have a higher 414 incidence are the ones with β ca = β ac = 0 which reduces the number of infected 415 individuals in total. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted November 4, 2020. ; https://doi.org/10.1101/2020.11.02.20224584 doi: medRxiv preprint and mental) in the population. Therefore, the optimization task of finding the minimal 423 portion of the population to be locked down such that the epidemic will be constrained 424 is important.

The optimization problem can be written formally as 426 min La,Lc

where L a and L c are the portion of adults and children in lockdown. We ran the spatial 427 model multiple times where each time a combination of (L a , L c ) = {(0.1i, 0.1i)} 10 i=0 of 428 the population is assumed to be in lockdown. Fig 8a shows the results of this calculation. 429 Each dot represents R 0 of each case L a , L c . The black grid shows the threshold R 0 = 1. 430 The behavior of R 0 as a function of L a , L c has been retrieved using the LMS method 431 and takes the form

obtained with a coefficient of determination R 2 = 0.971. Therefore, it is safe to 433 claim that function R 0 (a, c) is well fitting the data despite the stochastic noise of the 434 simulation and presents a fair approximation for the R 0 behavior as a function of L a , L c . 435 Using Eq (19) it is possible to find the constraints of L a , L c such that R 0 ≤ 1, by solving 436 R 0 (L a , L c ) ≤ 1.

Closing only schools without any reduction in adults going to work does not prevent 438 an epidemic outbreak while locking down half of the adult population will prevent an 439 outbreak. Lockdown of children (L c ) have a minor effect relatively to lockdown adults 440 (L a ) as shown in both Eq (19) and 

and obtained with a coefficient of determination R 2 = 0.840.

In addition, we repeat the analysis by numerically solving the system (Eq (9-10)). 445 The lockdown is reflected in the model as the infection rate between every two individuals. 446 L a affects the interactions between adults and either other adults or children β aa , β ac 447 and L c affects the interactions between children and either other children or adults 448 β cc , β ca . Mark the original values of β aa , β ac , β ca , β cc from Table 1 as β * aa , β * ac , β * ca , β * cc . 449 Each time, the system solved where (L a , L c ) = {(0.1i, 0.1i)} 10 i=0 such that β aa = β * aa · L a , 450 β ac = β * ac · L a , β ca = β * ca · L c , and β cc = β * cc · L c . The behavior of R 0 as a function of L a , L c has been derived using the LMS method 461 and takes the form

November 2, 2020 20/28 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted November 4, 2020. ; https://doi.org/10.1101/2020.11.02.20224584 doi: medRxiv preprint (a) Average R0. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted November 4, 2020. ; From [4] there was 90232 infected at August 15. Given that the average recover rate 482 of children is 2 days and adults is 14 days (see Table 1 ) and that children is 28% of 483 the population, we obtain that there was approximately 18899 infected adults and 868 484 infected children. We assume that the amount of recovered individuals is heterogeneous 485 in the population and that only adults die from the epidemic until this point. Eq (23) 486 shows the initial conditions for August 15 (t 0 ), Israel. Considering that, longer school day is reducing the infection rate (Fig 7a) , we solve 497 numerically the hybrid model and calculate R 0 for the case when the schools are open 498 (from September 1). We assume that children go to school every day (except weekends) 499 for nine hours while the adults go to work for the same period. We obtain that the 500 difference average R 0 of the case with the optimal NPI (R 0 = 1.45) and the historical 501 NPI (R 0 = 2.28) is ∆R 0 = 0.83 (see Fig 11) .

This study presents a model showing the effects of population age, time of day, and 504 gathering location on the spread of the epidemic. Based on the proposed model Eqs (9-505 10), the non-trivial equilibria of the system is presented in Table 3 . Based on Lemma 1, 506 the equilibria do not take the form of the asymptotic case (Eq (11)) of the model and 507 therefore are not stable. Although the proposed model is a simplification of the total 508 complexity of the COVID-19 epidemic, the asymptotic solution without any intervention 509 will result in a high percentage 39% children, 89% adults and 38.5% children, 100% 510 adults of infected individuals as shown in Figs 3 and 6, respectively.

In addition, the asymptotic solutions of the model are insignificant for themselves as 512 November 2, 2020 24/28 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted November 4, 2020. ; https://doi.org/10.1101/2020.11.02.20224584 doi: medRxiv preprint in practice these results occur after A days (13 years). The main epidemic dynamics as 513 obtained from both the SIRD model Eqs (9) (10) and the hybrid model is a few month 514 long, as shown in Figs 3 and 6. Besides, the asymptotic solution can be used to easily 515 obtain D a and D c as they accumulative sums of infected infected adults and children 516 multiplied by a factor ψ a and ψ c and do not change after the epidemic is over. As a 517 result, α = 2.1 · 10 −4 is not used as for 82 days (as shown in Fig 3) only 0.017% of the 518 children will become adults which effectually has only a minor affect on the epidemic 519 results. 520 We presented the proposed model as a stochastic process using the Markov chain 521 modeling approach as shown in Lemma 1. Using this representation, it is easy to find 522 the stationary state of the chain (Eq (11)) and therefore the possible asymptotic forms 523 of the SIRD model Eqs (9-10).

Furthermore, two policies for controlling the epidemic and preventing outbreak have 525 been investigated using the hybrid model. First, the duration of working and schooling 526 day has an influence (up to 10%) on the maximum of infected individuals but can prevent 527 outbreaks under the assumption that during this time adults do not contact children at 528 all as shown in Fig 7. These results match the conclusions reported by Keskinocak et 529 al. (2020) for Georgia state [35] . In addition, these results match the results reported 530 by Di Domenico [16] that school closure has a minor influence on the epidemic peak, 531 just delays it. Second, a partial lockdown of adults and children show that adults and 532 children lockdown has a similar first-order effect (linear coefficients of L a and L c ) but 533 the combination between the two (L a L c ) has a bigger weight on the average infection 534 rate as shown in Eqs (19) and (21) . These results match the conclusions reported by 535 Aglar et al. (2020) for the state of Georgia regarding the voluntary lockdown with 536 school closure [36] . By the same token, a lockdown of approximately half of the adult 537 population prevents an outbreak where children lockdown has a less significant influence 538 as shown in Figs 8 and 9. These results are based on values from Table 1 , which can vary 539 significantly between countries and depend on other hyperparameters such as population 540 density, age distribution of the population and others.

The model explains the dynamics that took place between August 15 and September 542 15 (2020) in Israel as presented in Fig 11 with mean square error (MSE) of 0.205. The 543 opening of the schools is equal to reducing the lockdown for children with some relaxation 544 (smaller L a ) of the lockdown for adults (because if the children in school they can go to 545 work) with less then half (50%) of the adult population voluntary in lockdown. Also, 546 the schooling hours were reduced to less than 6 hours each day. From Figs 7a and 8c 547 it is easy to see that this policy predicts high increase in the infection rate, as indeed 548 happened. Keeping the schools open while keeping the increase in the infection rate 549 from increasing significantly is possible if the schooling hours are longer (8-9 hours) as 550 shown in Fig 7a. The influence of this policy in Israel during the school opening which 551 take place in September 1 show that the R 0 can be reduced by 0.83 in comparison to 552 a policy in which children go to school every other day for five hours, as shown in Fig 553  11 . Also, if at least half of the adult population will be in lockdown the influence of the 554 schools on the infection rate will be relatively small as shown in Fig 8b. 

In the case of a future pandemic virus, researchers will be able to use our approach 556 to predict the exact consequences of choosing "Lockdown" strategies, especially those 557 needed for pandemics with a lack of immunity in the world's population. Understanding 558 the age-based dynamics in COVID-19 spread in the population will be crucial for the 559 development of an optimal NPI policy, as well as for optimizing current polices and 560 analyzing clinical data. The further extensions of the model will be used to learn how to 561 manage international travel between countries with a range of healthcare system abilities 562 in both the context of the COVID-19 epidemic and for other epidemics scenarios.

November 2, 2020 25/28 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted November 4, 2020. ; https://doi.org/10.1101/2020.11.02.20224584 doi: medRxiv preprint

Clinical progression of patients with COVID-19 in

A familial cluster of pneumonia associated with the 2019 novel coronavirus indicating person-toperson transmission: a study of a family cluster

Note from the editors: World Health Organization declares novel coronavirus (2019-nCoV) sixth public health emergency of international concern

Epidemic spread in human networks

Epidemiological Characteristics of 2143 Pediatric Patients With 2019 Coronavirus Disease in China

Clinical and epidemiological characteristics of 1420 European patients with mild-to-moderate coronavirus disease 2019

Effects of age and sex on recovery from COVID-19: Analysis of 5769 Israeli patients

Early antiviral treatment contributes to alleviate the severity and improve the prognosis of patients with novel coronavirus disease (COVID-19)

Mathematical models of SIR disease spread with combined non-sexual and sexual transmission routes. Infectious Disease Modelling

Mathematical modelling of COVID-19 transmission and mitigation strategies in the population of Ontario

Statistics-based Predictions of Coronavirus Epidemic Spreading in Mainland China

Mathematical modeling of the spread of the coronavirus disease 2019 (COVID-19) taking into account the undetected infections. The case of China. Communications in nonlinear science and numerical simulation

Artificial Intelligence (AI) Provided Early Detection of the Coronavirus (COVID-19) in China and Will Influence Future Urban Health Policy Internationally

Imitation dynamics in the mitigation of the novel coronavirus disease (COVID-19) outbreak in Wuhan

Impact of lockdown on COVID-19 epidemic in Ile-de-France and possible exit strategies

A Case Series of Children With 2019 Novel Coronavirus Infection: Clinical and Epidemiological Features

Estimation of the asymptomatic ratio of novel coronavirus infections (COVID-19)

COVID-19 epidemic: Disease characteristics in children

Simulating the spread of COVID-19 via a spatially-resolved susceptible-exposed-infected-recovered-deceased (SEIRD) model with heterogeneous diffusion

System inference for the spatio-temporal evolution of infectious diseases: Michigan in the time of COVID-19

A contribution to the mathematical theory of epidemics

Rational evaluation of various epidemic models based on the COVID-19 data of China. medRxiv

Modeling seasonal influenze in Israel

Births: Final Data for 2018

Deaths: Final Data for 2017

Cardiovascular Disease, Drug Therapy, and Mortality in Covid-19

Characteristics of and important lessons from the coronavirus disease 2019 (COVID-19) outbreak in China: summary of a report of 72314 cases from the Chinese Center for Disease Control and Prevention

COVID-19 in children: the link in the transmission chain

Non-homogeneous Markov Chains

Reproduction numbers and subthreshold endemic equilibria for compartmental models of disease transmission

Numerical Methods for Least Squares Problems

Polynomial Regression with Response Surface Analysis: A Powerful Approach for Examining Moderation and Overcoming Limitations of Difference Scores

Evaluating Scenarios for School Reopening under COVID19. medRxiv

Homebound by COVID19: The Benets and Consequences of Non-pharmaceutical Intervention Strategies

COVID-19 lockdown: if, when and how