Infectious Disease Modeling with Socio-Viral Behavioral Aspects—Lessons Learned from the Spread of SARS-CoV-2 in a University

When it comes to understanding the spread of COVID-19, recent studies have shown that pathogens can be transmitted in two ways: direct contact and airborne pathogens. While the former is strongly related to the distancing behavior of people in society, the latter are associated with the length of the period in which the airborne pathogens remain active. Considering those facts, we constructed a compartmental model with a time-dependent transmission rate that incorporates the two sources of infection. This paper provides an analytical and numerical study of the model that validates trivial insights related to disease spread in a responsive society. As a case study, we applied the model to the COVID-19 spread data from a university environment, namely, the Institut Teknologi Bandung, Indonesia, during its early reopening stage, with a constant number of students. The results show a significant fit between the rendered model and the recorded cases of infections. The extrapolated trajectories indicate the resurgence of cases as students’ interaction distance approaches its natural level. The assessment of several strategies is undertaken in this study in order to assist with the school reopening process.


Introduction
In epidemiology, compartmental models are general modeling techniques used to understand the spread of disease, and they commonly consider three variables: S for those who are susceptible, I for those who are infected, and R for individuals who have recovered. Variations of the generic SIR model are available: the SIS model accommodates temporal immunity [1], the SEIR model best represents the spread of disease with a significant latency period [2], and there are even combinations of the two [3]. The convenience of compartmental models in respect of adding more variables has resulted in their being widely used in infectious disease modeling [4]. Besides providing each state's estimated figure, this approach can also provide the reproductive ratio, which represents the expected number of secondary cases generated by one primary case [5][6][7]. In most of the constructed models, the reproductive ratio acts as a crucial threshold; above one indicates endemic, while below one indicates disease-free [8]. This is crucial for policymakers when regulating whether or not to ease restrictions amid disease spread.
However, generic compartmental models are sometimes based on assumptions that are not necessarily relevant; the population is considered closed in SIR models, whereas complete isolation was not followed in most regions, making them vulnerable to changes in the neighboring communities [9]. Another assumption that is commonly used in a generic model is that transmission and recovery rates are assumed to remain constant over time. Such a scenario will best represent disease that spreads in a population with no response to current disease prevalence, meaning that a high or low number of recorded cases will not affect the average socio-behavior of the population. The simplest case to consider is a disease spread within a closed population of sheep in a field [10]. When it comes to a human population, people's psychological behavior causes them to reduce their interaction intensity as the declared number of cases increases, which ought to vary the viral transmissibility [11]. Moreover, setting constant rates of transmission and recovery results in a high number of projected infected cases once the model is applied to a vast and highly populated community; this could be at the scale of entire nations [12]. According to recent studies, SIR-based predictions using early data for COVID-19 cases have shown an enormous figure for predicted cases, with the peak reaching up to 15-30% of the total population [13,14]. Nevertheless, an absence of the psychological behavior of the population could overestimate the prediction figure [15].
According to recent studies, there are so many studies that discuss the spread of the COVID-19 disease. Researchers developed various models and approaches from all over the world [16,17]. However, in this paper, we will discuss two major sources of transmission in some infectious diseases: direct contact and airborne transmission. In respect of the former, it is quite obvious that human-to-human transmission is mainly caused by direct contact such as talking at a close distance. The smaller the average interaction distance of people within a population, the greater the chance for pathogens to spread. By incorporating the effect of human psychological behaviors, it is natural to expect an increase in the average interaction distance given a high disease prevalence in a specific population, which will lead to a reduction in viral transmissibility. However, the latter source of transmission opens up possibilities for infections induced by the presence of airborne pathogens. This method of transmission is found in the spread of TB [18] and SARS-CoV-2 [19]. Although airborne pathogens can infect susceptible individuals, some studies have shown that most airborne pathogens can only last for a certain period. Mycobacterium tuberculosis, which attacks lungs and causes TB, can stay in the air for several hours depending on the environment [20], and SARS-CoV-2 can only last for hours in the air but can survive for up to a week on plastic [21]. In disease modeling, taking airborne pathogens into account is crucial, especially for those that have a significant period of viral survivability in the air.
The incorporation of the psychological behavior of society into responses to disease prevalence has been introduced in several works, such as Hua-Li et al. [11] and Oluyori et al. [22]. In practice, the authors define saturated transmission rates that are dependent on the figure of disease prevalence. The transmission rate is expected to increase for a low disease prevalence and start decreasing once the prevalence exceeds its critical point [11]. In 2021, Cabrera et al. [23] introduced a compartmental model that incorporates a sociobehavioral aspect in a slightly different way; they introduced the interaction distance to measure societal behaviors in response to disease prevalence. Hence, the nonlinear transmission rate integrates the interaction distances. However, the effect of airborne pathogens is rarely incorporated. One study conducted by Bazant and Bush in 2021 [24] demonstrates the significant effect of airborne transmissions on society regarding activities. Although airborne pathogens, especially SARS-CoV-2, can only last for hours, indoor transmission is crucial for infectious disease modeling, especially for school or office environments involving many indoor activities.
In this study, we constructed an SIR-based mathematical system that accommodates the two major causes of infection: direct contact and airborne transmission. The former source of infection, representing the socio-behavioral aspect, is based on the measure of the interaction distance of people in society. In 2021, Cabrera et al. [23] proposed adding a new variable that determines the interaction distance over time. The closer the interaction distance, the higher the chance of disease spread. The latter cause of infection, which represents viral characteristics, is incorporated by defining another variable that solely represents the concentration of pathogens in the air over time. We expect that the longer the pathogens can last in the air, the higher the concentrations over time, which leads to a higher chance of disease spread. Hence, the newly added variables will govern the transmission rate that eventually depends on the socio-viral behavioral aspects. In the analysis of the constructed model, we provide numerical results in respect of infections under different socio-viral behavioral aspects. The model performs well in depicting the spread of disease in societies under different rates of response, different rates of resistance to adopting new habits, and under different characteristics of the concerned diseases. As a case study, we applied the constructed models to the SARS-CoV-2 spread data that were collected in a university environment (Institut Teknologi Bandung College) in January 2022. The choice to use data from a university was made to ensure homogeneous sociobehavioral aspects for the whole society; no demographic is taken into account due to the homogeneity assumption [25][26][27]. The small scale of a university environment also ensures the involvement of pathogens in the air; the larger the scale of the observation, the smaller the effect of pathogens in the air. Lastly, we utilized the extrapolated figures to assess some strategic action plans related to SARS-CoV-2 infections in educational environments; school reopening schemes and vaccination implementation [28].

