Front. Appl. Math. Stat. Frontiers in Applied Mathematics and Statistics Front. Appl. Math. Stat. 2297-4687 Frontiers Media S.A. 10.3389/fams.2022.1096141 Applied Mathematics and Statistics Original Research Dynamic analysis and optimal control of COVID-19 with comorbidity: A modeling study of Indonesia Rois Muhammad Abdurrahman 1 Fatmawati 1 * Alfiniyah Cicik 1 Chukwu Chidozie W. 2 1Department of Mathematics, Faculty of Science and Technology, Universitas Airlangga, Surabaya, Indonesia 2Department of Mathematics, Wake Forest University, Winston-Salem, NC, United States

Edited by: Asep K. Supriatna, Universitas Padjadjaran, Indonesia

Reviewed by: Ebenezer Bonyah, University of Education, Winneba, Ghana; Pankaj Tiwari, University of Kalyani, India

 *Correspondence: Fatmawati ✉ fatmawati@fst.unair.ac.id

This article was submitted to Mathematical Biology, a section of the journal Frontiers in Applied Mathematics and Statistics

 06 01 2023 2022 8 1096141 11 11 2022 07 12 2022 Copyright © 2023 Rois, Fatmawati, Alfiniyah and Chukwu. 2023 Rois, Fatmawati, Alfiniyah and Chukwu

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Comorbidity is defined as the coexistence of two or more diseases in a person at the same time. The mathematical analysis of the COVID-19 model with comorbidities presented includes model validation of cumulative cases infected with COVID-19 from 1 November 2020 to 19 May 2021 in Indonesia, followed by positivity and boundedness solutions, equilibrium point, basic reproduction number (R0), and stability of the equilibrium point. A sensitivity analysis was carried out to determine how the parameters affect the spread. Disease-free equilibrium points are asymptotically stable locally and globally if R0 < 1 and endemic equilibrium points exist, locally and globally asymptotically stable if R0 > 1. In addition, this disease is endemic in Indonesia, with R0 = 1.47. Furthermore, two optimal controls, namely public education and increased medical care, are included in the model to determine the best strategy to reduce the spread of the disease. Overall, the two control measures were equally effective in suppressing the spread of the disease as the number of COVID-19 infections was significantly reduced. Thus, it was concluded that more attention should be paid to patients with COVID-19 with underlying comorbid conditions because the probability of being infected with COVID-19 is higher and mortality in this population is much higher. Finally, the combined control strategy is an optimal strategy that provides an effective guarantee to protect the public from the COVID-19 infection based on numerical simulations and cost evaluations.

 COVID-19 comorbidity stability sensitivity analysis optimal control cost evaluation 1. Introduction

The COVID-19 virus was reported in the Wuhan-Hubei Province of China in December 2019 and was spread rapidly to various parts of the world [1–6]. Symptoms are usually mild and appear gradually. In general, the symptoms of COVID-19 are fever, dry cough, and tiredness. In addition, there are other symptoms such as chest pain and tenderness, nasal congestion, headache, conjunctivitis, diarrhea, loss of sense of taste or smell, skin rash, or discoloration of the fingers or toes [6]. The symptoms experienced are usually mild and appear gradually. Furthermore, moderate and severe infection symptoms can occur in humans and appear gradually, such as having fever and cough accompanied by difficulty in breathing or shortness of breath, chest pain, and others [1, 6]. Individuals with previous comorbidity (such as diabetes, lung, and heart disease) are more likely to develop severe disease with stronger COVID-19 symptoms than individuals who do not have a comorbidity [7, 8]. In the case of COVID-19 comorbidity in Indonesia, 12 different diseases have been recorded, which range from the most at risk to the least at risk, namely hypertension, diabetes mellitus, heart, pregnancy, lung, kidney, immune disorders, cancer, other respiratory disorders, asthma, tuberculosis, and liver [9].

The first case in Indonesia was reported directly by President Jokowi Widodo on 2 March 2020 and there were as many as two people infected, namely a mother and child suspected of contracting it from a Japanese citizen [10]. Data from web [11], on 2 October 2020, to be precise, indicate that Indonesia was ranked 23 out of 215 countries reported being infected with 295,499 confirmed cases, 10,972 reported deaths, and 221,340 reported recoveries. Meanwhile, according to data on 14 June 2021, Indonesia was ranked 18 out of 222 countries reported being infected, with 1,919,547 confirmed cases, 53,116 reported deaths, and 1,751,234 reported recoveries.

The increasing number of COVID-19 cases requires a control strategy to control the COVID-19 outbreak. Control technique isolation and individual quarantine are the most efficient measures whenever a new outbreak occurs in a region without a vaccine or therapy [12, 13]. Several appeals or mitigations from WHO to control COVID-19 are social distancing, use of masks in public places, and intensive contact tracing (tracing) followed by quarantine of individuals who have the potential to contract the disease, and isolation of infected individuals in hospitals or independently [14]. Therefore, public education plays an important role in controlling the outbreak because it can convey information regarding how to prevent and reduce the transmission of COVID-19.

Furthermore, it is necessary to use mathematical modeling to determine the spread of COVID-19 infection and whether the control measures are effective. WHO also acknowledges that mathematical modeling can help health decision-makers (doctors or health professionals) and policymakers make decisions or find solutions (governments) [15]. The Susceptible-Infected-Removed (SIR) model is a mathematical representation of how diseases spread. The SIR model was first developed in 1927 by Kermack and McKendrick, who established it as a reference work and contributed significantly to the development of the mathematical theory of disease transmission [16, 17]. Several studies are related to the spread of disease, for example, research on the Coronavirus that caused SARS [18] and MERS [19, 20].

Soewono [21] applied the SEIR model, which has four subpopulations: susceptible (S), exposed (E), infected (I), and recovered (R), to simulate the spread of COVID-19. This model is an improvement on the SIR COVID-19 model. Furthermore, Das et al. [22] add a subpopulation of C (infected with comorbidity), so that the population is divided into five subpopulations, namely S, E, I, C, and R. The comorbidity referred to in this study is a general congenital disease, while research from Omame et al. [23] also proposed a comorbidity COVID-19 model, Omame et al. model coinfection with comorbidities (especially diabetes mellitus). So, Omame et al. built a model by dividing the population into eight subpopulations, namely susceptible (S), susceptible to comorbidity (Sc), individuals infected with COVID-19 without comorbidities (I), isolation and hospitalization for individuals infected with COVID-19 without comorbidity (H), recovered from COVID-19 but without comorbidity (R), infected with COVID-19 and comorbidity (Ic), isolation and hospitalization for those infected with COVID-19 and comorbidity (Hc), and recovered from COVID-19 but with comorbidity (Rc). In another study, Jia et al. [24] by incorporating subpopulations of isolation (H) and quarantine (Q), the model provided divides the population into seven subpopulations, namely S, E, I, A, Q, H, and R. The model is also based on the most recent data from the WHO, indicating that susceptible individuals must first be quarantined to stop the further spread. Research on COVID-19 was also conducted by Prathumwan et al. [25] by adding quarantine subpopulations (Q) and isolation (H) as well so that the model constructed has six subpopulations, namely S, E, I, Q, H, and R. The mathematical model that has been formed needs control to reduce the number of COVID-19 infections. Researchers discussing control issues include Deressa and Duress [26], Olaniyi et al. [27], and Das et al. [22]. Deressa and Duressa provide three controls, namely public education, protecting yourself from COVID-19 infection (such as wearing masks, washing hands, and maintaining distance), and treating individuals infected with COVID-19 in hospitals. In comparison, Olaniyi et al. provide two controls, public education and individual care management in hospitals. Other researchers, Das et al. [22], provide two controls to reduce the number of infected with comorbidity and without comorbidity, namely the control other than using drugs and the vaccination process. There are many studies related to COVID-19 besides those mentioned above, see for example the following literature studies [28–54].

