Mathematical models for assessing vaccination scenarios in several provinces in Indonesia

To mitigate casualties from the COVID-19 outbreak, this study aims at assessing the optimal vaccination scenarios, considering several existing healthcare conditions and 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 allows for different treatments among different age groups of the population. The simulation results show that the timing and period of the vaccination should be well planned and prioritizing particular age groups will give a significant impact on the total number of casualties.


Introduction
In March 2020, the World Health Organization 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 (COVID-19 coronavirus pandemic, 2020). From (Herawati and Setiyabudi, 2021) based on initial data of Indonesia, the infected people are mostly men (56.5%) and in the productive age, 31e59 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%). The study in (Zakianis et al., 2021) reported 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 the population's 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 border closing and a mass vaccination, several countries have succeeded in suppressing the number of active cases. On the contrary, in Indonesia the number keeps increasing, even though the vaccination program has been implemented since early 2021. The death toll in Indonesia is around 115 thousands people (Kawal informasi seputar, 2020), which is the second largest number among Asian countries, after India (Indonesia's coronavirus, 2020). Recently, the situation is even more frustrating due to the media news reporting some reinfection of COVID-19 (covid-19-coronavirus-vaccine, 2020), even though it has no clear evidence (Tang et al.,). Many presumed that it is due to the occurrence of new variants of the virus, that have been reported in some countries, including Brazil (Da Silva & Pena, 2021) and UK (Pan et al., 2021).
Inevitably, one way to mitigate the COVID-19 outbreaks is to develop effective vaccines and produce them massively for all affected countries. Vaccine could prevent a susceptible person from being infected at least for a time period or even for a lifetime. In early 2021, Indonesia has started a program of first jab with Sinovac, a COVID-19 vaccine developed by a Chinabased biopharmaceutical company, which has been tested in Indonesia for a Phase 3 clinical trial (Sinovac launches Phase 3, 2020). This program has brought a fresh hope for COVID-19 mitigation in Indonesia, even though resistance still occurs in some areas.
However, the quick need for vaccination is like a double-edged sword. Due to the history of other virus-induced illnesses, as well as issues of vaccine mismatch and suspected side effects on immunocompromised individuals, vaccine administration should be well observed. Instead of suppressing the morbidity and mortality of COVID-19, a lack of consideration on vaccination scenarios could also cause unwanted results, such as COVID-19 s outbreak, ineffective vaccination, etc. These pose challenges to the policy-makers and researchers to consider the best scenario of vaccination.
Recently, several studies have been conducted that model COVID-19 vaccination. (Bitsouni et al., 2020) discussed the influence of vaccine effectiveness using a SEIAR (Susceptible-Exposed-Infected-Asymptomatic-Recovered) model based on cases in Italy. The study focused 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. (Mukandavire et al., 2020) estimated the effect of social distancing implementation and explored vaccine efficacy scenarios based on cases in South Africa. (Lin et al., 2020), proposed a mathematical model to understand the spread of the virus as the response to the individual behavioral reaction and governmental action, such as travel restriction and quarantine.
In order to describe the behavior of COVID-19 spread in several provinces of Indonesia, i.e. Jakarta, Banten, and West Java, we propose mathematical models, based on non-age-structured and age-structured SIQRD (Susceptible-Infected-Quarantine-Recovery-Death) models. The level of effectiveness of the vaccination program is put into the model by considering 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 of the recovered people to be re-infected after a particular time. We assume that this re-infection can be caused by the same or different variant of viruses. By using these constructed models and simulation scenarios, we propose an estimation of optimum vaccination schedule 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 different age groups and evaluate the results.
The major questions we would like to address in this paper are: (Q1) Which mathematical model best represents the COVID-19 spread in Jakarta, Banten and West Java of Indonesia? Does the reinfection of recovered people cause the increase of spreading 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 optimal vaccination scenario considering the existing government's capacities, including the healthcare facility, and 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 some age groups over the others?
The answers to these questions may help the policymaker, especially in Indonesia, for planning a vaccination program in the near future.

