Managing bed capacity and timing of interventions: a COVID-19 model considering behavior and underreporting

1.   Introduction

Shortly after the first local transmissions of the coronavirus disease 2019 (COVID-19) in the Philippines were confirmed, the government imposed social distancing policies, termed community quarantines, which were largely implemented by the police and military ^[1,2]. By March 30, 2020, the country only had six laboratories that accommodated up to 1000 tests daily ^[3]. Contact tracing began slowly due to an insufficient number of contact tracers ^[4]. Testing capacity was increased to about 35000 tests daily by the end of September 2020 ^[5]. From April to September 2020, the weekly positivity rate ranged between 4.5 and 28%, higher than the 5% threshold set by the World Health Organization during this time ^[6,7].

The Philippines' capital region, Metro Manila, comprised around 12% of the country's population in 2020 ^[8]. Metro Manila was placed under Enhanced Community Quarantine (ECQ) on March 16, 2020 ^[9]. Under ECQ, movement was restricted to essential goods and services. Public transportation and mass gatherings were suspended ^[2]. People were encouraged to work from home and businesses were advised to do transactions online ^[1]. After two months, ECQ was replaced by the Modified Enhanced Community Quarantine (MECQ), a transition phase before easing further to General Community Quarantine (GCQ). On June 1, 2020, Metro Manila was placed under GCQ, wherein public transportation and other establishments, except those for leisure, were allowed to operate ^[2]. A surge in the number of cases occurred from July to August 2020 and consequently, Metro Manila was again placed under MECQ. During this time, the utilization of ICU, isolation, and ward beds in Metro Manila reached 77%, 74%, and 84%, respectively, placing most facilities on critical or high-risk, and prompting 80 medical societies, representing 80000 doctors and a million nurses, to demand a 'timeout' ^[10,11]. On August 19, 2020, Metro Manila returned to GCQ ^[12]. By the end of September 2020, 53% of the 309303 total confirmed cases in the Philippines belonged to Metro Manila ^[5].

Because of the lack of vaccines and limited antiviral therapies during the early phase of the COVID-19 pandemic, NPIs such as wearing masks, school and workplace closures, and travel restrictions were crucial disease control measures. In the Philippines, compliance with policies was not only prompted by public health campaigns, but also driven by uncertainty and anxiety about the disease, and fear of getting reprimanded by the authorities ^[13,14,15]. Some of those who got infected suffered stigma and were blamed for not following the protocols ^[16,17]. A study among low-income households in the Philippines done in the early phase of the pandemic reported that 66% of respondents who might experience symptoms considered staying at home instead of seeking medical attention ^[13].

Mathematical modeling has been extensively used to understand the dynamics of COVID-19 transmission ^{[18,19,20,21,22,23,24]}. Non-pharmaceutical interventions, behavior change, and underreporting of cases have been incorporated into mathematical models of COVID-19 ^{[25,26,27,28,29,30,31,32]}. In this study, we aim to investigate the effects of reporting and behavior on the timing and magnitude of the peak of COVID-19 infections in Metro Manila from March to September 2020. We extend the SEIQR model by Kim et al., which includes a compartment for behavior-changed susceptible individuals ^[28], by adding an unreported compartment to account for individuals who were undetected due to inadequate testing and tracing, or unwillingness to be detected. We provide a detailed mathematical analysis of the model showing the existence, nonnegativity, boundedness of solutions, computation of the threshold parameter, and two forms of disease-free equilibria (DFEs). We show the global stability of one of the DFEs by utilizing a suitable Lyapunov function and analyzing the asymptotic behavior of the states. Furthermore, a bi-objective optimization problem is formulated that allows the easing from ECQ to GCQ at an earlier time and minimizes the number of additional beds necessary to ensure sufficient capacity in healthcare facilities. The multiple optimal solutions (Pareto optimal solutions) of the bi-objective problem offer trade-off solutions in which the user, a policymaker or healthcare authority, can use in decision-making.

2.   Methods

2.1.   Data

The number of cumulative confirmed cases from March 8 to September 30, 2020, was obtained from the Philippine Department of Health (DOH) data drop ^[6]. These data were used to estimate the rates of transmission, behavior change, and reporting. The data on COVID-19 bed capacity and occupancy rate in Metro Manila were gathered from the number of occupied and available isolation, ward, and ICU beds from the weekly DOH bulletin from June 20 to September 26, 2020 ^[33]. The total population of Metro Manila was set to 13484462, based on the 2020 census data of the Philippine Statistics Authority ^[8].

2.2.   Mathematical model

The model we present is an extension of the model in ^[28] wherein an unreported compartment is added to represent the undetected or unreported COVID-19 cases in the early phase of the pandemic in the Philippines. The constant total population $N$ is divided into eight mutually exclusive compartments: susceptible ($S$), behavior-changed susceptible ($S_F$), exposed ($E$), reported infectious ($I$), unreported infectious ($I_u$), isolated ($Q$), recovered ($R$), and deaths ($D$). A schematic diagram of the model is shown in Figure 1.

