SCHOOL RE-OPENING SIMULATIONS WITH COVID-19 AGENT-BASED MODEL FOR QUEZON CITY, PHILIPPINES

The COVID-19 pandemic prompted schools to close in 2020, and with its negative impact on the quality of learning, the conversation shifts to the re-opening scenarios. In this study, we coupled a COVID-19 agent-based model (ABM) with social contact probabilities from pre-pandemic estimates. We then simulated school re-opening and vaccination scenarios in Quezon City, Philippines using our ABM. Our toy simulations suggest that the city could already re-open schools with 50% vaccination coverage. However, we suggest that students shall be vaccinated first, mask-wearing and physical distance shall be strictly observed, and schools shall only be reopened by 25% as a precaution. Policymakers may take insights from the study.


INTRODUCTION
In the Philippines, schools have closed since March 2020 when the COVID-19 pandemic hit. But this closure may result in longterm losses, warned the Asian Development Bank (2021) in an earlier report, projecting a $180-worth annual income decline with the quality of training the students are getting via distance learning. However, it also remains that the pandemic is a threat to life. With the threat of COVID-19 to life and of school closure to students' future, we ask: Should we already re-open schools? This study explores this question by simulating school reopening scenarios, vaccination scenarios, and the combination, using a COVID-19 agent-based model coupled with contact probabilities designed for Quezon City, Philippines, a coronavirus hotspot.

Agent-Based Modelling
Agent-based models are computer programs designed to describe a system. They are composed of three important elements: agents, environment, and time. The agents are independent computing entities; they compose the system. In the modeling of COVID-19, the agents represent human beings, the agents of disease transmission. Next, the environment is the physical setting of an agent. In this study, the environment is the geographic space of Quezon City, Philippines. Finally, time refers to the scale of a model step. In this study, a single model step is equivalent to one day.

COVID-19 Agent-Based Model for Quezon City, Philippines
We developed a COVID-19 agent-based model (ABM) for Quezon City, Philippines. We based the design of our new model on our previous works (Bongolan et al., 2021;Rayo et al., 2020;Minoza et al., 2020). In addition, we coupled the equivalent model with social contact probabilities from the pre-pandemic estimates of Prem et al. (2018). The contact probabilities were used to quantify school re-opening scenarios.

Mathematical Basis:
The Age-Stratified, Quarantine-Modified SEIR with Non-Linear Incidence Rates (ASQ-SEIR-NLIR) is a compartmental model describing COVID-19 transmission (Bongolan et al., 2021;Rayo et al., 2020;Minoza et al., 2020). It is defined by a system of differential equations: From the classic compartmental model, SEIR stands for susceptible-exposed-infectious-removed, four classes where human beings belong concerning the disease. In addition, parameters , , and  represent the transmission, incubation, and removal rates, respectively.
Three modifications were made to the classic SEIR to account for factors affecting the transmission of COVID-19. The timevarying quarantine variable Q(t) accounts for quarantine (Bongolan et al., 2021); the age-stratified infection expectation U accounts for age (Rayo et al., 2020); and the non-linear incidence rates  and  account for behavioral factors (e.g., wearing of masks, physical distance, etc.) and disease-resistance factors (e.g., natural immunity, healthy living, etc.) (Minoza et al., 2020).

Equivalent Agent-Based Model Parameters:
We transcribed ASQ-SEIR-NLIR as an agent-based model. To do this, we converted its parameters as agent-based model parameters.
First, the Age-Stratification Theory (age-strat) (Bongolan et al., 2021), accounted as U in ASQ-SEIR-NLIR (Rayo et al., 2020), was implicitly treated in the new model. This was done by agestratifying the parameters of the model. Parameters  (transmission rate) and  (removal rate) are now multi-valued parameters where each element represents the rate for each age group, i.e., 0 to 4, 5 to 9, …, and 80 and above. In addition, the removal rate  has been disaggregated into two separate parameters: the death and the recovery probabilities. The way to compute for the values of the three parameters also adopts the computation in age-strat: The age distributions for infections, deaths, and recoveries are treated as proxy values for the three parameters: infection, death, and recovery probabilities.
Next, parameter  (incubation rate) was converted as the average incubation period (days) of agents once exposed to an infectious agent. This is analogous to the COVID-19 incubation period. We also added a similar parameter for infection: the average infectious period, also in days.
Finally, the behavioral factors and disease-resistance factors  and  (Minoza et al., 2020) were disaggregated as follows: wearing of masks, the percentage of individuals wearing a mask; physical distancing, the percentage of individuals observing physical distance; protection from wearing of masks, the percentage of protection from contracting the disease when wearing a mask; and finally, protection from physical distancing, the percentage of protection from contracting the disease when observing physical distance For the time-varying quarantine variable Q(t), the conversion is accounted for through the social contact probabilities which is discussed in section 2.2.4.

