The model aims to emulate a closed human population in a village. We refer to this village as Barangay [29], the smallest political unit in the Philippines. In the model, two adjacent barangays, Barangay A and Barangay B, were considered and each of which consists of residents, their homes, and marketplaces where they can move in and out. Since the model replicates the state of barangays during the ECQ, it incorporates LGU guidelines as other factors affecting the simulation of the disease spread. These government guidelines, described in Table 1, are the implementation of FM-Pass, curfews, testing, and barangay lockdown.

The following points are the model assumptions:

The model considers three types of agents: a resident, a household, and a marketplace. A household is a representation of a house, while a marketplace is an abstraction of an actual public market where people buy their necessities. A patch, including its Moore neighborhood with r equal to 1, indicates a marketplace [30]. Both agents are immobile and do not affect the interactions of other agents. The only mobile type of agent would be the residents, which is defined as a person or an individual who lives in a barangay, resides in a household, and can travel to a marketplace. The behaviors of residents can change during the whole duration of the simulation (refer to the process overview and simulation section for further details). Table 2 tabulates and describes each of the agents' states and properties.

Table 3 summarizes the definitions, values, and references of each parameter in the model. Some of the parameters in the model are empirically defined. Since model parameters are sensitive to initial values, parameter variability is considered in the model. Furthermore, several parameters are allowed to be set differently for each barangay.

The model incorporated testing compliance, i.e., the voluntary-testing-compliance parameter, which is the probability of an infected resident being tested and consequently being isolated from the system. The parameters fm-pass and brgy-lockdown are based on actual government guidelines imposed during ECQ. For the relief system, the three relief types were mainly assumptions based on the experience during the food relief rollout in Davao City.

The range of the number of households is defined based on how dense the model space would be in its maximum range, and vice versa in its minimum. The same goes for the number of marketplaces for both barangays.

The global variables for the time, epidemiological values, and space are initialized first. For the globals, the day is first initialized to the value −1. The constant, ave-incubation-period, is set to the product 7 × 24, which is equivalent to 7 days. The value is described in a global meta-analysis and Chinese observational study as the mean incubation period [31]. The mean susceptibility for each barangay with r0_average, i.e., average-r0-brgy-A or average-r0-brgy-B parameter of Barangay A and Barangay B, respectively, is given by Equation 1. This equation was derived by empirically fitting the results of the proportion of susceptible individuals as a function of the reproduction number of the ABM. The details of the empirical fitting are found in Table 4.

By default, the model dimensions are set through the model settings. Specifically, each grid has a size of 13.4 × 13.4 pixels. The location of the origin is centered in the model space, and the minimum and maximum x-coordinates are −40 and 40 grids in dimension, while the minimum and maximum y-coordinates are −20 and 20 grids in dimension, respectively. A grid is assumed to be a scaled-down version of a 100 sq. m block. The model space will be modified during initialization if the parameter brgy-lockdown is set to TRUE. If this is the case, the lockdown barrier will be set up by indicating a unique different color to the origin along the y-axis, different from an empty space. This will create a clear separation between Barangay A and Barangay B.

As mentioned, a single marketplace in the model space is a single grid including its Moore neighborhood with r = 1. The marketplaces will be initialized first before all agents. By adjusting the parameters, market-count-brgy-A and market-count-brgy-B, the number of marketplaces for each barangay will be created. Each marketplace is assigned a unique random grid in the model space of every setup. The marketplaces of Barangay A are created first followed by the marketplaces of Barangay B.

A single household is a single grid in the model space or a hundred square meter space. After setting up the marketplaces, the households will be initialized next. Therefore, the marketplaces, the parameters, brgy-A-households, and brgy-B-households define the number of households for each barangay. The households are randomly assigned a position, ensuring that no two households share the same location in the model space. The occupied? state is set to FALSE by default initially. Finally, for each household, the value of persons-in-household is the maximum occupancy a household can hold. It is randomized through normal distribution with a mean and standard deviation of 4 and 1, respectively [32].

A single grid can be occupied by multiple residents. The initial total number of residents for each barangay depends on the total sum of the persons-in-household value of all households in a barangay. A condition for a resident to be assigned a permanent household is to move it to the location of a household with its occupied? state set to FALSE and the brgy-A? or brgy-B? state set to TRUE, depending on which barangay the resident is assigned. Once a household reaches its maximum occupancy, the occupied? state of the household is then set to TRUE.