Proposed model and analysis
In this section, we modify the standard SIRD model, that has been commonly used in modeling vaccination of other virusinduced illnesses, such as influenza (Feng et al., 2011;Hethcote & Waltman, 1973;Qian et al., 2020;Tasmi, 2016). The obtained models, SIQRD, representing the impacted populations due to COVID-19 spread are constructed with and without age structures. The modification includes the addition of the compartment Q representing the dynamics of quarantined people by reckoning the importance of governmental action (Lin et al., 2020). The relation among variables are necessarily constructed based on the observation of the real examined problem.

Non-age-structured model
Using the flow diagram in Fig. 1, the vaccination program transfers 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 z. Although there is no scientifically clear evidence of reinfected cases of COVID-19 according to (Tang et al.,), yet it is important to follow-up the recovered population for being reinfected to prevent a further spread. There is a possibility that reinfection could be caused by new variants of viruses. 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 that the susceptible people are infected only due to their contact with people in compartment I. The hospitalization or successful isolation of the quarantined ones in compartment Q will not cause infection. We define 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.
dSðtÞ dt (1) NðtÞ ¼ N À DðtÞ that represents the number of living people that will generate susceptible newborns, we have dN dt NðtÞ: 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 p ¼ m.
In this model, there are viral-related and intervention-related parameters and their values are assumed to be constant. The former ones consist of transmission rate, recovery and death rates due to COVID-19 and reinfection rate. Besides the natural birth and death rates, the remaining parameters are considered to be the intervention-related parameters. The values of some of the parameters are given by assumption that follow the existing findings. The natural birth and death rates are obtained from the statistical data of the provinces. The quarantine rate is based on the time needed for a person to be not-infectious anymore 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 that the method can get 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 the 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 be k periods from the  beginning of the scenario. Here simply we assume that 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, y(t) will be a piecewise function that has constant values during each period of time. The mathematical formulation of vaccination rate 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 y þ kD is the end of k-th vaccination period, and D (in days) is the length of each period that has the same rate of inoculation.
Firstly, the constructed model is simulated within the first interval when there is no vaccination, or v(t) ¼ 0, and the result will become the reference to make the comparison to the other periods of the vaccination program. In this paper, the rate of vaccination c i , i 2 {1, 2, …, k}, during the i-th time interval, is calculated by solving an optimization problem based on some healthcare aspects existing in Indonesia. The detailed description of the optimization problem is explained further in Section 5.1.