Agent Transition Rules:
In our new model, we disaggregated the removed (R) compartment into three separate compartments: recovered (R), dead (D), and vaccinated (V). We assume that once recovered from the disease, an agent is forever immune (a common assumption for viral infections). Note that we also assume that vaccinated individuals belong to the removed compartment. This implies that once vaccinated, an agent is also forever immune from the disease.
The S-E-I-R-D-V transition of agents is a probability game. A susceptible agent becomes exposed when it interacts with an infectious agent, is not protected by wearing a mask or observing physical distance, and the infection probability is hit. Once exposed, an agent becomes infectious after the average incubation period. An infectious agent then either dies or recovers after the average infectious period and when either of the death probability or the recovery probability is hit. For the vaccinated, we assume that only susceptible agents are eligible for the vaccination, and vaccination only happens at the beginning of the simulation. This implies that the vaccination curve is constant, essentially freezing vaccination throughout the simulation.
The transition rules assume that the total population N is constant, i.e., N = S + E + I + R + D + V; no deaths from other causes or additional agents from births, i.e., for now, we ignore the population's vital statistics.

Social Contact Probabilities for the Philippines:
The social contact probabilities define the probability of interaction among agents. These probabilities are computed from the prepandemic social contact estimates of Prem et al. (2018) for the Philippines. These are shown as contact matrices ( Figure 1). They define the average number of interactions individuals have in their day, in different locations, i.e., school, work, home, and others. Model. Green dots represent the susceptible, orange for exposed, red for infectious, and grey for removed (i.e., dead, recovered, and vaccinated).
The contact matrices allowed us to simulate school re-opening scenarios. The school matrix was scaled according to the reopening scenario of interest. For a no school opening scenario, the scale is 0, resulting in a zero matrix; for a 25% opening, the scale is 0.25, i.e., school matrix * 0.25; for a 50% opening, the scale is 0.50, i.e., school matrix * 0.50; and so on. The resulting matrices were then added (matrix addition) and max normalized so that the values inside the resulting matrix are between 0 and 1 (social contact probabilities). These probabilities define the probability of interaction among agents. For instance, when two agents, say agents 1 and 2, are in the same space, the probability that they will interact will depend on the age groups where agents 1 and 2 belong. The contact probabilities offer flexibility in controlling the interaction of the agents.
With the contact probabilities, quarantine scenarios could also be simulated. Adopting the design of the time-varying quarantine variable Q(t) of ASQ-SEIR-NLIR, the effectiveness of quarantine could be simulated, same as school re-opening, by scaling the contact probabilities: contact probabilities * 0.60 represents 40% effectiveness; contact probabilities * 0.40 represents 60% effectiveness; and so on. In addition, the scaling may be done by row and column to achieve age-group level granularity. This is especially useful in cases where quarantine is implemented by age, e.g., only individuals between ages 15 -65 are allowed to go outside their homes. This age-nuanced feature is one of the core strengths of using our contact probabilities.

Model Implementation using Python Mesa and MesaGeo:
The model was implemented using two application libraries for Python: Mesa (Project Mesa Team, 2021) and MesaGeo (Corvince, 2018). The Mesa library forms the core of the model, accounting for agent instantiation, agent calls, data collection, and model starting and stopping. Essentially, Mesa provided the scaffold for the model. The MesaGeo library on one hand is an extension library of Mesa. It allowed us to integrate and use a GIS component in the model (Figure 2).
The GIS component provided the environment component of the model, i.e., the six districts of Quezon City, Philippines, and the visualization for the environment and the agents inside it. In addition, it also let us to capture population density. Initially, agents are assigned at their home districts, resulting in higher density in districts with small geographic areas and big populations.