Context
Humans are mobile creatures who move in their part of an environment; they may meet an acquaintance or not. When the former scenario happens, they will likely move closer to reaching out to that acquaintance [29]. This phenomenon exemplifies the importance of interpersonal space (IPS) and peripersonal space (PPS) in which humans can perform body-environment interactions [30]. Although the dimensions of IPS and PPS include all directions, previous studies have only focused on a specific distance, i.e., the distance from the front of the person [31]. When it comes to understanding infectious diseases, the front-directed PPS is essential since most diseases, including SARS-CoV-2, are transmitted via the front parts of the human body. One unit that measures the intensity of PPS contact is the interaction distance, in which the closer the distance, the more intense the contact, which leads to an increase in the risk of disease transmission [32]. According to Sorokowska et al. [33], the preferred interpersonal distance of humans differs between different types of social relations (strangers, acquaintances, and partners). Table 1 provides a global comparison in respect of interaction distance.
Other than the interaction distance that causes direct transmissions, airborne transmission of some diseases is now widely recognized, especially for the spread of COVID-19 [34,35]. This approach accounts for the plausibility of infections caused by pathogenbearing aerosols that are fine enough to be continuously mixed through an indoor space. Every infected individual present will contribute to the production of droplets containing the virus. Bazant and Bush [24], in their COVID-19 study, estimated the concentration of pathogens produced by a single infected individual in a well-mixed room for every breath, and for whispering and talking indoors.
However, other studies have shown that pathogens can remain active on other media, such as copper, cardboard, and plastic [21], for a certain period. Hence, other than significant airborne transmission indoors, pathogens that are attached to other media can also infect susceptible individuals. A study by Doremalen et al. provides the estimated critical periods of SARS-CoV (1 and 2) before they become inactive; these are given in Table 2. The estimations show that SARS-CoV can last up to 12 h in the air but can last longer on other media. This fact should indicate the importance of airborne pathogens and their attachment to other media in respect of understanding viral transmission. In this study we construct a mathematical model that incorporates both socio-behavioral and airborne pathogen effects.

Proposed Model
In this study we used a generic model, but we separated those who had and had not received vaccines. This modification was based on the fact that the presence of immune titer in the human body can significantly prevent people from becoming infected, offering up to 90% protection [36]. Hence, there are three main state variables: susceptible (S), currently infected individuals (I), and removed individuals (R), with the total of six state variables created by adding subscripts v and u to each of the main states, representing the categories of being vaccinated and not, respectively. As shown in Figure 1, new infected individuals are generated from both S u and S v , caused by a direct interaction between susceptible and infectious individuals. After a specific period of infections, infected individuals will enter R, which represents being immune or deceased. We assume that there is no demographic change, which implies a constant population size: , for t ≥ 0, and N = N u + N v with a constant proportion of vaccinated and unvaccinated individuals. The model also assumes no significant difference in the recovery rates of vaccinated and unvaccinated individuals. individuals will enter , which represents being immune or deceased. We assume that there is no demographic change, which implies a constant population size: = ( ) + ( ) + ( ) and = ( ) + ( ) + ( ), for ≥ 0, and = + with a constant proportion of vaccinated and unvaccinated individuals. The model also assumes no significant difference in the recovery rates of vaccinated and unvaccinated individuals. As shown in Figure 1, there are three parameters involved: transmission rate ( ), recovery rate ( ), and vaccine effectiveness ( ). The last two parameters are observable, i.e., their values can be measured and estimated using relevant information. Vaccine effectiveness, which ranges from 0 to 1, represents the protection induced by the vaccine. The higher this value, the lower the chance of people becoming infected once they interact with infectious individuals. Limited to the COVID-19 vaccine, the vaccine efficacy should vary depending on the manufacturer and COVID-19 variants [36]. The value of represents the rate of recovery, which governs the speed of transition from to . To make this realistic, can be considered as the average infection period. In contrast, the rate of transmission is unclear in terms of its physical representation; it summarizes all factors that produce infections. Hence, the value of is considered unobservable. To incorporate the two major causes of infection as mentioned in Section 2, we added two additional lines to the system that represent the dynamics of the interaction distance and the pathogen concentrations . The final two variables dictate the dynamics of resulting in the transmission rate that depends on the socio-viral behavioral aspect. A mathematical representation of the constructed model is given in the following form: with a constant population size . The other two additional variables are (in meters) and (in quanta/m ), representing the average interaction distance and viral loads over time. The formulation of was first introduced by Cabrera et al. [23] along with the definition of the natural distancing habit * that could differ from one society to otherssymbol * denotes the average of natural distancing behavior of society. The complete additional lines are given in the following systems:  As shown in Figure 1, there are three parameters involved: transmission rate (β), recovery rate (γ), and vaccine effectiveness (ρ). The last two parameters are observable, i.e., their values can be measured and estimated using relevant information. Vaccine effectiveness, which ranges from 0 to 1, represents the protection induced by the vaccine. The higher this value, the lower the chance of people becoming infected once they interact with infectious individuals. Limited to the COVID-19 vaccine, the vaccine efficacy should vary depending on the manufacturer and COVID-19 variants [36]. The value of γ represents the rate of recovery, which governs the speed of transition from I to R. To make this realistic, γ −1 can be considered as the average infection period. In contrast, the rate of transmission β is unclear in terms of its physical representation; it summarizes all factors that produce infections. Hence, the value of β is considered unobservable. To incorporate the two major causes of infection as mentioned in Section 2, we added two additional lines to the system that represent the dynamics of the interaction distance D and the pathogen concentrations V. The final two variables dictate the dynamics of β resulting in the transmission rate that depends on the socio-viral behavioral aspect. A mathematical representation of the constructed model is given in the following form: with a constant population size N. The other two additional variables are D (in meters) and V (in quanta/m 3 ), representing the average interaction distance and viral loads over time. The formulation of D was first introduced by Cabrera et al. [23] along with the definition of the natural distancing habit D * that could differ from one society to others-symbol D * denotes the average of natural distancing behavior of society. The complete additional lines are given in the following systems: with non-negative initial conditions S 0 , V 0 that are evaluated at the initial point t = 0. It is natural to assume that I 0 u = (1 − α)I 0 and I 0 v = αI 0 , for I 0 = I 0 u + I 0 v , with α (in percent relative to the population size) representing the vaccine coverage. The addition of the two variables involves another four parameters. On one hand, the value of λ 2 (distance/time) represents how quickly people react to the current disease prevalence, i.e., the so-called rate of social response. By neglecting the first term, there are two scenarios that increase the interaction distance D: high values of the rate of response λ 2 or the disease prevalence I. Interestingly, setting λ 2 equals zero will lead to a situation where a society pays no attention to the current disease spread. Such a scenario drives the society to resort to their natural interaction distance D for λ 1 = 0. On the other hand, the rate λ 1 (1/time) measures the rate of resistance in society, per distance unit, to changing distancing behavior. It represents how quickly individuals return to their natural interaction distance D * or their natural distancing habits. This rate is strongly related to the distancing culture. When we set a high value of λ 1 , this results in a situation in which the society has a strong culture embedded, making it resistant to changes in behavior amid the current pandemic. In this study, we restrict the plausibility of λ 1 = 0 since we assume that every society has its own resistance in changing habits. When the disease prevalence approaches zero, then D approaches −λ 1 (D − D * ), which leads to the convergence of D to D * regardless of the initial condition D 0 . More detailed formal analysis of System (1) and (2) are given in Appendix A.
While the first equation of System (2) portrays the socio-behavioral aspects, the second equation portrays the concentration of the pathogens. The rate λ 3 (quanta/(time m 3 ·person)) denotes the average concentration of viral/pathogens emitted by one infected individual per unit time. Face coverings and the practice of other social and respiratory etiquette will likely reduce the value of λ 3 and hence reduce the number of pathogens emitted into the air. The discharged microbes will remain suspended in the air in dust particles, respiratory particles, and water droplets [37]. However, pathogens will not last forever in the air (or other media); they will decay due to natural and human intervention. On the other hand, parameter λ 4 (1/time) denotes the removal rate of viral quanta in the air. A higher intervention of humans in the community, including through air filtering and periodical sanitation, can increase λ 4 and hence allow more microbes to decay or be inactive [38]. However, in most cases, λ 4 will only account for the natural effect of pathogen removal (subject to ambient temperature, humidity [21,39], and sunlight [40]), while human intervention can be represented by another functional term added to the dynamic of V [41]. Eventually, λ 1 , λ 2, and λ 3 represent the socio-behavioral aspects in society while λ 4 represents the characteristics of the pathogens.
Since the model adopts a uni-flow, then there exists τ such as I(t) < ε, t > τ, for every ε > 0. In terms of epidemiology, the virus will always be eradicated to zero for large values of t since people will accumulate in the removed compartments. For the dynamics of D, the second and third terms approach zero as t approaches infinity, leaving only the first term that converges to D * . However, the presence of V is strongly related to the presence of infectious individuals, who will vanish once the disease vanishes, no matter how large the initial condition. It should be noted that the proposed models do not consider reinfection or susceptible newborns. Hence, multiple disease outbreaks (if any) are expected to be driven by the change in interaction distance in society.

