The Signature Features of COVID-19 Pandemic in a Hybrid Mathematical Model-Implications for Optimal Work-School Lockdown Policy.

The new COVID-19 pandemic has challenged policymakers on key issues. 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. Mathematical modeling in non-pharmaceutical intervention policy management has proven to be a major weapon in this fight due to the lack of an effective COVID-19 vaccine. A new hybrid model for COVID-19 dynamics using both an age-structured mathematical model based on the SIRD model and spatio-temporal model in silico is presented, analyzing the data of COVID-19 in Israel. Using the hybrid model, a method for estimating the reproduction number of an epidemic in real-time from the data of daily notification of cases is introduced. The results of the proposed model are confirmed by the Israeli Lockdown experience with a mean square error of 0.205 over 2 weeks. The use of mathematical models promises to reduce the uncertainty in the choice of "Lockdown" policies. The unique use of contact details from 2 classes (children and adults), the interaction of populations depending on the time of day, and several physical locations, allow a new look at the differential dynamics of the spread and control of infection.

Introduction and Related Work

At the beginning of 2020, the novel severe acute respiratory syndrome coronavirus 2 (SARS‐CoV‐2), also known as COVID‐19, reached Europe and the western world from China.[ ] Similar to other diseases from the coronavirus family, COVID‐19 is transmitted human‐to‐human, but it turned out that COVID‐19 is more infectious and transmissible than previous coronavirus.[ ]

The World Health Organization (WHO) has declared COVID‐2019 a public health emergency of international concern.[ , ] Currently, due to a lack of an efficient vaccine or clinical treatment to COVID‐19, policy‐makers are forced to rely on non‐pharmaceutical intervention (NPI) policies to reduce the infection rate and control the epidemic. A few examples of NPI policies are masks, social distancing, work capsules, and partial to a full lockdown of central locations (restaurants, malls, offices, etc.). These policies reduce the infection rate but have an influence on the economy, healthcare system, and social life.

Multiple studies have been conducted to study epidemics in general and the COVID‐19 epidemic in particular from both biological and epidemiological perspective providing the information about the epidemics behavior with relevant data for later modeling.[ , , , , , ] In addition, these studies examine the clinical characteristics of COVID‐19, enabling a deeper understanding of the disease spreads in the population.

Mathematical models are shown to be a useful tool for policy‐makers to make data‐driven decisions based on investigation of different scenarios and their outcomes in a controlled manner.[ , , ] These mathematical models can be divided into two main groups.

First, models aim to predict the different parameters such as the total COVID‐19 related death and peak in hospitalized individuals, given the historical data up to some point. For example, Nesteruk[ ] used the data from January 16 to February 9 (2020), from mainland China with the continuous SIR (S‐susceptible, I‐infected, R‐recovered) model. Nesteruk[ ] fitted the SIR model using the least mean square method, which resulted in poor fitting to later officially confirmed infected cases.[ ] This can be explained by the fact that the SIR model with no modification is too simplistic to correctly represent the dynamics of the COVID‐19 epidemic. In ref. [12], a simple calculation method was proposed that allows quick results which can predict the COVID‐19 spread with a fine accuracy if calculated using a very accurate data of infections and recoveries.

Model proposed by Tuite et al.[ ] which used data from January 25 to March 1 (2020) Ontario, Canada, is SEIR (E‐exposed) model, the extension of the SIR model. Authors of ref. [11] took into consideration four levels of infection severity, social isolation in the exposed and infected states, hospitalization dynamics, and death state.[ ] Estimated that 56% of the Ontario population would be infected over the course of the epidemic with a peak of 55 500 cases in intensive care units (ICU) and 107 000 total cases. At the same time, in all of Canada there are only half of this number of infections.[ ] This error in their prediction is associated with incomplete information due to lower frequency of tests and inaccurate ICU records. Nevertheless, Tuite et al.'s[ ] model presents a more detailed dynamic between the infected population and the healthcare response which policy‐makers can take into consideration.

Ivorra et al.[ ] proposed nine states model divides individuals into isolated, hospitalized, or dead. In addition, authors of ref. [13] take into consideration the fact that the recorded confirmed data is under‐sampled.