Simulation of School Re-opening and Vaccination
Scenarios in Quezon City, Philippines

School Re-opening Sensitivity Analysis:
We performed a school re-opening sensitivity analysis. Five scenarios were simulated using our agent-based model: 1) no opening (control scenario), 2) 25% opening, 3) 50% opening, 4) 75% opening, and 5) 100% opening. These scenarios could be interpreted as cross-generational school re-opening schemes, i.e., all school levels from pre-school to post-graduate re-open at the same time. Alternatively, the school re-opening percentages could be interpreted as the resulting interactions after reducing the contact of students in schools, which could be done in many ways, e.g., reducing the class size, cutting short class hours, etc. The bottom line is that the school interactions are reduced to a certain percentage. The contact matrices were used to capture school re-opening scenarios.

Model Parameters and SEIRDV Values:
The SEIRDV values for the model were computed from various sources. The infectious (active cases), recovered, and dead were counted from the COVID-19 data released on 24 July 2021 by the Department of Health Philippines (2021).
Next, we estimated exposed individuals based on the number of suspected COVID-19 cases in Quezon City, Philippines, and the positivity rate. The age distribution of the cases was estimated by multiplying the number of suspected COVID-19 cases with the positivity rate. The product was then distributed proportionally to the infection distribution in the city.
For the vaccinated, the age distribution was estimated from vaccination data (Rappler, 2021) in Quezon City as of July 2021. First, we counted half of those who only got the first dose. Second, we added to the count the entire number of fully vaccinated individuals. Third, we multiplied the sum with the average vaccine efficacy (70%) (Table 2, Appendix). Finally, we distributed the result by age proportional to the age distribution of the population of Quezon City, Philippines.
Finally, for the susceptible, the age-stratified sum (vector sum) of the exposed, infectious, dead, recovered, and vaccinated was subtracted (vector subtraction) from the age-stratified population of Quezon City, Philippines (QC Government, 2018). The resulting SEIRDV distribution was then applied to a population of N = 8060 agents.
With the resulting SEIRDV values and the parameter values in Table 1, we simulated school re-opening and vaccination scenarios. The model ran 20 times for each scenario with 200 model steps on each iteration. The averages of the results were then computed, compared, and validated.

Infection Probabilities
Age-stratified infection probabilities

School Re-opening Sensitivity Analysis:
A sensitivity analysis on school re-opening in Quezon City, Philippines was performed using our COVID-19 agent-based model. Figure 3 shows the changes in infections and deaths with respect to the control scenario (no re-opening). In an abrupt 100% school reopening, COVID-19 infections are projected to increase by up to 83% and deaths by up to 48%. This may be a catastrophic event in the context of a pandemic, i.e., 83 additional infections and 48 additional deaths for every 100 cases. In a 25% re-opening, on one hand, infections and deaths are projected to increase by 22% and 16%, respectively. Even in this modest re-opening scenario, the increase in infections and deaths are still substantial, i.e., 22 additional infections and 16 additional deaths for every 100 cases. These projections are alarming given that we already assume that everyone (100%) is wearing a mask and 60% are observing physical distance.  Figure 3 also suggests that as we increase the school re-opening coverage, the infection and death curves diverge. This behavior suggests that not everyone who will be infected will die. However, contrary to our previous findings (Bongolan et al., 2021;Rayo et al., 2020;Minoza et al., 2020), where the younger ones recover and the older ones die, we are now projecting that deaths will increase on the younger age groups if we re-open the schools. Figure 4 shows the projected infection and death age distributions for Quezon City, Philippines on various school re-opening scenarios. As expected, infections will increase in the age groups where the students belong, i.e., 5 to 9, 10 to 14, 15 to 19, and 20 to 24. However, inferring from previous findings, we did not expect that deaths will also increase in these age groups. The death age distribution in Figure 4 shows that deaths will occur even in the youngest age group, 0 to 4, and is the highest among the other age groups. We validated this behavior from the COVID-19 data released on 24 October 2021 by the Department of Health Philippines (2021), and indeed, infant death is higher than the other younger age groups ( Figure 5).
We speculate that this observed and projected change in the death distribution of Quezon City, Philippines is due to the acquired immunity of individuals from recovery and/or vaccination, especially the older ones. The senior citizens became the priority of the vaccination program in the Philippines, and the low projected deaths on the older age groups (60+) may be attributed to the protection given by the vaccines. With these findings, we recommend not to re-open schools in Quezon City, Philippines with the control condition, i.e., 16% vaccination. Moreover, we recommend that the Philippine government now focus on the vaccination of the younger age groups, i.e., the students, before reopening the schools. This is to afford the students the protection that they need once they enter the schools. But how much vaccination coverage is needed before schools could re-open?