By combining the research of Das et al. [22], Jia et al. [24], and Prathumwan et al. [25], the COVID-19 model will be constructed in this study. The discussion is divided into the following sections: The model formulation is presented in Section 2 followed by model validation and mathematical analysis in Section 3. A numerical simulation of the model without control is given in Section 4. Section 5 presents the model with controls and its simulation is given in Section 6. The last discussion on cost evaluation is presented in Section 7. The study is concluded with some key points in Section 8.

 2. Model formulation

We consider the new COVID-19 model with eight subpopulations, as shown in the compartment diagram in Figure 1.

Compartmental diagram of the COVID-19 model with comorbidity.

Based on the compartmental diagram in Figure 1 and the model assumptions, we have the following system of differential equations: dSdt = π-β1SIN-β2SCN-q1S-μS,dEdt = β1SIN+β2SCN-αE-μE,dIdt = ξαE-h1I-r1I-d1I-μI,dCdt = (1-ξ)αE-h2C-r2C-d2C-μC,dQdt =q1S-ρ1β1QIN-ρ2β2QCN-r3Q-μQ,dHdt = θh1I+δh2C+pρ1β1QIN+qρ2β2QCN   -r4H-d3H-μH,dJdt = (1-θ)h1I+(1-δ)h2C   +(1-p)ρ1β1QIN+(1-q)ρ2β2QCN-r5J-d4J-μJ,dRdt =r1I+r2C+r3Q+r4H+r5J-μR. In this model, the COVID-19 model is divided into susceptible (S), exposed (E), infected without comorbidity (I), infected with comorbidity (C), isolated (Q), treatment isolated (H), isolated without treatment (J), and recovered (R). Susceptible subpopulation increases with the recruitment or birth rate denoted by π and can be infected due to contact with infected individuals without comorbidity and with comorbidity denoted by β1 and β2, respectively. Susceptible individuals who are quarantined are denoted by q1 and cannot be returned to being susceptible due to the effects of public anxiety, which make some assumptions or opinions that susceptible individuals need to be quarantined, so that if quarantine is successful, then recovery is denoted by r3 and if not successful due to contact with infected individuals without comorbidity and with comorbidity, showing symptoms of being infected, then isolation is denoted by ρ1 and ρ2, respectively. Furthermore, p and q are the proportion of changes from quarantine to isolation. The progression from exposed to infection is denoted α, and ξ is the proportion of change from exposed to infection without comorbidity. From the infected subpopulation without comorbidity and with comorbidity, isolation is denoted by h1 and h2. The parameters r1, r2, r3, and r4 indicate the recovery rate of the subpopulations infected without comorbidity, infected with comorbidity, quarantine, isolated with treatment, and isolated without treatment. Furthermore, deaths from each subpopulation are denoted by μ and deaths from COVID-19 in subpopulations I, C, H, and J are denoted by d1, d2, d3, and d4.

 3. Mathematical analysis 3.1. Model validation

We calibrate our model (Equation 1) using cumulatively confirmed COVID-19 cases for Indonesia. We have retrieved COVID-19 case data from the Republic of Indonesia Task Force (SATGAS) situation report for the period 1 November 2020 to 19 May 2021 [9]. The parameter fitting uses the lsqcurvefit command, and the value of MAPE = 0.026022 is obtained. The results of the fitting parameters seem to match the infection case data as shown in Figure 2, and new parameter values are obtained according to conditions in Indonesia as follows in Table 1.

Parameter fitting results from the COVID-19 model.

Parameter values according to the fitting of the infected cases of COVID-19 in Indonesia.

Parameter Value Parameter Value Parameter Value
π 3783175.865 r5 0.088554 δ 0.00059843
β1 0.65799 h1 0.007884 p 0.090862
β2 0.79664 h2 0.00034162 q 0.28312
q1 0.16574 ρ1 0.99779 d1 0.00086579
r1 0.0068295 ρ2 0.9533 d2 0.022871
r2 0.0025349 α 0.25098 d3 0.36203
r3 0.030397 ξ 0.022219 d4 0.76233
r4 0.31851 θ 5.812 × 10−5 μ 0.0138
3.2. Positivity and boundedness of solutions

The change in the total population is given by dNdt=dSdt+dEdt+dIdt+dCdt+dQdt+dHdt+dJdt+dRdt,     =π-μN-d1I-d2C-d3H-d4J,          ≤π-μN, whose solutions give N(t)≤πμ+(N(0)-πμ)e-μt. Consequently as t → ∞, then limt→∞N(t)≤πμ . So, we can conclude that N is boundedness to N(t)≤πμ.

Considering the above solutions, we have that the model has a boundedness solution which is contained in a feasible region Ω, where Ω={(S, E, I, C, Q, H, J, R) | N(t)≤πμ}. Next, we show the positivity of solving the Equation (1) system by following Riyapan et al. [42] and Rois et al. [46], as follows:

Theorem 1. Let S, E, I, C, Q, H, J, and R be the system solutions (Equation 1). If S(0) ≥ 0, E(0) ≥ 0, I(0) ≥ 0, C(0) ≥ 0, Q(0) ≥ 0, H(0) ≥ 0, J(0) ≥ 0, and R(0) ≥ 0, then all solutions are positive for every t ≥ 0.

Proof.

Take the first equation of the system (Equation 1) as follows: dSdt=π-β1SIN-β2SCN-q1S-μS. Let η=β1IN+β2CN.                                     dSdt=π-S(η+q1+μ),    dSdt+S(η+q1+μ)=π,d(e(q1+μ)t+∫0tηdsS(t))dt=πe(q1+μ)t+∫0tηds, then a homogeneous solution is obtained d(e(q1+μ)t+∫0tηdsS(t))dt=0,                                S(t)=ke-(q1+μ)t-∫0tηds. Thus, let us assume that the solution is non-homogeneous S(t)=ke-(q1+μ)t-∫0tηds. Next, substituting the Equation (3) into the Equation (2) to get   dk(t)dt=πe(q1+μ)t+∫0tηds,    k(t)=∫0tπe(q1+μ)y+∫0yηdxdy+K. The Equation (4) is substituted into the Equation (3), we get

S(t)=∫0tπe(q1+μ)y+∫0yη dxdy×e-(q1+μ)t-∫0tη ds+S(0)e-(q1+μ)t-∫0tη ds. So, S(t) is positive for t ≥ 0.

Take the fifth equation of the system (Equation 1) as follows: dQdt=q1S-ρ1β1QIN-ρ2β2QCN-r3Q-μQ. Let ω=ρ1β1IN+ρ2β2CN.                                      dQdt=q1S-Q(ω+r3+μ),   dQdt+Q(ω+r3+μ)=q1S,d(e(r3+μ)t+∫0tωdsQ(t))dt=q1Se(r3+μ)t+∫0tωds, then a homogeneous solution is obtained d(e(r3+μ)t+∫0tωdsQ(t))dt=0,                                Q(t)=ke-(r3+μ)t-∫0tωds. Thus, let us assume that the solution is non-homogeneous Q(t)=ke-(r3+μ)t-∫0tωds. Next, substituting the Equation (6) into the Equation (5) to get dk(t)dt=q1Se(r3+μ)t+∫0tωds,k(t)=∫0tq1Se(r3+μ)y+∫0yωdxdy+K. The the Equation (7) is substituted into the Equation (6), we get

Q(t)=∫0tq1Se(r3+μ)y+∫0yω dxdy×e-(r3+μ)t-∫0tω ds+Q(0)e-(r3+μ)t-∫0tω ds. So, Q(t) is positive for t ≥ 0.

Take the second equation of the system (Equation 1) as follows: dEdt=β1SIN+β2SCN-αE-μE≥-αE-μE,

or     dE(t)dt≥-E(α+μ),∫dE(t)E≥∫-(α+μ)dt,     E(t)≥e-(α+μ)t+k,     E(t)≥E(0)e-(α+μ)t. Thus, E(t) is positive for t ≥ 0. Furthermore, in the same way as proof number 3, I(t), C(t), H(t), J(t), and R(t) can be shown respectively to be positive.     □

 3.3. Equilibrium point and basic reproduction number