Assuming a local, prevalence-based spread of fear of the disease and following the study of Perra et al. ^[34], the transition rate of a susceptible to a behavior-changed susceptible is given by $\beta_FQ/(N-Q-D)$. This means that a susceptible individual is more likely to change behavior as the number of confirmed cases among one's contacts increases. The movement back to $S$ is assumed to be influenced by the number of recoveries and susceptible individuals without behavior change ^[34]. As the recoveries and susceptible individuals increase among the contacts of a behavior-changed susceptible, the more likely the individual to exit the $S_F$ class and resume regular social behavior. The parameter $\mu$ represents the rate of the easing of behavior and its value is assumed to be 1/14 ^[28].

Susceptible individuals ($S$ and $S_F$) move to the exposed class upon contact with infectious individuals ($I$ and $I_u$) at a rate $\beta$. The transmission rate for the behavior-changed susceptible class is assumed to be reduced by a factor $\delta$. The reporting ratio $\rho$ partitions the exposed class to reported $I$ and unreported $I_u$ classes. Assuming that an individual becomes infectious 2 days before symptom onset ^[35], the mean incubation period of the original virus strain is 6 days ^[36], and the mean duration between symptom onset and the first medical consultation in the Philippines during this time was 6.75 days ^[37], we set the mean latent period ($1/\kappa$) to 4 days and the mean infectious period of the reported cases ($1/\alpha$) to 8.75 days. From isolation, individuals recover $1/\gamma$ days on average or die. The average fatality rate, which is the ratio of daily deaths to daily active cases and denoted by $f$, is set to 1.9% ^[38]. Those in the unreported class are assumed to have less severe symptoms and move to the recovered class $1/\eta$ days on average. The following system of differential equations describe the model used in the study:

where $S(0)\geq0, \ S_F(0)\geq0, \ E(0)\geq0, \ I(0)\geq0, \ I_u(0)\geq0, \ Q(0)\geq0, \ R(0)\geq0,$ and $D(0)\geq0$. We introduce the term $\tilde{N}$ because those in $Q$ and $D$ are assumed to be isolated and have no contact with the rest of the population. All the parameters are assumed to be positive and $\rho, f\in(0, 1]$. In the numerical simulations, we set the initial population of the infectious $I(0)$, exposed $E(0)$, and unreported $I_{u}(0)$ classes equal to the number of cases $1/\alpha$, $1/\alpha +1/\kappa$, and $10\times I_0$ days from March 8, 2020, respectively. The initial number of isolated individuals $Q(0)$ was the number of cases at the start of the estimation period. The initial susceptible population $S(0)$ was computed by getting the difference between the total population and $E(0)$, $I_{u}(0)$, $I(0)$, and $Q(0)$. The rest of the state variables were initially set to zero. The model parameters and initial values of the state variables are shown in Table 1.

2.3.   Least-squares fitting of parameters

The values of the transmission rate, reporting ratio, and reduction factor were estimated from the cumulative cases data in Metro Manila from March 8 to September 30, 2020. We divide the estimation period into two: Period 1 is from March 8 to May 31, and period 2 is from June 1 to September 30, 2020. Metro Manila was mostly under ECQ during period 1, while it was mostly under GCQ during period 2. It was during period 2 that the first major epidemic wave in the Philippines occurred. Since the intensity of NPIs and the behavior of the population during ECQ and GCQ vary, the values for the transmission rates ($\beta$ and $\beta_F$) and reporting ratio ($\rho$) in periods 1 and 2 are assumed to be different. The reduction in transmission ($\delta$) for the behavior-changed susceptible class is assumed to have the same value in the two periods. We denote the transmission rates and reporting ratio for periods 1 or 2 by the subscripts 1 or 2, respectively.

Denote by ${\bf p} = \left[\beta_1, \beta_2, \beta_{F, 1}, \beta_{F, 2}, \delta, \rho_1, \rho_2\right]$ the vector of parameters to be estimated and set the domain as $\Gamma = [0, 1]\times[0, 1]\times[0, 10^3]\times[0, 10^3]\times[0, 1]\times[0, 1]\times[0, 1]$. Define the objective functional $J:\Gamma\subset\mathbb{R}^7\rightarrow \mathbb{R}$ as

where $\tilde{C}(t_i)$ is the cumulative reported cases on day $t_i$ based on the data and $I(t_i; {\bf p})$ is the solution $I$ of (1) at $t = t_i$ using the parameter vector ${\bf p}$. To estimate the unknown parameters, we calculate the minimizer ${\bf p}^\star$ of the objective function (2.2) using a trust-region method based on the interior-reflective Newton method ^[40,41]. Each iteration of this approach solves a large linear system using preconditioned conjugate gradient method ^[42]. We utilize the Matlab built-in function $\texttt{lsqcurvefit}$ to implement this numerical optimization technique. The program requires the bounds, which we set based on the domain $\Gamma$, and initial guess, which we set to ${\bf p}_{init} = [0.3, 0.3,285,285, 0.2, 0.25, 0.25]$.

2.4.   LHS-PRCC and parameter bootstrapping

