The signature features of COVID-19 pandemic in a hybrid mathematical model - implications for optimal work-school lockdown policy

Two class age-structured epidemic model

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

Data from several epidemiological studies show that children and adults transmit the disease at different rates [17, 19]. In addition, adults, on average, have a much longer recovery duration compared to children [6, 8]. Therefore, a 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 [25].

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 N = |S| + |I| + |R|. When an individual in the susceptible group (S) exposed to the infection, it is transferred to the infected group (I). The individual stays in this group on average d_I→R 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 S_c, I_c, R_c, D_c and S_a, I_a, R_a, D_a 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 [26] while the mortality rate in 2017 was 8.6 for every 1000 individuals [27], 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.

Eq (1-10) describes the epidemic’s dynamics. In Eq (1), is the dynamical amount of susceptible individual children over time. It is affected by the following three terms. First, with rate β_cc each infected child infects susceptible children. Second, with rate β_ca 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. N_c is the size of the children population and used to take all variables as fixed proportions of the population N. In Eq (2), 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 β_ac each infected child infects the susceptible adult. Finally, with rate β_aa each infected adult infects a susceptible adult. N_a is the size of the adult population and used to take all variables as fixed proportions of the population N. In Eq (3), is the dynamical amount of infected individual children over time. It is affected by the following four terms. First, with rate β_ca each infected child infects the susceptible adult. Second, with rate β_cc each infected child infects a susceptible child. Third, individuals recover or die from the disease after a period γ_c. 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 Eq (4), is the dynamical amount of infected adult individuals over time. It is affected by the following four terms. First, with rate β_ac each infected child infects the susceptible adult. Second, with rate β_aa each infected adult infects a susceptible adult. Third, individuals recover or die from the disease after period γ_a. 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 Eq (5), 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 γ_c which is multiplied by the rate of children that do recover from the disease ρ_c. 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 Eq (6), 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 γ_a which is multiplied by the rate of adults that do recover from the disease ρ_a. 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 Eq (7), 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 γ_c which is multiplied by the rate of children that do not recover from the disease ψ_c. In Eq (8), 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 γ_a which is 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. 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 Fig 1 are modelled by the following system of eight coupled ordinary differential equations (ODE) Eq (9) with initial conditions at t = 0, Eq (10) 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 β_ca calculated from the data of the correlated source [6, 17]. The threshold of the children’s age to become adults in parameter A is set to 13 year 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. [6] (Table 2). The susceptible contacts in children who become infected due to direct disease transmission from an adult β_ca 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. [17]. Therefore, the infected rate is set to .

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. Eq (1-8) are ODEs, first-order, nonlinear, from ℝ to ℝ⁸, where ℝ is time (marked by t) and ℝ⁸ is the population distribution of all eight populations (marked by 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)). 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 (10).

The solutions of the system (9-10) are shown in Fig 3. In addition, the population size, adults, and children are set to be N = 8000000, N_c = 2240000, N_a = 5760000 to present the distribution of the Israeli population in 2017 as published by the Israeli central bureau of statistics. Fig 3 shows the population group sizes 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) 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 Fig 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.

• Download figure
• Open in new tab

Fig 3. Numerical simulation of trajectories of Eq (9 - 10) and the parameter values from 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.

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₀. In the standard SIR model [22], R₀ is defined to be the ratio between individual infection rate and recovery rate, where almost everyone is susceptible (namely, S_c + S_a ∼ N). For any start condition and parameters, at t →∞ the model’s state is an instance of Eq (10) according to Lemma 1.

where g is a function g: ℝ¹⁷ →ℝ³ = [g₁, g₂, g₃] such that g(ℙ), where Recall that N = |S_c| + |I_c | + |R_c | + |D_c | + |S_a | + |I_a | + |R_a | + |D_a |.

Proof. Define a non-homogeneous continuous-time Markov chain M with state space 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 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 states i₀, …i_t, j and h is an arbitrary small step in time [31]. M satisfies Eq (12) if both v₁(t), v₂(t) depend only on the last state X_t = i_t. 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) and (14), v₁(t + h) and v₂(t + h) depends only on the previous state X_t = i_t, respectively. Let π ∈ ℕ⁸ be the initial condition for the system. Therefore, the model’s state at some time t is defined by A^tπ. For any initial condition π, the vector lim_t→∞(A^tπ) takes the form where g₁, g₂, g₃ ∈ [0, …, N].

From Table 2 one can see that R_a, D_c, and D_a are stationary values because . By setting π = [N_c, N_a, 0, 0, 0, 0, 0, 0] one gets ∀t ∈ℕ: v₁(t) = v₂(t) = 0 and therefore the sub matrix of M with only the states S_c and S_a takes the form Therefore, M (t → ∞) = [0, N, 0, 0, 0, 0, 0, 0] because and for lim_t→∞(1 − α)^t = 0 meaning that asymptotically the population in S_c decreases 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 possible to use the next generation matrix approach [32]. 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 F_i be the rate of appearance of new infections in compartment i. Eq (1 - 8) takes the form: where In equilibrium, I_c and I_a are both equal to zero so the derivatives at equilibrium, focusing on I_c and I_a from Eqs. (3) and (4), are mapped to the third and forth elements in vectors F and V giving the matrices F and V. The next generation matrix is defined as FV⁻¹ [32]. Assume an initial condition marked as v = [v₁, v₂]. After one time unit, the amount of infected individuals can be calculated using v^′ = Gv. Therefore, for any start condition and model parameters one can find g by calculating where C ∈ℕ is the first index that satisfies I_c(C) + I_a(C) = 0.