Observability of Socio-Behavioral Parameters
As discussed in the previous section, the model has 3 parameters that are related to the socio-behavioral aspects of society: λ 1 , λ 2 , and λ 3 . It is clear from its definition that λ 3 is observable and that its value follows the estimations of the pathogen concentration per person per m 3 . Bazant and Bush [24] and Miller et al. [42] provided estimated concentrations for several expiratory activities. Calibrated normal speaking activity is estimated to produce 72 infections quanta/m 3 while superspreading activity can contribute up to 970 infections quanta/m 3 . However, the first two socio-behavioral parameters are not observable, i.e., the rate of social resistance λ 1 is not something that we can determine from the field. It combines all aspects that inhibit society in the change of behaviors.
The rate of social response, denoted by λ 2 , has a dimension of meters per unit of time. In the absence of λ 1 , the formula of D reduces to only D = λ 2 (I/N), with I = I u + I v . When I = N, then D = λ 2 , which is interpreted as the interaction distance increasing at the rate of λ 2 meters per unit time when the whole population is infected. Taking another scenario, I = 1 person results in D = λ 2 /N, which is considered as the λ 2 /N increment of the interaction distance per unit of time when the society contains 1 infected individual. Henceforth, λ 2 is related to the quantity of the change in D for a certain disease prevalence. To understand this parameter more, let us take the solution of Expecting the presence of an average prevalence of I in the length of time T 2 will drive people in society to interact at the distance of D, then λ 2 can be estimated using the following formula: Note that I/N represents the percentage of infections in society, i.e., the so-called point prevalence, denoted by a%. Therefore, by knowing that the society is practicing distancing habits of D = D once the point prevalence is roughly a%, we can estimate the expected value of λ 2 as the rate of social response amid the disease spread. Henceforth, We can also consider the dynamics of D in the given system. When (I u + I v )/N tends to zero, the effect of λ 2 is no longer significant; the whole second term will tend to zero, leav- The higher the value of λ 1 , the faster the dynamics of D to approach D * . It is easy to prove that lim t→∞ D(t) = D * , regardless of the value of D 0 . Hence, for an arbitrary small value ε > 0, there exists a value of T 1 that satisfies the following condition.
We can manipulate the solution of D(t) to reach D * + ε in t = T 1 by adjusting the value of λ 1 as given by: Henceforth, the rate of social resistance λ 1 can be evaluated using the estimated time for society to return to their natural interaction distance in the absence of disease spread, denoted by T 1 ; see Figure 2 for illustration. It should be noted that ε is an arbitrary small number ε > 0 divided by D 0 . According to Equation (7), λ 1 takes the log value of ε, which will be sensitive to the choice of ε. Hence, it is natural to assume the relative deviation from D * as ε = 1%, although the formula of λ 1 should clearly confirm that the value of λ 1 is dependent on the assumption.  Figure 2. Illustration of rate of social resistance by the given data for .

Contact and Airborne-Based Transmission Rate
The rate of transmission is defined to be related to the interaction distance ( ) an concentration of pathogens ( ). In this study we accommodate two methods of transmi sion: contact-based and airborne-based transmission. Contact-based transmission is a fected by the average interaction distance; the transmission rate decreases as the averag interaction-distance increases, as people practice social-distancing. However, a high co centration of airborne pathogens contributes to an increase in the transmission rate.
( , ) = * 2 * * + The definition of the contact-related transmission rate is adopted from [23], but w have added the effect of the current concentration of pathogens. The basic transmissio rate * (1/time) is defined as constant, representing the basic probability of transmissio per unit time. The second term (dimensionless) represents the effect of the average inte action-distance, which will decrease the overall as increases. The third term, how ever, represents how the concentration of pathogens affects the overall value at whic the risk of infection will rise as increases. To keep the effect of dimensionless, w divide by the standard number of quanta exhaled by infectors per individual per per unit time. The adjuster levels of and are added to be fitted to the data, represen ing the strength of each source of infection in society.

Recovery Rate
Recovery rates (1/time) denote the quantity representing how fast infected indivi uals recover from the disease and, hence, build their immunity [43]. For some infectiou diseases, the absence of healthcare might cause a longer infection period [44][45]27], sp cifically for COVID-19. Not limited to this disease, we define the implicitly time-depen ent recovery rate as follows: where denotes the state variable for infectious individuals and denotes the consta healthcare capacity (beds). In addition, and are both recovery rates but represe different situations: excessive beds and collapsing health systems. On one hand, when th number of beds is excessive, then each of the infected individuals receives proper trea ment and this leads to a shorter period of infections [45]. In other words, ( , ) w

Contact and Airborne-Based Transmission Rate
The rate of transmission is defined to be related to the interaction distance (D) and concentration of pathogens (V). In this study we accommodate two methods of transmission: contact-based and airborne-based transmission. Contact-based transmission is affected by the average interaction distance; the transmission rate decreases as the average interactiondistance increases, as people practice social-distancing. However, a high concentration of airborne pathogens contributes to an increase in the transmission rate.
The definition of the contact-related transmission rate is adopted from [23], but we have added the effect of the current concentration of pathogens. The basic transmission rate β * (1/time) is defined as constant, representing the basic probability of transmission per unit time. The second term (dimensionless) represents the effect of the average interactiondistance, which will decrease the overall β as D increases. The third term, however, represents how the concentration of pathogens affects the overall β value at which the risk of infection will rise as V increases. To keep the effect of V dimensionless, we divide V by the standard number of quanta exhaled by infectors per individual per m 3 per unit time. The adjuster levels of v and w are added to be fitted to the data, representing the strength of each source of infection in society.