Furthermore, a machine learning‐based model was proposed by Allam et al.[ ] using several machine learning algorithms (Knn, Random Forest, Support Vector Machine) on data of 53 clinical cases. This model was used to predict infection severity and spread. However, this approach suffers from an imbalanced sample of the distribution of severity in the population, as simple or asymptomatic cases are not reported, leading to poor representation of the real dynamic.

Second, models aim to analyze and optimize NPI policies. A model proposed by Zhao et al.[ ] is an extension to the SEIR model where the susceptible population is separated into two groups: individuals not taking infection‐prevention actions and people taking infection‐prevention actions as an NPI policy. In addition, authors of ref. [15] included a probability of willingness to take infection‐prevention actions with changes over time. They predicted 148.5 thousand infections by the end of May 2020 in Wuhan alone, while in all of China there were 84.5 thousand infections at the same time.[ ] The authors introduces a stochastic element to the SEIR model making it more robust for social changes happened during the epidemic.

Di Domenico et al.[ ] used data from March 17 to May 11 (2020) Île‐de‐France with a stochastic age‐structured transmission extension of the SEIR model integrating data on age profile and social contacts of four age‐based classes. In this model, hospitalization dynamics with ICU cases are taken into consideration in the model. Model shows that during full lockdown the reproductive number is estimated to be 0.68, due to an 81% reduction of the average number of contacts. These results show that dividing the population into several age‐based classes better represent the spread of COVID‐19 from an epidemiological perspective.[ , ]

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

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

Hybrid Model

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

In this section, we present a hybrid model (Figure  ) that is based on an SIRD model for two age classes using eight populations (dynamics are shown in Figure  ) and a spatial model where these populations are distributed in space and time between work, school, and home (Figure 5). We will define these sub‐models in the following sections.

Figure 1

Schematic view of the hybrid model's components and relationship between the ODE and the spatial models.

Figure 2

The Hybrid model's schematic view as a transition between disease stages, divided by age‐class.

Figure 5

The panel of the spatial model. From top to bottom: the time from the beginning of the simulation. Distribution of the population to susceptible, infected, recovered, and dead groups. The distribution of the population to susceptible, infected, recovered, and dead groups with separation to children and adults and their current location (home, work, school). The R 0 at a certain time and the average R 0 from the beginning of the simulation.

Table 1

The Hybrid model's parameter description, values, and sources

Parameter definition	Symbol	Value	Source
Children COVID‐19 threshold age in days [t]	A	4745 (13 years)	[6]
Children to adult each day transition rate [t^−1]	α:=1A	2.1 ·10^−4	[6]
Infected to recover average duration for children in days [t^−1]	γ_c	0.5	[17]
Infected to recover average duration for adults in days [t^−1]	γ_a	0.0714	[8]
Susceptible contacts in children which become infected due to direct disease transmission from an adult in a day [t^−1]	β_ca	0.266	[18]
Susceptible contacts in adults which become infected due to direct disease transmission from children in a day [t^−1]	β_ac	∼0	[20]
Susceptible contacts in children which become infected due to direct disease transmission from children in a day [t^−1]	β_cc	0.308	[19]
Susceptible contacts in adults which become infected due to direct disease transmission from an adult in a day [t^−1]	β_aa	0.308	[19]
The probability an infected adult will recover from the disease [1]	ρ_a	0.942, 0.98	[29, 30]
The probability an infected child will recover from the disease [1]	ρ_c	∼1	[31]
The probability an infected adult will not recover from the disease [1]	ψ_a	0.05, 0.02	[20, 30]
The probability an infected child will not recover from the disease [1]	ψ_c	∼0	[31]

Two Class Age‐Structured Epidemic Model

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

Data from several epidemiological studies show that children and adults transmit the disease at different rates.[ , ] In addition, adults, on average, have a much longer recovery duration compared to children .[ , ] Therefore, an SIR model which takes into consideration the different age groups better represents the epidemiological population dynamics. An extension of the SIR model to two age‐classes has been investigated for explanation of Polio outbreak by Bunimovich‐Mendrazitsky and Stone.[ ]

Model Definition