It is possible to retrieve R₀ from the next generation matrix G as it is the dominant eigenvalue of the matrix [25], 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 R₀ ≤ 1. 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 G can be calculated using the characteristic equation which is a second-order polynomial of R₀. The zero points of this polynomial are Both γ_c and γ_a are biological properties of the disease; on the other hand, β_aa, β_ac, β_ca, and β_cc can be managed using social distance, quarantine, masks, and other methods. Figs 4a - 4e presents five projections of stability from the parameter space {β_aa, β_ac, β_ca, β_cc} calculated using Eq (16). Fig 4a shows the β_aa − β_ac projection where β_cc = 0.308, β_ca = 0.266. There are no values such that R₀ < 1 and from the color gradient, it is possible to see that β_ac has slightly more influence on R₀ relatively to β_aa. Fig 4b shows the β_aa − β_ca projection where β_cc = 0.308, β_ac = 0. There are no values such that R₀ < 1 and from the color gradient, it is possible to see that β_aa has a minor to no influence on R₀ relatively to β_ca. Fig 4c shows the β_cc − β_ca projection where β_aa = 0.308, β_ca = 0.266. R₀ < 1 for any combination of β_cc, β_ac such that β_ac < 0.08 − 1.6β_cc. Fig 4d shows the β_cc − β_ac projection where β_aa = 0.308, β_ac = 0. R₀ < 1 for any combination of β_cc, β_ac such that β_cc ≤ 0.3. Fig 4e shows the β_aa − β_cc projection where β_ca = 0.266, β_ac = 0. R₀ < 1 for any combination of β_aa, β_cc such that β_cc < 0.9 − 0.155β_aa.

• Download figure
• Open in new tab

Fig 4. 2-dimensional projections of the {β_aa, β_ac, β_ca, β_cc} space and there influence on R₀.

The green section is where R₀ ≤ 1. The color gradient is from blue (lower R₀) to red (higher R₀). The model parameters used are α = 8.78 · 10⁻⁶, γ_c = 0.5, γ_a = 0.0714.

In order to obtain the equilibria points of the model and their stability properties, the Jacobian (J) is obtained by linearizing Eqs (3) and (4) in the system (Eqs (1 - 8)). The calculations show that The system (Eqs (9-10)) has four non-trivial equilibria which may be found by setting all rates in Eq (9) to zero; marked by EQ₁, EQ₂, EQ₃, and EQ₄. Thus equilibria are 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 unrealistic for real-world dynamics. Below we investigate the stability of four nontrivial equilibria.

EQ₁: is the case for I_c(t) = I_a(t) = 0 for every t ∈ℕ (see Table 3). 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 EQ₁ is derived. By settings J (EQ₁), all eigenvalues are negative for all parameter values except . Therefore, this equilibrium is stable if .

EQ₂: is the case for I_a(t) = 0 for every t ∈ℕ (see Table 3). This equilibrium can be achieved if β_ac = 0 because in any other case (assuming S_a(0) > 0) exists time t such that and therefore I_a(t + 1) ≠ 0.

EQ₃: is the case for I_c(t) = 0 for every t ∈ℕ (see Table 3). This equilibrium can be achieved if β_ca = 0 because in any other case (assuming S_c(0) > 0) exists time t such that and therefore I_c(t + 1) ≠ 0.

EQ₄: is the case for ∃t such that I_c(t) ≠ 0 and I_a(t) ≠ 0 (see Table 3). This equilibrium is the generic case for the model.

According to Lemma 1, EQ₂, EQ₃, and EQ₄ are not possible asymptotic forms of the system and are therefore not stable.

Spatial model

Fig 5 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 follows parameters are added to the hybrid model as part of the spatial model:

• Download figure
• Open in new tab

Fig 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₀ at a certain time and the average R₀ from the beginning of the simulation.

1. ϕ_ac the average number of meeting events between adults and children per hour.

2. ϕ_aa the average number of meeting events between adults and adults per hour.

3. ϕ_cc the average number of meeting events between children and children per hour.

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

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

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 one hour. The population size is constant during the simulation and initialized in the beginning of each iteration by setting children population size N_c and adult population size N_a. In each simulation step, the following three actions take place:

1. If a member is in the susceptible group and meets other members of the infected group there is a change of β_aa, β_ac, β_ca, or β_cc according to the age-class of the two members that the first will be infected.

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

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

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

Fig 6 shows the population group sizes 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) as an average of 10 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 .

• Download figure
• Open in new tab