Vaccination Sensitivity Analysis
We performed a vaccination sensitivity analysis to see the effects of vaccination on the dynamics of COVID-19 transmission. The vaccination coverage in Quezon City, Philippines as of July 2021 is 16%. We made this vaccination coverage the control scenario for our analysis, combined with the no school re-opening scenario. Figure 6 shows the change in cases as we increase the vaccination coverage. At 25% vaccination coverage, the projected change in infections is -1.9%, a decrease, and 1.6% in deaths, an increase. The seemingly inconsistent increase in deaths may be due to computational noise since the model is fundamentally stochastic. This noise could be confirmed in the distribution of cases in Figure 7. The control scenario (16%) and 25% vaccination scenario are overlapping at some points. This behavior is negligible since the difference between the control scenario (16%) and the 25% vaccination scenario is small (9%). We will only observe a substantial decrease in cases as we further increase the vaccination coverage.
The highest decrease in cases, as expected, is at 100% vaccination coverage where the decrease in infections and deaths are almost half (51% and 48%, respectively). However, due to various reasons including hesitancy and the lack of vaccine supply, a 100% vaccination coverage may be unrealistic. Moreover, our simulations suggest that even in a 100% vaccination coverage, a zero-transmission scenario is impossible to achieve. This can be explained by the fact that the assumed average efficacy is only 70%. There is also an observed decrease in the effect of vaccination to infections as we increase the coverage: ~20% change from 25% coverage to 50%, ~15% from 50% coverage to 75%, and ~13% from 75% coverage to 100%. Hence, a 100% vaccination coverage should not be the goal of the policymakers because it simply is impractical and unrealistic.
What they could instead do is to look for the vaccination coverage where a downtick in cases could already be observed. A downtick in cases could be thought of as a car that is stopping. There will still be infections, there will still be deaths, but we are sure that the spread is already stopping. We think that the downtick points, the points at which infections or deaths start to decrease, are indicators of a safe school reopening. To find these points, we compared two school re-opening scenarios under various vaccination coverages. Figure 8 shows that in a 100% school re-opening with around 50% vaccination coverage, infections will begin to downtick. Moreover, with around 60% vaccination coverage, deaths will also begin to downtick. At these downtick points, the infections and deaths begin to plunge even when schools are 100% open. We interpret this range (50% -60%) as the minimum required vaccination coverage before re-opening the schools in Quezon City, Philippines.
However, our toy simulations shall be cautiously taken. The assumptions of our model shall be noted. To reiterate, we assume the following: the vaccination freezes at the start of the simulation, the average vaccine efficacy is 70%, school reopening simulations are cross-generational, and everyone is wearing a mask and 60% are observing physical distance.
As of 22 October 2021, Quezon City already fully vaccinated 55% of its entire population (Rappler, 2021). Hence, we think that schools could already re-open in Quezon City, but we recommend not to re-open the schools at 100% but only at 25% as a precaution. Also, we reiterate the need to vaccinate the students before re-opening schools to prevent avoidable deaths. Finally, mask-wearing and physical distancing shall be strictly observed once the schools re-opened.

CONCLUSION
In this study, we asked: Should we already re-open schools? We developed a COVID-19 agent-based model coupled with contact probabilities to explore this question. Our toy simulations suggest that the minimum required vaccination coverage for Quezon City, Philippines is 50%. With its 55% vaccination coverage, Quezon City, Philippines may already re-open schools. However, we suggest that students shall be vaccinated first, mask-wearing and physical distance shall be strictly observed, and schools shall only be re-opened by 25% as a precaution. Policymakers may take insights from the study. Average Vaccine Efficacy = (∑ * )/ ∑ = .