Recovery Rate
Recovery rates (1/time) denote the quantity representing how fast infected individuals recover from the disease and, hence, build their immunity [43]. For some infectious diseases, the absence of healthcare might cause a longer infection period [27,44,45], specifically for COVID-19. Not limited to this disease, we define the implicitly time-dependent recovery rate as follows: where I denotes the state variable for infectious individuals and K denotes the constant healthcare capacity (beds). In addition, γ 1 and γ 0 are both recovery rates but represent different situations: excessive beds and collapsing health systems. On one hand, when the number of beds is excessive, then each of the infected individuals receives proper treatment and this leads to a shorter period of infections [45]. In other words, γ(I, K) will achieve its maximum rate of recovery as K approaches infinity. Otherwise, γ(I, K) will converge to γ 0 as the number of burdens is higher relative to the healthcare services [11]. Hence, the former denotes the maximum recovery rate given the proper treatment, while the latter denotes the lowest recovery rate achieved by patients treating themselves in order to recover. Figure 3 depicts the functional parameters and their dependent variables. Figure 3 (left) illustrates the effect of the average interaction distance that results in higher values of β(D, V) as D approaches D * . On the other hand, the rate of recovery follows Equation (9), which lessens the rate of the increase in the burden of cases down to γ 0 . For the case with excessive healthcare capacities, the rate of recovery can be maximized up to γ 1 , as shown in Figure 3 (right).
Trop. Med. Infect. Dis. 2022, 7, x FOR PEER REVIEW 9 of 23 achieve its maximum rate of recovery as approaches infinity. Otherwise, ( , ) will converge to as the number of burdens is higher relative to the healthcare services [11]. Hence, the former denotes the maximum recovery rate given the proper treatment, while the latter denotes the lowest recovery rate achieved by patients treating themselves in order to recover. Figure 3 depicts the functional parameters and their dependent variables. Figure 3 (left) illustrates the effect of the average interaction distance that results in higher values of ( , ) as approaches * . On the other hand, the rate of recovery follows Equation (9), which lessens the rate of the increase in the burden of cases down to . For the case with excessive healthcare capacities, the rate of recovery can be maximized up to , as shown in Figure 3 (right).

Numerical Results
In this study, the behavior of society that is being accommodated by the model is the rate of social resistance and social response. Socio-resistance rate, denoted by , represents the resistance of society to distancing their interactions due to the prevalence of people when is not significantly zero. When the prevalence of people is close to zero, the resistance rate represents how fast the society moves back towards their natural interaction-distance * . In contrast, the rate of societal response represents the increase in interaction distancing per increase in point prevalence, which inhibits the disease spread when this value is high. In this section, we provide the number of infected individuals (per one thousand members of the population) for several values of and .

Variations under Different Society Behaviors
The rate of social response is given in three scenarios (low, moderate, and high response), by taking values of = 0.20, 0.53, and 0.87, respectively. These are based on the physical distancing campaign: (i) low social response drives people to physical distancing limited to = 1 m only, (ii) moderate can reach = 1.5 m, and (iii) high social response can reach up to 2 m. Table 3 shows the diverse approaches of countries in campaigning for physical distancing. We also set the rates of social resistance to = 0.15, 0.07, and 0.05, which are based on = 30, 60, and 90 days, respectively. The ranges of and produced by Formulas (5) and (7) conform to those used in Cabrera et al. [23].

Numerical Results
In this study, the behavior of society that is being accommodated by the model is the rate of social resistance and social response. Socio-resistance rate, denoted by λ 1 , represents the resistance of society to distancing their interactions due to the prevalence of people when I is not significantly zero. When the prevalence of people is close to zero, the resistance rate represents how fast the society moves back towards their natural interaction-distance D * . In contrast, the rate of societal response represents the increase in interaction distancing per increase in point prevalence, which inhibits the disease spread when this value is high. In this section, we provide the number of infected individuals (per one thousand members of the population) for several values of λ 1 and λ 2 .

Variations under Different Society Behaviors
The rate of social response is given in three scenarios (low, moderate, and high response), by taking values of λ 2 = 0.20, 0.53, and 0.87, respectively. These are based on the physical distancing campaign: (i) low social response drives people to physical distancing limited to D = 1 m only, (ii) moderate can reach D = 1.5 m, and (iii) high social response can reach up to 2 m. Table 3 shows the diverse approaches of countries in campaigning for physical distancing. We also set the rates of social resistance to λ 1 = 0.15, 0.07, and 0.05, which are based on T 1 = 30, 60, and 90 days, respectively. The ranges of λ 1 and λ 2 produced by Formulas (5) and (7) conform to those used in Cabrera et al. [23].  Figure 4 provides the numerical simulations for I u + I v and D for the different pairwise scenarios of λ 1 and λ 2 . As expected, the value of D(t) will vary over time-increases as the disease prevalence increases. In all sub-figures, all dynamics for D(t) always start from D * as its natural distancing behavior when disease prevalence is around zero (no new cases recorded). However, as the disease prevalence rises, people in society build awareness to practice physical distancing which then increases the average distancing behavior D. As the new cases decrease to zero, it is natural that people in society return to their natural distancing D * .   Figure 4 provides the numerical simulations for + and for the different pairwise scenarios of and . As expected, the value of ( ) will vary over time-increases as the disease prevalence increases. In all sub-figures, all dynamics for ( ) always start from * as its natural distancing behavior when disease prevalence is around zero (no new cases recorded). However, as the disease prevalence rises, people in society build awareness to practice physical distancing which then increases the average distancing behavior . As the new cases decrease to zero, it is natural that people in society return to their natural distancing * . The sub-figure in the left upper corner depicts the simulation results for a society with a lower response yet a higher resistance rate. Such a scenario results in a higher peak of the burden of cases relative to other scenarios. This result shows that if the society does not have enough awareness about the disease's prevalence, and has a strong resistance that inhibits the practice of physical distancing, the dynamics of will be likely in around * , which results in a higher number of cases. Societies that campaign for close physical distancing (e.g., 1 m only), and have tendencies to always practice their natural habits, are The sub-figure in the left upper corner depicts the simulation results for a society with a lower response yet a higher resistance rate. Such a scenario results in a higher peak of the burden of cases relative to other scenarios. This result shows that if the society does not have enough awareness about the disease's prevalence, and has a strong resistance that inhibits the practice of physical distancing, the dynamics of D will be likely in around D * , which results in a higher number of cases. Societies that campaign for close physical distancing (e.g., 1 m only), and have tendencies to always practice their natural habits, are likely represented by the left upper corner sub-figure. The figure situated at the center of the nine depicts simulations that apply to a society that has a considerably moderate level of resistance and response rate. The right lower corner depicts shows societies with a higher rate of response but a lower rate of resistance. Due to higher values of λ 2 , every individual in the society moves further away relative to other scenarios and this results in a significant change in D relative to the value of D * . When it comes to the figure of the burden of cases, this scenario estimates the lowest number of cases relative to other scenarios. Societies that practice physical distancing and have a tendency to keep practicing it in a longer period, even after the disease is no longer present, are best represented by this scenario, resulting in a lower burden of cases relative to other scenarios.

Variations under Different Pathogen Characteristics
Different pathogens lead to different survivability periods in the air or other media. The longer the pathogens are active as airborne pathogens, the more they accumulate, which increases the risk of infections. Characteristics of the observed pathogens are governed by parameter λ 4 . In System (2), the term −λ 4 V represents the concentration of pathogens per unit of time to become inactive. Figure 5 shows the dynamics of the disease prevalence I u + I v under different periods of pathogens lasting in the air for the same parameters as used in Figure 4. In the lower-right picture, it is shown that pathogens that can last up to 48 h (red) can accumulate up to 300,000 quanta pathogens per m 3 and drive infections to as high as 23%. Figure 5a demonstrates how the dynamics of I u + I v precede V on reaching a peak for exactly 2 days (48 h). It is natural to accept that the longer the period, the wider the gap between the occurrences of the two peaks. By setting a smaller period (higher λ 4 ), the dynamics of V decrease and so does I u + I v . Moreover, the peak of I u + I v shifts to the right (see Figure 5b,c). More results on the model's sensitivity analysis are provided in Appendices B and C. likely represented by the left upper corner sub-figure. The figure situated at the center of the nine depicts simulations that apply to a society that has a considerably moderate level of resistance and response rate. The right lower corner depicts shows societies with a higher rate of response but a lower rate of resistance. Due to higher values of , every individual in the society moves further away relative to other scenarios and this results in a significant change in relative to the value of * . When it comes to the figure of the burden of cases, this scenario estimates the lowest number of cases relative to other scenarios. Societies that practice physical distancing and have a tendency to keep practicing it in a longer period, even after the disease is no longer present, are best represented by this scenario, resulting in a lower burden of cases relative to other scenarios.