The model considers a constant population with a fixed number of individuals N. Each individual belongs to one of the three groups: susceptible (S), infected (I), and recovered (R) such that

. When an individual in the susceptible group (S) is exposed to the infection, it is transferred to the infected group (I). The individual stays in this group on average days, after which it is transferred to the recovered group (R).

In addition to these groups, we define a death group (D) which is associated with individuals that are not able to fully recover from the disease or succumb to death. Individuals from the infected group (I) recover from the disease in some chance and move to (R), while the others do not and move to (D). Therefore, in each time unit, some rate of infected individuals recover while others die or remain seriously ill.

We divide the population into two classes based on their age: children and adults because these groups experience the disease in varying degrees of severity and have different infection rates. In addition, adults and children are present in various discrete locations throughout many hours of the day which affects the spread dynamics. Individuals below age A are associated with the “children” age‐class while individuals in the complementary group are associated with the “adult” age‐class. The specific threshold age (A) may differ in different locations but the main goal is to divide the population into two representative age‐classes. Since it takes A years from birth to move from a child to an adult age group, the conversion rate is set as .

By expanding the designation to two age‐classes, we let and represent susceptible, infected, recovered, and death groups for children and adults, respectively such that

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

In addition, Kelvin and Halperin[ ] conclude that children may be asymptomatic but still act as transmission vectors for the virus. Thus, refers to children with asymptomatic infection who have minimal clinical symptoms, but are still able to infect others. As a result, it is assumed that the infant mortality rate is 0, as shown in the Table  .

Equations (2)–(12) describes the epidemic's dynamics.

In Equation (2), is the dynamical amount of susceptible individual children over time. It is affected by the following three terms: First, with rate , each infected child infects susceptible children. Second, with rate , each infected adult infects the susceptible children. Finally, children grow and pass from the children's age‐class to the adult's age‐class with transition rate α, reduced from the children's age‐class. is the size of the children population and used to take all variables as fixed proportions of the population N.

In Equation (3), is the dynamical amount of susceptible adult individuals over time. It is affected by the following three terms: First, children grow and pass from the children's age‐class to the adult's age‐class with transition rate α, added to the adult age‐class. Second, with rate , each infected child infects the susceptible adult. Finally, with rate , each infected adult infects a susceptible adult. is the size of the adult population and used to take all variables as fixed proportions of the population N.

In Equation (4), is the dynamical amount of infected individual children over time. It is affected by the following four terms. First, with rate , each infected child infects the susceptible adult. Second, with rate , each infected child infects a susceptible child. Third, individuals recover or die from the disease after a period . Finally, children grow and pass from the children's age‐class to the adult's age‐class with transition rate α, reduced from the adult's age‐class. While the last process has a minor impact relative to the first three processes, we do not neglect it in order to count edge‐cases.

In Equation (5), is the dynamical amount of infected adult individuals over time. It is affected by the following four terms. First, with rate , each infected child infects the susceptible adult. Second, with rate , each infected adult infects a susceptible adult. Third, individuals recover or die from the disease after period . Finally, children grow and pass from the children's age‐class to the adult's age‐class with transition rate α, added to the adult's age‐class. While the last process has a minor impact relative to the first three processes, we do not neglect it to count edge‐cases.

In Equation (6), is the dynamical amount of recovered individual children over time. It is affected by the following two terms: First, in each point, a portion of the infected children recover after period , which is multiplied by the rate of children that do recover from the disease . Second, children grow from birth and pass from the children's age‐class to the adult age‐class with transition rate α, reduced from the children's age‐class.

In Equation (7), is the dynamical amount of recovered adult individuals over time. It is affected by the following two terms. First, in each point, a portion of the infected adults recover after period which is multiplied by the rate of adults that do recover from the disease . Second, children grow from birth and pass from the children's age‐class to the adult age‐class with transition rate α, added to the adult age‐class.

In Equation (8), is the dynamical amount of dead individual children over time. It is affected by the portion of the infected children that do not recover after period which is multiplied by the rate of children that do not recover from the disease .

In Equation (9), is the dynamical amount of dead adult individuals over time. It is affected by a portion of the infected adult that do not recover after period which is multiplied by the rate of adults that do not recover from the disease .