Sensitivity analysis is a numerical technique that is widely used in identifying and ranking critical parameters for a given model output ^[43]. A parameter is said to be influential to an output if small perturbations of its value lead to significant changes in the model output. In this work, we use the Partial Rank Correlation Coefficient (PRCC) method paired with the Latin Hypercube Sampling (LHS) technique. LHS-PRCC is a global sensitivity analysis technique and is widely used in infectious disease modeling ^{[44,45,46,47,48]}. We investigate the sensitivity of the ten parameters ${\bf \chi}: = \left[\beta, \rho, \alpha, \delta, \beta_F, \eta, \mu, \kappa, \gamma, f\right]\in[0.01, 1]\times[0.01, 1]$$\times[0.01, 0.5]\times[0.01, 1]\times[0.01, 10^3]$$\times[1/40, 1/7]\times[0.05, 1]\times[0.1, 0.5]\times[1/40, 1/7]\times[0.01, 0.1]$. To consider every infection, we use the cumulative number of infected individuals,

as the model output. For each parameter $\chi_j,$ where $j\in\{1, 2, \ldots, 10\}$, we sample 10000 values $\{{\bf X}_{ij}\}_{i = 1}^{10000}$, all following a uniform distribution, using LHS. Hence, ${\bf X}\in\mathbb{R}^{10000\times 10}$ with each row representing a sampled combination of ${\bf \chi}$. For brevity, we denote ${\bf X}_{i}$ as the $i$th row of ${\bf X}$ and ${\bf X}_{j}$ as the $j$th column of ${\bf X}$. Denote by ${\bf Y}\in\mathbb{R}^{10000}$ the output vector whose $i$th component is calculated by evaluating $\mathcal{F}\left(t; {\bf X}_{i}\right)$. The correlation coefficient (CC) is a metric to quantify the linear association between input and output. If the data is ranked-transformed first before calculating the CC, the result is a Spearman or rank correlation coefficient (RCC). Ranking is done by sorting both ${\bf X}_j$ and ${\bf Y}$ in descending order, and the integer rank values ${\bf X}_j^R$ and ${\bf Y}^R$, respectively, are stored. Partial correlation coefficient (PCC) is used to characterize the linear relationship between ${\bf X}_j$ and ${\bf Y}$ after the linear effects on ${\bf Y}$ of the other inputs are discounted ^[43]. The partial rank correlation coefficient (PRCC) between ${\bf X}_j$ and ${\bf Y}$ is the PCC between ${\bf X}_j^R$ and ${\bf Y}^R$. The PRCC is computed numerically using the Matlab built-in function $\texttt{partialcorr}$, with the $\texttt{'type'}$ set to $\texttt{'Spearman'}$. Note that the output ${\bf Y}$ is time-dependent as seen in (2.3). In our simulations, the PRCC values are calculated at five time points: April 19, May 31, August 2, October 4, and November 1, 2020.

Parameter bootstrapping is a statistical technique to quantify uncertainty and construct confidence intervals of estimated parameters. In this study, we utilize the algorithm introduced in ^[49], where large samples of synthetic data sets using the estimated model parameters are generated assuming a certain probability distribution structure. In our simulations, parameters are re-estimated from 10000 (number of realizations) synthetic data sets each generated by assuming a Poisson error structure. The mean, standard deviation, and 95% confidence intervals of the re-estimated parameters are determined. This technique is summarized in Algorithm 1.

2.5.   Optimization problem

Using the bed occupancy data from the DOH ^[33], we calculated that an average of 16% of the active cases $Q(t)$ occupied COVID-19 beds from June to September 2020. Fitting the weekly data on available beds, a linear function representing 75% of the bed capacity was obtained. The DOH categorizes a facility as high risk if the bed occupancy is 70% to 85%, and critical if bed occupancy is greater than 85%. From July 18 to August 8, 2020, most facilities in Metro Manila were on critical or high-risk, with a combined COVID-19 bed occupancy (isolation, ward, and ICU beds) exceeding 75% of the capacity. Here, we propose an optimization approach to determine the number of additional beds per day so that the number of cases requiring beds does not exceed 75% of the total bed capacity and if it is possible to transition from ECQ to GCQ earlier than June 1, 2020.

We denote by $Q_H(t)$ the number of active cases requiring beds and $H(t)$ the linear, time-dependent, data-fitted bed capacity. We consider two objectives:

$(1)$ Minimize the number of additional hospital beds needed per day ($\omega$) so that the 75% capacity is not reached;

$(2)$ Determine an earlier timing ($\tau$) of easing from ECQ to GCQ.

The problem can be formulated as a bi-objective constrained optimization problem expressed as follows:

such that

where $Q_H(t) = 0.16\cdot Q(t)$, $H_0$ is the baseline number of beds at time $\tau$ based on the data, and $Q(t)$ is solved from (2.1) by setting $\tau$ as the day when GCQ started. Since this is a bi-objective optimization problem, the solution is not unique but a pareto optimal set. We solve (2.4) and (2.5) using Genetic Algorithm, which has found a growing number of applications in various fields of science and engineering ^[50,51,52]. In particular, we implement the Matlab built-in function $\texttt{gamultiobj}$, which is based on a variant of Non-dominated Sorting Genetic Algorithm Ⅱ (NSGA-Ⅱ) ^[53,54].

3.   Results

3.1.   Model analysis