The equilibrium point of the system (Equation 1) is obtained by setting the right-hand side to zeros. Therefore, the first equilibrium point is obtained, namely the disease-free equilibrium point, as follows: X0=(S0, E0, I0, C0, Q0, H0, J0, R0)     =(πa1, 0, 0, 0, πq1a1a5,0, 0, πq1r3a1a5μ). Where a1 = q1 + μ and a5 = r3 + μ.

Furthermore, the basic reproduction number, denoted by R0, is obtained using the next-generation matrix method [55, 56]. The constituent components of the next-generation matrix method only consist of infected subpopulation groups, namely: f=[β1SIN+β2SCN00], and v=[a2E-ξαE+a3I-(1-ξ)αE+a4C]. The partial derivative evaluated at X0 gives         F(X0)=[0β1a5μ(a5μ+q1μ+q1r3)β2a5μ(a5μ+q1μ+q1r3)000000]and V(X0)=[a200-ξαa30-(1-ξ)α0a4]. The inverse of the V(X0) matrix is

V-1=[1a200ξαa2a31a30(1-ξ)αa2a401a4]. Based on the F(X0) and V−1(X0) matrices, the next-generation matrix FV−1 can be formed so that we can obtain

FV-1=[β1a5μξαa2a3(a5μ+q1μ+q1r3)+β2a5μα(1-ξ)a2a4(a5μ+q1μ+q1r3)β1a5μa3(a5μ+q1μ+q1r3)β1a5μa4(a5μ+q1μ+q1r3)000000].

So, the basic reproduction number is obtained based on the eigenvalues of the FV−1 matrix as follows: R0=ρ(M)=a5μα(β1ξa4+β2a3(1-ξ))a2a3a4(a5μ+q1μ+q1r3). Next, the second equilibrium point is obtained, namely the endemic equilibrium point X* = (S*, E*, I*, C*, Q*, H*, J*, R*) with S*=πa5A1R0,E*=πa2R0(R0-1), I*=πξαa2a3R0(R0-1),C*=(1-ξ)παa2a4R0(R0-1),Q*=a2a3a4a5q1πA3,H*=(R0-1)A4, J*=(R0-1)A5,R*=(R0-1)A6+r3a2a3a4a5q1πμA3, with A1 = a5μ + μq1 + q1r3, A2 = ρ1β1a4μαξ + ρ2β2a3μα(1 − ξ), A3 = A1A2(R0 − 1) + a2a3a4a5A1R0, A4=παa6R0(θh1ξa2a3+δh2(1-ξ)a2a4+pρ1β1q1a4a5μξA3+qρ2β2q1a3a5μ(1-ξ)A3), A5=παa7R0((1-θ)h1ξa2a3+(1-δ)h2(1-ξ)a2a4+(1-p)ρ1β1μξa4a5q1A3+(1-q)ρ2β2μa3a5q1(1-ξ)A3), and A6=r1παξμa2a3R0+r2πα(1-ξ)μa2a4R0+r4A4μ+r5A5μ.

The existence of the endemic equilibrium point X* depends on the value of R0. If the value of R0 < 1 is taken, then the endemic equilibrium point X* does not exist because it is clear that E*, I*, C*, H*, and J* are obtained negative. If R0 = 1, then we get the equilibrium point X* = X0, which causes the equilibrium point X* not to exist. Furthermore, if R0 > 1, then we get S*, E*, I*, C*, Q*, H*, J*, and R* are positive and the endemic equilibrium point X* exists.

 3.4. Local stability

The local stability of the equilibrium point is obtained by linearizing system (Equation 1), which yields the following jacobian matrix below:

J=[B1-B2-a1B1B1-β1SNB1-β2SNB1B1B1B1B2-B1-B2-a2β1SN-B1β2SN-B1-B1-B1-B1-B10ξα-a3000000(1-ξ)α0-a40000q1+B3B3B3-ρ1β1QNB3-ρ2β2QNB3-B4-a5B3B3B3-B5-B5B6-B5B7-B5B8-B5-B5-a6-B5-B5-B9-B9B10-B9B11-B9B12-B9-B9-B9-a7-B900r1r2r3r4r5-μ]

With B1=β1SIN2+β2SCN2, B2=β1IN+β2CN, B3=ρ1β1QIN2+ρ2β2QCN2, $B4=ρ1β1IN+ρ2β2CN, B5=pρ1β1QIN2+qρ2β2QCN2, B6=θh1+pρ1β1QN, B7=δh2+qρ2β2QN, B8=pρ1β1IN+qρ2β2CN, B9=(1-p)ρ1β1QIN2+(1-q)ρ2β2QCN2, B10=(1-θ)h1+(1-p)ρ1β1QN, B11=(1-δ)h2+(1-q)ρ2β2QN, and B12=(1-p)ρ1β1IN+(1-q)ρ2β2CN.

3.4.1. Local stability of disease-free equilibrium point

Evaluating (8) at X0 yields

J(X0)=[-a10-β1C1-β2C100000-a2β1C1β2C100000ξα-a3000000(1-ξ)α0-a40000q10-ρ1β1C2-ρ2β2C2-a500000θh1+pρ1β1C2δh2+qρ2β2C20-a60000(1-θ)h1+(1-p)ρ1β1C2(1-δ)h2+(1-q)ρ2β2C200-a7000r1r2r3r4r5-μ],

With C1=a5μ(a5+q1)μ+q1r3 and C2=q1μ(a5+q1)μ+q1r3. The characteristic equation for the |J(X0) − λI| = 0 is as follows: (-a1-λ)(-a5-λ)(-a6-λ)(-a7-λ)(-μ-λ)|-a2-λβ1C1β2C1ξα-a3-λ0(1-ξ)α0-a4-λ|=0. Based on the Equation (9), we obtain the eigenvalues λ1 = λ2 = λ3 = λ4 = λ5 < 0. Therefore, the stability of the disease-free equilibrium point depends on M1=|-a2-λβ1C1β2C1ξα-a3-λ0(1-ξ)α0-a4-λ|. From the Equation (10), we obtain the following characteristic equation: λ3+k1λ2+k2λ+k3=0, with

k1 =a2+a3+a4,k2 =a2a3(1-R0)+a2a4(1-R0)+a3a4+a4β1αC1ξa3+a3β2αC1(1-ξ)a4, andk3 =a2a3a4(1-R0).

In the Equation (11), it is clear that k1 > 0, and if R0 < 1, then k2 > 0 and k3 > 0. Therefore, the stability property of the equilibrium point X0 is established using the Routh–Hurwitz criterion. Furthermore, the equilibrium point X0 is asymptotically stable if and only if it satisfies the following criteria:

k1 > 0,

k3 > 0, and

k1k2 − k3 > 0.

Criteria Equations (1) and (2) have been met so that the disease-free equilibrium point X0 is locally asymptotically stable if it meets k1k2 − k3 > 0 where k1k2-k3>0,(a2+a3+a4)(a2a4(1-R0)+a3a4+a4β1αC1ξa3+a3β2αC1(1-ξ)a4)+a22a3(1-R0)+a2a32(1-R0)>0. It is clear that the Routh–Hurwitz criteria are satisfied; thus, the roots of the characteristic Equation (11) have negative real parts. Therefore, the disease-free equilibrium point is asymptotically locally stable if R0 < 1.

 3.4.2. Local stability of the endemic equilibrium point

Evaluating (8) at X* yields

J(X*)=[B1-B13B1B1-β1S*NB1-β2S*NB1B1B1B1B2-B1-B14β1S*N-B1β2S*N-B1-B1-B1-B1-B10ξα-a3000000(1-ξ)α0-a40000q1+B3B3B3-ρ1β1Q*NB3-ρ2β2Q*NB3-B15B3B3B3-B5-B5B6-B5B7-B5B8-B5-B16-B5-B5-B9-B9B10-B9B11-B9B12-B9-B9-B17-B900r1r2r3r4r5-μ],