It should be noted that both Equations (8) and (9) can be presented as one equation as and are accumulative populations which do not infect the dynamics separately. Nevertheless, to provide a better representation of the age‐based death in the epidemic, we chose to divide the dead individuals into the same age‐groups which allows later analysis of the death of children and adults separately. In summary, the interactions between disease stages presented in Figure 1 are modeled by the following system of eight coupled ordinary differential equations (ODE), Equation (11) with initial conditions at , Equation (12) as:

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

Table 2

The four equilibria solutions for Equations (11) and (12)

	EQ_1	EQ_2	EQ_3	EQ_4
S_c	0	(γ_c−α)N_cβ_cc	β_caNa2(α−γ_a)β_aaN_c	(α+γ_c)Ic∗β_ccI_c+β_caI_a
S_a	N	β_caNc2(γ_c−α)N_aβ_cc	(α−γ_a)N_aβ_aa	(γ_a−α)Ic∗β_acIc∗+β_aaIa∗
I_c	0	αN_cβ_cc	αN_cβ_ca	Ic∗
I_a	0	0	γ_aρ_aN_cβ_ca	Ia∗
R_c	0	γ_cρ_cN_cβ_cc	0	γ_cρ_cIc∗α
R_a	0	0	αN_cβ_ca	−γ_aρ_aIa∗α
D_c	0	αN_cβ_cc−−γ_cρ_cN_cβ_cc	0	Ic∗−γ_cρ_cIc∗α
D_a	0	0	αN_cβ_ca−γ_aρ_aN_cβ_ca	Ia∗+γ_aρ_aIa∗α

Numerical Solution

To obtain a better understanding of how different parameter values influence the system dynamics, we illustrate the behavior of the system using numerical analysis. Equations (2)–(9) are ODEs, first‐order, nonlinear, from to , where is time (marked by t) and is the population distribution of all eight populations (marked by , , , , , , , and ). We processed eight‐order Runge–Kutta integration on the differential system to enable numerical simulations. Runge–Kutta integration of the equations was implemented by Octave programming (version 5.2.0), using standard program lsode, for a set of initial conditions described in Eq (12).

The solutions of the system ((11), (12)) are shown in Figure  . In addition, the population size, adults, and children are set to be to present the distribution of the Israeli population in 2017 as published by the Israeli central bureau of statistics. Figure 3 shows the population group sizes of , , , , , , , and where the x‐axis in all eight graphs represents the time (in days) that has passed from the beginning of the dynamics and the y‐axis is the size of each population, respectively. A maximum in the percent of infected children and adults (39%, 89%) is reached in the 8th and 11th days as shown in Figure 3. Therefore, the maximum infected population was 89% of the whole population. Besides, all children infected and recovered after 17 days while all adults recovered after 82 days.

Figure 3

Numerical simulation of trajectories of Equations (11) and (12) and the parameter values from Table 1. The graphs show the evolution in time (days) of , , , , , , , and . 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.

Stability Analysis

In epidemiology, it is essential to quantify the severity of outbreaks of infectious diseases. The standard parameter indicates the severity called the basic reproduction number R 0. In the standard SIR model,[ ] R 0 is defined to be the ratio between individual infection rate and recovery rate, where almost everyone is susceptible (namely, ).

In order to find the asymptotic form g for given initial conditions and parameters, it is possible to use the next generation matrix approach.[ ] Let be the rate of transfer of individuals by all means except of appearance of new infections in compartment i, and be the rate of transfer of individuals out of compartment i. Let be the rate of appearance of new infections in compartment i. Equations (2)–(9) take the form: where

In equilibrium, and are both equal to zero so the derivatives at equilibrium, focusing on and from Equations (4) and (5), are mapped to the third and forth elements in vectors F and V giving the matrices and .

The next generation matrix is defined as .[ ]

Assume an initial condition marked as . After one time unit, the amount of infected individuals can be calculated using . Therefore, for any start condition and model parameters one can find the asymptotic form g by calculating where is the first index that satisfies .