All residents have the susceptible? state initially set to TRUE, then the rest of the EIR (exposed?, infected?, and recovered?) states are assigned FALSE. The susceptibility-rate value is randomized through normal distribution with the mean and standard deviation of the barangay's mean susceptibility (see Equation 1) and 0.5 as the value of sd-susceptibility or the standard deviation, respectively. The same is applied to the incubation-duration but with a mean of ave-incubation-period, or 7 days, and a standard deviation of 1 day or 24 ticks.

Once the aforementioned setup procedures are done, a number of infected agents for Barangay A and Barangay B are randomly chosen based on the values of the parameters, initial-infected-residents-Brgy-A, and initial-infected-residents-Brgy-B, respectively.

The simulation is an iterative process wherein each loop represents 1 h. Every timestep performs the behaviors of agents and updates all the states. Only the residents consist of complex processes. Table 5 lists and describes the behavior of the residents.

The model was formulated with the spatial perspective of Davao City and the early situation of the pandemic in mind. Therefore, to test this agent-based model, scenarios will be established with specific parameter settings matching the aforementioned circumstance. It is also assumed that different variants are in circulation to replicate the present conditions where several variants are already spreading in the community. For this study, three barangay scenarios were hypothetically defined namely: the worst-case scenario, compromise-case scenario, and best-case scenario. All three scenarios' parameter values are specified in Table 6.

The worst-case scenario expects the worst results. For one, there is less compliance with testing which means less detection and isolation of infected residents. The household densities in Barangay A and Barangay B are the highest in this scenario as well. Furthermore, there is only one marketplace for every barangay which means many residents would possibly converge to these individual marketplaces. A sample figure of what the Barangays in the worst scenario might look like upon initialization is given in Figure 3. The visualization shows dense households in both barangays.

The compromise scenario expects better results compared with the worst scenario. A sample visualization of the initialization in this scenario is also given in Figure 3. For instance, the residents have more compliance than in the worst-case scenario, but it is still not the ideal testing compliance rate. There is lesser household density in this scenario compared with the worst-case scenario, and there are more marketplaces within each of the barangays.

Finally, the best-case scenario expects the relatively best situation of a community epidemiologically. The parameters in this scenario have the highest voluntary compliance compared with the other two scenarios. The number of households for each barangay is also the lowest compared with the other two scenarios, and the number of available marketplaces can abundantly cater to the population. A sample initialization of this scenario is shown in Figure 3.

The agent-based model of this article was executed using the NetLogo version 6.0.4 [17]. The minimum system requirements to run the model were based on the system requirements set by the NetLogo software. Consequently, the model was run using a Windows 10, 64-bit laptop with 16 GB RAM, and a GTX 1070 graphics card.

A total of 100 simulation runs were performed with 8,760 ticks each (~1 year) for each of the scenarios. Specific parameters such as fm-pass and brgy-lockdown are set to TRUE and FALSE, respectively, by default. Two variants are in circulation within these two barangays—one with a high transmission rate (93%) and one with a low (30%) transmission rate. During runtime, the total infected cases per day from Barangays A and B were gathered in every simulation run per food relief system and per barangay scenario to generate a 1-year epidemic curve with a 95% confidence interval estimate. Such an epidemic curve will provide a quick overview of the possible effect of the food relief system. Furthermore, the simulated epidemic curve peak among all the possible information was obtained since the peak is a representation of the possible healthcare burden [33]. Consequently, the simulated maximum number of COVID-19 cases, i.e., epidemic curve peak, for each food relief system under a particular barangay scenario were compared and statistically tested at a 5% level of significance. The most recommended relief type for every scenario is then identified and further discussed in the next section.

The leftmost image in Figure 5 displays the distribution of the simulated maximum number of COVID-19 cases and it visually appears that at least one food relief system yielded a relatively lower number of COVID-19 cases. In the figure, the highest distribution of average infected cases is the “irregular” relief type, while the “regular and sufficient” relief type has the lowest distribution of a number of cases. The Kruskal–Wallis test for equal medians revealed that there is a significant difference between median COVID-19 cases. Furthermore, its post-hoc test revealed that all three relief systems are significantly different from each other in this scenario.

