key: cord-0735771-cgjgengq authors: Nuraini, N.; Sukandar, K.K.; Hadisoemarto, P.; Susanto, H.; Hasan, A.I.; Sumarti, N. title: Mathematical models for assessing vaccination scenarios in several provinces in Indonesia date: 2021-09-24 journal: Infect Dis Model DOI: 10.1016/j.idm.2021.09.002 sha: bbdeb2a5680d19a3306cc5f429fcb2f52c57a157 doc_id: 735771 cord_uid: cgjgengq To mitigate more casualties from the COVID-19 outbreak, this study aimed to assess optimal vaccination scenarios, considering some existing healthcare conditions and some assumptions, by developing SIQRD (Susceptible-Infected-Quarantine-Recovery-Death) models for Jakarta, West Java, and Banten, in Indonesia. The models include an age-structured dynamic transmission model that naturally could give different treatments among age groups of population. The simulation results show that the timing and period's length of the vaccination should be well planned and prioritizing particular age groups will give a significant impact on the total number of casualties. In March 2020, the World Health Organization has declared the COVID-19 outbreak as a global pandemic. As of August 2021, it was reported that at least 206 million people have been infected worldwide, with death toll more than four million [1] . From [2] based on initial data of Indonesia, the infected people are mostly men (56.5%) and in the productive age, 5 31-59 years old (57.5%). Most deaths occurred at aged greater than 60 years (43.6%). The most recurrent clinical symptom was cough (77.8%), the most recurrent co-morbidity was hypertension (52.4%), and the province with the highest COVID-19 incidence was Jakarta (34.3%). Authors in [3] reports that the average COVID-19 incident rate in Jakarta is 99.8 per 10,000 population. Risk factors for the spread of COVID-19 were associated with population's 10 high level of education, which can reflect a higher economic status to the population and a tendency to be more mobile. By implementing well-planned and systematic mitigation strategies, including the closing of country's border and a mass vaccination, some countries have been succeeded in suppress-ing the number of active cases. By contrast, the number of active cases in Indonesia keeps increasing, even though the vaccination program have been implemented since the early of 2021. The death toll in Indonesia is around 115 thousands people [4] , which is the second largest number among Asian countries, after India [5] . Recently, the situation is even more frustrating due to the media news reporting some reinfection of COVID- 19 [6] , even though it has no clear evidence regarding the significant reinfection of COVID-19 in the early of 20 pandemic [7] . Many people presume the occurrence of some new variant of virus, where it also has been reported in some countries, including Brazil [8] and UK [9] . Inevitably, one way to mitigate the COVID-19 outbreaks is to develop effective vaccines and produce them massively for all affected countries. A vaccine could prevent a susceptible person from being infected at least for a time period or even for a lifetime. In early 2021, 25 Indonesia have started the program of first jab with Sinovac, a COVID-19 vaccine developed by a China-based biopharmaceutical company, which has been tested in Indonesia for a Phase 3 clinical trial [10] . This program have brought fresh hope for dealing with COVID-19 mitigation in Indonesia, even though any resistance still occurs in some areas. However, the need of the vaccine existence hastily is like a double-edged sword. Due to 30 the history of other virus-induces illnesses, issues of vaccine mismatch and its suspected side effect on immunocompromised individuals, this existence of vaccine should be well observed. Instead of suppressing the morbidity and mortality of COVID-19, a lack of consideration on vaccination scenarios could also cause several unwanted results, i.e. COVID-19's second outbreak, ineffective vaccination, etc. This will challenge policy-makers and researchers 35 to consider the best scenario of vaccination for suppressing the mortality and morbidity of COVID-19. Recently, several studies have been conducted for discussing models on COVID-19 vaccination. Bitsouni et al. [11] discussed the vaccine effectiveness using the SEIAR (Susceptible-Exposed-Infected-Asymptomatic-Recovered ) model based on cases in Italy. The study focused 40 on the risk of infection spread, the peak prevalence of infection and the time at which the peak prevalence occurred. The paper by Zindoga et al. [12] estimated the effect of social distancing implementation and explored vaccine efficacy scenarios based on cases in South Africa. In [13] , Daihai He, et al. proposed a mathematical model to understand the spread of the virus as the response to the individual behavioral reaction and governmental action, 45 such as travel restriction, and quarantine. In order to describe the behavior of COVID-19 spread on several provinces of Indonesia, i.e. Jakarta, Banten, and West Java, we propose mathematical models, based on non-agestructured and age-structured of SIQRD (Susceptible-Infected-Quarantine-Recovery-Death) models. The level of effectiveness of the vaccination program is put into the model by con-50 sidering the proportion of those who likely to recover or be immune to the illness after getting vaccinated. Simultaneously, there is a chance of reinfection by allowing a portion number of recovered people to be re-infected after a particular time. We expect this re-infection occurrence could be caused by the same or other variant of viruses. By using these constructed models and simulation scenarios, we propose an estimation of optimum vaccination sched-55 ule considering the vaccination cost, healthcare capacity, and vaccination capacity per day. Based on the simulation, the effect of the timing to begin the vaccination on the mortality and the morbidity cases, is also observed. Eventually, we consider several scenarios in prioritizing some age groups and evaluate the results. 2 J o u r n a l P r e -p r o o f The major questions we would like to address in this paper are: (Q1) What is the mathematical model that best represents COVID-19 spread in provinces Jakarta, Banten and West Java, Indonesia? Does the reinfection of recovered people urge the increase of spread significantly? (Q2) What is the effect of vaccination on mortality and morbidity cases of all age groups? When is the optimal timing to begin the vaccination? (Q3) Can we find an optimum vaccination scenario due to the existing government's capacities, including the healthcare facility, estimated maximum capacity of vaccination per day? How many people should be vaccinated for mitigating the spread? (Q4) If we consider the age-structure of susceptible people in those provinces, should we prioritize on some age groups over others? The answers to these questions may help the policymaker, especially in Indonesia, for planning the vaccination program in the near future. In this section, we modifies the standard model SIRD, that has been commonly used in modeling of the vaccination of other virus-induced illness such as influenza [14, 15, 16, 17] . The obtained models, SIQRD, representing the impacted populations due to COVID-19 spread are constructed based on non-age and age structures consequently. The modification includes the addition of the compartment Q representing the dynamics of quarantined people by reckoning the importance of the governmental action [13] . The relation among variables are necessarily constructed based on the observation on the real examined problem. Using the flow diagram in Figure 1 , the vaccination program transfer the suspected people directly into compartment R. The vaccination scenarios will be accommodated by determining a particular function of v(t). Some people in R, the recovered and vaccinated people, possibly can be infected again, so they are transferred back to S with transfer rate ζ. Al-85 though there is no scientifically clear evidence of reinfected cases of COVID-19 according to [7] , yet it is important to follow-up the recovered population for being reinfected to prevent the further spread. There is a possibility that reinfection could be caused by new variants of the virus. On the other hand, the vaccine may not give the perfect protection from infection due to a mismatch problem of the type of virus being used for the vaccine. We can assume 90 the susceptible people are infected only due to their contact with people in compartment I. The the hospitalization or successful isolation of the quarantined ones in compartment Q will not cause the infection. We defined compartment R as the total number of immune people due to both full recovery from the illness and taking vaccination. Using the description of variables and parameters in Table 1 , the constructed mathematical model is given as follows. (1) in which system (1) satisfies We defineN (t) = N −D(t) represents the figure of living people that will generate susceptible newborns, so we have For the sake of simplicity, the number of population is assumed to be constant over the time horizon, or dN dt = 0, so we have π = µ. In this model, there are viral-related and intervention-related parameters and their values are assumed constant. The former ones consist of transmission rate, recovery and death rates 100 due to COVID-19 and reinfection rate. Besides the natural birth and death rates, the re-maining parameters are considered to be the intervention-related parameters. The values of some parameters are given by assumption that follow the existing findings. The natural birth and death rates are obtained from the statistical data of related province being observed. The quarantine rate is based on the time needed for a person to be not-infectious anymore 105 to other persons. The vaccine efficacy and reinfection rate are given, but the values can be simulated in order to see the dynamic of the variables' changes. Other parameter values are estimated/calculated as the results of the implemented numerical method, so that graphs of variables being observed will be the closest the method can gets to the real data provided for each province. Further explanation of these parameters' estimation is given in section 4.1. In the constructed models, the vaccination process impacts on a direct transfer of people from compartment S to compartment R due to the emerging of immune system inside the body of vaccinated people. To develop the vaccination process, we design a periodical schedule, where the vaccination was given several times at certain timings and its value rate will be constant between two timings. There will be k times of vaccination program so there will 115 be k periods from the beginning of scenario. Here simply we assume the jab needed is once per individual, whereas two-jab vaccination can be accommodated in the model further by only adjusting its parameter values . Having determined the above definition, v(t) will be a piecewise function that has constant values during each period of time. The mathematical formulation of vaccination rate 120 v(t) is as follows: where c j is a positive constant, t 0 is the initial time of the pandemic, t v is the time of the first vaccination shot, t V = t v + k∆ is the end of k-th vaccination period, and ∆ (in days) is the length of each period that haves the same rate of inoculation . Firstly, the constructed models is simulated inside the first period before t v when there is no the vaccination program, 125 or v(t) = 0 and the result will become the reference to make the comparison to other periods with the vaccination programs. In this paper, the rate of vaccination c i , i ∈ {1, 2, · · · , k}, during the i-th time interval, is calculated by solving the optimization problem based on some healthcare aspects existing in Indonesia. The detailed description of the optimization problem is explained further in section 5.1. whereN Vaccination rate v(t) (2) is also implemented into this age-structured model where the vaccination rate v i (t) refers to the vaccination rate of the age group-i. For example, v 2 (t) is defined the vaccination rate of the 2 nd age group, where values c ij are positive constants, for i = 1, 2, . . . , 5 and j = 1, 2, . . . , k. In the first and second equations of system (3), the transmission rate β ij = β i C ij gives 140 cross transmission between age groupsi and j, where β i represents the probability of infection level in age group i, and C ij is a matrix describing contacts between any pair of age group i and j, where i, j = 1, 2, . . . , 5. Here the unit of C ij and β ij is 1/day. Function v i (t) gives the potent vaccination rate of group-i, which shows the success of vaccination program based on its efficacy for age group-i. We define the values of C ij based on the results from the study by Data of COVID-19 victims of provinces Jakarta, Banten, and West Java is retrieved from https://kawalcovid19.id/ [4] . It consists of time-series data of active cases, total recovered 150 cases, and total deaths from late March until late October 2020, which is called as Dataset 1 following which will be used to extract the parameters of the systems (1) and (3). The chosen time period of the data is expected to represent the real figure of COVID-19 spread in each province before the vaccine inoculation. On the other hand, it is an unfortunate that the data of COVID-19 victims along the observed time period does not include details on their age 155 for those provinces. We define values of parameter β i for system 3 from comparisons of the population pyramid of a state in USA that has the desired data and has the similar portions of age groups with of the provinces being observed. The best choice is Connecticut, USA, where its data resembles enough the province's data. The COVID-19 data of Connecticut, USA, is retrieved from https://data.ct.gov/ [19] . This website provides time-series data on 160 total infections by age group, which is called as Dataset 2. The ratio between the number of COVID-19 victims and the total number of population per age group in each province is assumed to have same percentage for each age group in Connecticut. Data on population 6 J o u r n a l P r e -p r o o f by age group in provinces Jakarta, Banten, West Java, and Connecticut are retrieved from [20, 21] and [22] . 165 Firstly, we introduce two additional variables; CR(t), representing the total number of recovery due to COVID-19, and V (t), depicting the total figure of vaccine inoculation at time-t. The rate of first variable is proportional to the existing number of people in Q, or In this stage, the value of transmission rate, recovery rate, and death rate due to COVID-19 are estimated numerically using the Least Square Method (LSM), so that the values of Q(t), CR(t) and D(t) generated from system (1) are close enough to the respective real data of active cases (DAT A Q ), total recovered (DAT A CR ), and total death (DAT A D ) due to COVID-19. The initial value of compartment I(t) is unknown since it also includes the not-quarantined ones, so we estimate its value. Simultaneously, we define where their values are the solution of the following minimization problem. and D is a search domain. Remember that DAT A Q , DAT A QR , DAT A D are the real retrieved data. Basically, the best estimated value of T makes the objective func-170 tion (4) is near zero so it means the model dynamics is close to the real data. The second variable is defined following this differential equation. Note that the number of vaccinated people V (t) at time t is not always fully added up to the number of immune people R(t) because of the effect of vaccine efficacy. The values V (t) are computed as the results of solving its differential equation (5) together with SIQRD model, using the built-in function in MATLAB, i.e. ode45. In the fitting process, the natural recruitment and the natural death rates are π = µ = 1 365·70 where the life expectancy in Indonesian people is 70 years, based on [23]. The other assumed values of parameters are q = 0.4 (the quarantine rate), ζ = 1 150 (reinfection rate) according to [24] , and v(t) = 0 for t in a period when the vaccine is not available. Later, we set η = 80% (the vaccine efficacy) when the vaccine is available. Table 2 for systems (1) and (3), the estimated values of parameters and initial condition of I(t) in Jakarta are given. Based on [25] , there were only 2 cases reported in the first day of pandemic in Jakarta. Using the Jakarta's basic reproduction number that is slightly higher than 1, approximately there were more than 60 persons had been infected that day. The simulation of the model utilizing the estimated parameters and initial condition 185 is depicted by Figure 2 plotted together with the existing data. The dynamics of total recovered and deceased due to COVID-19 are well-fitted to the data, although the figure for total deceased is extremely smaller than that of recovered. However, particularly in Jakarta, the active case dynamic is close-fitted to the data only in the late simulation since it underestimates the fluctuating data at other times. Given the assumption of constant 190 parameters over time, the model is not be able to capture the data precisely, especially when it fluctuates. However, the overall fitting results are well-fitted in most of the observed time in Jakarta. The other two fitting results for Banten and West Java are given in Appendix A and both results show a well-fitted model most of the time. To construct the age structured model in system (3) for each province, we apply Con-195 necticut data by using Algorithm 1. Note that we use the same values of β, γ, δ and I(t 0 ) from the non-age-structured model for each province. Steps: • Estimate the recovery and death rate γ dan δ for Connecticut, USA using the non-age-structured data applied to the non-age-structured model, • Using Dataset 2 applied to system 3, estimate the probability of success transmission, denoted as x i , . Notice that this quantity represents the proportion of the probability of success transmission in Connecticut, AS. Assuming that the proportion of success transmission follows the Connecticut, assume β i = α ·x i as the probability of success transmission in the observed provinces, • Use Dataset 1 and estimated parameters, γ, δ, for the observed provinces to estimate the best value of α. The estimation is obtained by minimizing the square error of data, non-age-structured active cases data, and model, The obtained estimation values of β i of Jakarta, for i ∈ {1, 2, 3, 4, 5} are shown in Table 200 3. The largest values are owned by age-groups 3,4, and 5 using Algorithm 1. The numbers of active cases in Jakarta are dominated by three older age-groups as seen in Figure 3 . The results for other provinces are available in the Appendix. β 1 β 2 β 3 β 4 β 5 0.0251 0.0355 0.1122 0.1698 0.4488 Table 3 : Obtained values β i for the age-structured model of Jakarta Having developed the SIQRD model of both non-age-structured and age-structured, first 205 we analyze the dynamics of all obtained models for v(t) = 0 or without vaccination program, therefore we can project the peak of outbreak that may happen in the near future. Figure 4 (a) portrays the dynamics of Q(t), R(t) and D(t) for non age-structured model of Jakarta. It shows that the peak of active cases will be in February 2021. On the other hand, the number of immune people will largely increase and then decrease as some of them become 210 susceptible again. This condition is possibly causing the second outbreak of the disease later. Notice that for a certain region, the result given by the non-age-structured model is similar to the result given by the age-structured model, e.g. peak occurrence time. However, they are not precisely similar in numbers since we use different objective functions on estimating its parameters. For instance, both Figure 4 (a) and (b) are showing that the peak of outbreak will occur on February 2021 in Jakarta. However, they are numerically different since the non-age-structured model estimated that the number of active cases reaches nearly 80,000 cases but the age-structured model estimated 65,000. (a) (b) Figure 4 : SIQRD simulations without vaccination in Jakarta: (a) non-age-structured; (b) age-structured model In Figure 4 (a), the graph of total number of deaths D(t) increases because the existence of infectious people. Based on the age-structures model, the dynamics of active cases of Jakarta 220 are classified into 5 age groups in Figure 4 (b). Similar to Figure 4(a) , the age-structured model also illustrates the second outbreak due to the reinfection for each age-groups. By using the constructed models, this simulation answers (Q1), we can see that the reinfection played a vital role in producing the second outbreak. The peak of outbreak projections of SIQRD model in Banten and West Java are given by 225 Figure 5 . Different from the projection of Jakarta, the dynamics of active cases will increase until August 2021. This is possible because the populations of Banten and West Java are larger than the population of Jakarta. Nevertheless, the general behavior of models is similar where the reinfection factor could cause the second peak of outbreak. Now we discuss whether the existence of vaccination program will generally reduce the number of active cases and total deaths or not. Based on equation (2), we set k = 12 and the length of one period was 30 days, so the whole vaccination period is about one year. Set c j be any values randomly chosen, with 0 ≤ c j ≤ 10 −3 , which represents the vaccination rate in the jth-period, j = 1, 2, . . . , k. The vaccination program is set to begin in the third week 235 of October 2020. It is shown Figure 6 , the numbers of active cases and total deaths after 1-2 months of vaccination in Jakarta become less than half of the related numbers in the model without vaccination. Later in section 5, we will determine the optimal values of v i (t) that met some desired constraints. There will be several scenarios of vaccination being assessed on particular age 240 groups so we could analyse the urgency of vaccination program. We want to observe the dynamics on the numbers of active cases Q(t) using simulations with different values of vaccine efficacy η and quarantine rates q. The first simulation is done by varying the vaccine efficacy and keeping the quarantine rate constant, where the results are in Figure 7 (a). The second turn of simulation is using the opposite scheme, where the results are in Figure 7(b) . The values for v(t) are chosen as the same as the ones already used in the previous subsection. the number of active cases. Remember thatv(t) = ηv(t). When the vaccine efficacy is low, we require high values of v(t) to suppress the numbers of mortality and morbidity of COVID-19. In Figure 7 (b), shows that the higher the value of the quarantine rate, the lower the number of active cases is. As a conclusion, we need vaccine with as high efficacy as possible, and also need a quite high quarantine rate to suppress the number of active cases. The total number 255 of deaths will follow the same. Previously, the simulation of the models is carried out without any constraints in finding the solutions of SIQRD systems (1) and (3). Now we examine several scenarios for vaccination with the main objective of effectively reducing the numbers of active cases and death toll 260 at the minimum cost. It is obvious that the higher vaccination rate needs large number of vaccine provided by the government, and so it requires very high expenses. We establish an approximation problem dealing with vaccination function v(t) with some constraints based on the real situation. There are conditions that should be considered in this optimization process in order to have feasible solution for the real problem, which is as follows: • The maximum number of active cases does not exceed the maximum capacity of the existing healthcare facility, denoted by K 1 (persons), • The number of daily vaccine does not exceed the maximum shots of vaccination per day provided by Indonesian government, denoted by K 2 (per day). The model being observed first is the non age-structured model. Some of results of this 270 optimization problem using this model will be used to execute some scenarios using the age-structured model later. Since comorbidities contribute a high number of deaths, we assume that comorbidities exist in all age groups. Thus, applying vaccines to all age groups simultaneously can be one of the scenarios worth considering. The scenarios are implemented in finding answer to questions (Q1)-(Q3). In accordance with the required constraints, we define the optimization problem as follows: Minimize where t ∈ [t v , t V ] as defined in Section 2.1 and ω represents the expense needed to vaccine a single person. Besides that, the objective function f has to meet the following constraints: As ∆ represents the time period length of each vaccination, it is defined as 30 days which means the rate of vaccination will be kept to level off in certain values for every 30 days in a row, roughly a month. This assumption is based on the initial plan of government on providing the vaccine in the term of months. According to [26] considering the readiness of the vaccine, the first vaccine inoculations are set in several phases having different length 285 periods, e.g. 4 months (Jan-Apr 2021) for the first phase and 3 months later (May-Jul 2021) for the second phase of vaccination. Hence, choosing ∆ = 30 days is reasonable and it allows the rates of vaccine to vary even though they are in the same phase of vaccination. For the maximum number of active cases K 1 and the maximum number of vaccinated persons K 2 , we defined two possible conditions; limited and good facilities, and analysed whether there are 290 solutions for this approximation or not. Practically, the rate of vaccination v(t) is obtained by means of solving the defined optimization problem. The original constrained problem is first converted into the unconstrained one by adding the penalty value once the point violates the constraints following which it is solved numerically using the trust-region method. According to [27] , Indonesia government is able to provide about 31.000 COVID-19 testing 295 devices per day nationwide. However these testings kits are not evenly distributed to all provinces in Indonesia. Respectively, Jakarta, Banten and West Java get approximately 31%, 3.4%, and 10% of the total number of provided testing kits. We assume the limited facility condition based on this news, where the maximum numbers of vaccines received by Jakarta, Banten and West Java are 10.000 , 1.100, and 3.100 vaccines per day, respectively. On the 300 other hand, the government claimed there will be 1 million vaccines per day nationwide [28] . We assume the good facility condition based on this news, so the maximum availability of vaccines are 310.000 vaccines for Jakarta, 100.000 vaccines for West Java, and 34.000 vaccines for Banten. The simulation using the limited and good conditions of facility is examined when the 305 vaccination starts in October 2020. Even after solving the optimum process problem with limited facility condition, solution has not been found yet. So that the good facility condition is used for simulations from now onward. In showing the results, the output graph D(t) and Q(t) are always equipped with the one coming out from the model without vaccination so we can see their stark differences. We set the vaccination program starting in different timeline; October 2020, January 2021 and after the peak of outbreak for each province by setting different t v for each simulation. The three mentioned scenarios are aimed to compare the expected situation if the first vaccine jab is delivered at three different times. The value K 1 is adjusted for each province to be 315 higher than the existing number of active cases at the time the vaccination starts. The values of K 1 and K 2 are given in Table 4 . Times of vaccination Jakarta Banten West Java 30,000 310,000 20,000 34,000 110,000 100,000 January 2021 65,000 310,000 20,000 34,000 110,000 100,000 After the peak of outbreak 78,500 310,000 100.000 34,000 780,000 100,000 There are many possibilities of obtained values v(t) coming as the solutions of the optimization problem. Figure 8 ing for only after three periods of vaccination. In Table 5 , we can see that none of the constraints is violated. The total number of vaccinated people is about 5.14 million, which is only 49.6% of Jakarta's population. This scenario looks promising where the remaining active cases number is only a single person at the final time t V , or a year later. Unfortunately, this scenario is impossible to apply since in fact the vaccine is not ready yet. Some news 330 media predicted the readiness of the vaccine is not earlier than January 2021. If the starting time is in January 2020, Figure 10 (a)-(b) and Table 6 give the results for Jakarta. It shows the number of susceptible people rapidly decreased for only after second period of vaccination. The fluctuated number of vaccination per day scheduled for each 340 period, plotted in blue thick line, is one of solutions found from the approximation process. It is interesting to find out whether the rate of vaccination in constant value will be effective or not, because it will be simpler to operate in the real condition. We have this kind of scenario later in this paper. in Jakarta. The vaccination is successful because the remaining number of active cases is only one. The simulation of the vaccination in Banten and West Java are given in Figure 11 , where the number of active cases is limited to the value of K 1 . The number of vaccinated people is 2.1 million or 17% of the population in Banten, and 10.1 million or 22.7% of the population in West Java. Now the starting time of vaccination is after the number of active cases reaches its peak. So the starting time is different for each province, April 2021 in Jakarta, 1.25 million in Banten, and 8.5 million in West Java. Figure 12 (a) shows the steeper decrease of Q(t) than of the decrease without vaccination in Jakarta. The number of death D(t) seems to stabilize to about 4,000 due to the peak of outbreak happened before the vaccination. From Table 7 , 355 total number of people to be vaccinated is 3.19 million or 30.4% of population in Jakarta. If we compare among figures 8(a), 10(a) and 12(a) in sequence, total numbers of death are increasing, which means the later the starting time of vaccination, the higher the number of the pandemic casualties. hence lift the number of active cases down. The calculated rates of inoculation are varied each month which is based on the definition of v(t) in section 2.1 which is too complicated if it comes to the practical implementation. The extremely fluctuating rates of jab can cause extreme addition and subtraction of medical personnel each month. In a response to that, the constant rate of vaccination is further discussed in section 5.1.3. In Figure 10 (a), the vaccination for 12 months could make the number of active cases starting to decrease significantly in the first month of vaccination period. It is interesting to see whether one time vaccination in the first month could really work. Most of the time when the first attempt showed good result, vaccination program could be immediately stopped to 370 reduce cost. Figure 13 depicts a simulation of this one-time vaccination for Jakarta, where v(t) = 0.0031 for time t in January 2020, and v(t) = 0 for other months. It shows a significant decrease of the active cases at the beginning of vaccination, but it will start to increase from 375 August 2021 onward, so this may cause another outbreak in the future. This prediction is clearly proven when we plot the graphs in longer time-span in Figure 14 . On the left, the vaccination consistently given for 12 months from January 2021 can make the non-existence of active cases happen until the beginning of year 2024. On the right, the one-time vaccination potentially impacts on another outbreak in July 2022. The increase of death toll seems to be slowed down about one year, then it will increase to the same number as the obtained number from the non-vaccine model with few-month delay (a) (b) Figure 14 : (a) 12-month vaccination; (b) One-month vaccination in Jakarta The scenarios from the previous subsection assess the timing and frequency of vaccination. The scheme of vaccination obtained from the approximation process gives fluctuating values of vaccination rates making it hard to be implemented . It means there will be fluctuating needs of the vaccination workers per period. Now we simulate the vaccination with constant rates and analyze the dynamics of the obtained numbers of active cases. There are 3 scenarios with different values of constant rates starting from January 2021 390 based on the optimal vaccination rates obtained on Table 6 . The constant vaccination rate for the first, second and third scenarios are the average of rates (c = 0.0014), the maximum rate (c max = 0.0031), and the minimum rate (c min = 0.0005), respectively. The number of vaccinated people of those scenarios are plotted in Figure 15 , where the original optimum values, the average, the maximum, the average and the minimum rates are respectively rep-395 resented by blue stripes, upper magenta stripes, middle magenta stripes, and lower magenta lines. Due to the decreasing number of susceptible people, the number of vaccinated ones is also decreasing, even though the vaccination rates are constant. Let the optimum fluctuated rates be the benchmark of the uses the healthcare facility that requires an expense calculated from equation (6). The expenses calculated from implementing 400 those scenarios will be compared to this benchmark expense. Figure 16 shows the active cases Q(t) from all scenarios where the optimum vaccination rate is plotted in solid blue line. The average rate scenario gives values Q(t) exceeding the maximum healthcare capacity (the green line) and it requires 99.1% of the benchmark cost. The maximum rate scenario gives plot of Q(t) that resembles the plot from the optimum rates, but it requires almost twice of the 405 benchmark cost, i.e. 192%. Finally, the minimum rate scenario gives insignificant reduction of Q(t) compared to the plot of Q(t) without vaccination. It concludes the optimum fluctuated rates are the best scheme among other constant rate scenarios based on the dynamic of active cases and the required cost. Considering the results of the simulation of starting time scenarios, the vaccination program ideally should starts in the middle of the pandemic, which is October 2020. The total number of deaths can be suppressed significantly. Unfortunately, this scenario is impossible to apply since the vaccine was not ready yet at that time. If the beginning of vaccination program is delayed as in the scenarios, the required amount of vaccine will decrease, so it 415 seems better to begin vaccination program after the peak of outbreak with the least expense. However, delaying starting time will result in increasing total deaths. In fact, the decrease of total number of deaths is not significant if we choose the latest scenario. Having simulated on the consistency on the frequency of vaccination, it shows in Figure 13 (a) that only a single time of vaccination gives insignificant reduction of the number of 420 active cases at the end of the intended vaccination program. Moreover, the final number of infectious people in the simulation, I(t V ), is about 800 people in Jakarta, which is still too high so it could trigger another peak of outbreak through reinfection. It is concluded that a consistent schedule of vaccination will significantly reduce the number of active cases at the end of vaccination program. The pattern of the vaccination rate is also interesting to observe. Having seen the results from the first two scenarios, where the vaccine was applied before the number of active cases reaches its peak, the vaccination rate obtained from the optimization procedure is high at the beginning of vaccination schedule. It seems the healthcare facility is not yet at maximum capacity K 1 , so much effort can be used to reduce the number of active cases by maximizing 430 the vaccination shots per day. On the other hand, if the vaccine is applied after the peak of outbreak, the vaccination rate will be low at the beginning and then high at the end of the vaccination. In this scenario, the healthcare facility is already at the maximum capacity, so the main focus is preventing another peak of outbreak to come. Now four scenarios of prioritizing particular age groups are implemented on the agestructured model, where the priority is respectively given to groups of the active and older people, groups of 20 years old and younger, group of active people only, or alternating target's group in each period. The starting time is January 2020 and we use the scheme of vaccination rates shown in Figure 10 . During the first six months, the vaccine shots are given to workers and high-risk people, which are 20 and older, then they are given to younger groups in the remaining months. As seen in Figure 17 (a), the darker color of the cells for certain periods means higher rate of vaccination. Figure 17(b) illustrates the simulation result of the total numbers Q(t) and 445 D(t) of all age-groups using this scenario in Jakarta. The simulation results from Banten and West Java are given in Figure 18 . The vaccination performs well in the reduction of the numbers of active cases and death toll. In this scenario, the first interval of six months is the vaccination time only for the younger 450 people with ages 0-19, and the remaining time is for the active and older people's groups, as scheduled in Figure 19 (a). The dynamics of total numbers Q(t) and D(t) of all age-groups using this scenario in Jakarta, Banten and West Java are described in Figure 19 (b) and Figure 20 , respectively. It seems this scenario gives insignificant decrease on the total numbers of active cases and death toll. In this scenario, we change the target of age-groups in certain ways that is shown in Figure 21 (a). High rate vaccination in the first month is given to active and older people, ages 20 and older. Figure 21 (b) and Figure 22 describe the dynamics of total number Q(t) and D(t) of all age-groups obtained using this scenario in Jakarta, Banten, and West Java. Compared to those obtained by prioritizing active and older people, this scenario yields the dynamics of active cases having thicker tail. This argument can be seen in every region we observed. Thus, this scenario is not likely preferable than the first scenario prioritizing the active and older people. If the Indonesian government plans to conduct the vaccination only to the group of active people only, ages 20-49, we do the simulation in this scenario. Figure 23 depicts the dynamics of Q i (t) once the third age-group is vaccinated. On the other hand, we also provide the dynamics of active cases that results from the other three scenarios as a comparison. Figure 24 shows that the vaccination prioritizing the active and older people is more preferable due 470 to its significant reduction of cases. This argument is also valid to the results given in Banten and West Java as depicted by Figure 25 . The simulations show that the vaccination scenario by first targeting active and older adults (above 50) (group 3-5) is better than other scenarios. The numbers of active cases 475 and also total deaths are decreased significantly. In the result of scenario 2, the decrease is insignificant, because the transmission rate of virus among 20-year-old and younger people is much lower than of other age groups, as shown in Table 2 . After conducting the simulation by targeting the active people only, the dynamics from other age-groups on the number of active cases, Q i (t), i = 3, are also decreasing once people 480 in the third age-group was vaccinated. It is logic because the key of the transmission is contact among people. The active people have larger access to people of other age-groups so the transmission across age groups is highly possible. Once the biggest source of infectious people is vaccinated, it affects the number of active cases on the other group of age. Having constructed the SIQRD for non-age and age structures models and developed several scenarios on the implementation of vaccination program, we concluded our work by the following findings related to the proposed questions: • Modifying the SIRD into the SIQRD model by adding quarantine, reinfection, and even vaccination aspects are considered capable to represent how COVID-19 spreads 490 in several provinces in Indonesia. Utilizing the existing data and information related to those provinces, the graph of the simulated model is well to resemble the figure that of happening in the real situation. • The proper timing of vaccine implementation is in the early stage of the pandemic. This implementation suppresses the number of active cases immediately, and consequently the total deaths. After the number of active cases reaches its peak, the implementation of the vaccination program is not reduce the total deaths significantly. • The vaccination should be implemented consistently following the planned schedule for a certain period. The vaccination's implementation for only one or two months does not reduce the number of infectious persons, and eventually it fails to prevent the second 500 and more peaks of outbreak. • The prioritization of the active and older adults (above 50) for vaccination over others and prioritization of active people only reduce significantly the total deaths. Appendix A. Graphs of non age-structured model Those figures above are depicted by the estimation results of β, γ, δ, and I(t 0 ) for Banten and West Java. The obtained parameters of the non age-structured model are given by the following table. Region β γ δ I(t 0 ) R 0 Banten 0.5968 0.0520 0.0024 9 1.0447 West Java 0.6016 0.0354 0.0012 20 1.0528 Sociodemographic and environmental health risk factor of COVID-19 in Jakarta, Indonesia: An ecological study Kawal informasi seputar COVID-19 (in Bahasa) Indonesia's coronavirus fatalities are the highest in Southeast Asia Can you get infected again after recovering from coronavirus? Reinfection or reactivation of severe acute respiratory syndrome coronavirus 2: a systematic review Collapse of the public health system and the emergence of 545 new variants during the second wave of the COVID-19 pandemic in Brazil Sinovac launches Phase 3 trial for COVID-19 vaccine in Indonesia, reports Phase 2 details A model for the outbreak of COVID-19: Vaccine effectiveness in a case study of Italy Quantifying early COVID-19 outbreak transmission in South Africa and exploring vaccine efficacy scenarios A conceptual model for the coronavirus disease 2019 (COVID-19) outbreak in Wuhan, China with individual reaction and governmental action Optimal Vaccination and Treatment Schedules in a Deterministic Modeling the impact of mass influenza vaccination and public health interventions on COVID-19 epidemics with limited detection capability Optimal Vaccination Schedules in a Deterministic Epidemic Model Modeling the Effects of Vaccination and Treatment on Pandemic Influenza Social Contacts and mixing patterns relevant to the spread of infectious disease COVID-19 Cases and Deaths by Age Group COVID-19-Cases-and-Deaths-by-Age-Group Jumlah Penduduk Menurut Kelompok Umur di Provinsi (in Bahasa) Pusat Data Ekonomi dan Bisnis Indonesia Connecticut's Official State Website Indonesia: Life Expectancy Humoral Immune Response to SARS-CoV-2 in Iceland, The NEW 590 ENGLAND JOURNAL of MEDICINE Kapan Sebenarnya Corona Pertama Kali Masuk RI? Covid-19 Phase-2 Vaccine Raw Materials Have Arrived in Indonesia Pemerintah RI Targetkan Suntik Vaksin Corona ke 1 Juta Orang Per Hari (in Bahasa) Indonesia masih kejar target 30.000 tes PCR per hari (in Bahasa) We provide tables representing the estimation results of vaccination scheme for Banten and West Java with several starting times of vaccination.