It is possible to retrieve R 0 from the next generation matrix as it is the dominant eigenvalue of the matrix,[ ] as this matrix describes the total amount of infections caused by each class over the lifetime of the infection. The epidemic is assumed to be stable if . This means that for each infected individual there is less than one infected individual in any of the groups in the next time unit. The dominant eigenvalue of can be calculated using the characteristic equation which is a second‐order polynomial of R 0. The zero points of this polynomial are

Both and are biological properties of the disease; on the other hand, , and can be managed using social distance, quarantine, masks, and other methods. Figure  presents five projections of stability from the parameter space calculated using Equation (21). Figure 4a shows the projection where . There are no values such that and from the color gradient, it is possible to see that has slightly more influence on R 0 relatively to . Figure 4b shows the projection where . There are no values such that and from the color gradient, it is possible to see that has a minor to no influence on R 0 relatively to . Figure 4c shows the projection where . for any combination of such that . Figure 4d shows the projection where . for any combination of such that . Figure 4e shows the projection where . for any combination of such that .

Figure 4

2D projections of the space and their influence on R 0, as obtained from Equation (21). The green section is where . The color gradient is from blue (lower R 0) to red (higher R 0). The model parameters used are .

In order to obtain the equilibria points of the model and their stability properties, the Jacobian (J) is obtained by linearizing Equations (4) and (5) in the system (Equations (2)–(9)). The calculations show that

The system (Equations (11) and (12)) has four non‐trivial equilibria which may be found by setting all rates in Equation (11) to zero; marked by , and . Thus, equilibria are provided as the state vector of the eight population states in Equation (11), where

There are other equilibria for trivial cases such as or which are unrealistic for real‐world dynamics. Below, we investigate the stability of four nontrivial equilibria.

: is the case for for every (see Table 2). This is a trivial equilibrium where there is no epidemic at all. Because the model does not take into consideration the birth of new children, after A days all children become adults and is derived. By settings , all eigenvalues are negative for all parameter values except . Therefore, this equilibrium is stable if .

: is the case for for every (see Table 2). This equilibrium can be achieved if because in any other case (assuming ) exists time t such that and therefore . This equilibrium corresponding to the case where adults do not infect children at all which is improbable as a real world scenario.

: is the case for for every (see Table 2). This equilibrium can be achieved if because in any other case (assuming ) exists time t such that and therefore . This equilibrium corresponding to the case where children do not infect adults at all which is similarly to , improbable as a real world scenario.

: is the case for such that and (see Table 2). This equilibrium is the generic case for the model. This equilibrium does not have a unique epidemiological properties.

The equilibria , and are not stable because in each of the cases, either or . Assuming, the equilibria is kept until some time such that is defined as the first time t such that if > 0 or if . At time , or which is different then the values of and for equilibria , and .

Spatial Model

Figure  shows the spatial model schema presenting the populations distribution and locations (work/school, and home) at some time of the day. In addition to the ODE model's parameters (shown in Table 1), the following parameters are added to the hybrid model as part of the spatial model: We assume the transition from home to either work or school and back is immediate and that everybody is following the same clock. Each simulation step simulates 1 h. The population size is constant during the simulation and initialized in the beginning of each iteration by setting children population size and adult population size . In each simulation step, the following three actions take place: The spatial model adds day–night circle and three main locations to the dynamics of the hybrid model. A description of the whole spatio‐temporal dynamics of the Hybrid model as a system of ODEs can be found in Equations (S1)– (S16), Supporting Information.

, the average number of meeting events between adults and children per hour.

, the average number of meeting events between adults and adults per hour.

, the average number of meeting events between children and children per hour.

, hours of the day that children are at home and hours of the day that children are at school.

, hours of the day that adults are at home and hours of the day that adults are at work.

If a member is in the susceptible group and meets other members of the infected group, there is a change of , or according to the age‐class of the two members that the first will be infected.

Each infected child or adult that was infected for simulation steps becomes either recovered or deceased, respectively.

According to the hour of the day, the adults transition to home or work and the children to home or school.

Figure  shows the population group sizes of , , , , , , , and as an average of ten iterations of the hybrid model. The x‐axis is the time (in days) that has passed from the beginning of the epidemic and the y‐axis is the normalized size of each population, respectively. The parameters used in the simulations are .