When both barangays are densely populated, disease transmission is the highest with the “irregular” food relief system as more contact between residents can be expected due to their basic needs than with the other food relief system [34]. The disease transmission is further expedited with limited available marketplaces as agents tend to cluster, making the virus easily transmitted [35]. Meanwhile, it is noteworthy to mention that despite all the factors that constitute the worst scenario, with the “regular and sufficient” food relief type, the virus cannot propagate as shown in the worst scenario epidemic curve in Figure 4 since the residents rarely go outside their households. Hence, the best relief system for a locality with similar circumstances as the worst parameter is a “regular and sufficient” food relief system.

The distribution of the simulated maximum number of COVID-19 cases for the three relief systems in the compromise-case scenario is visualized in the middle image in Figure 5. The “irregular” relief system still showed the highest maximum number of COVID-19 cases. On the contrary, the “regular but insufficient” and “regular and sufficient” relief types showed little difference in the simulated maximum number of COVID-19 cases. This is also observed in the epidemic curve in Figure 4.

As in the worst-case scenario, there is a significant difference in the maximum number of COVID-19 cases among the three relief systems. Moreover, the three food relief systems are significantly different from each other in the compromise scenario. While the result of the statistical analysis points out that the relief systems are significantly different, the plot shown in Figure 5 would suggest almost no difference between “regular but insufficient” and “regular and sufficient” reliefs. However, it could be argued that having a single infected case could still be a cause for a lot of transmissions [36]; therefore, the best relief system to opt for in this scenario would be “regular and sufficient.” However, since the Philippines is subject to certain bureaucratic and government procurement complexity [37] that could limit a barangay's resources, then the least relief system the locality could implement would be “regular but insufficient.”

Finally, the best-case scenario expects the best situation of a community, epidemiologically. The parameters in this scenario have the highest voluntary compliance compared with the other scenarios. The number of households for each barangay is also lower, and there are more marketplaces within the simulated space.

The rightmost image in Figure 5 shows similar averages among all relief systems. Among the three food relief systems, the “regular and sufficient” relief type yielded no variability. Variability is essential in comparing scenarios so we cannot include it in the analysis. The relief systems left to compare are “regular but insufficient” and “insufficient.” Due to less density, as well as the increase in the number of marketplaces that results in less convergence of individuals, the virus could no longer sustain its spread which, in the simulation runs, resulted in an average of only one maximum case [35]. This barangay scenario, therefore, has favorable conditions for the virus to die out. These favorable conditions are compounded due to the movement of individuals which is affected depending on the relief type. Hence, based on the results, the regularity and sufficiency of food relief do not matter. In other words, an “irregular” relief type would suffice in similar barangay scenarios such as this.

This model aimed to describe the early improvised interventions imposed by the local government unit of Davao City during the onset of the pandemic in the Philippines: FM-PASS coding, curfews, and lockdowns. The concept of the model revolves around movement restriction through incentives, i.e., the food relief system as an added non-pharmaceutical intervention. To our knowledge, there have been no other models so far that have conceptualized similar interventions in their models. Distributing food relief has also been done in different countries, and common problems persisted [38, 39]. In Uganda, for example, problems such as delays, poor quality of food aid, and arrests due to strict quarantine policies are rampant and are being criticized [39]. Upon further examination of Uganda's situation, one recommendation by Nathan et al. [39] is decentralized food distribution. In this study, individuals going to decentralized food sources (e.g., marketplaces), in a way, limit the spread of the virus between villages and the distribution of the food relief further controls the possible viral spread within the village.

Most of the COVID-19 models are based on an SEIR compartmental framework, which is also the basis of our epidemiological model. Since the model's transmission dynamics heavily depend on the mobility of the agents, such results are best compared to models studying the effect of mobility on COVID-19 incidences. Oka et al. [40] used an SIR stochastic model to describe the spatiotemporal evolution of COVID-19 across several provincial regions of China. A critical component of their model is capturing the impact of human movement as part of disease transmission. Notably, Oka et al. [40] took into consideration external infections from other regions and examined scenarios wherein movement restrictions were not in effect across regions. The simulated number of infected individuals expectedly showed an increasing trend, and delaying further restrictions predicted an even higher growth. Such results of Oka et al. [40] have been reflected by our model, especially in the “irregular” relief distribution for the worst-case and compromise-case scenarios.