In the following theorems, we establish the existence, nonnegativity, and boundedness of the solutions to model (2.1) for $t\geq0.$ Moreover, we derive the disease-free equilibria and a threshold parameter for the given model.

Theorem 1 (Existence and uniqueness of solutions). There exists a unique solution to $(2.1)$ for $t\geq0$.

Proof. Model (2.1) can be written as

where the column vector $\mathit{\boldsymbol{x}}(t) = (S(t), S_F(t), E(t), I(t), I_u(t), Q(t), R(t), D(t))^T\in\mathbb{R}^8$ defines a mapping from $[0, +\infty]$ to $\mathbb{R}^8$. Moreover, let $\mathit{\boldsymbol{F}}(t, x(t)) = (F_1(t, x(t)), \ldots, F_8(t, x(t)))^T$ be a vector-valued function from $\mathbb{R}^8$ to $\mathbb{R}^8$ equal to the right-hand side of the equations in (2.1). It can be verified that $\mathit{\boldsymbol{F}}$ is continuous and each of the first-order partial derivatives of $\mathit{\boldsymbol{F}}$ with respect to $\mathit{\boldsymbol{x}}$ is continuous in its domain. By the existence and uniqueness theorem in ^[55], system (2.1) has a unique solution for $t\geq0$.

Theorem 2 (Nonnegativity and boundedness of solutions). Solutions to system $(2.1)$ with nonnegative initial conditions will remain nonnegative for all $t > 0$. Moreover, the solutions are bounded.

Proof. Using the notation in (3.1), assume that $\mathit{\boldsymbol{x}}_0\geq{\bf{0}}$. First, we show by contradiction that $S(t)\geq0$ for all $t > 0$. Following the approach in ^[56,57], suppose that there exists $t_1 > 0$ such that $S(t_1) = 0, \ S'(t_1) < 0, \ S_F(t) > 0, \ E(t) > 0, \ I(t) > 0, \ I_u(t) > 0, \ Q(t) > 0, \ R(t) > 0$, and $D(t) > 0$ for $t\in(0, t_1)$. Then, from the first equation of system (2.1),

which is a contradiction to the assumption $S'(t_1) < 0$. Here we choose $t_1$ such that $S_F(t_1) > 0$ and $R(t_1) > 0$. Hence, $S(t)\geq0$ for all $t\geq0$. Similarly, we can show that $S_F, \ I$, and $I_u$ remain nonnegative for $t\geq0$.

Next, since $I(t)\geq0$ for all $t\geq0$ then $Q'(t) = \alpha I-\gamma Q\geq -\gamma Q$ and hence, $Q(t)\geq Q(0)e^{-\gamma t}\geq0$. It also follows that $D'(t) = f\gamma Q\geq0$ and thus, $D(t)\geq0$ for all $t\geq0$. In a similar approach, it can be shown that $E(t)$ and $R(t)$ are nonnegative for all $t\geq0$.

Since $N$ is constant and $S, S_F, E, I, I_u, Q, R, D\geq 0$, then from the last equation of (2.1), $S, S_F, E, I, I_u, Q, R, D\leq N$.

By equating the right-hand side of (2.1) to zero, the model has DFEs of the following forms:

for any $0 \le \theta \le 1$ and $0 \le \phi \le 1$ such that $\theta + \phi \le 1$. Using the next-generation method and following the notations in ^[58], the matrix $F$ of new infections and the matrix $V$ representing the flow of individuals between compartments are given by

The basic reproduction number $\mathcal{R}_0$ is calculated by finding the spectral radius of the matrix $FV^{-1}$ and is given by

Moreover, the reproductive number $\mathcal{R}(t)$, which is the average number of secondary infections from an individual during one's infectious period, is expressed as

In the following theorems, we analyze the stability of the DFEs of the model.

Theorem 3 (Local asymptotic stability). $\mathcal{E}_1$ is locally asymptotically stable if $\mathcal{R}_0 < 1$ and $\mathcal{E}_2$ is unstable.

Proof. We calculate the eigenvalues $\lambda$ of the Jacobian matrix derived from system (2.1) at the DFEs. For $\mathcal{E}_1$, the corresponding characteristic polynomial is given by

where

and $\epsilon = 1/(\phi + \theta)$. It is clear from $\mathcal{C}_1(\lambda)$ that five of its roots are nonpositive. Thus, it is left for us to show that all of the roots of the polynomial $\mathcal{P}_1(\lambda)$ are negative or have negative real parts. By the Routh-Hurwitz criteria, it is sufficient to show that $a_1 > 0$, $a_3 > 0$, and $a_1 a_2 - a_3 > 0$ ^[59].

Suppose $\mathcal{R}_0 < 1$. Since all of the model parameters are positive, $a_1 > 0$. Next, we can rewrite $a_3$ as

Since $\mathcal{R}_0 < 1$, then $a_3 > 0$. Note that

where

It suffices to show that $A_0, A_1, A_2, A_3, A_4 > 0$. By expanding $A_0$, its last term will be cancelled out and so $A_0 > 0$. Clearly, $A_1 > 0$. Define