Figure 6

Average of ten iterations of numerical simulation of the hybrid model where the parameter values are taken from 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 on the 9th and 12th day, as shown in Figure 6. In addition, all children were infected and recovered after 12 days while all adults recovered after 28 days. The simulation shows similar results to the results derived by solving the model (Equations (11) and (12)) as presented in Figure 3. The main difference between the two is that the simulation predicts 4.43 times shorter duration from the time of maximum infected adults to the time that all adults are either dead or recovered, as the infected adult population equals zero on the and 28th days. In addition, the maximum of infected adults is reached on the 11th and 12th days resulting in 71 and 16 days for full recovery.

Results

Outbreak Analysis

The proposed hybrid model provide an in silico environment allowing relatively fast, cheap, and accurate analysis of different policies on the COVID‐19 spread dynamics.

One of the major hopes of politicians is for an NPI policy in which the epidemic does not reach an outbreak at any time after they initialized a given decision. We define this condition mathematically as follows:

Definition 2. For given initial condition and model's parameter . A solution of Equations (11) and (12) is defined as outbreak dynamics if .

We will examine two policies, based on this condition, to determine if each one is possible to fulfill the condition. If so, the optimal NPI policy is based on the parameter‐space criteria:

The influence of the duration of the work/school day.

Lockdown in homes with partial to full separation between individuals.

Duration of Working and Schooling Day

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

Based on the proposed spatial model, we change the number of hours children and adults spend in school and work each day, respectively. Figure  shows the average R 0 as a function of the duration in hours of the working and schooling () day. The dots are the calculated values from the simulator and the surface is a fitting function . The fitting function is calculated using the least mean square (LMS) method.[ ] In order to use the LMS method, one needs to define the family function approximating a function. The family function has been chosen to balance between the accuracy of the sampled data on the one hand and simplicity of usage on the other,[ ] and was obtained with a coefficient of determination . Similarly, Figure 7c shows the average as a function of the duration in hours of the working and schooling () day. The dots in Figure 7c are the calculated values from the spatial model and the surface is the fitting function obtained with a coefficient of determination .

Figure 7

Analysis of the epidemic spread as a function of working () and schooling () duration in hours each day.

The duration of either working or schooling day has a local minimum at with , as shown in Figure 7a. The optimal point from Equation (25) is with obtained by setting the gradient of to zero. The inconsistency in the values is resulted by the error in the fitting function but both show the same behavior in which longer working hours reduce significantly the infection rate but too many hours increase the infection rate back as the lack of circulation of contagious adults and children during the day helps to reduce the infection rate. Furthermore, Figure 7b presents a binary classification where red cells are cases with an outbreak and green cells are cases without outbreak during the simulation as a function of the working () and schooling () duration in hours each day.

On the other hand, the duration of either working or schooling day has an effect of 10% (as the maximal value is 82.5 while the minimal is 72.5) on the maximal infected percent of individuals from the population as shown in Figure 7c and by setting the limits in Equation (26). There is a sharp decline between a shorter working–schooling duration of 10 h or less and longer than that. This is associated with the fact that with more than 12 h (half‐day) the dynamics that have a higher incidence are the ones with which reduces the number of infected individuals in total.

Lockdown Policy

Partial or full lockdown of several locations or entire countries were broadly used at the beginning of the COVID‐19 outbreak as a policy to reduce the number of infected individuals and control the spread dynamics. The lockdown policy yields social distance which reduces the ratio of infection. On the other hand, this policy negatively affects the economy and mental health, and increases the presence of other diseases (both physical and mental) in the population. Therefore, the optimization task of finding the minimal portion of the population to be locked down such that the epidemic will be constrained is important.

The lockdown policy is similar to the schooling‐working hours policy in the manner that both modify the spatial dynamics of the population. Nevertheless, the schooling‐working hours policy defined the number of hours all the children and working adults populations go to school and work, respectively, while the lockdown policy keeps part (or all) the population at home all day long alongside the remain part of the population keeps the regular working and schooling hours. In addition, the lockdown policy isolates individuals at home, which is expressed by the fact that individuals can contact with them but they can not initial an contact with other individuals while this constraint does not take place in the working‐schooling hours policy.