With B1=β1S*I*N2+β2S*C*N2, B2=β1I*N+β2C*N, B3=ρ1β1Q*I*N2+ρ2β2Q*C*N2, B4=ρ1β1I*N+ρ2β2C*N, B5=pρ1β1Q*I*N2+qρ2β2Q*C*N2, B6=θh1+pρ1β1Q*N, B7=δh2+qρ2β2Q*N, B8=pρ1β1I*N+qρ2β2C*N, B9=(1-p)ρ1β1Q*I*N2+(1-q)ρ2β2Q*C*N2, B10=(1-θ)h1+(1-p)ρ1β1Q*N, B11=(1-δ)h2+(1-q)ρ2β2Q*N, B12=(1-p)ρ1β1IN+(1-q)ρ2β2CN, B13 = B2 + a1, B14 = B2 + a2, B15 = B4 + a5, B16 = B5 + a6, and B17 = B9 + a7.

The characteristic equation of the |J(X*) − λI| = 0 is λ8+k1λ7+k2λ6+k3λ5+k4λ4+k5λ3+k6λ2+k7λ+k8=0, It is difficult to prove analytically that all eigenvalues of J have negative real parts for R0 > 1. However, from our numerical simulations (case R0 > 1), all eigenvalues have negative real parts.

Figure 3 gives the projection of three orbits of three different initial conditions when R0 > 1 on the I − C plane. The component (I, C) of the equilibrium X* is not (0, 0). This simulation indicates that the endemic equilibrium X* is locally asymptotically stable when R0 > 1.

Projection of three orbits of the model on the I − C plane.

 3.5. Global stability analysis

In this study, we prove the global stability of disease-free and endemic equilibrium points by constructing the suitable Lyapunov function and following the theorem from Alligood et al. [57].

3.5.1. Global stability of the disease-free equilibrium point

Theorem 2. Disease-free equilibrium point X0 is globally asymptotically stable if R0 < 1 and unstable if R0 > 1.

Proof. Defined Lyapunov function L=κ1E+κ2I+κ3C, where κ1=a1a2a3a4,κ2=πβ1a2a4N,andκ3=πβ2a2a3N. The function L needs to be proven to determine whether Lyapunov is strong or weak for X0. L(x⃗*)=L(S0, E0, I0, C0, Q0, H0, J0, R0),L(x⃗*)=κ1E0+κ2I0+κ3C0=0. It is proven that L(x⃗*)=0. Next, L(x⃗)=κ1E+κ2I+κ3C, Because ∀(S, E, I, C, Q, H, J, R) ≠ (S0, E0, I0, C0, Q0, H0, J0, R0), so it is proved that L(x⃗)>0.

Thus, the Equation (13) can be reduced to ∂L∂t=κ1dEdt+κ2dIdt+κ3dCdt     =κ1(β1SI+β2SCN-a2E)         +κ2(ξαE-a3I)+κ3((1-ξ)αE-a4C), so, we obtain =a1a2a3a4(β1SI+β2SCN-a2E+πβ1a2a4ξαEa1a2a3a4N+πβ2a2a3(1-ξ)αEa1a2a3a4N-πβ1Ia1N-πβ2Ca1N). Let S=πa1, so we get ∂L∂t=a1a22a3a4E(R0-1). Based on the description above, it can be concluded that ∂L∂t<0 if R0 < 1 and ∂L∂t=0 if E = 0. Hence, by Lasalle's invariance principle, the disease-free equilibrium point in the spread of COVID-19 (X0) is globally asymptotically stable if R0 < 1.     □

 3.5.2. Global stability of the endemic equilibrium point

Theorem 3. If R0 > 1, then the endemic equilibrium point X* is said to be globally asymptotically stable.

Proof. The Lyapunov function is defined as follows: L=12[SS+EE+II+CC+QQ+HH+JJ+RR]2, Where SS=(S-S*), EE=(E-E*), II=(I-I*), CC=(C-C*), QQ=(Q-Q*), HH=(H-H*), JJ=(J-J*), and dan RR=(R-R*).

The function L needs to be proven to determine whether the Lyapunov function is strong or weak for X*. L(x⃗*)=L(S*, E*, I*, C*, Q*, H*, J*, R*),L(x⃗*)=12[SS*+EE*+II*+CC*+QQ*+HH*+JJ*+RR*]2=0. Where SS*=(S*-S*), EE*=(E*-E*), II*=(I*-I*), CC*=(C*-C*), QQ*=(Q*-Q*), HH*=(H*-H*), JJ*=(J*-J*), and RR*=(R*-R*). It is proven that L(x⃗*)=0.

Next, L(x⃗)=12[SS+EE+II+CC+QQ+HH+JJ+RR]2. Because ∀(S, E, I, C, Q, H, J, R) ≠ (S0, E0, I0, C0, Q0, H0, J0, R0), so that it is proven that L(x⃗)>0. Next, we check that the Equation (14) is reduced to ∂L∂t=[SS+EE+II+CC+QQ+HH+JJ+RR]          ddt[S+E+I+C+Q+H+J+R],    =[SS+EE+II+CC+QQ+HH+JJ+RR]          [π-μ(S+E+I+C+Q+H+J+R)          -d1I-d2C-d3H-d4J], Let π=μ(S*+E*+I*+C*+Q*+H*+J*+R*)+d1I*+d2C*+d3H*+d4J*. So that it gives

=[SS+EE+II+CC+QQ+HH+JJ+RR]×[-[μ(SS+EE+II+CC+QQ+HH+JJ+RR)+d1(I-I*)+d2(C-C*)+d3(H-H*)+d4(J-J*)]],=-[SS+EE+II+CC+QQ+HH+JJ+RR]×[μ(SS+EE+II+CC+QQ+HH+JJ+RR)+d1(I-I*)+d2(C-C*)+d3(H-H*)+d4(J-J*)]]. Based on the description above, it can be concluded that ∂L∂t<0 if R0 > 1 and ∂L∂t=0 if S = S*, E = E*, I = I*, C = C*, Q = Q*, H = H*, J = J*, and R = R*. Hence, by Lasalle's invariance principle, it means that the endemic equilibrium point in the spread of COVID-19 (X*) is globally asymptotically stable if R0 > 1.     □

The terms of existence and the type of stability of the equilibrium point of the system of equations are summarized in Table 2.

Stability conditions.

Equilibrium point Existence requirement Global stability type Stability condition
X0 None Asymptotically stable R0 < 1
X* R0 > 1 Asymptotically stable R0 > 1 and strong Lyapunov function
3.6. Sensitivity analysis

The sensitivity analysis aims to determine the parameters that cause the spread of the COVID-19 virus. The sensitivity index of the basic reproduction number depends on the differentiation of the parameters contained in the basic reproduction number [58, 59]. Sensitivity index R0 to the parameters is as follows:

IαR0=∂R0∂ααR0=μα+μ,Iβ1R0=∂R0∂β1β1R0=β1ξa4β1ξa4+β2a3(1-ξ),Iβ2R0=∂R0∂β2β2R0=β2a3(1-ξ)β1ξa4+β2a3(1-ξ),IξR0=∂R0∂ξξR0=β1ξa4-β2ξa3β1ξa4+β2a3(1-ξ),Iq1R0=∂R0∂q1q1R0=-q1q1+μ,Ih1R0=∂R0∂h1h1R0=a4β1h1ξ-a3(β2a3(1-ξ)+β1a4ξ),Ir1R0=∂R0∂r1r1R0=β1r1a4ξ-a3(β2a3(1-ξ)+β1a4ξ),Id1R0=∂R0∂d1d1R0=β1d1a4ξ-a3(β2a3(1-ξ)+β1a4ξ),Ih2R0=∂R0∂h2h2R0=β2h2a3(1-ξ)-a4(β2a3(1-ξ)+β1a4ξ),Ir2R0=∂R0∂r2r2R0=β2r2a3(1-ξ)-a4(β2a3(1-ξ)+β1a4ξ),Id2R0=∂R0∂d2d2R0=β2r2a3(1-ξ)-a4(β2a3(1-ξ)+β1a4ξ),IμR0=∂R0∂μμR0=-(α((-2ξ+2)β2+2β1ξ)μ5+a8(1-ξ)β2+β1a9ξ)μ4+2a10((1-ξ)a3β2+β1ξa4)μ3a12a22a32a42-(((-a12)q1+a112)α+a112(a12+q1)(ξ-1)β2-(((a11)q1-a122)α-a122(a11+q1))β1ξ)μ2a12a22a32a42-2(a11αq1(β1ξ+β2(1-ξ))a11μ+a12α((ξ-1)a11β2-β1ξa12)q1a11)a12a22a32a42,