Age-structured model
The age-structured model of SIQRD is similar to the previous system of Eq. (1). For each compartment, we define five age groups, i.e. 0-9, 10e19, 20e49, 50e59, and 60 or higher, so all age groups have a compartment assigned to them. For age group i 2 {1, 2, 3, 4, 5}, the obtained model is given as.
Vaccination rate v(t) given by Eq.
(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, y 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 b ij ¼ b i C ij gives cross transmission between age groups i and j, where b i represents the probability of infection level in age group i, and C ij is a matrix describing contacts Natural death rate 1/day assumed between any pair of age group i and j, where i, j ¼ 1, 2, …, 5. Here the unit of C 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 the age group i. We define the values of C ij based on the results from the study in (Mossong et al., 2008), while the values of b i are estimated based on the real data provided in Indonesia.

Datasets
The data of COVID-19 in Jakarta, Banten, and West Java are retrieved from https://kawalcovid19.id/(Kawal informasi seputar, 2020). It consists of time-series data of active cases, total recovered cases, and total deaths from late March until late October 2020, which we call 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 unfortunate that the data of COVID-19 victims along the observed time period do not include details of their age. We define values of the parameter b i for system of Eq. (3) from comparing the population pyramid of a state in the USA that has the desired data and similar portions of age groups with that of the provinces being observed. The best choice is Connecticut, USA, where its data resembles enough the provinces data. The COVID-19 data of Connecticut, USA, is retrieved from https://data.ct.gov/(COVID-19 cases and deaths by age group, 2020). This website provides time-series data on total infections by age group, which we will call 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 be the same with that in Connecticut. Data on population by age group in Jakarta, Banten, West Java, and Connecticut are retrieved from (Jumlah Penduduk Menurut, 2020; Pusat Data Ekonomi, 2020) and (Connecticut's official state website, 2020).

Parameters estimation
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 the first variable is proportional to the existing number of people in Q, or dCRðtÞ dt ¼ gQðtÞ: In this stage, the transmission, recovery, and death rates 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 (DATA Q ), total recovered (DATA CR ), and total death (DATA D ) due to COVID-19. The initial value of compartment I(t) is unknown since it also includes the non-quarantined ones. We therefore estimate its value. Simultaneously, we define where their values are the solution of the following minimization problem.
Remember that DATA Q , DATA QR , DATA D are the actual retrieved data. Basically, the best estimated value of T makes the objective function given by Eq. (4) close to zero, which means that the model dynamics is close to the real data. The second variable is defined following this differential equation dVðtÞ dt ¼ yðtÞSðtÞ: 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 Eq. (5) together with the SIQRD model, using the built-in function ode45 in MATLAB.
In the fitting process, the natural recruitment and the natural death rates are p ¼ m ¼ 1 365,70 where the life expectancy in Indonesia is 70 years, based on (Indonesia: Life Expectancy, 2020). The other assumed values of parameters are q ¼ 0.4 (the quarantine rate), z ¼ 1 150 (reinfection rate) according to (Gudbjartsson, 2020), and v(t) ¼ 0 for t in a period when the vaccine is not available. Later, we set h ¼ 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 (Kapan Sebenarnya Corona, 2020), 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 infected that day. The simulation of the model utilizing the estimated parameters and initial condition is depicted in Fig. 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 parameters over time, the model is not able to capture the data precisely, especially when they fluctuate. However, the data 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, where the model also fit the data most of the time.

As in
To construct the age structured model in system of Eq. (3) for each province, we apply the Connecticut data by using Algorithm 4.1. Note that we use the same values of b, g, d and I(t 0 ) from the non-age-structured model for each province.

Algorithm 1. Estimating the values of b i
The obtained estimation of b i of Jakarta, for i 2 {1, 2, 3, 4, 5} are shown in Table 3. The largest values are owned by agegroups 3,4, and 5 using Algorithm 4.1. The numbers of active cases in Jakarta are dominated by three older age-groups as seen in Fig. 3. The results for other provinces are available in the Appendix.

Dynamics of models with and without vaccination
Having developed the SIQRD model with and without age classes, first we analyze the dynamics of the models for y(t) ¼ 0 or without vaccination to project the peak of the outbreak that may happen in the near future. Fig. 4(a) portrays the dynamics of Q(t), R(t) and D(t) for the 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 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 at given by the age-structured one, e.g. in the peak occurrence time. However, they are not precisely the same in numbers since we use different objective functions on estimating its parameters. For instance, both Fig. 4(a) and (b) are showing that the peak of outbreak will occur in 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 while the age-structured one estimated 65,000.
In Fig. 4(a), the graph of the total number of deaths D(t) increases because of the existence of infectious people. Based on the age-structured model, the dynamics of active cases in Jakarta are classified into 5 age groups in Fig. 4(b). Similar to  Table 3 Obtained values b i for the age-structured model of Jakarta.   Fig. 4(a), the age-structured model also illustrates the second outbreak due to the reinfection for each age-groups. By using the constructed model, this simulation shows that 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 in Fig. 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 that of Jakarta. Nevertheless, the general behavior of the model is similar where the reinfection factor could cause the second peak of outbreak. Now we discuss whether or not the presence of vaccination will generally reduce the number of active cases and total deaths. Based on Eq.
(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 to 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 of October 2020. It is shown in Fig. 6, that the numbers of active cases and total deaths after 1e2 months of vaccination in Jakarta become less than half of the 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 in particular age groups so we could analyze the urgency of a vaccination program.

Sensitivity analysis on vaccine efficacy and quarantine rate
We want to observe the dynamics of the active cases Q(t) using simulations with different values of vaccine efficacy h and quarantine rate q. The first simulation is done by varying the vaccine efficacy and keeping the quarantine rate constant, where the results are in Fig. 7(a). The second simulation is using the opposite scheme, where the results are in Fig. 7(b). The values for y(t) are chosen to be the same as the ones already used in the previous subsection. Fig. 7(a) depicts the change in the numbers of active cases as the vaccine efficacy varies. The greater the vaccine efficacy, the more effective the vaccination is in suppressing the number of active cases. Remember that b vðtÞ ¼ hvð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. Fig. 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 a vaccine  with as high efficacy as possible, and a quite high quarantine rate to suppress the number of active cases. The total number of deaths will follow the same pattern.

Several scenarios of vaccination
Previously, the simulations were 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 at the minimum cost. It is obvious that high vaccination rate needs a large number of vaccines provided by the government, and as such it requires high expenses. We establish an approximation problem dealing with the 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 a feasible solution for the real problem, which are 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 the Indonesian government, denoted by K 2 (per day).
The model being observed first is the non age-structured one. Some of the results of this optimization problem using this model will be used to execute scenarios using the age-structured model later. Since comorbidities contribute to a high number of deaths, we assume that they 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 to find the answer to questions (Q1)-(Q3).

Non age-structured model
In accordance with the required constraints, we define the optimization problem as follows: Minimize where t 2 [t v , t V ] is defined in Section 2.1 and u represents the cost needed to vaccine a single person. Besides that, the objective function f has to meet the following constraints: As D represents the time interval of each vaccination, it is defined as 30 days, which means that the rate of vaccination will be kept to level off in certain values for every 30 days in a row. This assumption is based on the initial plan of the government on providing the vaccine in the term of months. According to (Biofarma, 2021) considering the readiness of the vaccine, the first vaccine inoculations are set in several phases having different length periods, e.g. 4 months (JaneApr 2021) for the first phase and 3 months later (MayeJul 2021) for the second phase of vaccination. Hence, choosing D ¼ 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 define two possible conditions; limited and good facilities, and analysed whether there are 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 an unconstrained one by adding a penalty value once the point violates the constraints following which it is solved numerically using the trust-region method. According to (Pemerintah RI Targetkan, 2020), the Indonesian government is able to provide about 31,000 COVID-19 testing 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 other hand, the government claimed that there will be 1 million vaccines per day nationwide (Indonesia masih kejar, 2020). 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 model using the assumption of limited and good conditions of facility has been examined with the vaccination starting in October 2020. Solving the optimum process problem with limited facility condition, we could not find any solution. Therefore, 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 compared with that coming from the model without vaccination so we can see their stark differences.

Starting times of vaccination
We set the vaccination program to start at different times: October 2020, January 2021 and after the peak of outbreak for each province by setting different t v for each simulation. The three scenarios are aimed to compare the expected situation if the first vaccine jab is delivered at different times. The value K 1 is adjusted for each province to be higher than the existing number of active cases at the time when the vaccination starts. The values of K 1 and K 2 are given in Table 4.
There are many possible v(t) as the solutions of the optimization problem. Fig. 8(a) and (b) are an example of v(t) that has a minimum value of Eq. (6) for Jakarta. Fig. 8(a) shows the graph of Q(t) that is successfully maintained under the constraint K 1 ¼ 30, 000. Fig. 8(b) shows the number of vaccination per day in Jakarta, represented by the blue thick lines. The number of susceptible people plotted in orange was drastically decreasing 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 Table 4 Values of K 1 and K 2 for each province.
Times of vaccination Jakarta Banten West Java October 2020 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 population. This scenario looks promising where the remaining active cases 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 media predicted the readiness of the vaccine is not earlier than January 2021. Fig. 9(a) and (b) show the number of active cases for Banten and West Java. It shows that the number of people needed to be vaccinated is about 2.9 million people or 23.3% of the population in Banten, and 12.5 million people or 26% of the population in West Java. Tables representing the estimation results of vaccination scheme for Banten and West Java are given in the Appendices.
If the starting time is in January 2021, Figs. 10(a)-(b) and Table 6 give the results for Jakarta. It shows that the number of susceptible people rapidly decreases only after the second period of vaccination. The fluctuating number of vaccination per   day scheduled for each period, plotted as the blue thick line, is one of the solutions found from the approximation process. It is interesting to see whether a constant rate of vaccination 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. From Table 6, the number of vaccinated people is 4.48 million or 43.2% of the population in Jakarta. The vaccination is successful because the remaining active case is only one. The simulations of the vaccination for Banten and West Java are given in Fig. 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. The starting time is therefore different for each province, i.e., April 2021 in Jakarta, October 2021 in Banten, and September 2021 in West Java. Fig. 12(a)   shows a steeper decrease of Q(t) than that of without vaccination in Jakarta. The number of death D(t) seems to stabilize to about 4,000 as the peak of outbreak happened before the vaccination. From Table 7, the total number of people to be vaccinated is 3.19 million or 30.4% of population in Jakarta. If we compare between Figs. 8(a), 10(a) and 12(a), the total numbers of death are increasing, which means that the later the starting time of vaccination, the higher the number of the pandemic casualties. Overall, the presence of vaccine jabs can cause a significant drop of susceptible people and 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.

Frequency of the vaccination
In Fig. 10(a), vaccination for 12 months can make the number of active cases decrease significantly in the first month of the vaccination period. It is interesting to see whether a one time vaccination in the first month could really work. Most of the time when the first attempt shows a good result, the vaccination program could be immediately stopped to reduce cost. Fig. 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 the other months. It shows a significant decrease of the active cases at the beginning of vaccination, but it will start to increase from August 2021 onward, which may cause another outbreak in the future. This prediction is clearly seen when we plot the graphs in a longer time-span in Fig. 14. On the left, a vaccination consistently given for 12 months from January 2021 can make active cases disappear until the beginning of year 2024. On the right, a one-time vaccination potentially creates 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 that obtained from the non-vaccine model with a few month delay.

Constant vaccination rates
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 that make it hard to be implemented. It means that there will be fluctuating needs of the vaccination workers per period. Now we simulate a the vaccination program with constant rates and analyze the dynamics of the active cases.  There are 3 scenarios with different values of constant rates starting from January 2021 based on the optimal vaccination rates obtained in 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 Fig. 15, where the original optimum values, the maximum, the average and the minimum rates are respectively represented 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 of the healthcare facility that requires a cost calculated from Eq. (6). The expenses calculated from implementing those scenarios will be compared to this benchmark expense. Fig. 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 a plot of Q(t) that resembles the plot from the optimum rates, but it requires almost twice of the 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 that the optimum fluctuated rates are the best scheme among other constant rate scenarios based on the dynamic of active cases and the required cost.

Discussion
Considering the simulations for different starting time scenarios, the vaccination program ideally should start 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 the vaccination program is delayed, the required amount of vaccines will decrease. It therefore seems better to begin a vaccination program after the peak of the outbreak. However, delaying the starting time will result in increasing total deaths. In fact, the decrease of the total number of deaths is not significant if we choose the latest scenario.
Having simulated different frequency of vaccination, it shows in Fig. 13(a) that only a single time of vaccination gives insignificant reduction of the number of 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 that could trigger another peak of outbreak through reinfection. It is concluded that a consistent schedule of vaccination is necessary to significantly reduce the number of active cases.
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 the vaccination schedule. It seems that 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 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.

Age-structured model
Now four scenarios of prioritizing particular age groups are implemented on the age-structured 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 targeted group in each period. The starting time is January 2020 and we use the scheme of vaccination rates shown in Fig. 10.

Groups of active and older adults (above 50)
During the first six months, vaccine shots are given to workers and high-risk people, which are 20 and older, then in the remaining months, they are given to younger groups. As seen in Fig. 17(a), the darker color of the cells for certain periods means higher rate of vaccination. Fig. 17(b) illustrates the simulation result of the total numbers Q(t) and D(t) of all age-groups N. Nuraini, K.K. Sukandar, P. Hadisoemarto et al. Infectious Disease Modelling 6 (2021 using this scenario in Jakarta. The simulation results from Banten and West Java are given in Fig. 18. The vaccination performs well in reducing the numbers of active cases and death toll.

Groups of 20 years old and younger persons
In this scenario, the first interval of six months is the vaccination time only for the younger people of age 0e19, and the remaining time is for the active and older people's groups, as scheduled in Fig. 19(a). The dynamics of the total numbers Q(t) and D(t) of all age-groups using this scenario in Jakarta, Banten and West Java are described in Figs. 19(b) and 20, respectively. It seems that this scenario gives insignificant decrease on the total numbers of active cases and death toll.

Alternating targeted group
In this scenario, we change the target of age-groups in certain ways that is shown in Fig. 21(a). High rate vaccination in the first month is given to the active and older people, i.e., age 20 and older. Figs. 21(b) and 22 describe the dynamics of the 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 a thicker tail in the active cases. This observation can be seen in every region we observed. Thus, this scenario is not likely preferable than the first one, i.e., prioritizing the active and older people.

Only the active people
Assuming that the Indonesian government plans to conduct vaccination only to the group of active people, i.e., age 20e49, we perform simulations in this scenario. Fig. 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. Fig. 24 shows that vaccination prioritizing the active and older people is more preferable due to its significant reduction of cases. This argument is also valid to the results given in Banten and West Java as depicted by Fig. 25.

Discussion
Our simulations show that the vaccination scenario by first targeting the active and older adults (above 50) (group 3e5) is better than the other considered scenarios. The numbers of active cases and total deaths decreased significantly. In the result of scenario 2, the decrease is insignificant, because the transmission rate of virus among the 20-year-old and younger people is much lower than that of the other age groups, as shown in Table 2.  After conducting simulations by targeting the active people only, we obtain that the dynamics of the number of active cases from the other age groups, Q i (t), i s 3, also decreased once people in the third age-group has been vaccinated. It is understandable because the key of the disease transmission is the contact among people. As the active people have larger access to people from the other age-groups, the transmission across age groups becomes highly possible. Once the biggest source of infectious people has been vaccinated, the number of active cases in the other groups of age can be expected to decrease significantly.

Conclusions and recommendations
Having constructed the SIQRD model with and without age structures and developed several scenarios on the implementation of a vaccination program, we concluded our work in the following findings that are related to the proposed questions: Modifying the SIRD into the SIQRD model by adding quarantine, reinfection, and even vaccination aspects has been considered to be capable in representing how COVID-19 spreads in several provinces in Indonesia. Utilizing the existing data and information related to those provinces, the graphs of the simulations well resemble the corresponding figure in the real situation. Vaccination should be implemented in the early stage of the pandemic. This is to suppress the number of active cases immediately, and consequently the total deaths. After the active cases reach their peaks, vaccination does not reduce the total deaths significantly. Vaccination should be implemented consistently following a schedule for a certain period. An implementation for only one or two months does not reduce the number of infectious persons, and eventually it will fail to prevent another outbreak. Prioritizing the active and older adults (above 50) for vaccination over others and prioritizing the active people only reduce significantly the total deaths.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgement
NN and KKS gratefully acknowledges the financial support provided by PDUPT 2021.

Appendix A. Graphs of non age-structured model
For Banten and West Java, the comparisons of simulated and real data are following. Those figures above are depicted by the estimation results of b, g, d, and I(t 0 ) for Banten and West Java. The obtained parameters of the non age-structured model are given by the following table.

Appendix B. Estimated values of transmission rate
We provide tables representing Obtained values b i for the age-structured model of Banten and West Java.

Appendix C. Vaccination scheme
We provide tables representing the estimation results of the vaccination scheme for Banten and West Java with several different starting times.

Banten
Vaccination in October 2020;  After the peak of outbreak. 2. West Java Vaccination in October 2020;  After the peak of outbreak.