The optimization problem can be written formally as where and are the portion of adults and children in lockdown. We ran the spatial model multiple times where each time a combination of of the population is assumed to be in lockdown. Figure  shows the results of this calculation. Each dot represents R 0 of each case . The black grid shows the threshold .

Figure 8

Analysis of the epidemic spread as a function of adult lockdown () and children lockdown () using the computer simulation of the hybrid model.

The behavior of R 0 as a function of has been retrieved using the LMS method and takes the form

obtained with a coefficient of determination . Therefore, it is safe to claim that function is well fitting the data despite the stochastic noise of the simulation and presents a fair approximation for the R 0 behavior as a function of . Using Equation (28), it is possible to find the constraints of such that , by solving .

Closing only schools without any reduction in adults going to work does not prevent an epidemic outbreak, while locking down half of the adult population will prevent an outbreak. Lockdown of children have a minor effect relatively to lockdown adults as shown in both Equation (28) and Figure 8b. In addition, the same phenomena repeat in the max infected individuals, as shown in Figure 8c. The surface calculated using the LMS method: and obtained with a coefficient of determination .

In addition, we repeat the analysis by numerically solving the system (Equations (11) and (12)). The lockdown is reflected in the model as the infection rate between every two individuals. affects the interactions between adults and either other adults or children , and affect the interactions between children and either other children or adults . Mark the original values of from Table 1 as . Each time, the system solved where such that , , , and .

Figure  shows the results of this calculation based on Equation (21). Each dot represents R 0 of each case . The black grid shows the threshold . Figure 9b presents a binary classification of the cases with outbreak dynamics (red) or without (green). Figure 9b is similar to Figure 8b such that adult lockdown has a higher influence on the outbreak and predicts that a much higher percentage of the population will be in lockdown. Again, the difference in the percent of adult's lockdown such that can be associated with the slower infection rate in hours of each day because part of the adults change their location which in it turns influence the infection rate. In addition, the recovery process is faster in the spatial model (as shown in Figures 3 and 6).

Figure 9

Analysis of the epidemic spread as a function of adult lockdown () and children lockdown () by the hybrid model.

The behavior of R 0 as a function of has been derived using the LMS method and takes the form obtained with a coefficient of determination .

In addition, Figure 9c presents the max infected individuals as a function of adult and children lockdown. The surface has been calculated using the LMS method and takes the form: obtained with a coefficient of determination . Figures 8c and 9c are presenting the same behavior while Figure 8c predicts 5% more infected individuals on average. The influence of the adult lockdown is one order of magnitude more significant than that of the children's lockdown in a first and second order approximation as shown in Equations (29) and (31). That is, where there is no lockdown for children , there is a local maximum in infected individuals regardless of the lockdown of adults .

Validation of the Hybrid Model on the Dynamics of the COVID‐19 Pandemic in Israel

In Israel, in February 21, the first COVID‐19 infected individual was detected. On March 25, a national quarantine was implemented that continued for 2 months with local relief and restriction on several occasions and in several cities. Figure  presents the daily new confirmed cases from August 15 to September 28, respectively. Where the blue (solid) line marks the date when the schools returned to almost full capacity, the green (solid‐dotted) lines are linear regression for before and after school opening, respectively. It is easy to see a significant increment in the number of new daily cases.

Figure 10

Daily confirmed cases in Israel between August 15 and September 29 (2020).[ ] The blue (solid) line is the full opening of schools. The green (solid‐dotted) lines are a linear regression for before and after school opening, respectively.

From ref. [4], there were 90 232 infected at August 15. Given that the average recovery rate of children is 2 days and adults is 14 days (see Table 1) and that children is 28% of the population, we obtain that there was approximately 18 899 infected adults and 868 infected children. We assume that the amount of recovered individuals is heterogeneous in the population and that only adults die from the epidemic until this point. Equation (32) shows the initial conditions for August 15 (t 0), Israel.