With a8 = (4d1 + α + d2 + h1 + h2 + q1 + 4r1 + r2), a9 = d1 + α + 4d2 + h1 + 4h2 + q1 + r1 + 4r2, a10 = (d1 + α + d2 + h1 + h2 + q1 + r1 + r2), a11 = h1 + r1 + d1, and a12 = h2 + r2 + d2.

The parameter sensitivity index is shown in Table 3.

Sensitivity analysis.

No Parameter Sensitivity index
1 q1 −0.9761
2 β2 0.9754
3 μ 0.5715
4 d2 −0.5636
5 r2 −0.0625
6 α 0.0522
7 β1 0.0246
8 h2 −0.0084
9 h1 −0.0066
10 r1 −0.0057
11 ξ 0.0025
12 d1 −0.0007

From Table 3, we can see that the most sensitive parameters are q1 and β2. A positive index means that if we reduce the parameter by almost 10%, then the value of the basic reproduction number can decrease by 10%.

 4. Simulation of the model without control

This section presents a numerical solution of system (1) using the Fourth-order Runge–Kutta method. The parameter values used in this simulation are shown in Table 1 and three different initial values. In this simulation, the stability of the disease-free equilibrium point is shown from the parameter values given in Table 1, except for the parameter q1 = 0.56574, the value is R0 = 0.456 < 1. Based on the value of these parameters, the disease-free equilibrium point is obtained, namely X0 = (6527542, 0, 0, 0, 83496205, 0, 0, 183499868).

Next, the graph of the solution R0 < 1 is obtained in Figure 4.

COVID-19 model solution graph for R0 < 1.

The analysis results from 3.4.1, and Theorem 2 are illustrated numerically. Numerical simulations with some initial values show that the graph solution is toward and close to the disease-free equilibrium point X0 (converging toward the disease-free equilibrium point X0). Based on the graph, this means that after all this time, no individual has been infected with COVID-19. The numerical results support the analysis that if R0 < 1, then the disease-free equilibrium point X0 is asymptotically stable locally and globally given different initial values.

We show the stability of the endemic equilibrium point, the parameter values in Table 1 are used, and three different values, so we get R0 = 1.47 > 1. Similarly, the disease-free equilibrium point is obtained X0 = (21071493, 0, 0, 0, 79018695, 0, 0, 174052991) and the endemic equilibrium point is given as

X*=(S*, E*, I*, C*, Q*, H*, J*, R*),     =(6031017, 10196616, 1933377, 63222893,           2119885, 362896, 798735, 30671327). The solution graph for the case of R0 > 1 is obtained as follows in Figure 5.

The COVID-19 model solution graph for R0 > 1.

The numerical simulation results support the analysis from 3.4.2 and Theorem 3 that some of the initial values given are obtained by the graph of the solution leading to the endemic equilibrium point X* (converging to the endemic equilibrium point X*), which means there is a spread of disease due to COVID-19. The numerical simulation results follow the analysis that if R0 > 1, the endemic equilibrium point X* is asymptotically stable locally and globally with different initial values. Based on the given parameter values, we obtain R0 > 1. This means that there is an outbreak of disease due to COVID-19. Therefore, it is necessary to take control measures to reduce the outbreak.

4.1. Effect parameters

The effect of parameters on R0 was analyzed using contour plots. We choose two significant parameters, q1 and β2, and provide a contour plot as a function of R0. The impact of some R0 parameters is further investigated in Figure 6. Figure 6 shows that increasing the parameter q1 decreases the value of R0. This implies that increasing quarantines has the effect of reducing the spread of COVID-19. Meanwhile, an increase in the β2 parameter resulted in an increase in the R0 value and implied that an increase in contacts of individuals with comorbidities would increase the spread of COVID-19, especially individuals with comorbidities. Therefore, increasing quarantine and reducing contact with individuals with comorbidities are important.

Effect parameters of R0 (β2 × q1 ∈ [0.1 : 0.79664] × [0.1 : 0.16574]).

 5. Optimal control problem 5.1. Comorbidity COVID-19 model with optimal control

The control variable given to the COVID-19 model consists of preventive measures through education (u1) and individual treatment efforts for infected (u2). So the model with control is given as follows: dSdt=π-(1-u1)(β1SI+β2SC)N-q1S-μS,dEdt=(1-u1)(β1SI+β2SC)N-αE-μE,dIdt=ξαE-(h1+u2)I-r1I-d1I-μI,dCdt=(1-ξ)αE-(h2+u2)C-r2C-d2C-μC,dQdt=q1S-ρ1β1QIN-ρ2β2QCN-r3Q-μQ,dHdt=(θh1+u2)I+pρ1β1QIN+qρ2β2QCN+(δh2+u2)C-r4H-d3H-μH,dJdt=((1-θ)h1+u2)I+((1-δ)h2+u2)C+(1-p)ρ1β1QIN+(1-q)ρ2β2QCN-a7J,dRdt=r1I+r2C+r3Q+r4H+r5J-μR.

The function that minimizes the number of infected cases without comorbidity (I) and the number of infected cases with comorbidity (C) over a time interval [0, T] can be defined as J(u1,u2)=∫0T(I(t)+C(t)+12(A1u12+A2u22)), Where A1 and A2 are the relative cost associated with the controls u1 and u2, and T is the final time. The aim of the control is to minimize the cost function. J(u1*,u2*)=min J(u1,u2), Subject to the system (Equation 15), where 0 ≤ (u1, u2) ≤ 1 and t ∈ (0, T).

 5.2. Optimal control analysis

The Hamilton function can be defined as follows: H=I+C+12(A1u12+A2u22)+λ1(π-(1-u1)(β1SI+β2SC)N-q1S-μS)+λ2((1-u1)(β1SI+β2SC)N-αE-μE)+λ3(ξαE-(h1+u2)I-r1I-d1I-μI)+λ4((1-ξ)αE-(h2+u2)C-r2C-d2C-μC)+λ5(q1S-ρ1β1QIN-ρ2β2QCN-r3Q-μQ)+λ6((θh1+u2)I+pρ1β1QIN+qρ2β2QCN+(δh2+u2)C-r4H-d3H-μH)+λ7(((1-θ)h1+u2)I+((1-δ)h2+u2)C+(1-p)ρ1β1QIN+(1-q)ρ2β2QCN-a7J)+λ8(r1I+r2C+r3Q+r4H+r5J-μR). Based on Pontryagin's principle, the Hamilton function will reach an optimal solution if it satisfies the state equation and the costate equation, and the condition is stationary. The equation state is obtained by deriving the Hamilton function (Equation 17) for each variable costate as Equation (15). Next, the equation costate is the negative value of the derivative of the Hamilton function (Equation 17) for each variable state as follows: dλ1dt=-∂H∂S=(λ1-λ2)((1-u1)(β1I+β2C)N+(1-u1)(β1SI+β2SC)N2)+q1(λ1-λ5)+λ1μ+ρ1β1QIN2(pλ6+(1-p)λ7-λ5)+ρ2β2QCN2(qλ6+(1-q)λ7-λ5),dλ2dt=-∂H∂E=(λ2-λ1)(1-u1)(CSβ2+ISβ1)N2+λ2α+λ2μ-λ3αξ-λ4α(1-ξ)+ρ1β1QIN2(λ6p+(1-p)λ7-λ5)+ρ2β2QCN2(λ6q+(1-q)λ7-λ5),dλ3dt=-∂H∂I=(λ1-λ2)((1-u1)Sβ1N-(1-u1)Sβ1N)+h1(λ3-θλ6-(1-θ)λ7)+u2(λ3-λ6-λ7)+r1(λ3-λ8)-λ3(-d1-μ)+ρ1β1QN(λ5-λ6p-λ7(1-p))-1,dλ4dt=-∂H∂C=(λ2-λ1)(-(1-u1)Sβ2N+(1-u1)(CSβ2+ISβ1)N2)+h2(λ4-λ6δ-λ7(1-δ))+r2(λ4-λ8)+u2(λ4-λ6-λ7)-λ4(-d2-μ)+ρ1β1QIN2(-λ5+λ6p+λ7(1-p))+(ρ2β2QN-ρ2β2QCN2)(λ5-λ6q-λ7(1-q))-1,dλ5dt=-∂H∂Q=(λ2-λ1)((1-u1)(CSβ2+ISβ1)N2)+r3(λ5-λ8)+(ρ1β1IN-ρ1β1QIN2)(λ5-pλ6-(1-p)λ7)+(ρ2β2CN-ρ2β2QCN2)(λ5-λ6q-λ7(1-q))+λ5μ,dλ6dt=-∂H∂H=(λ2-λ1)((1-u1)(CSβ2+ISβ1)N2)+r4(λ6-λ8)+ρ1β1QIN2(-λ5+λ6p+λ7(1-p))+ρ2β2QCN2(-λ5+λ6q+λ7(1-q)),dλ7dt=-∂H∂J=(λ2-λ1)((1-u1)(CSβ2+ISβ1)N2)+r5(λ7-λ8)+ρ1β1QIN2(-λ5+λ6p+λ7(1-p))+ρ2β2QCN2(-λ5+λ6q+λ7(1-q))+λ7(d4+μ),dλ8dt=-∂H∂R=(λ2-λ1)((1-u1)(CSβ2+ISβ1)N2)+ρ1β1QIN2(-λ5+λ6p+λ7(1-p))+ρ2β2QCN2(-λ5+λ6q+λ7(1-q))+λ8μ. With transverse condition λ1(T) = λ2(T) = λ3(T) = λ4(T) = λ5(T) = λ6(T) = λ7(T) = λ8(T) = 0.