so that $\mathcal{R}_0 = \mathcal{R}_{0, I} + \mathcal{R}_{0, I_u}$. Since $0 < \mathcal{R}_0 < 1$, then $\mathcal{R}_{0, I} < 1$ and $\mathcal{R}_{0, I_u} < 1$. Hence,

Using the above-mentioned inequalities, it follows that $A_2$ and $A_3$ are both positive. Finally,

Therefore, $\mathcal{E}_1$ is stable.

For $\mathcal{E}_2$, the characteristic polynomial is given by

for some third-degree polynomial $\mathcal{P}_2(\lambda)$. Observe that one of the eigenvalues is $\mu > 0$. Thus, $\mathcal{E}_2$ is unstable.

Remark. If $\theta = 0$, then $a_3$ in $(3.5)$ and $a_1a_2-a_3$ in $(3.6)$ are both positive regardless of the value of $\mathcal{R}_0$, which means that $[0, 0, 0, 0, 0, 0, \phi N, (1-\phi)N)]$ is stable. Meanwhile if $\phi = 0$, then $[\theta N, 0, 0, 0, 0, 0, 0, (1-\theta)N)]$ is unstable if $\mathcal{R}_0 > 1$ since $a_3 < 0$. By some algebraic manipulation, $\mathcal{E}_1$ is unstable if $\phi < \theta(\mathcal{R}_0-1)$.

In the next theorem, the global stability of the DFE $\mathcal{E}_1$ was studied. Lyapunov stability was established using Lyapunov's direct method and asymptotic stability was shown following the definition in ^[60].

Theorem 4 (Global asymptotic stability of $\mathcal{E}_1$). The DFE $\mathcal{E}_1$ of the model $(2.1)$ is globally asymptotically stable in the domain

whenever $\mathcal{R}_0 < 1$.

Proof. Consider the Lyapunov function

where $\mathit{\boldsymbol{x}}$ is the vector of state variables as in (3.1) and $\mathcal{R}_0 = \frac{\beta\rho}{\alpha}+\frac{\beta(1-\rho)}{\eta}$ is the basic reproduction number. Note that $V$ is continuous and differentiable in $\Omega$, $V(\mathcal{E}_1) = 0$, and $V(\mathit{\boldsymbol{x}}) > 0$ for $\mathit{\boldsymbol{x}}\in\Omega, \ \mathit{\boldsymbol{x}}\neq\mathcal{E}_1$.

Next, we obtain $\dot{V}(\mathit{\boldsymbol{x}}) = \mathcal{R}_0E'(t) + \dfrac{\beta}{\alpha}I'(t) + \dfrac{\beta}{\eta}I'_u(t)$. Note that $\dot{V}(\mathcal{E}_1) = 0$. Substituting the expressions for $E'(t), I'(t),$ and $I'_u(t)$ from the model and simplifying, we get

since $\frac{\delta S_F+S}{\tilde{N}}\leq1$. It follows that $\dot{V}(\mathit{\boldsymbol{x}})\leq0$ for all $\mathit{\boldsymbol{x}}\in\Omega$ if $\mathcal{R}_0 < 1$. Therefore, $\mathcal{E}_1$ is Lyapunov stable ^[60]. To show that the DFE $\mathcal{E}_1$ is globally asymptotically stable, we need to show that $\mathit{\boldsymbol{x}}\rightarrow \mathcal{E}_1$ as $t\rightarrow \infty$ ^[60].

First, we show that $E, I, I_u\rightarrow0$ by showing that $V\rightarrow0$ as $t\rightarrow \infty$. Since $V$ is decreasing and nonnegative whenever $\mathcal{R}_0 < 1$, then by the monotone convergence theorem, $V$ converges to, say, $\xi$ as $t\rightarrow \infty$. We claim that $\xi = 0$. Suppose otherwise, that is, $\xi > 0$. Then $\xi\leq V(\mathit{\boldsymbol{x}}(t))\leq V(\mathit{\boldsymbol{x}}(0))$. Define the set $U: = \{\mathit{\boldsymbol{x}}\in\Omega: 0 < \xi\leq V(\mathit{\boldsymbol{x}}(t))\leq V(\mathit{\boldsymbol{x}}(0))\}$. Clearly, $U$ is compact and $\mathcal{E}_1\notin U$. Since $\dot{V}$ is continuous, then by the extreme value theorem, there exists $\zeta > 0$ such that

Because $\dot{V}(\mathit{\boldsymbol{x}}(t))\leq-\zeta$ for all $t$, we get

If we take $T > V(\mathit{\boldsymbol{x}}(0))/\zeta$, then $V(\mathit{\boldsymbol{x}}(T)) < 0$, which is a contradiction. Therefore, $V\rightarrow 0$ as $t\rightarrow \infty$, which implies that $E(t), I(t), I_u(t)\rightarrow 0$.

Next, since $R'(t), D'(t)\geq0$ then $R(t)$ and $D(t)$ are both increasing. Moreover, since $R(t)$ and $D(t)$ are bounded, then by the monotone convergence theorem, there exist $R^\star, D^\star \in[0, N]$ such that $R(t)\rightarrow R^\star$ and $D(t)\rightarrow D^\star$.