Fig 6. Average of 10 iteration 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 9_th and 12_th day, as shown in Fig 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 (Eqs (9-10)) as presented in Fig 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 82_nd and 28_th days. In addition, the maximum of infected adults is reached on the 11_th and 12_th days resulting in 71 and 16 days for full recovery.

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 a 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:

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:

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

2. 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. Fig 7a shows the average R₀ as a function of the duration in hours of the working (d_W) and schooling (d_S) day. The dots are the calculated values from the simulator and the surface is a fitting function R₀(d_W, d_S). The fitting function is calculated using the least mean square (LMS) method [33]. In order to use the LMS method, one needs to define the family function approximating a function. The family function

• Download figure
• Open in new tab

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

has been chosen to balance between the accuracy of the sampled data on the one hand and simplicity of usage on the other [34], and was obtained with a coefficient of determination R² = 0.815. Similarly, Fig 7c shows the average max_t(I_c(t) + I_a(t)) as a function of the duration in hours of the working (d_W) and schooling (d_S) day. The dots in Fig 7c are the calculated values from the spatial model and the surface is the fitting function obtained with a coefficient of determination R² = 0.748.

The duration of either working or schooling day has a local minimum at (d_W, d_S) = (18, 14) with R₀ = 0.808, as shown in Fig 7a. The optimal point from Eq (17) is (d_W, d_S) = (9, 17.5) with R₀ = 0.967 obtained by setting the gradient of R₀(d_W, d_S) 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, Fig 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 (d_W) and schooling (d_S) 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 Fig 7c and by setting the limits d_W = {0, 24}, d_S = 0, 24 in Eq (18). There is a sharp decline between a shorter working-schooling duration of 10 hours or less and longer than that. This is associated with the fact that with more than 12 hours (half-day) the dynamics that have a higher incidence are the ones with β_ca = β_ac = 0 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 optimization problem can be written formally as where L_a and L_c 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. Fig 8a shows the results of this calculation. Each dot represents R₀ of each case L_a, L_c. The black grid shows the threshold R₀ = 1.

• Download figure
• Open in new tab

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

The behavior of R₀ as a function of L_a, L_c has been retrieved using the LMS method and takes the form obtained with a coefficient of determination R² = 0.971. Therefore, it is safe to claim that function R₀(a, c) is well fitting the data despite the stochastic noise of the simulation and presents a fair approximation for the R₀ behavior as a function of L_a, L_c. Using Eq (19) it is possible to find the constraints of L_a, L_c such that R₀ ≤ 1, by solving R₀(L_a, L_c) ≤ 1.

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 (L_c) have a minor effect relatively to lockdown adults (L_a) as shown in both Eq (19) and Fig 8b. In addition, the same phenomena repeat in the max infected individuals as shown in Fig 8c. The surface calculated using the LMS method: and obtained with a coefficient of determination R² = 0.840.

In addition, we repeat the analysis by numerically solving the system (Eq (9-10)). The lockdown is reflected in the model as the infection rate between every two individuals. L_a affects the interactions between adults and either other adults or children β_aa, β_ac and L_c affects the interactions between children and either other children or adults β_cc, β_ca. Mark the original values of β_aa, β_ac, β_ca, β_cc from Table 1 Each time, the system solved where such and .

Fig 9a shows the results of this calculation based on Eq (16). Each dot represents R₀ of each case L_a, L_c. The black grid shows the threshold R₀ = 1. Fig 9b presents a binary classification of the cases with outbreak dynamics (red) or without (green). Fig 9b is similar to Fig 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 (L_a) such that R₀ < 1 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 Figs 3 and 6).

• Download figure
• Open in new tab

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

The behavior of R₀ as a function of L_a, L_c has been derived using the LMS method and takes the form obtained with a coefficient of determination R² = 0.952.

In addition, Fig 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 R² = 0.702. Figures 8c and 9c are presenting the same behavior while Fig 8c predicts 5% more infected individuals on average. The influence of the adult lockdown is one order of magnitude more significant that of the children’s lockdown in a first and second order approximation as shown in Eqs (20) and (22). That is, where there is no lockdown for children (L_c = 0) there is a local maximum in infected individuals regardless of the lockdown of adults (L_a).

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 two months with local relief and restriction on several occasions and in several cities. Fig 10 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.

• Download figure
• Open in new tab

Fig 10. Daily confirmed cases in Israel between August 15 and September 29 (2020) [4].

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 [4] there was 90232 infected at August 15. Given that the average recover 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 18899 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. Eq (23) shows the initial conditions for August 15 (t₀), Israel. To test the hybrid model, we calculated the R₀ 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 [4]. We solve numerically the hybrid model and calculate R₀ for the case when the schools are open versus the case when the schools are closed. If schools are open, children go to school three times a week, every other day and do not go on weekends for five hours. If schools are closed, children stayed at home all the time. In both cases, we assume that adults will work eight hours each day, except on weekends when they stay home. The obtained values of R₀ for both cases are shown in Fig 11 with mean square error (MSE) of 0.205.

• Download figure
• Open in new tab

Fig 11. Daily R₀ between August 15 and September 15 (2020) comparison between the historical data and the hybrid models predictions.

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