The stationary condition for the optimal control problem (16) is obtained by deriving the Hamilton function (17) on the control variables u1 and u2 (∂H∂u1=0, ∂H∂u2=0), successively obtained

u1=(β1SI+β2SC)(λ2-λ1)NA1 andu2=I(λ3-λ6-λ7)+C(λ4-λ6-λ7)A2. The control variables in the COVID-19 model with preventive measures through education and treatment efforts for infected individuals are defined as 0 ≤ u1 ≤ 1 and 0 ≤ u2 ≤ 1. So, the optimal control u1* and u2* can be expressed as

u1*=maks{0,min((β1S*I*+β2S*C*)(λ2-λ1)NA1,1)  } andu2*=maks{0,min(I*(λ3-λ6-λ7)+C*(λ4-λ6-λ7)A2,1)}. The optimal system is obtained by substituting the optimal control variables u1* and u2* into the system of state (Equation 15) and costate (Equation 18) equations.

 6. Simulation of the model with control

The method used in solving this optimal control problem is the forward–backward sweep method. In this numerical simulation, the parameter values used are presented in Table 1 according to the state of the COVID-19 case in Indonesia. Next, the initial values given are as follows S0 = 270, 911, 990, E0 = 1, 000, 000, I0 = 412, 784, C0 = 500, 000, Q0 = 100, 000, H0 = 56, 899, J0 = 200, 000, and R0 = 341, 942, with simulation intervals t ∈ [0, 100]. The results of the optimal numerical control simulation are presented as follows:

6.1. Using the control strategy <italic>u</italic><sub>1</sub> ≠ 0 and <italic>u</italic><sub>2</sub> = 0

Figure 7 shows the control strategy u1 ≠ 0 and u2 = 0. This is the result of education that makes all individuals always be careful, such as interacting outside the home. Implementing this strategy on subpopulations exposed, infected with comorbidity, infected without comorbidity, and isolation is significantly reduced. Implementing the control strategy u1 ≠ 0 and u2 = 0 can also increase the subpopulation of quarantine. Furthermore, the control strategy profiles u1 ≠ 0 and u2 = 0 to reduce the number of COVID-19 cases during t = 100 are presented in Figure 8.

The optimal control simulation results with u1 ≠ 0.

The optimal control profile with u1 ≠ 0.

The control strategy u1 ≠ 0 and u2 = 0 is given by one (maximum) from the beginning of the period to t = 99.9 and decreases significantly to zero at the end of the period. Control is terminated at the period's end, meaning no more control is given.

 6.2. Using the control strategy <italic>u</italic><sub>1</sub> = 0 and <italic>u</italic><sub>2</sub> ≠ 0

Figure 9 shows the control strategy u1 = 0 and u2 ≠ 0. This strategy can reduce the subpopulation infected with comorbid and without comorbid because there is an increase in the care of infected individuals. Implementing this strategy on subpopulations exposed, infected with comorbidity, infected without comorbidity, and isolation is significantly reduced. Implementing the control strategies u1 = 0 and u2 ≠ 0 can also increase the subpopulation of quarantine. Furthermore, the profiles of u1 = 0 and u2 ≠ 0 control strategies to reduce the number of COVID-19 cases for t = 100 are presented in Figure 10.

The optimal control simulation results with u2 ≠ 0.

Optimal control profile with u2 ≠ 0.

The control strategy u1 = 0 and u2 ≠ 0 is given by one (maximum) from the beginning to t = 48.6 and decreases slowly until t = 100 reaches zero which means control is stopped at the end of the period.

 6.3. Using the control strategy <italic>u</italic><sub>1</sub> ≠ 0 and <italic>u</italic><sub>2</sub> ≠ 0

Figure 11 indicates a combined control strategy. This results from education that makes all individuals always be careful (such as interacting outside the home) and increased care for infected individuals. Combined control strategies can control or reduce deployment significantly. Furthermore, the combined control strategy profile to reduce the number of COVID-19 cases during t = 100 is presented in Figure 12.

Combined optimal control simulation results.

Combined control profile.

The combined control strategy consists of two concurrent administrations of control. The control u1 is given by one (maximum) up to t = 65.4 and decreases significantly until the end of the period reaches zero. Then, control u2 is assigned one (maximum) until t = 39.5 and then decreases until t = 100 slowly reaches zero. Both controls are terminated at the end of the period, which means they are no longer given control of u1 and u2.

 6.4. Comparison of total infections using all strategic control scenarios

The varying initial values of the exposed subpopulations are given. Total infected subpopulations for different initial conditions from exposed subpopulations are E(0) = 200000, E(0) = 1000000, E(0) = 10000000, and E(0) = 100000000 successively from left to right using the three control strategies shown in Figure 13.

(A–D) Total subpopulation infected using various strategies.

From Figure 13, it can be seen that the number of infected subpopulations was reduced by applying the third strategy compared to other strategies. Based on strategies 1–3, it can be concluded that strategy 3 is the best strategy to minimize the number of people infected with COVID-19 in the community.

 7. Cost evaluation

The cost evaluation aims to determine the most minimal cost-effectiveness strategy of COVID-19 spread control measures. Cost evaluation in this study uses average cost-effectiveness ratio (ACER) and incremental cost-effectiveness ratio (ICER). According to the approach to cost-effectiveness analysis, ACER is defined mathematically as follows: ACER=Objective function (J)Total number of infections averted. The strategy with the smallest ACER value is the most cost-effective and is obtained in Table 4 as follows:

The incremental cost-effectiveness ratio, which compares two intervention options vying for the same scarce resources, typically tracks costs and health benefits changes. ICER is defined as follows when considering strategies p and q as two competing control intervention techniques: ICER=Change in total costs in strategies p and qChange in control benefits in strategies p and q. Next, ICER was calculated to determine the most cost-effective strategy out of all the control strategies. First, the competition for strategies 1 and 2 is calculated as follows: ICER(1)=49,247,000,000-053,846,000,000-0=0.9146,ICER(2)=2,991,000,000-49,247,000,00053,893,000,000-53,846,000,000=-46,256,000,00047,000,000=-984.1702, The ICER results in strategy 1 were greater than in strategy 2, so educational controls alone were more expensive and ineffective than medical care enhancement controls. Therefore, strategy 1 is removed from the possible control strategies. Next, the ICER for strategies 2 and 3 is recalculated as follows: ICER(2)=2,991,000,000-053,893,000,000-0=0.5130,ICER(3)=1,794,700,000-2,991,000,00053,894,000,000-53,893,000,000=-1,196,300,0001,000,000=-1,196.3. Strategy 2 has a higher ICER value than strategy 3. So, strategy 3 (combined control) is the best control strategy of all options because of its cost-effectiveness and prevention of the spread of infectious diseases.