Variations under Different Pathogen Characteristics
Different pathogens lead to different survivability periods in the air or other media. The longer the pathogens are active as airborne pathogens, the more they accumulate, which increases the risk of infections. Characteristics of the observed pathogens are governed by parameter . In System (2), the term − represents the concentration of pathogens per unit of time to become inactive. Figure 5 shows the dynamics of the disease prevalence + under different periods of pathogens lasting in the air for the same parameters as used in Figure 4. In the lower-right picture, it is shown that pathogens that can last up to 48 h (red) can accumulate up to 300,000 quanta pathogens per m and drive infections to as high as 23%. Figure 5a demonstrates how the dynamics of + precede on reaching a peak for exactly 2 days (48 h). It is natural to accept that the longer the period, the wider the gap between the occurrences of the two peaks. By setting a smaller period (higher ), the dynamics of decrease and so does + . Moreover, the peak of + shifts to the right (see Figure 5b,c). More results on the model's sensitivity analysis are provided in Appendices B and C.

Case Study: SARS-CoV-2 Spread in School
As mentioned in the previous section, the proposed model incorporates the sociobehavioral aspects of people in the society combined with the effect of airborne transmission. When it comes to socio-behavioral aspects, we included social resistance and social response amid the disease spread, which limit the scope of the implementation. At the scale of nations, people in society comprise all levels of education, culture, habits, or even wealth [48], which leads to a variety in perceptions when dealing with disease spread; some may have higher awareness but some may not. This fact challenges the modeler regarding how to estimate λ 1 and λ 2 that will accurately portray the society. Hence, we designed the model to be applied to the scale of an educational or office environment. It is natural to expect the homogeneity of socio-behaviors, even homogeneity in age, in schools or offices. These limitations also support the involvement of airborne transmission due to the indoor activity in schools or offices [24]. Henceforth, this section provides the applications of the proposed model to understand the disease spread in a university environment.

Dataset and Parameters' Estimation
We collected the data in respect of the SARS-CoV-2 spread in a college environment, namely the Institut Teknologi Bandung (ITB), and data range from early January until late April 2022. The data comprise record daily cases, current active cases, and the total number of recovered individuals out of all enrolled students, lecturers, and college staff. Although students and staff do not stay at the college 24/7, it is reasonable to assume that they spend most of the time in the college environment. Here, we exclude the enrolled students that were infected in other cities due to the hybrid (online-offline) learning practice. The data are privately available at https://covidtrak.itb.ac.id/ (accessed on 1 April 2022), which is only accessible by ITB COVID-19 task-force members.
In terms of the parameter estimation, we only used data for the daily new cases from early January 2022 until late April 2022, which will be later denoted as D a . Given in Table 4, System 204 and 210 involve 11 parameters, with only three of them being estimated by the integration of data D a , namely β * , v, and w, while other assumptions are as follows: the population size N equals 4000 (according to the report of the initial school reopening), the average vaccine efficacy ρ = 0.37 for SinoVac [36], γ 0 = 1/14, and γ 1 = 1/6, which represent the rate of recovery under lack of and excessive healthcare, respectively. In order to obtain the estimations of β * , v, and w, we used a Markov Chain Monte Carlo (MCMC) method to estimate the whole distribution. The complete Bayesian hierarchy for the MCMC method is provided in Appendix D. Figure 6 shows the estimated posterior distribution of β * , v, and w that was implemented to the data that resemble the recorded daily new cases.  Removal rate of airborne pathogens 2 1/day [21] Figure 6. (a-c) Posterior distribution for * , , and , estimated using MCMC method with prior of normal distributions: ( , ), with being the estimated single point and being a higher value to acquire the possibility of achieving the global minima. We generated three independent samples to portray the posterior density to ensure its consistency. (d) The comparison between the data (daily new cases) and the model with its 90% prediction interval.

Projected Number of Cases
Assuming no further changes in all parameters, the estimated posterior distribution of * , , and can be used to sample their values and generate the extrapolated trajectories for all states of the proposed model. Figure 7a and b show the projections of the disease prevalence in the university from early 2022 until mid-2023. Both consistently predict a significant decrease in the number of cases from May 2022, which implies a decrease in the average interaction distance , approaching the social natural distancing * . Figure 7a -c clearly show that the figures are estimated with a relatively narrow prediction interval, which leads to high confidence in the results under the hold assumptions. As the average interaction distance is around * , or, in other words, people in society behave as if there is no disease, the expected number of cases shown in Figure 7a and b increases in August 2022 and peaks in around October 2022, though the prediction interval is relatively wider compared to the previous period. These simulations show that the number of cases is expected to increase as approaches * , without even considering reinfection due to immunity waning. Figure 6. (a-c) Posterior distribution for β * , w and w, estimated using MCMC method with prior of normal distributions: N (µ, σ), with µ being the estimated single point and σ being a higher value to acquire the possibility of achieving the global minima. We generated three independent samples to portray the posterior density to ensure its consistency. (d) The comparison between the data (daily new cases) and the model with its 90% prediction interval.

Projected Number of Cases
Assuming no further changes in all parameters, the estimated posterior distribution of β * , v, and w can be used to sample their values and generate the extrapolated trajectories for all states of the proposed model. Figure 7a and b show the projections of the disease prevalence in the university from early 2022 until mid-2023. Both consistently predict a significant decrease in the number of cases from May 2022, which implies a decrease in the average interaction distance D, approaching the social natural distancing D * . Figure 7a-c clearly show that the figures are estimated with a relatively narrow prediction interval, which leads to high confidence in the results under the hold assumptions. As the average interaction distance D is around D * , or, in other words, people in society behave as if there is no disease, the expected number of cases shown in Figure 7a,b increases in August 2022 and peaks in around October 2022, though the prediction interval is relatively wider compared to the previous period. These simulations show that the number of cases is expected to increase as D approaches D * , without even considering reinfection due to immunity waning. were evaluated using the estimated parameters obtained in the previous section, with the assumption of no significant change in parameters for the 365-day prediction intervals: (a,b) extrapolated number of active cases and daily new cases that levels off in around May to July 2022, but starts to increase on August and peaks on October 2022, (c) average interaction distance of society that will approach its natural distancing of * = 1.2 m as the number of cases decreases, but increases when the case resurgence is identified, (d) time-dependent transmission rate that gradually increasing as approaches * , (e) dynamic for viral load over time that resembles that for active cases, and (f) time-dependent recovery rate that its values are bounded by and . Figure 7e depicts the dynamics of the pathogen concentrations in the observed area per m , which resembles the dynamics of the active cases over time. As stated in the previous section, the longer the pathogens can last in the air (or other media), the further the shift to the right relative to the dynamics of the active cases. In other words, the presence of + contributes to the presence of pathogens that govern the rate of transmissibility. Figure 7d and f show the dynamics of = ( , ) and = ( , ). Although they are not directly dependent on time, they are time-dependent due to the dependency of , , and to the unit of time. During the training time (initial time until the dashed lines), the rate of transmission decreases due to the significant deviation of the average interaction distance relative to * . It is expected that the trajectories of will increase during the prediction interval due to the decreasing values of . For the rate of transmission, its value is always bounded within the − ribbon. The rate is expected to approach the maximum value of as the burden of cases approaches zero; otherwise, it approaches the minimum . For Figure 7f, we set = 100, which represents the ability of the university hospital to accommodate only 100 patients at one time. This assumption causes a significant decrease in as the expected + exceeds the value of , depicting the ineffective health service as the burden exceeds its capacity.