To test the hybrid model, we calculated the R 0 on average per day for the period from August 15 to September 1 and from September 1 to September 15, dividing the number of new cases by the total number of cases on that day (obtained from ref. [4]). We solve numerically the hybrid model and calculate R 0 for the case in which schools are open versus the case in which schools are closed. If schools are open, children go to school three times a week (every other day, except weekends) for 5 h each day. If schools are closed, children stay at home all the time. In both cases, we assume that adults work 8 h each day, except on weekends when they stay home. The obtained values of R 0 for both cases are shown in Figure  , with a mean square error (MSE) of 0.205.

Figure 11

R 0 between August 15 and September 15 (2020) comparison between the historical data and the hybrid model predictions.

Considering that, longer school day is reducing the infection rate (Figure 7a), we solve numerically the hybrid model and calculate R 0 for the case when the schools are open (from September 1). We assume that children go to school every day (except weekends) for 9 h while the adults go to work for the same period. We obtain that the difference in the average R 0 of the case with the optimal NPI () and the historical NPI () is (see Figure 11).

Discussion

This study presents a model showing the effects of population age, time of day, and gathering location on the spread of the epidemic. Based on the proposed model Equations (11) and (12), the non‐trivial equilibria of the system is presented in Table 2. Although the proposed model is a simplification of the total complexity of the COVID‐19 epidemic, the asymptotic solution without any intervention will result in a high percentage 39% children, 89% adults and 38.5% children, 100% adults of infected individuals as shown in Figures 3 and 6, respectively.

In addition, the asymptotic solutions of the model are insignificant for themselves as in practice these results occur after A days (13 years). The main epidemic dynamics as obtained from both the SIRD model Equations (11) and (12) and the hybrid model is a few month long, as shown in Figures 3 and 6. Besides, the asymptotic solution can be used to easily obtain and as they accumulative sums of infected infected adults and children multiplied by a factor and and do not change after the epidemic is over. As a result, is not used as for 82 days (as shown in Figure 3) only 0.017% of the children will become adults which effectually has only a minor affect on the epidemic results.

Furthermore, two policies for controlling the epidemic and preventing outbreak have been investigated using the hybrid model. First, the duration of working and schooling day has an influence (up to 10%) on the maximum of infected individuals but can prevent outbreaks under the assumption that during this time adults do not contact children at all as shown in Figure 7. These results match the conclusions reported by Keskinocak et al. (2020) for Georgia state.[ ] In addition, these results match the results reported by Di Domenico[ ] that school closure has a minor influence on the epidemic peak, just delays it. Second, a partial lockdown of adults and children shows that adults' and children's lockdowns have a similar first‐order effect (linear coefficients of and ) but the combination between the two ( has a bigger weight on the average infection rate as shown in Equations (28) and (30). These results match the conclusions reported by Aglar et al. (2020) for the state of Georgia regarding the voluntary lockdown with school closure.[ ] By the same token, a lockdown of approximately half of the adult population prevents an outbreak where children lockdown has a less significant influence as shown in Figures 8 and 9. These results are based on values from Table 1, which can vary significantly between countries and depend on other hyperparameters such as population density, age distribution of the population and others.

The model explains the dynamics that took place between August 15 and September 15 (2020) in Israel as presented in Figure 11, with MSE of 0.205. The opening of the schools is equal to reducing the lockdown for children with some relaxation (smaller ) of the lockdown for adults (because if the children are in school they can go to work) with less than half (50%) of the adult population voluntarily in lockdown. Also, the schooling hours were reduced to less than 6 h each day. From Figures 7a and 8c, it is easy to see that this policy predicts high increase in the infection rate, as indeed happened. Keeping the schools open while keeping the increase in the infection rate from increasing significantly is possible if the schooling hours are longer (8–9 h) as shown in Figure 7a. The influence of this policy in Israel during the school opening which take place in September 1 shows that the R 0 can be reduced by 0.83 in comparison to a policy in which children go to school every other day for 5 h, as shown in Figure 11. Also, if at least half of the adult population will be in lockdown, the influence of the schools on the infection rate will be relatively small as shown in Figure 8b.

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

Conflict of Interest

The authors declare no conflict of interest.

Supporting Information

Click here for additional data file.