Total infections prevented, total costs, and ACER for strategies 1, 2, and 3.

Strategy Infections prevented Total cost ACER
No strategy 0 0 0
Strategy 1 53,846,000,000 49,247,000 0.9146
Strategy 2 53,893,000,000 2,991,000 0.0555
Strategy 3 53,894,000,000 1,794,700 0.0333
8. Conclusion

In this study, we have proposed a mathematical model of COVID-19 with comorbidities and added control of community education and improvement of medical care. The proposed model has been calibrated using cumulative confirmed infection cases in Indonesia. The basic reproduction number has been calculated by the next-generation matrix method. The model has an asymptotically stable disease-free equilibrium, provided that the basic reproduction number is <1. Furthermore, the model has an asymptotically stable endemic equilibrium, provided that the basic reproduction number is more than one. Individuals with comorbidity have a greater risk of infection, so there is a need for more supervision and preventive measures such as wearing masks, maintaining distance, and proper sanitation. Public education can be through social media, TV, radio, print media, and others to control the COVID-19 pandemic in Indonesia. Based on the model analysis, it is found that the COVID-19 pandemic can be controlled and eradicated if the value of R0 < 1 by providing public education control and improving medical care. The sensitivity analysis results show that the most influential parameters are quarantine and contact with infected individuals, so educating the public to reduce disease transmission is important. After public education was given, the community became aware of the COVID-19 outbreak and began to reduce contact with other people. Likewise, the Indonesian government imposed large-scale social restrictions (PSBB) and enforced restrictions on community activities (PPKM) with four levels aiming to reduce infection and reduce social contact, educational institutions conducted online classes, webinars, etc. In addition to public education, increased medical care also need to be given to individuals who are already infected so that they recover quickly and that the epidemic is resolved soon. Furthermore, from the numerical results and cost-effectiveness analysis on the optimal control problem, it is found that applying a combination of controls can give the best results compared to a single control. This study can be extended in various ways, including considering the stochastic, time delay, and fractional derivative versions of this model. In addition, providing control variations (such as the presence of vaccination) combines the dynamics of two COVID-19 strains with another comorbidity and considers the COVID-19 vaccination model.

 Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

 Author contributions

