A Model-Based Strategy on COVID-19 Vaccine Roll-out in the Philippines ====================================================================== * Rey Audie S. Escosio * Olive R. Cawiding * Bryan S. Hernandez * Renier G. Mendoza * Victoria May P. Mendoza * Rhudaina Z. Mohammad * Carlene P.C. Pilar-Arceo * Pamela Kim N. Salonga * Fatima Lois E. Suarez * Polly W. Sy * Thomas Herald M. Vergara * Aurelio A. de los Reyes V ## Abstract Coronavirus disease 2019 (COVID-19) is an infectious disease caused by severe acute respiratory syndrome coronavirus 2. Millions of people have fallen sick, and some have died due to this affliction that has spread across the globe. The current pandemic has disrupted normal day-to-day human life, causing a profound social and economic burden. Vaccination is an important control measure that could significantly reduce the incidence of cases and mortality if properly and efficiently distributed. In this work, an age-structured model of COVID-19 transmission, incorporating an unreported infectious compartment, is developed. Three age groups are considered, namely: *young* (0-19 years), *adult* (20-64 years), and *elderly* (65+ years). The transmission and reporting rates are determined for each group by utilizing the number of COVID-19 cases in the National Capital Region in the Philippines. Optimal control theory is employed to identify the best vaccine allocation to different age groups. Further, three different vaccination periods are considered to reflect phases of vaccination priority groups: the first, second, and third account for the inoculation of the elderly, adult and elderly, and all three age groups, respectively. This study could guide in making informed decisions in mitigating a population-structured disease transmission under limited resources. Keywords * COVID-19 * mathematical modeling * aged-structured model * optimal control theory * vaccination ## 1 Introduction Coronavirus disease 2019 (COVID-19) is an infectious disease due to severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). Two years have passed since the World Health Organization characterized the COVID-19 outbreak as a pandemic [11], with over 468 million confirmed cases and just over 6 million reported deaths globally at the time of writing [10]. Its emergence has disrupted not only public health but also the global economy. Unfortunately, with the highly contagious Omicron variant driving the current increase of reported cases around the world [10], this pandemic is still not over. To analyze the transmission dynamics of this infectious disease, mathematical modeling has extensively been used and has played an essential role in predicting, assessing, and controlling outbreaks [44]. Numerous compartmental epidemiological models, such as SIR (susceptible-infectious-recovered), SIRS (susceptible-infectious-recovered-susceptible), and SEIR (susceptible-exposed-infectious-recovered) models, have been developed in the past decades, to represent the transmission of infectious diseases [16, 17, 25, 32]. Moreover, to properly identify interventions that can prevent outbreaks, one needs to take into account the social structure and mixing patterns, which vary across countries in different stages of development and with different demographics, in the epidemic model [39]. In the Philippines, early studies on COVID-19 revealed that while only eight percent of the population are in the age group 60 and over, they account for the majority of all COVID-19 deaths so far; and while individuals below age 60 dominate the reported cases, they account for only a third of total deaths [2, 12, 43]. Vaccine is a biological preparation that stimulates one’s immune response against a particular disease. Vaccine administration known as vaccination is an effective method designed to prevent an infectious disease [5]. To reduce the risk of contracting COVID-19, several vaccines have been approved worldwide for public use [1]. However, due to global demand and limited supply of vaccines against the virus, analysis of priority in the vaccine roll-out is vital. Several studies use optimal control theory to investigate compartmental models in epidemiology. In particular, Gaff and Schaefer (2009) focused on finding the optimal response balancing vaccination and treatment strategies that will minimize the number of infected individuals with variations of standard SIR, SIRS, and SEIR epidemiological models [24]. In addition, Hansen and Day (2011) extended their work on the optimal control of a basic SIR epidemic model and obtained policies for an isolation-only model, a vaccination-only model, and a combined isolation–vaccination model [26]. For the COVID-19-specific study, Kantner and Koprucki (2020) provided a study for the epidemics with purely non-pharmaceutical interventions based on an extended SEIR model and optimal control theory [29]. Moreover, Djidjou-Demasse et al. (2020) investigated optimal COVID-19 epidemic control via social distancing and lockdown measures until the deployment of vaccines [20]. Additionally, Obsu and Balcha (2020) used optimal control to minimize the number of exposed and infected populations that took into account the cost of implementation [36]. In the Philippine setting, Estadilla et al. (2021) considered the impact of vaccine supplies and delays on the optimal control of the COVID-19 pandemic [23]. Furthermore, Olivares and Staffetti (2021) considered optimal control to suggest a suitable schedule for vaccination and testing policies to control the spread of COVID-19 [37]. Finally, age-structured non-pharmaceutical interventions for optimal control of the COVID-19 epidemic were established by Richard et al. in the Burkina-Faso, France, and Vietnam settings [42]. In this work, we focus on the National Capital Region (NCR), which is the most densely populated region and the center of the economy of the Philippines [8]. Our study generally aims to develop an age-structured model that could describe the spread of COVID-19 in the NCR and identify effective control strategies based on optimal control theory. Specifically, this work aims to consider a dynamic model incorporating the effects of age-dependent susceptibility to COVID-19 infection in the NCR; obtain time-dependent optimal vaccination policies based on priority groups in different periods of vaccine roll-out; and assess the corresponding allocation of vaccines to each group. To the best of our knowledge, this is the first study that applies optimal control theory to a COVID-19 age-structured model in the Philippines. Since testing and contact tracing of suspected COVID-19 cases in the Philippines were lagging [27], our model incorporates possible unreporting in the different age groups. The capacity of daily vaccine administration is also considered to identify the prioritization of vaccination in different age groups so that the number of infections is minimized. ## 2 Methodology ### 2.1 Age-Structured Model In this work, a modified deterministic SEIR compartmental framework is used to model the progression of COVID-19 in the National Capital Region (NCR), Philippines. We consider the total population (*N*) as a sum of the four epidemiological classes: susceptibles (*S*), exposed (*E*), infected (*I*), and recovered (*R*). Individuals from the susceptible population can be exposed to the disease from interaction with infected individuals at a transmission rate *λ*. Upon acquiring infection, they move to the exposed class and stay in this stage for an average period of, ![Graphic][1] where *α* is the latency rate. In our model, we assume that the *incubation period*, which is the time period from acquiring the virus until developing symptoms, has the same duration as the *latent period*, the period from acquiring the virus until becoming infectious [13, 15, 31, 33, 40]. After this period, the infection worsens, the exposed individual develops infectivity, and is then moved to the infectious class. In this model, the infectious class is subdivided into two. First, the reported subclass *I*r is composed of individuals who are tested and recorded as virus-positive with reporting rate *ρ*. The second is the unreported subclass *I*u, which considers infectious individuals who are not reported due to the absence of testing brought about by socio-economic burden, under-reporting, and/or misreporting. Taking into account the possible transmission contributions of these two subclasses, we consider the force of infec-tion ![Graphic][2], where *β* is the transmission rate and *σ* is a multiplicative factor denoting relative infectiousness of *I*. Note that *I* includes both asymptomatic and symptomatic infectious individuals who are not in quarantine or isolation. After being infectious, an individual may either progress to the recovered class with a recovery rate *γ*, or die due to the disease with a disease-induced death rate *µ*. The model limits the progression of the disease to a single infection, i.e., recovered individuals are no longer infectious and are immune to reinfection. It has been observed that COVID-19 cases and risk severity increases with age [19,21]. Since population structure plays a significant role on the disease transmission dynamics [14], each epidemiological class is further divided into three age groups: *young* (0-19 years), *adult* (20-64 years), and *elderly* (65+ years). The grouping is consequential to the following protocols implemented by the Philippine government during the community quarantine placed into effect since March 2020: (1) the group of individuals below 20 years old is comprised mainly of students whose social mixing is impacted by the suspension of face-to-face classes, (2) the majority of the individuals aged 20-64 years old are considered as the working population who are authorized to attend to the basic needs of their families, and (3) the persons who are 65 years old and older are prohibited to go out of their residences, except for unavoidable circumstances [6]. A contact matrix *C* = (*C**ij*) denoting the interactions between an age group *i* with other age groups *j* is incorporated to account for the effect of mixing patterns in the disease transmission [39–41]. Using this, we modify the force of infection *λ**i* for a susceptible individual in age group *i* as follows ![Formula][3] where *β**i* is the transmission rate of age group *i*, the notations ![Graphic][4] are the reported infectious, unreported infectious, and total population of age group *j*, respectively, *σ**j* is a multiplicative factor relative to the transmission contribution of ![Graphic][5] in group *j*, and *i, j* = 1, 2, 3, correspond to each age group. Figure 1 schematically depicts the disease progression for an age group *i*. ![Figure 1:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2022/05/27/2022.05.27.22275675/F1.medium.gif) [Figure 1:](http://medrxiv.org/content/early/2022/05/27/2022.05.27.22275675/F1) Figure 1: Schematic diagram of the COVID-19 disease transmission model. Compartments in the model are classified into distinct epidemiological classes: Susceptible (*S*), Exposed (*E*), Reported Infectious (*I*r), Unreported Infectious (*I*u), and Recovered *R*. Each compartment is then stratified into three age groups: Young (0-19), Adult (20-64), and Elderly (65+). The subscript *i* denotes a particular age group. The disease transmission dynamics can be described by the following set of coupled ordinary differential equations ![Formula][6] where *λ**i* is the force of infection, *α**i* is the latency rate, *ρ**i* is the reporting rate, *µ**i* is the death rate, and *γ**i* is the recovery rate, for a particular age group *i*. Table 2.1 lists the set of parameters, their descriptions and units, including initial values for the state variables of the three different age groups. View this table: [Table 2.1:](http://medrxiv.org/content/early/2022/05/27/2022.05.27.22275675/T1) Table 2.1: Parameters and initial values for the different age groups. ### 2.2 Data Acquisition Age-specific case numbers compiled and reported by the Department of Health (DOH) in the Philippines are used in this study. Information from January 1, 2021 until February 28, 2021, on the incidence cases for young, adult, and elderly groups are extracted from the DOH Data Drop [3]. The study considers only the National Capital Region (NCR) as the area of interest due to its soaring number of cases, its high mobility and population density, and most importantly, its significant share of the vaccine allocation. The daily incidence cases for the three different age groups considered are shown in Figure 2(A-C). The biggest share of reported cases is among the adult group with 82.80%, followed by the young group with 9.27%, while the elderly group had the least contribution of only 7.93% of the total cumulative cases since January 1, 2021, as depicted in Figure 2(D). An increasing trend of cumulative cases can also be observed across the three age groups. Note that the adult group constitutes the biggest population share of about 59.16%, and the relative population shares of young and elderly groups are, 36.91% and 3.93%, respectively. However, in terms of infection rate, the elderly group has the highest susceptibility to infection with 3.4%, followed by the adult and young groups, with 2.3% and 0.3%, respectively (refer to Figure 2(E)). ![Figure 2:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2022/05/27/2022.05.27.22275675/F2.medium.gif) [Figure 2:](http://medrxiv.org/content/early/2022/05/27/2022.05.27.22275675/F2) Figure 2: (A-C) Daily reported cases for young (blue), adult (red), and elderly (yellow) group with the corresponding 7-day moving average (black curve). (D) Cumulative cases for the three different age groups from January 1 to February 28, 2021. Adult group comprises the bulk of COVID-19 cases. (E) Infection rate for elderly, adult, and young populations based on cumulative cases until February 28, 2021. The elderly group shows the highest susceptibility to infection. ### 2.3 Contact Matrices Transmission of an infectious disease such as COVID-19 is influenced by the social structure involving person-to-person interactions which vary by age and locations – at home, at the workplace, at school, or in the community [33, 41]. Our model considers the effect and contribution of the various social mixing in these different locations. Matrices for household (*D**h*), workplace (*D**w*), school (*D**s*), and community (*D**c*) contact can be obtained from [41] and are shown in Figure 3. The data consists of four 16 *×* 16 matrices binned in 15 five-year age levels (0 *−* 4, 5 *−* 9, …, 70 *−* 74) and a single age group with ages 75 and above. The total of the four matrices when added can represent the whole contact matrix *D* = (*D**kl*), i.e., ![Figure 3:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2022/05/27/2022.05.27.22275675/F3.medium.gif) [Figure 3:](http://medrxiv.org/content/early/2022/05/27/2022.05.27.22275675/F3) Figure 3: Baseline contact matrices in the Philippines during pre-COVID-19 pandemic showing strong interactions within the same age group in different locations. Figure is reproduced from [41]. For visualization, the shown matrices *D**†* are flipped upside down, i.e., *D**†* = 𝕁*D* where 𝕁 is the backward identity. ![Formula][7] The matrix entries *D**kl* show a pre-pandemic view of the number of contacts in the *k*th age group of a certain *j*th age group population in the Philippines. We reduce the size of these matrices from 16 *×* 16 to a mere 3 *×* 3 corresponding to our clumped age groups which are denoted by *C* = (*C**ij*). The process is detailed in Appendix A. Additionally, this 3 *×* 3 total contact matrix *C* = (*C**ij*) can be rewritten accordingly based on [30] to integrate the non-pharmaceutical interventions during the pandemic, i.e., ![Formula][8] where the coefficients *v**h*, *v**w*, *v**s*, and *v**c* correspond to the implemented contact interventions in house-holds, workplaces, schools, and communities, respectively. Note that these equate to one if no intervention is implemented, and equate to zero if interactions in such location are entirely halted. One obvious value is *v**s* = 0 for schools have been closed since the start of the community quarantine. We assign *v**h* = 1 since household members are maintained to be close contacts. On the other hand, the regulated workforce and community interactions can be assumed to be at 50% or *v**w* = *v**c* = 0.50. Applying these coefficients, we get the contact matrices in Figure 4. ![Figure 4:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2022/05/27/2022.05.27.22275675/F4.medium.gif) [Figure 4:](http://medrxiv.org/content/early/2022/05/27/2022.05.27.22275675/F4) Figure 4: Assumed mitigated 3 *×* 3 contact matrices in the Philippines considering the young, adult, and elderly population groups in total, household, workplace, school, and community context. For visualization, the shown matrices *C**†* are flipped upside down, i.e., *C**†* = 𝕁*C* where 𝕁 is the backward identity. The final contact matrix *C* = (*C**ij*) is shown in the left of Figure 4. For the young compartment, the first row of the contact matrix is [*C*1*J*] = [4.086, 4.135, 0.261] where *J* = 1, 2, 3. Even with the interventions brought about by government policies such as the closure of schools, this row contains high values due to continuous contact among the juvenile population in households and communities. As for adults, we have the second row [*C*2*J*] = [2.000, 5.888, 0.200] where *J* = 1, 2, 3. For this age group, close contacts are largely observed in the workplace and communities, even though only 50% of the population is allowed to interact. On the other hand, the elderly population has the least contact with all age groups, including itself, as seen in row [*C*3*J*] = [0.902, 2.145, 0.372] where *J* = 1, 2, 3. This can be attributed to prohibitions of elderly citizens from activities beyond their homes. ### 2.4 Parameter Estimation We used the population projection of DOH for 2021 as the initial population for the system (i.e., *N* = 13, 966, 223). This can be divided accordingly into the young (*N*1,0 = 5, 155, 544), adult (*N*2,0 = 8, 261, 835), and elderly (*N*3,0 = 548, 844) compartments using the available population fractions data. The remaining initial values for the compartments are seen in Table 2.1. From Table 2.1, the total *E* is approximated from the result of the 2020 COVID-19 model [22]. On the other hand, the initial reported infectious ![Graphic][9], deceased *D*, and recovered *R* populations are derived from the tallied data in [3] and divided according to the size of the age groups. The initial unreported compartment is approximated as ![Graphic][10], immediately after parameter estimation. As for the initial susceptible population *S*, it can be calculated as ![Graphic][11]. The multiplicative factor *σ* for unreported infections is held at 2 for all age groups. This presumes that undetected infectious interactions are twice those of the detected. The latency rate *α* is assumed to be 1*/*5, corresponding to 5 days of disease progression from infection to the infectious state. The death rate *µ* and recovery rate *γ* are calculated from the data in [3] since dates of confirmation, death, and recovery are available. Lastly, the transmission rate *β* and reporting rate *ρ* will be identified via parameter estimation. The process of parameter estimation aims to identify selected parameters as a result of minimizing an objective function. In this paper, we use the sum of the squared difference between the selected data and the solution of the model as the objective function to identify the transmission rate *β* and reporting rate *ρ* of each age group. The data used are the daily reported cases from [3] during the dates of January 1 until February 28, 2021, dubbed the pre-vaccination period. For the objective function, we use the following three arbitrary differential equations ![Formula][12] where ![Graphic][13] denotes the cumulative new reported cases for age group *i*. The differences of these solutions are the reported daily new cases. To incorporate the probable late reporting of cases, we consider the 7-day moving average as the base data. A derivative-free simplex method is used which is capable of constrained and multivariate optimization. Latin hypercube sampling [28, 35] is also done to better sample the initial guesses of the method. Upon choosing the best fit among the estimated parameters, bootstrapping [18] is performed by redoing the parameter estimation multiple times with the best model fit as the data used including Poisson noise. By increasing the number of realizations, the estimated parameters may change and a confidence interval will be determined. ### 2.5 Optimal Control Problem For the control problem in this study, we pay particular attention to vaccination efforts. Three controls are considered, each of which corresponds to the vaccination efforts for each age group. Specifically, the controls *u*1(*t*), *u*2(*t*), *u*3(*t*) denote the vaccination effort for the young, adult, and elderly age group, respectively. Incorporating these controls in Model (2.2) yields the following *disease-controlled model* ![Formula][14] where *θ**i* denote vaccine efficacy for age group *i*. Note that an assumption of this system is that only susceptible individuals are candidates for vaccination. Here, the Pontryagin Maximum Principle (PMP) is employed to formulate the optimal control problem [34]. The objective is to minimize the number of infected individuals in ![Graphic][15] and ![Graphic][16] across different age groups with minimum vaccination administrative cost. This yields the following objective functional ![Formula][17] where [*t*, *t**f*] refers to the period of vaccination effort, and for *i* = 1, 2, 3, *B**i*’s are weight constants representing the costs required to administer the corresponding control measures. Note that the controls are expressed as quadratic functions to take into account the nonlinear costs for the implementation of each control. Utilizing PMP leads to minimizing the objective functional in Equation (2.5) subject to the constraints given by Model (2.4). We want to identify the optimal controls ![Graphic][18], and ![Graphic][19] such that ![Formula][20] where ![Formula][21] The upper (*b**i*) and lower (*a**i*) bounds for each control represent the maximum and minimum implementation efforts, respectively. If *u**i* = 0, then it suggests that no vaccination effort is employed, while *u**i* = 1 denotes that maximum effort is exerted for vaccine administration in the corresponding age group *i*. A theorem stating the existence of the optimal controls and corresponding states is given in Appendix B. ## 3 Results and Discussion ### 3.1 Estimated Transmission and Reporting Rates To complete the parameters in Table 2.1, we use parameter estimation using new cases from January 1 up to February 28, 2021, to get the transmission and reporting rates for each of the three age groups. Figure 5 shows the best fit of the model compared to the data on cumulative cases for the young, adult, and elderly groups. From the parameter bootstrapping, the 95% confidence interval is also shown as bands around the best model fit. We close up as well on the cumulative data and model solution at the endpoint on February 28, 2021. ![Figure 5:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2022/05/27/2022.05.27.22275675/F5.medium.gif) [Figure 5:](http://medrxiv.org/content/early/2022/05/27/2022.05.27.22275675/F5) Figure 5: Comparison between raw data (black dots), 7-day moving average data (black line), and model output with 95% confidence interval for young (blue), adult (red), and elderly (yellow) compartments for cumulative data. Close-up (violet) of the cumulative data on the last day, February 28, 2021, is also shown for comparison between the data and the model. Meanwhile, the estimated parameters and their mean and median values from the 1000 bootstrap realizations are presented in Table 3.1 and shown as distributions in Figure 6. Transmission is highest among the elderly with *β*3 = 0.0288 and least among the young age group with *β*1 = 0.0057. The mean reporting rate for the young, adult, and elderly are 22.36%, 70.14%, and 96.06%, respectively. Since the elderly are more likely to experience severe symptoms than the younger age groups, then we observe a higher reporting rate in the elderly group. Meanwhile, there is possible under-reporting in children and teenagers since testing resources are primarily given to priority groups such as front-line workers and symptomatic patients. View this table: [Table 3.1:](http://medrxiv.org/content/early/2022/05/27/2022.05.27.22275675/T2) Table 3.1: Best parameter estimates for the transmission and reporting rates using the least-squares method, and the mean and median of the 1000 re-estimates in bootstrapping ![Figure 6:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2022/05/27/2022.05.27.22275675/F6.medium.gif) [Figure 6:](http://medrxiv.org/content/early/2022/05/27/2022.05.27.22275675/F6) Figure 6: Histogram of parameter estimates for each age group after 1000 realizations from bootstrapping with their respective mean (red, dashed) and median (black, solid) values. ### 3.2 Optimal Vaccination Strategies The main strategy used in this study is adding vaccination efforts according to the age-dependent vaccine prioritization imposed by the government. We use optimal controls for each age group to represent this pharmaceutical intervention. These optimal vaccination strategies are obtained by solving Model (2.4) with vaccine efficacy set to *θ**i* = 90% for all age groups, *i* = 1, 2, 3, and with weight parameters *B*1 = 10*−*6, *B*2 = 10*−*6, and *B*3 = 10*−*5, shown in Table 3.2. View this table: [Table 3.2:](http://medrxiv.org/content/early/2022/05/27/2022.05.27.22275675/T3) Table 3.2: Values used for the numerical simulations of the optimal control problem for each time period and specific to each age group. A one-year period from March 1, 2021, the official starting date of vaccination in the Philippines, is chosen as the time window for the optimal vaccination strategies. We divide this vaccination period into three: first, March 1, 2021, to June 6, 2021, where the vaccine inoculation is only administered to the elderly [9]; second, June 7, 2021, to October 14, 2021, where the adult population, most of whom are in the working class, started to be inoculated [4]; and last, October 15, 2021, to February 28, 2022, where the young population was allowed to be vaccinated [7]. Furthermore, we assume that once vaccination becomes available to a particular age group, the maximum vaccination rate for the first period of implementation is low due to the adjustments in implementation, and speeds up in the next vaccination period. It is also assumed that in the third period, vaccine supply is sufficient so that all control parameters for vaccination can be increased to 0.1. Specifically, for the elderly population, the maximum vaccination rate in the first period is set at 0.01, increased to 0.03 in the second period, and increased further in the third period to 0.1. For the adult population, it is set at zero in the first period as no vaccination is available yet, then it is set at 0.01 in the second period, and then increased to 0.1 in the third period. Lastly, for the young population, the maximum vaccination rate is set at zero in the first two periods and is increased to 0.1 in the last period. Once vaccination becomes available to a certain age group, implementation is continuous even if they are not the priority group. Then, the lower bounds for the vaccination efforts are all set at a relatively lower value of 0.001. These additional values used for the optimal control problem are summarized in Table 3.2. Figure 7 shows the optimal vaccination strategy obtained for all the age groups. The daily scheduled vaccination effort for the elderly (*u*3) shall be maintained at the upper bound for Periods 1 and 2. In particular, the optimal control chooses the maximum vaccination effort to prevent the elderly and the adult groups from getting infected. It suggests that efficient and effective implementation of vaccine administration, exhausting all possible limited resources. In Period 3, more vaccine supply is available, vaccine administration is easier, and vaccination hesitancy is lower than the previous periods. However, the results depict that the vaccination efforts for all the groups should be at the maximum for the first few months and must be sustained for the adult and young groups throughout the considered inoculation duration. It is more likely that the target number of elderly to be inoculated is achieved around the midway of Period 3. Thus, the control effort for the elderly can be decreased and then re-allocated to the other groups. ![Figure 7:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2022/05/27/2022.05.27.22275675/F7.medium.gif) [Figure 7:](http://medrxiv.org/content/early/2022/05/27/2022.05.27.22275675/F7) Figure 7: Control profiles for the vaccination of the elderly, adult, and young populations for Period 1 (March 1, 2021 - June 6, 2021), Period 2 (June 7, 2021 - October 14, 2021), and Period 3 (October 15, 2021 - February 28, 2022). For the daily scheduled vaccination effort for the adults (*u*2), once they are allowed to be vaccinated (from June onward), a maximum daily vaccination effort shall be scheduled in Period 2. The maximum effort shall remain for Period 3 until before the end of February when vaccination can be lowered. Finally, the daily scheduled vaccination effort for the younger population (*u*1) shall also be at maximum, once they are allowed to be vaccinated (from October 15 onward). Vaccination can be decreased towards the end of February, a little later than the lowering of efforts for the adult population. Generally, throughout the course of the first two periods, in order to reach optimal results, the vaccination efforts must be set to the maximum value. For the third period, the control efforts can be lowered towards the end with a relatively earlier lowering of vaccination efforts for the elderly population which can be attributed to their comparatively smaller population. Figure 8 shows that without intervention, the projected number of daily infected cases can reach up to 21, 208 for the elderly population, 317, 640 for the adult population, and 128, 509 for the young population. Meanwhile, the optimal vaccination strategies result in a significant decrease in the number of infected individuals across all age groups. The controls result in a peak of 63, 753 cases for the adult age group, 28, 177 for the young age group, and 1, 537 for the elderly age group. It can also be observed that even without vaccination in the period June to October, the age group 0 *−* 19 exhibits a significant decrease in the infected population. This means that the vaccination of the adult and elderly populations have a great impact on reducing the number of infections in the younger population. It is observed from Figure 8 (left panel) that the number of daily incidence cases does not decrease immediately within the initial phase of vaccination roll-out, e.g., first and second period for elderly and adult, respectively. This can be attributed to the delayed immunization effect of the vaccines. Though the young individuals are vaccinated in the third period, the incidence cases are already decreasing due to the vaccination of elderly and adult population in the previous periods. ![Figure 8:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2022/05/27/2022.05.27.22275675/F8.medium.gif) [Figure 8:](http://medrxiv.org/content/early/2022/05/27/2022.05.27.22275675/F8) Figure 8: Total daily active cases without (dashed) and with (solid) optimal control strategies for both reported and unreported compartments for young, adult, and elderly populations. Figure 9 shows a comparison of the infected populations with and without control. Even with control, the largest contributor to the total infections is the adult age group which is followed by the young population and then the elderly population. This is intuitive since the adult population comprises the bulk of the total population, and also the most mobile and interacting group. The inoculation of the elderly population during the first period led to a decrease in cumulative cases for all age groups, ranging from 41% for the young, 43% for the adult, and 49% for the elderly. When adults were allowed to get vaccines during the second period, a massive decrease in infections can be seen from uncontrolled to the controlled system. The percent decreases are 397% for the elderly, 251% for the adult, and 254% for the young. However, inoculation for the young population in the third vaccination period decreased this even further wherein a 392% drop can be observed for the young, 440% for the adult, and 897% for the elderly. ![Figure 9:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2022/05/27/2022.05.27.22275675/F9.medium.gif) [Figure 9:](http://medrxiv.org/content/early/2022/05/27/2022.05.27.22275675/F9) Figure 9: Upper panel shows the cumulative new reported COVID-19 cases without control (left panel, dashed) and with control (right panel, solid) for the young (blue), adult (red), and elderly (yellow). Lower panel depicts the cumulative cases at the end of each vaccination period and their respective percent decreases. ## 4 Conclusion In this paper, we developed an age-structured compartmental model for COVID-19 with an unreported infectious class. The epidemiological compartments are divided into three groups: *young* (0-19 years), *adult* (20-64 years), and *elderly* (65+ years), according to vaccine prioritization. We denote these groups by *i* = 1, 2, and 3, respectively. A contact matrix is used to incorporate the interaction among individuals across various age groups. Specifically, each entry of the matrix takes into account the contributions from the number of daily contacts between the considered age group *i* with other age groups *j* in household, workplace, school, and community. Furthermore, we introduced the force of infection denoted by *λ**i* for a susceptible individual for the *i*th age group. It is composed of two factors. The first one is the transmission rate for this particular age group. The second one is the sum of the product of the associated entry of the contact matrix and the proportion of the infected cases over all three age groups *j* = 1, 2, 3, for the given specific age group *i*. Moreover, the unreported cases are multiplied by a factor denoted by *σ**i*, which is the relative contribution to transmission of the unreported class for the associated age group. Most of the model parameters are obtained from the literature or aggregated from available data. Reporting and transmission rates are estimated because they are country/region-specific and are not readily available. Estimation was performed using a least-squares formulation by fitting the model to the cumulative case data. Results show that both reporting and transmission rates are higher in the older age group. Our simulation suggests that thousands of cases in young and adult populations were undetected. Uncertainty analysis showed the reliability of the estimates. Given the limited availability and administration of the vaccines in the first half of the year 2021, we analyzed the effect of age-targeted vaccine distribution by comparing the model output when no intervention is imposed, against the output when vaccination is introduced in the model. The results of the numerical simulations for the controls show that efforts from the government in implementing its vaccination drive in the country should always be at maximum until towards the end of the last period when vaccination can be lowered initially for the elderly population first, then for the adult population, and lastly, for the young population. Moreover, implementing the vaccine strategies as early as possible can reduce the cumulative reported cases in the young, adult, and elderly population by 392%, 440%, and 897%, respectively. The results from the optimal control problem can be attributed to the proportion of each age group to the total population, the transmission rate specific to each age group, and the order of prioritization in vaccination. Our model took into account age structure and unreporting. The model was applied when the vaccination program started, assuming that one vaccine dose can already provide immunity. Further model iterations are needed to include the impact of the succeeding shots and the effect of waning immunity. At the beginning of the vaccination period, only one variant was dominant in the Philippines so variants were also not incorporated into the model. Nevertheless, integrating these into the model is an exciting research direction, but would require other vaccination and breakthrough infection data and would demand a separate study. Ultimately, the aim of this study is to guide the policymakers in crafting strategies that can mitigate the spread of an infectious disease considering age-dependent transmission and limited vaccine supply. Although the model is tested using Philippine data, the model is general enough that it can be applied to other regions with similar situations. ## Data Availability All data produced in the present study are available upon request to the authors. ## Acknowledgement This project was supported by a grant from the University of the Philippines Computational Research Laboratory. ## Appendix A. Computation of the 3 *×* 3 Contact Matrices The four 16 *×* 16 pre-pandemic contact matrices for different interactions and their total can be reduced to a size that clumps together the five-year age levels into three main age groups for the young, adult, and elderly populations. Let *C* be the corresponding 3 *×*3 total contact matrix from *D*. The superscripts also follow for reduced contact matrices for different interactions, e.g., *C**h* is the reduced *D**h*. The young age group contains ages from 0 to 19 that correspond to the first four age levels. The adult age group covers the most age levels which are from 20 to 64. Lastly, ages 65 and above will be classified under the elderly age group. From here, we can determine sets containing the indices of each of these age groups, i.e., ![Formula][22] where *X*1 is for the young age group, *X*2 is for the adult age group, and *X*3 is for the elderly age group. This can be observed on Table .1. View this table: [Table 1:](http://medrxiv.org/content/early/2022/05/27/2022.05.27.22275675/T4) Table 1: Age groups and their respective population fractions based on the 2015 NCR Population. Using the latest available data for the National Capital Region from the Philippine Statistics Authority [8], we can create a 16 *×* 1 vector *d* = (*d**k*) that contains the population fraction *d**k* of each of the age level *k* for *D*. We determine its elements by dividing the age group population over the total population. The values of population fraction *d**k* for *k* = 1, 2, …, 16 are shown in Table .1. The formula of obtaining *C* = (*C**ij*) is ![Formula][23] where *i, j ∈ {*1, 2, 3*}* correspond to the young, adult, and elderly populations, respectively. For example, if we want to get *C*13 from the total contact matrix *D*, which is the contact value of the young population relative to the elderly population, we consider the indices *i ∈ X*1 and *j ∈ X*3. The denominator simply adds up the population fraction of the young population to get 0.3959. Meanwhile, the numerator considers the values for the 4 *×* 3 submatrix of *D* with row indices from *X*1 and column indices from *X*3. Each row is summed and multiplied by its respective population fraction. For this example, we will get the numerator as 0.2110. Hence, the value of *C*13 is 0.533. This can be observed as the ![Graphic][24] of the first matrix plot in Figure 10. ![Figure 10:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2022/05/27/2022.05.27.22275675/F10.medium.gif) [Figure 10:](http://medrxiv.org/content/early/2022/05/27/2022.05.27.22275675/F10) Figure 10: The obtained 3× 3 contact matrices considering the young, adult, and elderly population groups in total, household, workplace, school, and community context. For visualization, the shown matrices *C**†* are flipped upside down, i.e., *C**†* = 𝕁*C* where𝕁is the backward identity. ## Appendix B. Existence of Optimal Control and Adjoint Equations Theorem .1. *There exist optimal control* ![Graphic][25] *and corresponding optimal state* ![Graphic][26] *for i* = 1, 2, 3, *which minimize the objective functional* (2.5) *over all controls in* (2.7). *Given this optimal pair of solution* (**u, X**), *there exist adjoint functions ψ**k*, *k* = 1, 2, …, 15 *such that for the Hamiltonian given by* ![Formula][27] *we have* ![Formula][28] *with transversality conditions* ![Formula][29] *Furthermore*, ![Formula][30] *Proof*. The optimal control ![Graphic][31] satisfying Equation (2.6) exists due to the convexity of the integrand in Equation (2.5). We use Pontryagin’s Maximum Principle [38] to obtain the adjoint equations and transversality conditions. Differentiating the Hamiltonian *H* with respect to the state variables gives us the following system of equations: ![Formula][32] resulting to System (.2) with *ψ**k*(*t**f*) = 0 for *k* = 1, 2, …, 15. Moreover, we derive ![Graphic][33] on *U* by solving the optimality conditions ![Graphic][34] for *i* = 1, 2, 3, i.e., ![Formula][35] Evaluating at ![Graphic][36], and ![Graphic][37] on 𝒰, we get ![Formula][38] and by considering the upper (*b**i*) and lower (*a**i*) bounds for each control, we end up with Equations (.3). * Received May 27, 2022. * Revision received May 27, 2022. * Accepted May 27, 2022. * © 2022, Posted by Cold Spring Harbor Laboratory This pre-print is available under a Creative Commons License (Attribution-NonCommercial-NoDerivs 4.0 International), CC BY-NC-ND 4.0, as described at [http://creativecommons.org/licenses/by-nc-nd/4.0/](http://creativecommons.org/licenses/by-nc-nd/4.0/) ## References 1. [1].10 vaccines granted emergency use listing (EUL) by WHO. [https://covid19.trackvaccines.org/agency/who](https://covid19.trackvaccines.org/agency/who). 2. [2].The changing demographics of COVID-19 infections and deaths in the Philippines: how age-sex structure, living arrangement, and family ties intersect. [https://www.drdf.org.ph/research/covid-19/rb7](https://www.drdf.org.ph/research/covid-19/rb7). 3. [3].Department of Health (DOH) COVID-19 tracker. Available online: [https://doh.gov.ph/covid-19/case-tracker](https://doh.gov.ph/covid-19/case-tracker). 4. [4].Gov’t to start COVID-19 vaccination on A4 priority group on June 7. Available online:. [https://doh.gov.ph/node/29928](https://doh.gov.ph/node/29928). 5. [5].Immunization: The basics.Available online:. [https://www.cdc.gov/vaccines/vac-gen/imz-basics.htm](https://www.cdc.gov/vaccines/vac-gen/imz-basics.htm). 6. [6].OMNIBUS GUIDELINES ON THE IMPLEMENTATION OF COMMUNITY QUAR-ANTINE IN THE PHILIPPINES with amendments as of December 14, 2020. [https://doh.gov.ph/sites/default/files/health-update/20201214-omnibus-guidelines-on-the-implementation-of-community-quarantine-in-the-Philippines.pdf](https://doh.gov.ph/sites/default/files/health-update/20201214-omnibus-guidelines-on-the-implementation-of-community-quarantine-in-the-Philippines.pdf). 7. [7].Pediatric vaccination roll-out for 12 to 17 y.o. children with comorbidities kicks off. Available online:. [https://doh.gov.ph/press-release/PEDIATRIC-VACCINATION-ROLL-OUT-FOR-12-TO-17-Y.O.-CHILDREN-WITH-COMORBIDITIES-KICKS-OFF](https://doh.gov.ph/press-release/PEDIATRIC-VACCINATION-ROLL-OUT-FOR-12-TO-17-Y.O.-CHILDREN-WITH-COMORBIDITIES-KICKS-OFF). 8. [8].Philippine Statistics Authority. Available online: [https://psa.gov.ph](https://psa.gov.ph). 9. [9].TIMELINE: The Philippines’ 2021 COVID-19 vaccine plan. Available online:. [https://www.rappler.com/newsbreak/iq/timeline-philippines-2021-covid-19-vaccination-plan](https://www.rappler.com/newsbreak/iq/timeline-philippines-2021-covid-19-vaccination-plan). 10. [10].Weekly epidemiological update on COVID-19 - 22 March 2022. [https://www.who.int/publications/m/item/weekly-epidemiological-update-on-covid-19\---|22-march-2022](https://www.who.int/publications/m/item/weekly-epidemiological-update-on-covid-19\---|22-march-2022). Accessed: 2022-03-25. 11. [11].World Health Organization declares COVID-19 a ‘pandemic.’ here’s what that means. Available online: [https://time.com/5791661/who-coronavirus-pandemic-declaration](https://time.com/5791661/who-coronavirus-pandemic-declaration). 12. [12]. K. A. Agrupis, C. Smith, S. Suzuki, A. M. Villanueva, K. Ariyoshi, R. Solante, E. F. Telan, K. A. Estrada, A. C. Uichanco, J. Sagurit, J. Calayo, D. Umipig, Z. dela Merced, F. Villarama, E. Dimaano, J. B. Villarama, E. Lopez, and A. R. Sayo, Epidemiological and clinical characteristics of the first 500 confirmed COVID-19 inpatients in a tertiary infectious disease referral hospital in Manila, Philippines, Tropical Medicine and Health, 49 (2021). 13. [13]. M. Aguiar, E. M. Ortuondo, J. B. Van-Dierdonck, J. Mar, and N. Stollenwerk, Modelling COVID 19 in the Basque Country from introduction to control measure response, Scientific Reports, 10 (2020). 14. [14]. R. M. Anderson and R. May, Infectious Diseases of Humans: Dynamics and Control, Oxford University Press, 1991. 15. [15]. K. Blyuss and Y. Kyrychko, Effects of latency and age structure on the dynamics and containment of COVID-19, Journal of Theoretical Biology, 513 (2021), p. 110587. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.jtbi.2021.110587&link_type=DOI) 16. [16]. S. Bowong and J. Kurths, Modelling tuberculosis and hepatitis B co-infections, Mathematical Modelling of Natural Phenomena, 5 (2010), pp. 196–242. 17. [17]. F. Chamchod and N. F. Britton, On the dynamics of a two-strain influenza model with isolation, Mathematical Modelling of Natural Phenomena, 7 (2012), pp. 49–61. 18. [18]. G. Chowell, Fitting dynamic models to epidemic outbreaks with quantified uncertainty: A primer for parameter uncertainty, identifiability, and forecasts, Infectious Disease Modelling, 2 (2017), pp. 379–398. 19. [19]. N. G. Davies, P. Klepac, Y. Liu, K. Prem, M. Jit, C. C.-. working group, and R. M. Eggo, Age-dependent effects in the transmission and control of COVID-19 epidemics, Nature Medicine, 26 (2020), p. 1205–1211. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41591-020-0962-9&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=32546824&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F05%2F27%2F2022.05.27.22275675.atom) 20. [20]. R. Djidjou-Demasse, Y. Michalakis, M. Choisy, M. T. Sofonea, and S. Alizon, Optimal COVID-19 epidemic control until vaccine deployment, Medrxiv, (2020). 21. [21]. C. Eastin and T. Eastin, Epidemiological characteristics of 2143 pediatric patients with 2019 coronavirus disease in China, The Journal of Emergency Medicine, 58 (2020), pp. 712–713. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.jemermed.2020.04.006&link_type=DOI) 22. [22]. R. A. Escosio, E. V. G. Rosero, J. P. Viernes, R. R. Oryan, R. A. Joel Addawe, C. Arceo, V. M. Mendoza, M. A. Mata, and A. de los Reyes V, Mathematical modelling of COVID-19 spread in various regions in the Philippines. Work in progress. 23. [23]. C. D. S. Estadilla, J. Uyheng, E. P. de Lara-Tuprio, T. R. Teng, J. M. R. Macalalag, and M. R. J. E. Estuar, Impact of vaccine supplies and delays on optimal control of the COVID-19 pandemic: mapping interventions for the Philippines, Infect. Dis. Poverty, 10 (2021), p. 107. 24. [24]. H. Gaff and E. Schaefer, Optimal control applied to vaccination and treatment strategies for various epidemiological models, Math. Biosci. Eng., 6 (2009), pp. 469–492. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.3934/mbe.2009.6.469&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=19566121&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F05%2F27%2F2022.05.27.22275675.atom) 25. [25]. S. M. Garba and U. A. Danbaba, Modeling the effect of temperature variability on malaria control strategies, Mathematical Modelling of Natural Phenomena, 15 (2020), p. 65. 26. [26]. E. Hansen and T. Day, Optimal control of epidemics with limited resources, J. Math. Biol., 62 (2011), p. 423–451. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1007/s00285-010-0341-0&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=20407775&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F05%2F27%2F2022.05.27.22275675.atom) 27. [27]. K. Hapal, The Philippines’ COVID-19 response: Securitising the pandemic and disciplining the pasaway, Journal of Current Southeast Asian Affairs, 40 (2021), pp. 224–244. 28. [28]. R. L. Iman and J. C. Helton, An investigation of uncertainty and sensitivity analysis techniques for computer models, Risk Analysis, 8 (1988), pp. 71–90. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1111/j.1539-6924.1988.tb01155.x&link_type=DOI) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=A1988N161900012&link_type=ISI) 29. [29]. M. Kantner and T. Koprucki, Beyond just “flattening the curve”: Optimal control of epidemics with purely non-pharmaceutical interventions, Journal of Mathematics in Industry, 10 (2020), p. 23. 30. [30]. M. Kimathi, S. Mwalili, V. Ojiambo, and D. K. Gathungu, Age-structured model for COVID-19: Effectiveness of social distancing and contact reduction in Kenya, Infectious Disease Modelling, 6 (2021), pp. 15–23. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.idm.2020.10.012&link_type=DOI) 31. [31]. A. J. Kucharski, T. W. Russell, C. Diamond, Y. Liu, J. Edmunds, S. Funk, R. M. Eggo, F. Sun, M. Jit, J. D. Munday, N. Davies, A. Gimma, K. van Zandvoort, H. Gibbs, J. Hellewell, C. I. Jarvis, S. Clifford, B. J. Quilty, N. I. Bosse, S. Abbott, P. Klepac, and S. Flasche, Early dynamics of transmission and control of COVID-19: a mathematical modelling study, The Lancet Infectious Diseases, 20 (2020), pp. 553–558. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/s1473-3099(20)30144-4&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F05%2F27%2F2022.05.27.22275675.atom) 32. [32]. S. Kumar, R. Chauhan, J. Singh, and D. Kumar, A computational study of transmission dynamics for dengue fever with a fractional approach, Mathematical Modelling of Natural Phenomena, 16 (2021), p. 48. 33. [33]. Y. N. Kyrychko, K. B. Blyuss, and I. Brovchenko, Mathematical modelling of the dynamics and containment of COVID-19 in Ukraine, Scientific Reports, 10 (2020). 34. [34]. S. Lenhart and J. T. Workman, Optimal control applied to biological models, Chapman and Hall/CRC, 2007. 35. [35]. M. D. McKay, R. J. Beckman, and W. J. Conover, A comparison of three methods for selecting values of input variables in the analysis of output from a computer code, Technometrics, 21 (1979), pp. 239–245. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.2307/1268522&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=WOS:A1979GW1&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F05%2F27%2F2022.05.27.22275675.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=A1979GW16700012&link_type=ISI) 36. [36]. L. L. Obsu and S. F. Balcha, Optimal control strategies for the transmission risk of COVID-19, Journal of Biological Dynamics, 14 (2020), pp. 590–607. [PubMed](http://medrxiv.org/lookup/external-ref?access_num=32696723&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F05%2F27%2F2022.05.27.22275675.atom) 37. [37]. A. Olivares and E. Staffetti, Optimal control-based vaccination and testing strategies for COVID-19, Computer Methods and Programs in Biomedicine, 211 (2021), p. 106411. 38. [38]. L. S. Pontryagin, Mathematical theory of optimal processes, CRC press, 1987. 39. [39]. K. Prem, A. R. Cook, and M. Jit, Projecting social contact matrices in 152 countries using contact surveys and demographic data, PLoS computational biology, 13 (2017), p. e1005697. 40. [40]. K. Prem, Y. Liu, T. W. Russell, A. J. Kucharski, R. M. Eggo, N. Davies, S. Flasche, S. Clifford, C. A. B. Pearson, J. D. Munday, S. Abbott, H. Gibbs, A. Rosello, B. J. Quilty, T. Jombart, F. Sun, C. Diamond, A. Gimma, K. van Zandvoort, S. Funk, C. I. Jarvis, W. J. Edmunds, N. I. Bosse, J. Hellewell, M. Jit, and P. Klepac, The effect of control strategies to reduce social mixing on outcomes of the COVID-19 epidemic in Wuhan, China: a modelling study, The Lancet Public Health, 5 (2020), pp. e261–e270. 41. [41]. K. Prem, K. v. Zandvoort, P. Klepac, R. M. Eggo, N. G. Davies, C. for the Mathematical Modelling of Infectious Diseases COVID-19 Working Group, A. R. Cook, and M. Jit, Projecting contact matrices in 177 geographical regions: an update and comparison with empirical data for the COVID-19 era, PLoS computational biology, 17 (2021), p. e1009098. 42. [42]. Q. Richard, S. Alizon, M. Choisy, M. T. Sofonea, and R. Djidjou-Demasse, Age-structured non-pharmaceutical interventions for optimal control of COVID-19 epidemic, PLoS Comput. Biol., 17 (2021), p. e1008776. 43. [43]. E. P. Salva, J. B. Villarama, E. B. Lopez, A. R. Sayo, A. M. G. Villanueva, T. Edwards, S. M. Han, S. Suzuki, X. Seposo, K. Ariyoshi, and C. Smith, Epidemiological and clinical characteristics of patients with suspected COVID-19 admitted in Metro Manila, Philippines, Tropical Medicine and Health, 48 (2020). 44. [44]. C. Siettos and L. Russo, Mathematical modeling of infectious disease dynamics, Virulence, 4 (2013), pp. 295–306. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.4161/viru.24041&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F05%2F27%2F2022.05.27.22275675.atom) [1]: /embed/inline-graphic-1.gif [2]: /embed/inline-graphic-2.gif [3]: /embed/graphic-1.gif [4]: /embed/inline-graphic-3.gif [5]: /embed/inline-graphic-4.gif [6]: /embed/graphic-3.gif [7]: /embed/graphic-7.gif [8]: /embed/graphic-8.gif [9]: /embed/inline-graphic-5.gif [10]: /embed/inline-graphic-6.gif [11]: /embed/inline-graphic-7.gif [12]: /embed/graphic-10.gif [13]: /embed/inline-graphic-8.gif [14]: /embed/graphic-11.gif [15]: /embed/inline-graphic-9.gif [16]: /embed/inline-graphic-10.gif [17]: /embed/graphic-12.gif [18]: /embed/inline-graphic-11.gif [19]: /embed/inline-graphic-12.gif [20]: /embed/graphic-13.gif [21]: /embed/graphic-14.gif [22]: /embed/graphic-22.gif [23]: /embed/graphic-24.gif [24]: /embed/inline-graphic-13.gif [25]: /embed/inline-graphic-14.gif [26]: /embed/inline-graphic-15.gif [27]: /embed/graphic-26.gif [28]: /embed/graphic-27.gif [29]: /embed/graphic-28.gif [30]: /embed/graphic-29.gif [31]: /embed/inline-graphic-16.gif [32]: /embed/graphic-30.gif [33]: /embed/inline-graphic-17.gif [34]: /embed/inline-graphic-18.gif [35]: /embed/graphic-31.gif [36]: /embed/inline-graphic-19.gif [37]: /embed/inline-graphic-20.gif [38]: /embed/graphic-32.gif