Prospective Action Plans
Other than providing the extrapolated trajectories for all states, we are also interested in supplying numerical simulations related to prospective action plans for preventing the expected surge of COVID-19 in schools. This section provides the numerical assessment of three action plans: school reopening management, disinfection, and vaccine-related improvement. were evaluated using the estimated parameters obtained in the previous section, with the assumption of no significant change in parameters for the 365-day prediction intervals: (a,b) extrapolated number of active cases and daily new cases that levels off in around May to July 2022, but starts to increase on August and peaks on October 2022, (c) average interaction distance of society that will approach its natural distancing of D * = 1.2 m as the number of cases decreases, but increases when the case resurgence is identified, (d) time-dependent transmission rate that gradually increasing as D approaches D * , (e) dynamic for viral load over time that resembles that for active cases, and (f) time-dependent recovery rate that its values are bounded by γ 0 and γ 1 . Figure 7e depicts the dynamics of the pathogen concentrations in the observed area per m 3 , which resembles the dynamics of the active cases over time. As stated in the previous section, the longer the pathogens can last in the air (or other media), the further the shift to the right relative to the dynamics of the active cases. In other words, the presence of I u + I v contributes to the presence of pathogens that govern the rate of transmissibility. Figure 7d and f show the dynamics of β = β(D, V) and γ = γ(I, K). Although they are not directly dependent on time, they are time-dependent due to the dependency of D, V, and I to the unit of time. During the training time (initial time until the dashed lines), the rate of transmission β decreases due to the significant deviation of the average interaction distance relative to D * . It is expected that the trajectories of β will increase during the prediction interval due to the decreasing values of D. For the rate of transmission, its value is always bounded within the γ 0 − γ 1 ribbon. The rate is expected to approach the maximum value of γ as the burden of cases approaches zero; otherwise, it approaches the minimum γ. For Figure 7f, we set K = 100, which represents the ability of the university hospital to accommodate only 100 patients at one time. This assumption causes a significant decrease in γ as the expected I u + I v exceeds the value of K, depicting the ineffective health service as the burden exceeds its capacity.

Prospective Action Plans
Other than providing the extrapolated trajectories for all states, we are also interested in supplying numerical simulations related to prospective action plans for preventing the expected surge of COVID-19 in schools. This section provides the numerical assessment of three action plans: school reopening management, disinfection, and vaccinerelated improvement.

School Reopening Management
Although technologies support students in attending online classes, the practice of in-person classes is still preferable. This fact should be the main reason for the massive reopening of most Indonesian schools, regardless of the level of education. However, this should challenge the previous simulations due to a significant change in the number of individuals in a school as it is reopened. Henceforth, we provide numerical simulations of all states, more concerned with I u + I v , as the number of individuals in a school varies due to the school reopening. In practice, we assume that all individuals (students, lecturers, and staff) can be considered vulnerable to the disease. The higher the number of susceptible individuals, the more individuals can be infected. Hence, it is natural to analyze the effect of the increase in N on the dynamics of I u + I v .
Mathematically speaking, Substituting the derivatives of all states as stated in System (1), we have N = 0, meaning that the population size remains unchanged. However, we modified the model to accommodate the change in the population size due to the school reopening. Since we assume that all new individuals enter compartments S u and S v (with the proportion governed by the vaccine coverage), we add recruitment terms f and g for S u and S v as given by Equation (10).
This gives us N = f + g, for f = f (t) and g = g(t). Integrating both sides gives us N(t) = t 0 ( f (s) + g(s))ds. If we choose f (t) = (1 − α)N obj (t) and f (t) = αN obj (t), for a continuous and differentiable function N obj (t), then f (t) + g(t) = N obj (t) and we expect that N(t) = N obj . The simulation is conducted numerically, which includes the discretization of the time domain, and hence the condition of the differentiability of N obj is no longer relevant. The subscript objwhich stands for 'objective', denotes the preferred dynamics of N(t) that represent a certain school opening scheme. Hence, we can assess the effect of a specific school reopening scheme by choosing the appropriate function N obj that depicts the expected dynamics of the total individuals at any time t. Then, we choose three different N obj (t) values that represent three interesting school reopening schemes:

1.
No school reopening (benchmark) We preserve the size of the population as it was used to generate simulations in the previous section. We set N = 4000 for all t > 0, which leads to the constant population size for all time. This scenario is a benchmark for the other two scenarios.

2.
Gradual school reopening A gradual school reopening is a scheme that admits students and academical staff gradually until, at some point, the total number of students and staff is reached. In the Institut Teknologi Bandung (ITB), there are approximately 20,000 students and academical staff at any time for a non-pandemic era, which starts with only 4000 individuals in a pandemic era (January until April 2022). Hence, we choose a simple-bounded linearly increasing function N obj as given by: 20, 000 f or t > 201 t ∈ Z + , with t ∈ [0, 58], is the training data interval, which uses N = 4000. However, t ∈ [200, end] represents the total school reopening that starts in September 1, 2022, with N = 20, 000. The middle period of t ∈ [59, 200] represents a linearly gradual reopening from 4000 to 20,000. In practice, it is easy to add that f = 114(1 − α) and g = 114α during the period of reopening t ∈ [59, 200], and f = g = 0 otherwise.