MAR and F: conceptualization. MAR: methodology, software, writing—original draft preparation, and visualization. MAR, F, CA, and CWC: validation, investigation, and data curation. F, CA, and CWC: formal analysis. MAR, F, and CA: writing—review and editing. F and CA: supervision. All authors have read and agreed to the published version of the manuscript.

 Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

 References WHO. Novel Coronavirus. (2020). Available online at: https://www.who.int/indonesia/news/novel-coronavirus Wu Z McGoogan JM. Characteristics of and important lessons from the coronavirus disease 2019 (COVID-19) outbreak in China: summary of a report of 72 314 cases from the Chinese center for disease control and prevention. JAMA. (2020) 323:1239–42. 10.1001/jama.2020.264832091533 Kemkes. Pertanyaan dan Jawaban Terkait COVID-19 Kementerian Kesehatan. (2020). Available online at: https://www.kemkes.go.id/article/view/20030400008/FAQ-Coronavirus.html WHO. Novel Coronavirus (COVID-19). (2021). Available online at: https://www.moh.gov.sa/en/HealthAwareness/EducationalContent/Corona/Pages/corona.aspx WHO. Novel Coronavirus (2019-nCoV): Situation Report 1. (2020). Available online at: https://apps.who.int/iris/handle/10665/330760 Wu YC Chen CS Chan YJ. The outbreak of COVID-19: an overview. J Chin Med Assoc. (2020) 83:217–20. 10.1097/JCMA.000000000000027032134861 Guan W Liang W Zhao Y Liang H Chen Z Li Y . Comorbidityity and its impact on 1590 patients with COVID-19 in China: a nationwide analysis. Eur Respir J. (2020) 55:2000547. 10.1183/13993003.00547-202032217650 Yang J Zheng Y Gou X Pu K Chen Z Guo Q. International journal of infectious diseases prevalence of comorbidities and its effects in patients infected with SARS-CoV-2: a systematic review and meta-analysis. Int J Infect Dis. (2020) 94:91–5. 10.1016/j.ijid.2020.03.01732173574 Satgas. Peta Sebaran COVID-19. (2021). Available online at: https://covid19.go.id/peta-sebaran Adisasmito W Suwandono A Trihono Gani A Aisyah DN Solikha DA . Studi komparasi pembelajaran penanganan COVID-19 Indonesia-Korea Selatan. In: Direktorat Kesehatan dan Gizi Masyarakat Kementerian PPN/BAPPENAS. (2021). Available online at: https://infeksiemerging.kemkes.go.id/document/download/3ax61Bxrn5 Worldometers. Report Coronavirus Cases (2021). Available online at: https://www.worldometers.info/coronavirus/ Lemon SM Hamburg MA Frederick P Choffnes ER Mack A. Ethical and Legal Considerations in Mitigating Pandemic Disease: Workshop Summary. Washington, DC: National Academic Press (2007). 21595110 WHO. COVID-19 Strategy Update. (2020). Available online at: https://www.who.int/publications/m/item/covid-19-strategy-update WHO. Pertimbangan-Pertimbangan Untuk Karantina Individu Dalam Konteks Penanggulangan Penyakit Coronavirus (COVID-19). (2020). Available online at: https://www.who.int/docs/default-source/searo/indonesia/covid19/who-2019-covid19-ihr-quarantine-2020-indonesian.pdf?sfvrsn=31d7cbd8_2 Tang B Wang X Li Q Bragazzi NL Tang S Xiao Y . Estimation of the transmission risk of the 2019-nCoV and its implication for public health interventions. J Clin Med. (2020) 9:462. 10.3390/jcm902046232046137 Müller J Kuttler C. Methods and Models in Mathematical Biology. Berlin; Heidelberg: Springer-Verlag (2015). Murray JD. Mathematical Biology I: An Introduction. Berlin; Heidelberg: Springer-Verlag (2002). Feng Z. Final and peak epidemic sizes for SEIR models with quarantine and isolation. Math Biosci Eng. (2007). 4:675–86. 10.3934/mbe.2007.4.67517924718 Tahir M Shah SIA Zaman G Khan T. Stability behaviour of mathematical model MERS corona virus spread in population. Filomat. (2019) 33:3947–60. 10.2298/FIL1912947T Usaini S Hassan AS Garba SM Lubuma JMS. Modeling the transmission dynamics of the middle east respiratory syndrome coronavirus (MERS-CoV) with latent immigrants. J Interdisc Math. (2019) 22:903–30. 10.1080/09720502.2019.1692429 Soewono E. On the analysis of COVID-19 transmission in wuhan, diamond princess and jakarta-cluster. Commun Biomath Sci. (2020) 3:9–18. 10.5614/cbms.2020.3.1.2 Das P Upadhyay RK Misra AK Rihan FA Das P Ghosh D. Mathematical model of COVID-19 with comorbidity and controlling using non-pharmaceutical interventions and vaccination. Nonlinear Dyn. (2021) 106:1213–27. 10.1007/s11071-021-06517-w34031622 Omame A Sene N Nometa I Nwakanma CI Nwafor EU Iheonu NO . Analysis of COVID-19 and comorbidity co-infection model with optimal control. Optimal Control Appl Methods. (2021) 42:1568–90. 10.1002/oca.274834226774 Jia J Ding J Liu S Liao G Li J Duan B . Modeling the control of COVID-19: impact of policy interventions and meteorological factors. Electron J Diff Equ. (2020) 2020:23. Prathumwan D Trachoo K Chaiya I. Mathematical modeling for prediction dynamics of the coronavirus disease 2019 (COVID-19) pandemic, quarantine control measures. Symmetry. (2020) 12:1404. 10.3390/sym12091404 Deressa CT Duressa GF. Modeling and optimal control analysis of transmission dynamics of COVID-19: The case of Ethiopia. Alexandria Eng J. (2020) 60:719–32. 10.1016/j.aej.2020.10.004 Olaniyi S Obabiyi OS Okosun KO Oladipo AT Adewale SO. Mathematical modelling and optimal cost-effective control of COVID-19 transmission dynamics. Eur Phys J Plus. (2020) 135:938. 10.1140/epjp/s13360-020-00954-z33262923 Aldila D Khoshnaw SHA Safitri E Rais YA Bakry ARQ Samiadji BM . A mathematical study on the spread of COVID-19 considering social distancing and rapid assessment: the case of Jakarta, Indonesia. Chaos Solit Fractals. (2020) 139:110042. 10.1016/j.chaos.2020.11004232834600 Aldila D Ndii MZ Samiadji BM. Optimal control on COVID-19 eradication program in Indonesia under the effect of community awareness. Math Biosci Eng. (2020) 17:6355–89. 10.3934/mbe.202033533378859 Ali M Shah STH Imran M Khan A. The role of asymptomatic class, quarantine and isolation in the transmission of COVID-19. J Biol Dyn. (2020) 14:389–408. 10.1080/17513758.2020.177300032498655 Rois MA Tafrikan M Norasia Y Anggriani I Ghani M. SEIHR model on spread of COVID-19 and its simulation. Telematika. (2022) 15. Annas S Pratama MI Rifandi M Sanusi W Side S. Stability analysis and numerical simulation of SEIR model for pandemic COVID-19 spread in Indonesia. Chaos Solit Fractals. (2020) 139:110072. 10.1016/j.chaos.2020.11007232834616 Baba IA Nasidi BA Baleanu D. Optimal control model for the transmission of novel COVID-19. Comput Mater Continua. (2021) 66:3089–106. 10.32604/cmc.2021.012301 Li XP Darassi MH Khan MA Chukwu CW Alshahrani MY Shahrani MA . Assessing the potential impact of COVID-19 Omicron variant: insight through a fractional piecewise model. Results Phys. (2022) 38:105652. 10.1016/j.rinp.2022.10565235663799 Bajiya VP Bugalia S Tripathi JP. Mathematical modeling of COVID-19: impact of non-pharmaceutical interventions in india mathematical modeling of COVID-19: impact of non-pharmaceutical interventions in India. Chaos. (2020) 30:113143. 10.1063/5.002135333261327 Das P Nadim SS Das S Das P. Dynamics of COVID-19 transmission with comorbidity: a data driven modelling based approach. Nonlinear Dyn. (2021) 106:1197–211. 10.1007/s11071-021-06324-333716405 Diagne ML Rwezaura H Tchoumi SY Tchuenche JM. Mathematical model of COVID-19 with vaccination and treatment. Comput Math Methods Med. (2021) 2021:1250129. 10.1155/2021/125012934497662 Ghostine R Gharamti M Hassrouny S Hoteit I. An extended SEIR model with vaccination for forecasting the COVID-19 pandemic in Saudi Arabia using an ensemble Kalman filter. Mathematics. (2021) 9:636. 10.3390/math9060636 Kouidere A Khajji B Bhih AE Balatif O Rachik MA. Mathematical modeling with optimal control strategy of transmission of COVID-19 pandemic virus. Commun Math Biol Neurosci. (2020) 2020:24. 10.28919/cmbn/4599 Madubueze CE Onwubuya IO Dachollom S. Controlling the spread of COVID-19: optimal control analysis. Comput Math Methods Med. (2020) 2020:6862516. 10.1101/2020.06.08.2012539332963585 Naveed M Baleanu D Rafiq M Raza A Soori A Hasan . Dynamical behavior and sensitivity analysis of a delayed coronavirus epidemic model. Comput Mater Continua. (2020) 65:225–41. 10.32604/cmc.2020.011534 Riyapan P Shuaib SE Intarasit A. A mathematical model of COVID-19 pandemic: a case study of Bangkok, Thailand. Comput Math Methods Med. (2021) 2021:6664483. 10.1155/2021/666448333815565 Yang C Wang J. A mathematical model for the novel coronavirus epidemic in Wuhan, China. Math Biosci Eng. (2020) 17:2708–24. 10.3934/mbe.202014832233562 Youssef HM Alghamdi NA Ezzat MA El-bary AA Shawky AM. Study on the SEIQR model and applying the epidemiological rates of COVID-19 epidemic spread in Saudi Arabia. Infect Dis Model. (2021) 6:678–92. 10.1016/j.idm.2021.04.00533898884 Youssef HM Alghamdi NA Ezzat MA El-Bary AA Shawky AM. A new dynamical modeling SEIR with global analysis applied to the real data of spreading COVID-19 in Saudi Arabia. Math Biosci Eng. (2020) 17:7018–44. 10.3934/mbe.202036233378886 Rois MA Trisilowati T Habibah U. Dynamic analysis of COVID-19 model with quarantine and isolation. JTAM. (2021) 5:418–33. Rois MA Trisilowati Habibah U. Optimal control of mathematical model for COVID-19 with quarantine and isolation. Int J Eng Trends Technol. (2021) 69:154–60. 10.14445/22315381/IJETT-V69I6P22332498655 Musa SS Qureshi S Zhao S Yusuf A Mustapha UT He D. Mathematical modeling of COVID-19 epidemic with effect of awareness programs. Infect Dis Model. (2021) 6:448–60. 10.1016/j.idm.2021.01.01233619461 Chukwu CW Fatmawati. Modelling fractional-order dynamics of COVID-19 with environmental transmission and vaccination: a case study of Indonesia. AIMS Math. (2022) 7:4416–38. 10.3934/math.2022246 Zeb A Alzahrani E Erturk VS Zaman G. Mathematical Model for coronavirus disease 2019 (COVID-19) containing isolation class. BioMed Res Int. (2020) 2020:3452402. 10.1155/2020/345240232685469 Bonyah E Juga M Fatmawati. Fractional dynamics of coronavirus with comorbidity via Caputo-Fabrizio derivative. Commun Math Biol Neurosci. (2022) 2022:12. 10.28919/cmbn/6964 Fatmawati Yuliani E Alfiniyah C Juga ML Chukwu CW. On the modeling of COVID-19 transmission dynamics with two strains: insight through caputo fractional derivative. Fractal Fract. (2022) 6:346. 10.3390/fractalfract6070346 Majumder M Tiwari PK Pal S. Impact of saturated treatments on HIV-TB dual epidemic as a consequence of COVID-19: optimal control with awareness and treatment. Nonlinear Dyn. (2022) 109:143–76. 10.1007/s11071-022-07395-635431455 Rai RK Khajanchi S Tiwari PK Venturino E Misra AK. Impact of social media advertisements on the transmission dynamics of COVID-19 pandemic in India. J Appl Math Comput. (2022) 68:19–44. 10.1007/s12190-021-01507-y33679275 Brauer F Castillo-Chavez C. Mathematical Models in Population Biology and Epidemiology. New York, NY: Springer-Verlag New York (2012). Driessche PVD Watmough J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math Biosci. (2002) 180:29–48. 10.1016/S0025-5564(02)00108-612387915 Alligood KT Sauer TD Yorke JA. CHAOS: an Introduction to dynamical systems. In: Introduction To Computational Modeling Using C and Open-Source Tools. New York, NY; Berlin; Heidelberg: Springer-Verlag (2000). Chitnis N Hyman JM Cushing JM. Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model. Bull Math Biol. (2008) 70:1272. 10.1007/s11538-008-9299-018293044 Rois MA Trisilowati T Habibah U. Local sensitivity analysis of COVID-19 epidemic with quarantine and isolation using normalized index. Telematika. (2021) 14:13–24. 10.35671/telematika.v14i1.1191