Let $\varepsilon_1 > 0$. Since $I(t)\rightarrow 0$, then there exists $T_1 > 0$ such that $Q'(t)+\gamma Q(t) = \alpha I < \varepsilon_1$ for all $t > T_1$. Letting $\varepsilon_1\rightarrow 0$, we have $Q'(t)\leq -\gamma Q(t)$. Equivalently, $Q(t)\leq Q(T_1)e^{-\gamma t}$. Therefore, as $t\rightarrow \infty$, $Q(t)\rightarrow 0$.

By removing some negative terms and noting that $S\leq\tilde{N}$, it follows from the second equation in (2.1) that

Let $\varepsilon_2 > 0$. Because $Q(t)\rightarrow 0$, then there exists $T_2 > 0$ such that

for all $t > T_2$. Also, since $R(t)\rightarrow R^\star$ and $D(t)\rightarrow D^\star$, then there exists $T_3 > 0$ so that

for all $t > T_3$ and $\psi = \frac{\mu R^\star}{N-D^\star}$. It follows that

Taking $T^\star = \max\{T_2, T_3\}$, (3.9)–(3.11) imply that

for all $t > T^\star$. Because $S_F(t)\leq N$, (3.12) becomes

Letting $\varepsilon_2\rightarrow 0$, we get $S_F'(t)+\psi S_F(t)\leq 0$ as $t\rightarrow \infty$. Thus, $S_F(t)\leq S_F(T^\star)e^{-\psi t}$ and $S_F(t)\rightarrow 0$.

Finally, $S = N-(S_F+E+I+I_u+Q+R+D)$ implies $S(t)\rightarrow N-R^\star-D^\star = :S^\star$.

Because $S^\star+R^\star+D^\star = N$, then there exist $\theta, \phi\in[0, 1]$ with $\theta+\phi\leq 1$ so that $\mathcal{E}_1 = [\theta N, 0, 0, 0, 0, 0, \phi N, (1-\theta - \phi)N]$. Thus, if $\mathcal{R}_0 < 1$, all solutions of (2.1) converge to the DFE of the form $\mathcal{E}_1$.

3.2.   Parameter estimation

The best model fit for the cumulative and daily cases is shown as the black curves in Figure 2. The red circles represent the data points. The estimated transmission rates $\beta_1$ and $\beta_{F, 1}$ in period 1 were 0.199 and 471.057, respectively. In period 2, the transmission rate of the disease $\beta_{2}$ increased to 0.361, while the transmission rate of awareness or fear of the disease $\beta_{F, 2}$ dropped to 67.783. The reporting ratio $\rho_1$ in period 1 was estimated at 28%, which increased to 86% in period 2. The reduction in transmission induced by behavior change $\delta$ was estimated at 0.202. The parameter estimates are given in Table 1.

The reproductive number $\mathcal{R}(t)$ is shown as the blue curve in Figure 2. Initially, $\mathcal{R}(t)$ was at 1.9 then decreased to 0.9 by the end of period 1. During period 2, $\mathcal{R}(t)$ remained above 1 from early June until mid-August, with a peak of up to 1.5 in mid-July 2020.

3.3.   Sensitivity and uncertainty analyses

Using the cumulative number of infected individuals as the model output, the results of the sensitivity analysis showed that $\beta$ (range: $0.82$ to $0.92$) and $\delta$ (range: $0.52$ to $0.68$) have the highest PRCC values, followed by $\rho$ (range: $-0.49$ to $-0.44$) and the parameters for the infectious periods, $\alpha$ (range: $-0.48$ to $-0.46$) and $\eta$ (range: $-0.45$ to $-0.28$). PRCC values of $\kappa$ declined over time, with its highest value at 0.436. The rest of the parameters have small magnitudes of PRCC. Figure 3 illustrates the results. Moreover, results of the parameter bootstrapping showed that the re-estimated values of $\beta, \beta_F, \delta$ and $\rho$ follow a normal distribution and the mean values of the estimates all fall within their respective 95% confidence intervals. Figure 4 shows the distributions of the parameter re-estimates. The mean, standard deviation, and confidence interval are also indicated in the figure.

3.4.   Effects of behavior change and reporting

The solid curves in the upper panel of Figure 5 are the plots of the susceptible class $S,$ while the dashed lines are the behavior-changed susceptible class $S_F$. Here we investigate what happens if people changed their behavior one (orange), two (yellow), three (purple), or four (green) weeks earlier. At the beginning of period 2, we calculated that the proportion of $S_F$ with respect to the total susceptible population was 88% ($S_F$: 11849000; $S$: 1592900). To incorporate early behavior change, we scale the value of $\beta_F$ to yield the same proportion of $S_F$ one to four weeks before the start of GCQ. The black curves in Figure 5 represent the plots of the model using the parameters in Table 1.

During period 1, we observe a switch in the populations of $S_F$ and $S$. By the end of period 1, $S_F$ comprises the majority (88%) of the susceptible classes. The impact of early behavior change is seen in the daily and cumulative cases in period 2, shown in the bottom panels of Figure 5. As behavior changed one, two, three, or four weeks earlier, the cumulative cases decreased to 140468, 115573, 93041, or 73328 from 163191 (model, black), respectively. These translate to reductions of 14%, 29%, 43%, or 55% in cumulative cases by September 30, 2020. The peak of the daily cases also reduced to 2148, 1870, 1618, and 1392 from 2408 (model, black), and the timing was delayed from one to four weeks.