3.
Prevalence-tuned school opening The last scenario accommodates the response of the school officials to reduce the school capacity as the disease prevalence level increases. Hence, we assume that the number of N should be related to the number of I. We chose a negative exponential to represent the relation between N and I as follows: N obj (t) = (20, 000 − 4000)e −kI(t) + 4000.
This formulation suggests that as I is around zero, then the school officials are about to totally open the school, and N = 20, 000. The opposite condition with a large number of I forces the school restriction and allows only 4000 individuals. This formulation of N obj is not explicitly time-dependent; instead, it depends on the varying values of I(t). In practice, we can set f = (1 − α)(20, 000 − 4000)ke −kt I (t) and g = α(20, 000 − 4000)ke −kt I (t). Figure 8 shows the numerical assessment of the three school reopening schemes. The simulations in red are the results that act as a benchmark for the other scenarios. This benchmark scenario gives the constant population size that drives the resurgence of the active cases around October as the average interaction distance increases. However, adding more people into the school through the gradual reopening scheme leads to more infections recorded, which reach a peak around July-August 2022. The surge is expected to happen since we add more people as N increases from 4000 to 20,000 in early September 2022. However, the infection-tuned scheme allows more people to enter the school relative to the other two schemes, yet results in lower cases compared to the second scenario. This is caused by the response of the school officials to reducing the school participants as the cases start to increase. This is the reason why cases increase in the same period as the second scheme but are significantly lower. By this simulation, all scenarios of reopening drive more people to enter the school, leading to more infections. The next section shows how the vaccine-related improvement can solve the problem of reopening without expecting any surge in infections. The last scenario accommodates the response of the school officials to reduce the school capacity as the disease prevalence level increases. Hence, we assume that the number of should be related to the number of . We chose a negative exponential to represent the relation between and as follows: ( ) = (20,000 − 4000) ( ) + 4,000.
This formulation suggests that as is around zero, then the school officials are about to totally open the school, and = 20,000. The opposite condition with a large number of forces the school restriction and allows only 4000 individuals. This formulation of is not explicitly time-dependent; instead, it depends on the varying values of ( ). In practice, we can set = (1 − )(20,000 − 4000) ( ) and = (20,000 − 4000) ( ). Figure 8 shows the numerical assessment of the three school reopening schemes. The simulations in red are the results that act as a benchmark for the other scenarios. This benchmark scenario gives the constant population size that drives the resurgence of the active cases around October as the average interaction distance increases. However, adding more people into the school through the gradual reopening scheme leads to more infections recorded, which reach a peak around July-August 2022. The surge is expected to happen since we add more people as increases from 4000 to 20,000 in early September 2022. However, the infection-tuned scheme allows more people to enter the school relative to the other two schemes, yet results in lower cases compared to the second scenario. This is caused by the response of the school officials to reducing the school participants as the cases start to increase. This is the reason why cases increase in the same period as the second scheme but are significantly lower. By this simulation, all scenarios of reopening drive more people to enter the school, leading to more infections. The next section shows how the vaccine-related improvement can solve the problem of reopening without expecting any surge in infections.

Vaccine Coverage and Effectiveness Improvement
Other than the physical distancing campaign, vaccination is one of the control measures in the spread of COVID-19, especially in a school environment. It has been shown that any school reopening leads to more infections recorded within the society. This section provides a simulation of the three reopening schemes whilst also varying the vaccine efficacy. By April 2022, the current average of vaccine efficacy is around 37%, as most Indonesians have been inoculated twice with SinoVac, which has 37% effectiveness in response to the Omicron variant. The improvement of the vaccine efficacy can be achieved by campaigning for a third vaccine dose with higher efficacy, such as Moderna, Pfizer, or Oxford AstraZeneca. Figure 9 shows the numerical simulations for different values of vaccine efficacy: ρ = 37%, 50%, 65%, and 80%.

Vaccine Coverage and Effectiveness Improvement
Other than the physical distancing campaign, vaccination is one of the control measures in the spread of COVID-19, especially in a school environment. It has been shown that any school reopening leads to more infections recorded within the society. This section provides a simulation of the three reopening schemes whilst also varying the vaccine efficacy. By April 2022, the current average of vaccine efficacy is around 37%, as most Indonesians have been inoculated twice with SinoVac, which has 37% effectiveness in response to the Omicron variant. The improvement of the vaccine efficacy can be achieved by campaigning for a third vaccine dose with higher efficacy, such as Moderna, Pfizer, or Oxford AstraZeneca. Figure 9 shows the numerical simulations for different values of vaccine efficacy: = 37%, 50%, 65%, and 80%. The figure illustrates the effect of the improvement of the vaccine efficacy, by assuming that 80% of school attendees have received a full-dose vaccine with such efficacy. When the efficacy is improved from 37% to 50%, this affects the first scenario that has only 4000 and cases that are expected to occur in October 2022 vanish. However, = 50% is not enough to reduce the other two scenarios (blue and black) significantly as the expected cases remain high for such scenarios. For = 65%, most of the expected cases are reduced significantly. Lastly, = 80% or higher is expected to reduce a whole surge of cases, at least for the 365-day prediction interval. It can be seen that as the number of cases reduces to zero, the average interaction distance of people approaches its natural level, yet there is no trigger for more infections due to the acquisition of vaccine-induced immunity. These results suggest that any reopening scheme, up to a maximal school capacity of The figure illustrates the effect of the improvement of the vaccine efficacy, by assuming that 80% of school attendees have received a full-dose vaccine with such efficacy. When the efficacy is improved from 37% to 50%, this affects the first scenario that has only 4000 and cases that are expected to occur in October 2022 vanish. However, ρ = 50% is not enough to reduce the other two scenarios (blue and black) significantly as the expected cases remain high for such scenarios. For ρ = 65%, most of the expected cases are reduced significantly. Lastly, ρ = 80% or higher is expected to reduce a whole surge of cases, at least for the 365-day prediction interval. It can be seen that as the number of cases reduces to zero, the average interaction distance of people approaches its natural level, yet there is no trigger for more infections due to the acquisition of vaccine-induced immunity. These results suggest that any reopening scheme, up to a maximal school capacity of 20,000 individuals, will not lead to any surge in COVID-19 cases as long as 80% of the population has received a vaccine with a minimum 80% efficacy.

Conclusions
In this study, we provided a modified SIR-type model that incorporates socio-viral behavioral aspects. The first aspect (socio-behavioral) was added to the model by integrating the average interaction distance in society, while the other was added by integrating the critical period in which airborne pathogens remain active. In a general case, a society with a higher resistance rate λ 1 but a lower response rate λ 2 will record more total cases compared to other plausible scenarios. In other words, the mentioned scenario applies to society with people that are hardly accepting new distancing habits and that do not have the awareness of disease prevalence. In contrast, a society with people that easily adapt to new distancing behaviors due to disease transmission, representing a society with a higher λ 2 but a lower λ 1 , will result in the least total cases compared to other scenarios. Furthermore, varying the critical period for active airborne pathogens also influence the model behaviors. The higher the critical period, the longer the airborne pathogens actively contribute to the increase in transmission rate.
As a case study, we implemented the proposed data on the spread of COVID-19 in a school environment to preserve the assumption of homogeneity in the population. Using the data on infections, we inferred the unknown parameters using the Bayesian approach. We have shown that the rendered model is well-depicting the training data. Using the inferred parameters, we extrapolated the model and came up with the evidence for a resurgence of cases in around August 2022. The resurgence of the case is purely implied by the return of society to its natural distancing behavior Dˆ* when no new COVID-19 cases are recorded. The dynamics of airborne pathogens load V seem not to influence that significantly due to the short critical period of SARS-CoV-2 to remain active in the air.
In response to the resurgence of cases, we used the model to numerically assess some strategic actions, applicable to the school context, to prevent the resurgence. First, we define some reasonable school reopening schemes that influence the population size: no reopening, gradual reopening, and infection-tuned reopening. While the first has a constant population size N, the second is gradually increasing the population size until it reaches the maximum capacity. Different from the other two, the infection-tuned reopening is a scheme that increases N to its maximum capacity when no diseases are identified, but also allows for decreasing N as the number of cases increases. Though the third scheme seems not practical, an infection-tuned scheme is proven to be the most effective strategy to reopen the school and minimize the risk of the rerise of COVID-19. Second, since we have demonstrated that all school reopening schemes lead to the resurgence of COVID-19 cases, we provide a numerical simulation that justifies the importance of vaccine quality; coverage, and efficacy. With constant vaccination coverage, increasing the vaccine efficacy will reduce the risk of COVID-19 resurgence-a vaccine with an efficacy of more than 80% has been proven to effectively prevent the COVID-19 resurgence, regardless of how society behaves towards the disease spread. Acknowledgments: Authors would like to deliver a big appreciation for the ITB team for covidtrak that has voluntarily collecting data that was used for this research. We also should appreciate Yayasan ITB for funding this research. This research is funded by Riset Unggulan 2022 with research grant number of 1V/IT1.C02/TA.00/2022 and ITB.