Zhou et al. [41], on the contrary, built an SEIR model for COVID-19 using ordinary differential equations which explored the effect of varying intensity of movement limitations on COVID-19 transmission. In their model, they used mobile phone mobility data to establish the number of people moving to or within regions. The way they reduced movement is by setting the level of mobility by 20, 40, and 60% of infectious people only, whereas in this model, all agents are affected by movement reduction based on the level of regularity and sufficiency of food relief. A 60% reduction in movement simulated the flattest epidemic curve across all scenarios, while the 20% reduction in movement simulated the highest among the restriction scenarios. Our simulation results are also in agreement with Zhou et al. [41]—our highest epidemic curve was observed when agents produced more movement and interactions in situations when there is irregularity and insufficiency of food relief, while the opposite saw the immediate flattening of the curve.

One of the main differences between this model as compared with the other models is the behavioral assumptions that are taken into account when people are subjected to restrictions as well as the uncertainty of food security. Some work sectors have been resistant to the community lockdowns driven by their physiological needs as the work-from-home (WFH) setup during the pandemic is not amenable for them. Consequently, human mobility persists which underscores the importance of a food relief system as an additional layer of the implemented NPIs. For instance, the impact of community quarantine is affected by how agreeable the tertiary workforce is to a work-from-home setup [42]. It is then important for the government to take into consideration an assistance program or to allocate resources for the affected constituents within their jurisdiction.

This section serves as a demonstration of the possible application of the results. Davao City's barangays are chosen for an actual scenario comparison in this study. To generally describe the population characteristics of these barangays, they are grouped together into 11 districts. The boundaries of the districts and population densities can be visualized in Figure 6.

Along with the ECQ implementation in the city, the local government imposed a “clustering” method to curb the pandemic [43]. This method divided the barangays by districts into “clusters” where residents were only allowed to move within their assigned district cluster. These clusters were also strictly given designated marketplaces within their district where they are only permitted to purchase necessities. However, there is an uneven proportion of marketplaces assigned for each cluster, that is, some clusters are assigned a few marketplaces in proportion to their high population density barangays. This further solidifies the scenario comparison wherein the population is restricted to moving to marketplaces and households only.

From the aforementioned information, the Barangays can be categorized under the three scenarios as follows: barangays under districts Agdao, Bunawan, Buhangin, Poblacion, and Talomo are mostly categorized under the worst scenario due to their high population density and tightly clustered barangays to low marketplaces count ratio. As an example, Agdao district with 11 barangays and a high population density of 16,676.50 people per sq. km only has a single marketplace to cater to this many people. This discourages distancing measures inside the market and encourages the transmission of the virus, hence the worst case. For the compromise scenario, this could be best distinguished by the barangays under districts Baguio, Calinan, Tugbok, and Toril because of their average population density and several marketplaces proportionally provided. Finally, districts with barangays with more than sufficient marketplaces to accommodate a low population density are found in Paquibato and Marilog for the best case. Pertinent information for each district, its clustering, and recommendations are shown in Table 7.

To provide proper recommendations on which relief system can be applicable, the number of barangays with cases and the average daily attack rate per 100,000 (ADAR) for each of the barangays are determined and calculated. The duration computed for the ADAR is from the imposition of ECQ and clustering back on 04 April 2020, up to a month. It is important to note that at this time, the number of cases was still minuscule compared with the number of cases today. Figure 6 summarizes the COVID-19 ADAR results into a choropleth map.

The summary for each district in Table 7 shows that the districts under the worst scenario all had barangays with confirmed cases. It can also be observed that despite the relatively high number of available marketplaces for the Talomo district, half of its barangays had confirmed cases. Reiterating the results of the simulation runs, the districts with the worst scenario are where it is strongly recommended to reinforce a regular and sufficient food relief system. The data further validate that worst scenario barangays are prone to COVID-19 transmission and therefore recommended regular and sufficient relief system distribution.

Districts under the compromise scenario have fewer barangays with confirmed cases. Under the recommendations in Table 7, some of these districts are given regular but insufficient distribution and some are given regular and sufficient. This is based on the data that some of these districts are starting to have barangays with confirmed cases, while some continue to have no cases at the time. But as stated in the previous section, regularity still matters provided a sufficient supply of resources, and as a fallback due to logistical limitations, the least food relief distribution is regular but insufficient for these barangays.

Finally, the districts categorized as the best scenario are recommended for the irregular distribution of food relief due to having no cases during the whole duration of the ECQ.