In December 2019, a novel coronavirus disease (COVID-19) emerged in Wuhan, China [1] , and has since spread out to the rest of the world, evolving into a pandemic [2] . Due to its high infection rate, the virus SARS-CoV-2 is not only causing mortality, but it is also stressing national health systems due to the large number of infected people that need hospitalization, also causing a profound economic impact [3] [4] [5] . In the absence of a vaccine or an effective treatment for COVID-19, the role of public health measures mostly comprehends the so-called non-pharmaceutical interventions (NPIs) [6] , including social distancing and isolation, reducing economic activities, enforcing home office, moving to remote classes and closing social venues where people tend to agglomerate in close proximity. Initially, sanitary authorities opted for strict NPIs. As the pandemic evolved, other factors emerged and started an interplay with the health crisis. Therefore, measures were revisited and, oftentimes, relaxed. However, a critical question is still unresolved: What is the impact of reopening the school system for in-person classes [7, 8] ?
When to start reopening schools and the reopening policy are particularly important questions for multiple reasons: i) education is a basic human necessity; ii) in-person classes enable other economic activities because parents can go to work while their children are at school [9] , besides other markets associated with schools, such as transportation, food, clothing, among others; iii) in-person classes are often hard to replace with online learning, especially in developing countries and for impoverished families, due to lower availability of broadband Internet and fast laptop computers with cameras; iv) due to the relatively long time children spend in proximity to each other in schools, the potential for exposure is much higher, requiring several workarounds to reduce the transmission rate [10] , especially because young people tend to be asymptomatic while carrying high viral loads [11] , which increases the chances of an exponen-tial transmission growth, thereby requiring additional quarantines to prevent a health system collapse; and v) a well-crafted policy for returning to school may help preventing new waves of the pandemic, save lives, and reduce the amount of psychological stress due to the health crisis [12] . To understand the impact of the pandemic on schools and develop policies to deal with it, mathematical and computer simulation modeling is crucial [13] .
In this work, we set out to evaluate a range of strategies for opening the school system of the São Paulo Metropolitan Area (SPMA) using computer simulation. We chose this region as a case study because Brazil is one of the current epicenters of the pandemic [14, 15] , and the SPMA is the most relevant of the country due its huge population and economical importance.
The region contains 39 cities and 21.7 million inhabitants, representing about 10% of Brazil's population [16] . We analyze 3 different strategies of school reopening: (i) reopen schools with all students at once; (ii) reopen schools following a strategy based on the official plan of the São Paulo government [17] , which consists of 3 stages, carefully increasing the amount of students in each stage; and (iii) reopen schools only when a vaccine becomes available. We show results regarding the total amount of COVID-19 cases, critical cases and number of deaths.
Most models to analyze pandemics employ compartmental models, such as SIS, SIR, SEIR or SEIHR [18] . However, the assumptions of homogeneous random mixing in those models can lead to wrong projections. Also, the effects of non-pharmaceutical interventions (NPIs) [6] are difficult to implement in those models. To overcome these issues, we developed a heterogeneous and dynamic network stochastic model [19] [20] [21] . The model was implemented as a simulator written in the C++ language, which was chosen due its high performance, low overhead and powerful object-oriented features. The network was modeled as an undirected graph G = (V, E), in which the vertices V represent the people and the edges E the relations between people, which have an infection rate. Given two people v 1 and v 2 , such that v 1 ∈ V and v 2 ∈ V , if we have an edge e ∈ E such that e = (v 1 , v 2 ), the person v 1 can infect the person v 2 . To support heterogeneous infection rates, we color the edges of the graph, such that we can attribute different infection rates for each color. The color represents the type of relation. Although it is possible to have any kind of relations modeled, for this work we modeled the following relation types: home, community, workplace, schools and Inter-city.
Since our model supports heterogeneous networks, it can provide a better representation of the population [22] . Each person can have an arbitrary number of edges in the graph, and the relation type (color) of each edge is also arbitrary. Therefore, as in real life, each person can have different infection rates depending on the number of relations they have, as well as the type of each relation. Since each person has a different infection rate, each person also has a different basic reproduction rate (R 0 ). Therefore, to set the overall R 0 , we calibrate the infection rates over the entire network such that the average R 0 among all people is equal to the target value. As our model supports dynamic networks, we can change relations at any time during the simulation. We are able to modify the infection rates, as well as add and remove edges from the graph. This is useful to implement sanitary interventions. We based our model in the SEIR compartmental model [18] , but extended it to support more compartments. Each person can be in one of the following states:
Susceptible The person can be infected.
Infected The person is currently infected. An infected person can have different sub-states, as we will explain later.
Immune The person is immune, either from recovering from the disease or from taking a vaccine. The person will not contract the disease anymore.
Dead The person died.
As for infected people, each one can be in the following sub-states:
Incubation The person was infected but has not developed any symptoms yet and is not contagious.
Unreported The person is contagious but either is asymptomatic or has such light symptoms that it does not receive any health care. This person, in real life, would not be part of any official statistics, hence a sub-notified case.
Pre-symptomatic The person is contagious, but did not develop any symptoms yet.
Mild The person has mild symptoms and can recover at home.
Severe The person has severe symptoms and requires health care in a hospital.
Critical The person has critical symptoms and requires health care in an Intensive Care Unit (ICU). After recovering from critical symptoms, the state is set to severe.
The disease spreads from infected people to susceptible people according to the infection rate of each edge in the graph, which is evaluated using a Monte Carlo based approach. The probability of a susceptible person p to get infected in a simulation cycle is given by Equations 1 to 4. Equation 3 defines the number of neighbors of each vertex (person) that are infected -N (p). In Equation 4, the function Inf ectRate(p, n) returns the infection rate between p and n in the network.
In the simulation, we consider all the cities and population that comprehend the SPMA, including the demographic distribution. Data used to set the simulation parameters were ex- 
In the first strategy, we evaluate what could happen if all students went back to school at once, but following sanitary measures, such as wearing masks and frequently washing their hands. Figure 10 , following the same organization of the previous strategy ( Figure 9 ). For this experiment, we consider an isolation level of 50%. Compared to when all students go to school, we can observe that the total number of reported infected people, in the worst case, reduces by 96,216. This translates to 4,437 fewer deaths, as well as reduction of the critical patients peak to 4,890 (after schools opened), which represents an ICU occupancy of only 97.8%. In the best case, the maximum ICU occupancy is 48.8% after schools open.
Results only considering the school sub-population follow the same trend of the general results. Considering the worst case, we can reduce the number of reported cases by 29,200
people and the number of deaths by 54 people compared to when all students go to school.
Regarding the deaths expected for the worst case, there would be a total of 541 deaths inside the school sub-population, in which 356 deaths correspond to teachers, representing 65.8% of school-related deaths. Although there are many more students than teachers, teachers are older than the students, and hence have a higher mortality rate. If students strictly follow the sanitary measures, as considered by the best case, the difference between opening the schools with all students or with São Paulo's plan is lower, but still very significant. Overall, results show that, although the difference in the number of cases is low compared to when all students come back at once, using an opening strategy with several phases avoids the health system collapse, which would have a big impact in the actual number of deaths.
In Figure 11 , we analyze several isolation levels for the São Paulo strategy. An isolation level of 100% means that students inside the same classroom have no interaction between each other, while an isolation level of 0% means that each student is able to infect any other student inside the same classroom. As in previous experiments, we evaluate two different rates. If students follow the social distancing and sanitary measures very strictly, with a 100% isolation, there would be up to 678,729 less reported cases than students have a 0% of isolation, reflecting in 32,817 fewer deaths. Although 0% and 100% isolation levels may be infeasible, there is no way to predict which isolation level would be achieved. If we compare an 80% isolation to a 20% isolation, we reduce the number of deaths by up to 16,844. While comparing a range from 60% isolation to a 40% isolation, we still reduce the number of deaths significantly, by up to 4,748. Even with a 0% isolation, although the amount of deaths is higher, the maximum number of required ICUs was 5,212 after schools opened, only a little higher than the capacity of SPMA. As expected, results show that student behavior is a key factor for the success of this strategy. The major problem is that there is no way to predict how students are going to behave.
The last strategy we evaluate, shown in Figure 
This research was partially supported by a machine allocation on the Kabré supercomputer at the Costa Rica National High Technology Center. We would like to thank Marco A. Amato for reviewing an earlier version of this manuscript.
This work received no direct funding.
Conflicts of Interest: None.
In case the paper gets accepted for publication, we will fully disclose all data underlying the study in a public repository. 
The supplementary Material includes:
Figs. 5 to 11 Tables 1 to 4 References Other Supplementary Materials for this manuscript include the following:
Repository address of the simulator source code and input data: Will be disclosed upon paper approval.
The initial network was created in 5 stages. Note that the initial network considers a normal scenario with no pandemic and interventions. We explain how we handle the interventions in Section 1.2. The initial network was created as follows:
Home relations The number of people in the same home follows a distribution described in Table 1 and Figure 5a , which was derived from the microdata of the last official census from the Brazilian Institute of Geography and Statistics (IBGE) [24] . Inside the same home, every person is connected to all other people within the same home.
Community relations The number of relations of each person follows a distribution described in Table 1 and Figure 5b .
Workplace relations The number of people per company follows a distribution described in Table 1 and Figure 5c . We only model as workplace relations jobs that involve agglomeration of people, such as people that work in the commerce, offices, industry and public sectors. Other jobs, such as people that are autonomous workers or people that work with agriculture, are indirectly modeled in the community relations. The data regarding the number of employed people and per job type was extracted from IBGE [24] [25] [26] . In summary, 64% of people with 18 years old or above work. Among people that work, 45%
were modeled using workplace relations. Also, inside the same workplace, we consider that the probability of two people to be connected in the network as 50%.
School relations The number of students inside the same class and the number of students within each school follow the distributions described in Table 1 , Figure 5d and Figure 5e .
For the relations between students from different classrooms, two students have a 0.5% of probability of having a relation in the network. For the relations between students within the same classroom, two students have a 50% of probability of having relation in the network. We add one teacher for each classroom, in which the teacher has a relation to all the students of the same class in the network. According to the census [24], 77.8%
of people with 18 years old or below frequent the school.
Inter-city relations The number of people that routinely travel between the different cities relations was gathered from the census [24], which indicates a 17.8% of mobility between cities in the metropolitan region of São Paulo. However, the data does not show the cities, only if the person goes to a different city. To overcome this issue, we assume that 90% of the traffic from satellite cities go to the capital, and the rest is divided between the other cities in proportion to their population size. Regarding the traffic from the capital to satellite cities, the traffic is divided proportionally to their population size.
To handle interventions, we can modify the infection rates, as well as add and remove edges from the graph. Different types of intervention require different modifications. We handled the interventions evaluated in the paper as follows:
Quarantine To implement the quarantine, we first reduce the infection rates inside schools to zero. Afterwards, we remove the workplace relations corresponding to 50% of the workplaces. The 50% of the relations that are left correspond to essential services that can not be interrupted. Finally, we reduce the infection rate of community, workplace and inter-city relations to match the reduced reproduction rate, which also simulates other sanitary measures such as wearing masks. Home relations are not affected by quarantines. Re-opening schools with all students We only increase the infection rate of school relations, which was zeroed by the quarantine.
Re-opening schools with the São Paulo government strategy We first remove all school relations from the graph, backing up the data regarding which students and teachers belonged to the same classroom. We create 3 more relation types to represent each group of students. Afterwards, we divide each classroom in 3 and re-create the edges using the newly created relation types. To select which students go to school in a day, we just need to adjust the infection rates of the new relations, zeroing the relations corresponding to the students that will not go to school that day.
Vaccine To simulate the vaccines, we only need to change the state of the person being vaccinated to the immune state, carefully considering that it takes some time to the vaccine to make effect, and that not everyone that gets vaccinated will be immunized.
Several simulation parameters can be found in Table 2 [4, 16, 24-28]. The frequency distributions found in the table can be seen in Figure 6 .
In Table 3 , we present the probabilities of an infected person to develop symptoms per-age.
The biggest challenge to generate this information is that there is no official data regarding the number of unreported cases per age, as expected by the own definition of unreported. Therefore, we need to make several transformations with the reported data to take into account the unreported cases in the simulation. We now describe how we generated these information of Table 3 .
We first define some functions and variables that return the input data of the algorithm:
• Function people(age) returns the number of people of age age in the entire population, which we extracted from the census [16].
• Functions reportedM ild(age), reportedSevere(age) and reportedCritical(age) should return the number of people with mild, severe and critical symptoms per age, respectively.
We extracted this information from the microdata published by the Health Ministry [28].
• Variable globalRatioU nreported is the ratio of unreported cases in the entire population. We set it to 0.85 in this paper, since the vast majority of COVID infections are unreported [29].
Next, we define some variable and functions to be used later.
reportedM ild(age)
reportedSevere(age)
reportedCritical(age)
reported(age) = reportedM ild(age) + reportedSevere(age) + reportedCritical(age)
reported(age)
To calculate the probability of a person to be unreported, we first calculate the ratio of reported cases:
ratioReportedM ild(age) = reportedM ild(age) reported(age)
ratioReportedSevere(age) = reportedSevere(age) reported(age)
ratioReportedCritical(age) = reportedCritical(age) reported(age)
Then, we use Equation 14, which finds a value u that forces the mean number of unreported people to the desired value. In the equation, weight(age) returns the proportion of people of age age that should be unreported. Since there is no data for unreported (as expected by the own definition of unreported), we assume that weight has the same ratio of mild cases per age, as it is the softer symptom recorded by official data.
Find u, such that:
u · 9 age=0 weight(age) · people(age) 9 age=0 people(age)
After finding u, we calculate a temporary ratio of unreported cases by age, just used to calculate the symptomatic ratios:
tmpRatioU nreported(age) = weight(age) · u
Now that we know the ratio of unreported cases by age, we can calculate the relative proportion (weight) of symptomatic people:
weightM ild(age) = ratioReportedM ild(age) · (1 − tmpRatioU nreported(age)) (16)
weightSevere(age) = ratioReportedSevere(age) · (1 − tmpRatioU nreported(age)) (17) weightCritical(age) = ratioCriticalM ild(age) · (1 − tmpRatioU nreported(age)) (18)
Now we define our target averages for the ratio of mild, severe and critical cases relative to the entire population (not only relative to reported cases).
Then, we use Equations 22, 23 and 24, which find values m, s and c that force the mean ratio of people to the target value.
Find m, such that:
age=0 weightM ild(age) · people(age) 9 age=0 people(age)
Find s, such that:
s · 9 age=0 weightSevere(age) · people(age) 9 age=0 people(age)
Find c, such that:
age=0 weightCritical(age) · people(age) we zero the infection rates to home relations. Table 4 shows the death probability of each person depending on the symptoms and age.
As we can observe, the age of each person has a deep impact on the symptoms and death rate. Due to that, we use real demographic data from the census Vaccine -time to be effective 14 days. Vaccine -number of people vaccinated per day 300,000 people per day. Vaccine -immunity rate 90%.
Infection rate -home ×3.0 the community rate considering the beginning of the pandemic (not changed by any intervention). Infection rate -workplace ×1.5 the community rate. Infection rate -school
We evaluated both ×2 and ×4 the community rate. Infection rate -inter-city Equal to the community rate. Infection rate -community
We calculate it such that the average reproduction rate among all people matches the target reproduction rate. Table 4 : Probability of a person to die depending on their symptoms and age.