When the community quarantine was relaxed to GCQ, we observe that the size of $S_F$ declined, while $S$ increased. Around mid-July, when the number of daily cases were increasing ($\mathcal{R}(t) > 1$), the behavior of the susceptible classes switched again. From that point on until the peak of the first big wave in Metro Manila ($\sim$August 14 according to the model), the proportion of $S_F$ among the susceptible classes increased from 59% to 89%.

The upper panels in Figure 6 show the effect of varying the values of $\mu$ and $\beta_F$ on the cumulative cases and timing of the peak of infections. Higher values of $\beta_{F}$ and lower factors of the behavior change ease rate $\mu$ in period 2 result in notable reductions in the number of cases and delay in the occurrence of the peak. For instance, if in period 2 we set $\beta_{F, 2}$ as 3 times $\beta_{F, 1}$ and $\mu$ reduced by 90%, then the cumulative cases by the end of September 2020 would have been approximately 30000 and the peak would have occurred around July 25 (140 days from March 8). The blue area on the heatmap for peak timing indicates that the big wave in period 2 did not occur until September 2020.

Finally, the bottom panel in Figure 6 shows the effect of increasing the reporting ratio $\rho_1$ in period 1 on the cumulative cases by September 30, 2020. Only slight differences in the reported cumulative cases (blue bars; range: 11817 to 13421) were observed if $\rho_1$ was increased by factors of 1.5, 2, 2.5, or 3. On the other hand, cumulative cases, including the unreported (red bars), can be reduced by 29%, 45%, 54%, or 61%, if $\rho_1$ was increased by factors of 1.5, 2, 2.5, or 3, respectively.

3.5.   Optimal bed capacity and timing of policy change

Panel (A) in Figure 7 shows the fitted bed capacity $H(t)$ (dashed curve) from the data (red circles) and the number of cases requiring beds $Q_H(t)$ from the model. Note that the red circles depict 75% of the total COVID-19 bed capacity during this time and $Q_H(t)$ is 16% of $Q(t)$, representing the average number of active cases that occupy beds. By calculating the slope of $H(t)$, denoted by $\omega_\text{data}$, an estimated 33 beds were added per day from June 21 to September 30, 2022, in Metro Manila. Notably, $Q_H(t) > H(t)$ during the second MECQ, when healthcare workers demanded a 'timeout'. The peak of $Q_H(t)$ was calculated at 5307 cases.

The optimal solutions of the bi-objective optimization problem in (2.4) form a Pareto optimal set illustrated in Figure 7 Panel (B). The circles are all optimal solutions depicting different policies. The blue-colored optimal solution corresponds to the earliest easing to GCQ on May 2, 2020, and requires 131 additional beds per day to ensure that the bed capacity is adequate (up to 75% occupancy) during the surge in cases following the lifting of restrictions (blue curves in Panel (C)). On the other hand, the yellow-colored optimal solution has the least number of additional beds per day (47 beds) and the latest start of GCQ (on May 28, 2020). This solution has a delayed and lower peak of infections compared to the blue Pareto optimal solution (see Panel (C)). Compared to the curves in Panel (A), the number of cases requiring beds $Q_H(t; \omega^\star, \tau^\star)$ shown in Panel (C) is below the optimal bed capacity $H(t; \omega^\star, \tau^\star)$. Hence, constraint (2.5) of the optimization problem is satisfied. Panel (D) shows the peak sizes of the epidemic waves corresponding to the various Pareto optimal solutions. The smallest peak size (4807 cases, purple) is the optimal solution with GCQ starting on May 20, 2020, and with at least 56 additional beds.

4.   Discussion

Using the estimated parameter values, we observe that the model captures the trend of the daily and cumulative data from March until November 2020. The model shows a small peak in the number of daily cases (204 cases) and a slow increase in cumulative cases during period 1. A much higher peak (2408 cases) around mid-August 2020 and a sharp rise in cumulative cases are seen from the model during period 2. A delay of about one month between the drop in $\mathcal{R}_t$ and the decline in daily cases was also observed. Results of the parameter bootstrapping suggest good reliability of the estimated parameters. Moreover, sensitivity analysis showed that the transmission rate $\beta$ was the most sensitive parameter with respect to the number of cumulative infections. A higher reporting ratio $\rho$ or shorter mean infectious period of reported cases $1/\alpha$ reduces the cumulative infections. These results suggest that intensifying testing and tracing efforts can effectively reduce new infections. The average latent period ($1/\kappa$), which has the effect of delaying infection, becomes less sensitive to the cumulative number of infections as the epidemic progresses, while the mean infectious period of unreported cases ($1/\eta$) becomes more sensitive as the epidemic progressed.