Conflicts of Interest:
The authors have no conflicts of interest to declare. All co-authors have seen and agree with the contents of the manuscript and there is no financial interest to report. We certify that the submission is original work and is not under review at any other publication.

Appendix A. Model Analysis and Threshold Number
For the sake of simplicity, we drop the vaccination effect and hence merge all compartments having indices u and v together. To generate both the disease-free and endemic equilibria, we added the natural disease in all state compartment (SIR), with µ represents the natural death rate. By substituting β(D, V) and γ(I, K) with those given by Equations (3) and (4), we have our system be rewritten as follows.
All plausible equilibria of System (5) are obtained by solving this nonlinear system, that is modified by adding the recruitment and natural death rate in order to get the Endemic Equilibrium.
The disease-free equilibrium (DFE) can be obtained by pluggingÎ = 0 to the system, which leads to DFE = {S * , 0, R * , D * , 0}. Using the next generation matrix method, the formula of the basic reproductive ratio is given by vanishes to 1 as all states approach the DFE. This quantity will determine whether the state will approach the DFE (when R 0 < 1) or otherwise approach the other equilibrium point, namely the EE.
contrast, the rate of society response depicts the increase in interaction distancing per increase in point prevalence, which inhibits the disease spread when this value is set high. In this section, we provide the figure of infected individuals (per one thousand of population) as λ 1 and λ 2 is set varied.
Given in Figure A1, the setting of λ 1 = 0 while λ 2 = 0.5 depicts the situation that people in society tend to response to the figure prevalence by distancing interaction but have no intention to return on their previous interaction habit. The blue line in the bottomleft corner represents the average interaction-distance in the mentioned scenario that levels off in D = 4.2 without approaching the natural interaction-distance. By fixing λ 2 = 0.5, the increment of λ 1 will lower the average interaction-distance and fasten the dynamics of D to approach D * . Poor society in certain regions tend to hasten on returning back on their previous habits (to work, study, etc) once the decrease in the point prevalence are declared, which in this model is considered to have relatively higher λ 1 . In contrast, some regions in developed countries may have lower value of λ 1 . Hence, the model requires the homogeneity of society; the smaller the society, the easier to assume, such as: students and academic staff in closed school, workers. increase in point prevalence, which inhibits the disease spread when this value is set high.
In this section, we provide the figure of infected individuals (per one thousand of population) as and is set varied. Given in Figure A1, the setting of = 0 while = 0.5 depicts the situation that people in society tend to response to the figure prevalence by distancing interaction but have no intention to return on their previous interaction habit. The blue line in the bottomleft corner represents the average interaction-distance in the mentioned scenario that levels off in = 4.2 without approaching the natural interaction-distance. By fixing = 0.5, the increment of will lower the average interaction-distance and fasten the dynamics of to approach * . Poor society in certain regions tend to hasten on returning back on their previous habits (to work, study, etc) once the decrease in the point prevalence are declared, which in this model is considered to have relatively higher . In contrast, some regions in developed countries may have lower value of . Hence, the model requires the homogeneity of society; the smaller the society, the easier to assume, such as: students and academic staff in closed school, workers. Right-hand side figures in Figure A1 depicts how social response to point prevalence affects the dynamics of both infected individuals and interaction-distance. The scenario setting = 0 portrays the society with no attention to current disease spread, implying to a steady interaction-distance to its natural habit. The higher the setting of , the more responsive the people in society to the current point prevalence.  Figure A1 depicts how social response to point prevalence affects the dynamics of both infected individuals and interaction-distance. The scenario setting λ 1 = 0 portrays the society with no attention to current disease spread, implying to a steady interaction-distance to its natural habit. The higher the setting of λ 2 , the more responsive the people in society to the current point prevalence.

Appendix C. Numerical Simulations under Different Healthcare Capacity
Trop. Med. Infect. Dis. 2022, 7, x FOR PEER REVIEW 21 of 23 Figure A2. Numerical simulations under different values of healthcare capacities: (a) comparison between the dynamics of active cases and the rate of recovery using the same parameters as used in Figure A1, (b) and (c) sensitivity analysis of + and the rate of recovery under different value of healthcare capacity .

Appendix D. Bayesian Hierarchical for Parameters' Estimation
There are three parameters that are estimated by the integration of the provided data . We assume that those three parameters are each a realization of a complete posterior distribution. Let us assume that ( ), = 1,2, . . . . , denotes the daily cases of COVID-19 at day , while ; represents the dynamics of the daily cases evaluated from the proposed model given that = [ * , , ] . We expect that observation ( ) (experimental data or observed data) to be equal as the model response , plus the independent and identically distributed error , with mean zero and variance . Hence, we can write: ( ) = ; + with (0, ). The goal is to estimate the posteriod distribution of | ( ) , which quantify the probability of parameter given the set of observational data. | ( ) = ℒ( ( ) ) ( ) with ℒ( ( ) ) represent the likelihood and is the prior knowledge of the estimated parameters. In terms of the data fitting, it is common to define the likelihood function as given in the following formula: . We will use the MCMC method (with a Metropolis-Hastings Algorithm) to get samples from the posterior distributions given in 860. Since we have no prior knowledge for each of the estimated parameters, we choose the prior of normal distribution with zero mean, but high in variance. In MATLAB, we used the builtin function provided by Grinsted, A. available in https://www.mathworks.com/matlabcentral/fileexchange/49820-ensemble-mcmc-sampler (accessed on 17 June 2022). between the dynamics of active cases and the rate of recovery using the same parameters as used in Figure A1, (b,c) sensitivity analysis of I u + I v and the rate of recovery γ under different value of healthcare capacity K.

Appendix D. Bayesian Hierarchical for Parameters' Estimation
There are three parameters that are estimated by the integration of the provided data D α . We assume that those three parameters are each a realization of a complete posterior distribution. Let us assume that D α (i), i = 1, 2, . . . , n denotes the daily cases of COVID-19 at day i, while f i;θ represents the dynamics of the daily cases evaluated from the proposed model given thatθ = [β * , v, w] T . We expect that observation D α (i) (experimental data or observed data) to be equal as the model response f i,θ plus the independent and identically distributed error ε i , with mean zero and variance σ 2 . Hence, we can write: D α (i) = f i;θ + ε i with ε i N 0, σ 2 . The goal is to estimate the posteriod distribution of π θ |D α (i) , which quantify the probability of parameterθ given the set of observational data. π θ |D α (i) = L(D α (i) θ )π 0 θ N L with L(D α (i) θ ) represent the likelihood and π 0 θ is the prior knowledge of the estimated parameters. In terms of the data fitting, it is common to define the likelihood function as given in the following formula: (D α (i) θ ) = exp − SS q Var with SS s = ∑ n i D α (i) − f i;θ 2 . We will use the MCMC method (with a Metropolis-Hastings Algorithm) to get samples from the posterior distributions given in 860. Since we have no prior knowledge for each of the estimated parameters, we choose the prior of normal distribution with zero mean, but high in variance. In MATLAB, we used the built-in function provided by Grinsted, A. available in https://www.mathworks.com/ matlabcentral/fileexchange/49820-ensemble-mcmc-sampler (accessed on 17 June 2022).