Reporting was low in the early pandemic phase, possibly resulting from low testing capacity, slow contact tracing, uncertainty, and lack of knowledge about the disease and protocols, or fear of social stigma. This changed during period 2, where the estimated reporting ratio went up three times. These results are consistent with a study on time-varying under-reporting estimates in various countries, including the Philippines, during the same period ^[61]. The impact of the community quarantines imposed by the government is reflected in the reduction parameter $\delta$ and transmission rates $\beta$ and $\beta_F$. The estimated 20% reduction in transmission for the behavior-changed class compares with a mathematical model of COVID-19 transmission in the Philippines which showed that the minimum health standards reduced the probability of transmission per contact by $13-27\%$ ^[62]. As the community quarantine was relaxed from period 1 to 2, the transmission rate of the disease ($\beta$) increased from period 1 to 2, while the rate of behavior change or transmission rate of awareness or fear of COVID-19 ($\beta_F$) decreased from period 1 to 2.

Results in Figures 5 and 6 emphasize the importance of early public health campaigns and positive behavior changes (e.g., mask-wearing, improved hygiene practices, and social distancing) on reducing and delaying the peak of infections. Although fear and stigma can influence behavior changes ^[63,64], these can also affect reporting and negatively impact disease control. We see in Figure 6 that as reporting increased, total cumulative cases including the unreported decreased significantly.

Without vaccines and antiviral therapy, control of epidemic diseases relies on effective NPIs and the management of healthcare systems. The bi-objective optimization approach can be used as a decision support tool because of the multiple optimal solutions provided by the method, wherein policymakers can choose depending on how much priority is given to minimizing the number of additional beds or earlier easing of restrictions. Although the method is applied to the Philippines, the optimization approach can also be used by other cities or countries by adapting location-specific epidemiological parameters. For countries with limited resources, the solutions corresponding to later easing of restrictions and a smaller number of additional beds may be better options. On the other hand, for those that can provide sufficient additional beds, the approach can be used to identify the optimal timing of adjusting NPIs. Results in Figure 7 suggest that if Metro Manila eased to GCQ on June 1, 2020, at least 47 beds per day should have been prepared so that the bed occupancy in the capital did not reach critical or high-risk, and MECQ was not needed to be reimposed. The blue solutions in Figure 7 prioritizes the earlier timing ($\tau$) of easing protocols over the number of additional beds ($\omega$). With this policy, GCQ could have been started 30 days earlier. However, this requires 131 additional beds per day, which is 4 times $\omega_\text{data}$. On the other hand, the yellow solutions correspond to implementing GCQ on May 29, 2020. This would require 47 additional beds per day, which is still more than $\omega_\text{data}$. In Panel (D), the policy in purple has the lowest peak among all the Pareto solutions. For this policy, even though GCQ starts on May 20 (12 days earlier than what happened), the peak of cases (purple curve in Panel (C)) was 500 less than the peak from the model (black curve in (A)). This policy would have required 56 beds per day, which is almost double than $\omega_\text{data}$.

5.   Conclusions

In this work, we used an SEIQR model that considers behavior change and underreporting to study the spread of COVID-19 during the early phase of the pandemic in Metro Manila, Philippines. Behavior change can be influenced by awareness or fear of the disease, and willingness to observe NPIs such as social distancing and mask-wearing. It was incorporated into the model by introducing a two-compartment susceptible population: one for the behavior-changed population and the other for those with regular behavior. The probability of getting infected is reduced for the behavior-changed susceptible class. Due to limited testing and tracing, or negative attitudes of people towards seeking healthcare, a compartment for the unreported cases was also added.

The results of this study highlight the importance of early behavior change and a high reporting rate in reducing the number of cases and delaying the peak of infections. These can be done by intensifying case surveillance and public health campaigns promoting compliance with NPIs, seeking healthcare, and discouraging social stigma. Moreover, this study provides an optimization approach that quantifies the additional bed requirement when policies are eased. The approach can be helpful in planning strategies that address strengthening or easing of policies, especially during the early phase when NPIs were the only control measures. Although the study focused on the transmission of COVID-19 in the Philippines, the proposed model is general enough that it can be applied to any city or country. The optimization problem can also be applied to other disease outbreaks by adjusting key epidemiological parameters.

The model is best suited to represent the early phase of an epidemic. So, for a future work, the model can be extended to describe succeeding epidemic waves and incorporate vaccination, variants, and NPIs. Since the simulations did not consider the economic impact of NPIs, a model that maximizes the economic output of a city or country while minimizing the number of infections is another work that can be pursued. Moreover, one can also modify the model using fractional derivatives to account for memory effects ^[65,66]. Recently, numerous works on the theoretical analysis and applications of fractional differential equations have been done ^{[67,68,69,70,71,72]}. To the best of our knowledge, this approach has not yet been applied to a behavior change model similar to the one utilized in this study. Although a compelling research direction, this generalization requires an intensive investigation and demands a separate study.

Acknowledgments

This paper is supported by the Korea National Research Foundation (NRF) grant funded by the Korean government (MEST) (NRF-2021M3E5E308120711). This paper is also supported by the Korea National Research Foundation (NRF) grant funded by the Korean government (MEST) (NRF-2021R1A2C100448711). The authors acknowledge Dr. Peter Julian Cayton of the UP Resilience Institute for his assistance with the data collection.

Conflict of interest

All authors declare no conflicts of interest in this paper.