Optimal control problem arising from COVID-19 transmission model with rapid-test.

The world health organization (WHO) has declared the Coronavirus (COVID-19) a pandemic. In light of this ongoing global issue, different health and safety measure has been recommended by the WHO to ensure the proactive, comprehensive, and coordinated steps to bring back the whole world into a normal situation. This is an infectious disease and can be modeled as a system of non-linear differential equations with reaction rates which consider the rapid-test as the intervention program. Therefore, we have developed the biologically feasible region, i.e., positively invariant for the model and boundedness solution of the system. Our system becomes well-posed mathematically and epidemiologically for sensitive analysis and our analytical result shows an occurrence of a forward bifurcation when the basic reproduction number is equal to unity. Further, the local sensitivities for each model state concerning the model parameters are computed using three different techniques: non-normalizations, half-normalizations, and full normalizations. The numerical approximations have been measured by using System Biology Toolbox (SBedit) with MATLAB, and the model is analyzed graphically. Our result on the sensitivity analysis shows a potential of rapid-test for the eradication program of COVID-19. Therefore, we continue our result by reconstructing our model as an optimal control problem. Our numerical simulation shows a time-dependent rapid test intervention succeeded in suppressing the spread of COVID-19 effectively with a low cost of the intervention. Finally, we forecast three COVID-19 incidence data from China, Italy, and Pakistan. Our result suggests that Italy already shows a decreasing trend of cases, while Pakistan is getting closer to the peak of COVID-19.

Introduction

It is likely that this infectious disease (COVID-19) originated in an animal species, and then spread to humans. Person to person spread of the novel coronavirus reported daily throughout the world. This virus involves serious respiratory tract infections [1], [2]. Therefore, all the countries are making all-out efforts to deal with a rapidly evolving situation which is a challenge for the whole world. An emergency has been declared in infected areas of the world and a serious public health concern has been paid at a global level. Now, it is important to understand how to get control by monitoring the spreading of this disease. Within this urgency, doctors and paramedical staff are on the front line for treating the COVID-19 patients. While to stop the impact of this infection and to avoid further spreading some mathematical estimations are also being performed at each level. Some method has been proposed by the researches based on the mathematical modeling for calculating the basic reproduction number [3], [4].

Various ways have been conducted by the government all over the world to suppress the spread of COVID-19 on their countries, such as with social distancing, international travel restrictions, rapid-test, and even lockdown [5]. Many mathematical models conclude that lockdown is the best way to reduce the spread of COVID-19 effectively among all the aforementioned control strategies [6]. However, lockdown interventions are very risky for a country’s economic stability, Pakistan, India, Iran. Therefore, as a step to prevent the increasing number of infections, social distancing interventions to minimize the successful contact of infections and rapid-test to map the spread of infection into options in various countries [5], [7], instead of implementing lockdown in their countries. Computational results give an essential way to identify the key critical elements based on the modern decomposition techniques [8], [9], [10] in different available reaction routes [11], [12] that allow us to discuss the dynamical properties of the suggested models of the COVID-19. Recently, some models of the COVID-19 have suggested, they provide a good step forward to understand the dynamics of this disease [13], [14], [15], [16]. Accordingly, some suggested mathematical models were reviewed and some computational simulations investigated for the confirmed cases in China [17]. More recently, we developed an updated model of the COVID-19, we have also identified some key critical parameters with sensitivity analysis [18].

Although some mathematical models have been projected so far for new coronavirus disease prediction, a lot can still be improved. Defining such models based on mass action law with reaction rate constants and calculating the sensitivities for each model state with respect to model parameters could improve the outcomes. In a complicated modeling case like new coronavirus dynamics, it is necessary to pay more attention to the optimal control problem and sensitivity analysis more accurately and widely.

Here in this article, we further developed our previous model, some transmission paths and parameters are added. We focused on analyzing the effect of COVID-19 rapid-test as an alternative to suppress the spread of COVID-19. Another novelty of the paper is the identification of the critical model parameters, which makes it easy for the biologists to be used with less knowledge of mathematical modeling and also facilitates the improvement of the model for future development. Consequently, here we measure the effect of rapid-test infection identification on the COVID-19 free equilibrium point and the reproduction number for local stability. Interestingly, the optimal control problem applied to the established model shows that the time-dependent interventions which adapt to the number of infections are able to reduce the number of COVID-19 infections well and at a much lower cost. Finally, we give some short time forecasting of three countries (China, Italy and Pakistan) using our proposed model.

A mathematical model of COVID-19

Let assume the human population can be separated depending on their health status respected to infection status on COVID-19 disease, both visually (symptoms) or through a medical test. Next, let us consider that there is a random test to check whether someone is infected with COVID-19 or not. Then, we split the human class into 5 different classes.

Susceptible class (): This class presents a healthy individual.

Asymptomatic class (): This class presents an infected individual in the early stage of infection. They do not show any symptoms, but capable to infect through droplets or direct contact with the susceptible individual. Because they do not show any symptoms, these individuals are still easy to perform social contact with everyone since they do not realize that they are infected by COVID-19.

Symptomatic Unreported class (): This class presents an individual who gets infected, had symptoms of COVID-19, but did not detect by the government as a COVID-19 suspect.

Symptomatic reported class (): This class presents an individual who gets infected, shows a symptoms of COVID-19, and detected by the government, either it from a rapid test, or from voluntary action to report to the hospital. We assume that all individuals in this class will get a specific treatment and supervision, whether it is through monitored isolation or treatment in the hospital.

Recovered class (): This class present individual who get recovered from COVID-19, and had a temporal immunity.

The transmission diagram which illustrates the interaction between each class described in Fig. 1.

Fig. 1

Transmission diagram of COVID-19 with the effect of rapid-test.

The explanation about model construction is as follows. Susceptible individual increased caused by natural birth rate , and infection from and with effective contact rate and , respectively. Note that since undetected infected individuals still have full access to perform a contact social with randomly. On the other hand, has the smallest infection rate caused by an infected individual already detected by the government, which isolated in the hospital or monitored by the government to conduct self-isolation in their home. Individuals in increased cause by infection from , and decreased caused by recovery to with a rate of , progression in symptomatic to with a rate of and the infectious detected individual with the rate of . present both human awareness to report their health status to the government about their symptoms, so they can get treatment by the government, or detected by rapid test intervention by the government. Please note that since we assume that the government has more concerned to bring the symptomatic individual into the hospital. Undetected symptomatic individual increased by progression from and decreased by recovery with constant rate , rapid test and death rate induced by COVID-19 . Detected symptomatic infectious individual increased by progression from , rapid test from , and decreased by recovery rate and death rate induced by COVID-19. Last, recovered compartment increased by recovery rate from all infected individuals. Each compartment decreased by natural death rate .

Based on the transmission diagram in Fig. 1 and aforementioned explanation, the model equation to describe the effect of rapid test in the spread of COVID-19 is as follows: This model supplemented with the non-negative initial condition, and note that all parameters are positive.

Model analysis

Basic properties

Epidemiological meaningfulness of system (1) is one of the paramount analyzes in this section since it describes the human population. We perform the following theorem about the positiveness solution of the system (1).

Let the initial conditions:exist in the interval, then the solutionsandof system (1) are positive for all .

From system (1), we obtain: The above rates are all non-negative over their boundary planes of the non-negative cone . Therefore, we have the direction of vector fields intended inward from their boundaries. Consequently, we are starting from the non-negative initial conditions so that all the solutions of the system (1) remains positive for all the time . Hence, the following theorem implies the boundedness solution of the system (1).

The biologically feasible region is positively invariant for the model (1) .

Since we have that . Basically, with respect to the condition . Therefore, we have that to be positively invariant and attracting which suffices system (1) can be considered in . Hence, system (1) considered being well-posed mathematically and epidemiologically.

COVID-19 free equilibrium point and the reproduction number

The COVID-19 free equilibrium point of system (1) is given by: To analyze the local stability of , first, we construct the valued basic reproduction number of system (1) using the next-generation matrix approach (Please see [19] for further detail, and more example in [20], [21], [22], [23], [24], [25]). The basic reproduction number of system (1) is given by: where Note that is the summation of three types of “local” basic reproduction numbers depending on the origin of the infection, whether it from asymptomatic , undetected symptomatic or detected symptomatic individuals.

Having the basic reproduction number in hand, and using [29], we have the following theorem regarding the local stability of .

The COVID-19 free equilibriumof the system (1) is locally asymptotically stable if , and unstable if .

presents the expected number of secondary cases of COVID-19 which generated by a single infection introduced into a community of totally susceptible individuals. The results in Theorem 3 shows that COVID-19 can be eliminated in the community when the basic reproduction number is less than unity.

Endemic equilibrium point

The COVID-19 endemic equilibrium point of system (1) is given by where and is taken from the positive root of linear equation given by: where

It can be seen from (8) that will be positive if . The existence of the endemic equilibrium is given in the following theorem.

The COVID-19 endemic equilibrium pointof the system (1) exist if condition holds.

Bifurcation analysis

On the system (1), we assumed Therefore, system (1) can be re-written as Next, let as the bifurcation parameter. To do this, we solve respect to to yield which is given by

Next, substitute and to the Jacobian matrix of system (9) which will gave us: This matrix has one simple zero eigenvalues, while the other four are negative. Therefore, we can use a center-manifold approach to analyze the stability of the endemic equilibrium when close to one.

Firstly, we look for the right eigenvector and left eigenvector. Let vector as the right eigenvector of simple zero eigenvalue of . The right eigenvector is given by Similarly, let the left eigenvector of is given by . Therefore, the left eigenvector is obtained as follows: Since the eigenvector and , so there is no need to look for a partial derivative of and . Therefore, we find the derivatives of , , and to get the values and in the Castillo-Song bifurcation theorem [30]. From non-zero , and , the derivatives are as follows: So that and are obtained as follows: As and , the become unstable when , but close to one. At the same time, it appears an endemic equilibrium, which is locally asymptotically stable. The following results are stated in the form of the following theorem.

System (1) undergoes a forward bifurcation at .

To illustrate the result of Theorems 3, 4, and 5, we give a bifurcation diagram of system 1 presented by detecting symptomatic variables with in Fig. 2. To derive Fig. 2, we use parameter value as shown in Table 1 except which used as the bifurcation parameter. When , then we have that , which gives a stable COVID-19 free equilibrium. When , zero eigenvalues appear, and change of stability appears when we can see that the COVID-19 free equilibrium becomes unstable, while the endemic equilibrium arises and locally stable.

Fig. 2

Transcritical bifurcation diagram of system (1) using as the bifurcation parameter.

Table 1

The model parameters and initial model populations for COVID-19 epidemic outbreak with their biological definitions.

Symbols	Biological definitions	Estimated values	Sources
S(0)	Initial susceptible individuals	11.081×10^6	[26]
A(0)	Initial asymptomatic infected individuals	3.62	[26]
I(0)	Initial recorded symptomatic infected individuals	1	[26]
U(0)	Initial unrecorded symptomatic infected individuals	4.13	[26]
R(0)	Initial recovered individuals	0	[26]
β_1	Transition rate between susceptible and asymptomatic infected individuals	4.5×10^−7	Fixed
β_2	Transition rate between susceptible and unreported symptomatic infected individuals	1.45×10^−6	Fixed
β_3	Transition rate between susceptible and reported symptomatic infected individuals	8.68×10^−8	[27]
δ	Transition rate between asymptomatic infected and unreported symptomatic infected	0.0285	[26]
γ_1	Transition rate between asymptomatic infected and reported symptomatic infected	0.1142	[26]
γ_2	Transition rate between unreported symptomatic and reported symptomatic infected individuals	0.35	Fixed
η_1	The recovery rate of unreported asymptomatic infected case	0.13978	[28]
η_2	The recovery rate of unreported symptomatic infected case	0.2	Fixed
η_3	The recovery rate of reported symptomatic infected case	0.33029	[28]
ξ	The reported and unreported symptomatic death rate	1.7826×10^−5	[28]

Discussion on

Based on Theorem 3, Theorem 4, the basic reproduction number in (3) determines the existence or disappearance of COVID-19. This formed by three-component of a specific infection. The first component is , which described a number of new infection from direct contact between a susceptible individual with the asymptomatic individual during its infection period. To reduce this local basic reproduction number, various ways can be adopted, i.e., such as reducing to avoid infection with the asymptomatic individual. Unfortunately, the asymptomatic individual does not show any symptoms; therefore, avoid contact with this type of infected individual is difficult to be estimated. Therefore, the only way is to reduce random contact between each individual, whether it is a healthy or infected individual. In several countries [31], the government is campaigning for the use of medical masks to all people, regardless of whether they are infected or healthy individuals. The other way to reduce is to enlarging the rate of rapid-test . As mentioned before, the purpose of this rapid-test intervention is to map the infected individuals (asymptomatic and symptomatic), so several actions can be done to these individuals to prevent further infection. Enlarging will help the government to be able to focus all health measures on this individual group, such as monitoring self-quarantine, isolation in hospitals, and so forth. The long-term effect is that the government can reduce the number of infections in the field from this group of individuals. Another way to reduce is by increasing the recovery rate . Since these individuals do not show any symptoms, no medical intervention cannot be given to asymptomatic individuals. Therefore, encourage a healthy lifestyle to enhance self immunity through a media campaign is a reasonable option.

The second component in is , which describes new infection from the unreported symptomatic cases during its infection period. Since , reducing the progression rate will reduce . In COVID-19, the incubation period estimated between 2–14 days [32]. Smaller means that the infection needs a longer time to show the symptoms. The other way except increasing rapid-test rate , to reduce can be done by improving the quality and quantity of health services in hospitals, which our model represents, i.e., by reducing the death race induced by COVID-19 and increasing recovery rate . By increasing the capacity of the hospital, more infected individuals can get proper medical treatment, which is expected to reduce the rate of death due to disease, and shorten the duration of infection.

The third component in is which describes a new infection caused by symptomatic reported cases during its infection period. Here, note that this class is monitored by the government, whether it in the hospital, or monitored by the government tracking system to conduct self-quarantine in their home. This intervention will reduce . Discipline implementation of this contact reduction will reduce significantly, and more easy to be monitored compared with minimizing and .

To finalize our discussion on , we perform a sensitivity analysis of respect to and . Using parameters value in Table 1 except and , we plot a condition of on plane in Fig. 3. Since and are negative, the 2nd and rd(a) area represents a combination of and and gives , while 1st and rd(b) area representing a condition when . The first information that we can take from Fig. 3 is that is more sensitive in determining rather than . It means that rapid-test intervention into asymptomatic individuals is more advisable to encourage for the implementation in the field. The reason is this, by bringing/identify the infected individuals in the field will give the government flexibility and more focus in the intervention. Another information that can be taken from Fig. 3 is that the existence of “useless” intervention, that is the 1st area. The 1st area () representing a condition of is always larger than one, no matter the value of . On the other hand, the 2nd area (), represents an area when is always less than one even though the government does not take a rapid-test in unreported symptomatic individual. If the government only implements in the area between 0.0173 and 0.228, then implementation of must be taken carefully, since it can end up in an unsuccessful intervention (rd(b) area) or the successful intervention (rd(a) area). Therefore, implementing rapid-test need a very careful justification for the implementation, since a large rapid-test intervention leads to a very costly intervention, but small rapid-test could end up in an unsuccessful intervention.

Fig. 3

A condition of respect to and .

Sensitivity analysis

Consider a system of differential equations where and . The functions are often non-linear therefore model differential equations may not solve analytically. An important technique to analyze system (13) is the idea of sensitivity analysis. According to this approach, the sensitivity of each variable concerning parameters can be calculated. The main equation of sensitivity is given below The first order derivatives given in Eq. (14) represent the time-dependent sensitivities of all variables with respect to each parameter value . Furthermore, the differential equations can be solved for sensitivity coefficients as below Using the chain rule of differentiation, Eq. (15) can be further driven and the sensitivity equations take the Jacobian matrix as follows where the matrices and are defined by For more details and applications of sensitivity analysis in the field of systems biology, the readers are refereed to [33], [34], [35], [36], [37], [38], [39], [40]. The local sensitivity values are given in Eq. (16) can be computed using SimBiology Toolbox in MATLAB with three different techniques: non-normalizations, half normalizations, and full normalizations. Accordingly, in a complicated modeling case like new coronavirus dynamics, it is necessary to pay attention to sensitivity analysis more accurately and widely. This helps us to identify the key critical model parameters and to improve model dynamics.

Computational results

Mathematical models and computational simulations help and provide a good environment to analyze high dimensional models of infectious diseases. Such models may help biologists to predict future model dynamics and identify critical model parameters. The suggested mathematical models of COVID-19 are effective tools that give estimations and suggestions about controlling the virus and further preventions more effectively and widely. The values of parameters and initial populations in this study are obtained from the WHO situation reports (the National Health Commission of the Republic of China) presented in [26], [27], [28].

There are some numerical approximate solutions of the model equations (1) for different parameters and initial populations using System Biology Toolbox(SBedit) for MATLAB; see Fig. 4, Fig. 5, Fig. 6, Fig. 7, Fig. 8, Fig. 9, Fig. 10, Fig. 11. Accordingly, different model dynamics for initial populations are obtained based on changing the value of model parameters. Results in this study provide a good step forward in predicting the model dynamics in the future for development programs, interventions and health care strategies.

Fig. 4

Computational simulations for the model states given in system (1) of the COVID-19 using MATLAB; there are model dynamics of (a) susceptible individuals, (b) asymptomatic infected individuals, (c) unreported symptomatic infected individuals, (d) reported symptomatic infected individuals, (e) recovered individuals.

Fig. 5

Computational simulations for the model states given in system (1) of the COVID-19 using MATLAB; there is the relationship between the asymptomatic infected people and (a) susceptible individuals, (b) unreported symptomatic infected individuals, (c) reported symptomatic infected individuals, (d) recovered individuals.

Fig. 6

The effect of transition rate on (a) susceptible individuals (b) asymptomatic infected individuals, (c) unreported symptomatic infected individuals, (d) reported symptomatic infected individuals (e) recovered individuals, in computational simulations using MATLAB parameters used for blue lines, for green lines and for red lines. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 7

The effect of transition rate on (a) asymptomatic infected individuals, (b) unreported symptomatic infected individuals, (c) reported symptomatic infected individuals (d) recovered individuals, in computational simulations using MATLAB parameters used for blue lines, for green lines and for red lines. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 8

The effect of transition rate on (a) asymptomatic infected individuals, (b) unreported symptomatic infected individuals, (c) reported symptomatic infected individuals (d) recovered individuals, in computational simulations using MATLAB parameters used for blue lines, for green lines and for red lines. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 9

The effect of transition rate on (a) unreported symptomatic infected individuals, parameter used (blue line), (green line), (red line), (b) unreported symptomatic infected individuals, parameter used (red line), (green line), (blue line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 10

The effect of transition rate on (a) reported symptomatic infected individuals (b) recovered individuals, in computational simulations using MATLAB parameters used for blue lines, for green lines and for red lines. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 11

The effect of transition rate and on unreported symptomatic infected individuals in computational simulations using MATLAB parameters used (a) for blue lines, for green lines and for red lines, (b) for blue lines, for green lines and for red lines. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

To perform our numerical experiments in this section, due to a short time interval of simulation, we ignore new-born and natural death rate in our model. Therefore, we have that and . With this assumption, our model now read as:

The model dynamics of susceptible, asymptomatic infected, reported symptomatic infected, unreported symptomatic infected and recovered individuals are shown in Fig. 4. The number of susceptible individuals decreases dramatically and becomes stable after four days while the dynamics of recovered people increase gradually and get flat after 15 days. Interestingly, the number of asymptomatic infected individuals reaches a high level after 5 days while the number of infected people in both reported and unreported symptomatic are dramatically changed between 3 days to 15 days. Furthermore, Fig. 5 explains the relationship between asymptomatic infected people with the other groups in the COVID-19. There are almost the same model dynamics for reported and unreported symptomatic states whereas there are slightly different model dynamics for susceptible and recovered groups.

Fig. 6 shows that the impact of the transition rate on all model variables. The effect of this parameter can easily occur on in the dynamics of the model states. For example, if the value of is increased then the number of asymptomatic, unreported symptomatic and reported symptomatic infected people are also increased, see Fig. 6(b,c,d). Also, the dynamics of susceptible and recovered people become stable when the value of this parameter becomes larger and larger, see Fig. 6(a) and (e).

Fig. 7 shows that the impact of transition rate on asymptomatic infected people, reported symptomatic infected, unreported symptomatic infected, and recovered people. The effect of this parameter can easily occur on the variables and . It can be seen that the model dynamics for such states become flattered when the value of is increased, see Fig. 7(a,b,d). On the other hand, the number of reported symptomatic people becomes larger when the value of this parameter is increased, see Fig. 7(c). This is an important key element for controlling this disease.

The impact of transition rate on asymptomatic infected, reported symptomatic infected, unreported symptomatic infected, and recovered people is shown in Fig. 8. Intestinally, the number of infected people in groups is dramatically increased when the value of gets smaller, this is illustrated in Fig. 8(a,b,c). Furthermore, the dynamic of recovered individuals reaches stable very quickly when this parameter becomes large, see Fig. 8(d).

Fig. 9 shows that the impact of parameter on the dynamics of unreported symptomatic infected people. The number of infected people is significantly changed for different values of . The effect of transition rate on reported symptomatic infected individuals and recovered individuals are computed and given in Fig. 10. Moreover, the effect of transition rate and on unreported symptomatic infected individuals are shown in Fig. 11.

The idea of sensitivity analysis has an important role in identifying the model critical element. The main equation of local sensitivity is presented in Eq. (16). We use SimBiology Toolbox for MATLAB to calculate the local sensitivity of each model state concerning model parameters for the model equations (1). We compute the model sensitivities using three different techniques: non-normalizations, half normalizations and full normalizations; see Fig. 12. Interestingly, results provide us further understanding of the model and helps us to identify the key critical model parameters. For example, it generally seems that the susceptible, asymptomatic infected, recovered individuals are more sensitive to the set of model parameters compared to the reported and unreported symptomatic individuals, this result is based on non-normalization approach, see Fig. 12(a). Another interesting result is that the group of asymptomatic, reported symptomatic, unreported symptomatic people are very sensitive to almost all model parameters, see Fig. 12(b). This is an effective step to identify the model critical parameters for controlling the spread of COVID-19. Accordingly, the parameters are the key critical elements for understanding and to prevent this disease because they are very sensitive according to full normalization method presented in Fig. 12(c). As a result, identifying critical model parameters in this study based on computational simulations is an effective way to further study the model practically and theoretically and give some suggestions for future improvements of the novel coronavirus vaccination programs, interventions and controlling the spread of disease.

Fig. 12

The sensitivity of each model state with respect to model parameters in computational simulations for the coronavirus disease (COVID-19); (a) non-normalization sensitivity, (b) half normalization sensitivity, (c) full normalization sensitivity.

Optimal control problem

Characterization of the optimal control problem

In this section, we analyze the optimal control problem related to the model (1). This optimal control approach aims to minimize the number of an infected individual using rapid-test intervention. As explained in Section “A Mathematical Model of COVID-19” that with the rapid-test, policymakers could map and detect the infected individual; hence, the controlled isolation could be implemented to these individuals to avoid contact with a susceptible individual. Therefore, we re-define and in system (1) as a time-dependent parameter , respectively. Due to a short time of interval simulation and short term of COVID-19 pandemic, we neglect natural new-born and the natural death rate from our model, similarly with Section “Computational Results”. Also, since does not appear in another part of the equation in (17) except in , we may ignore from our optimal control problem. Therefore, system (17) now read as: Mathematically, our goal is to find an optimal rate of rapid test which minimize the following cost functional: The first three-component in the integrand is a cost related to a high number of COVID-19 infections in the community, while the last to component related to rapid-test intervention. Note that , and are the weight constant that will balance each component in this cost function.

Here we seek an optimal solution such that where be the set of admissible control. To investigate the existence of optimal control, we use Pontryagin’s Maximum Principle to govern the necessary condition of the optimal control problem. The Lagrangian for the optimal system (18) can be defined as: where for are the costate variable related to , respectively. The costate variable satisfies the following system of ordinary differential equations: with the transversality condition for . Using the optimality condition , we get: Taking into account the lower and upper bound for and , we get the characterization:

Numerical experiment on the optimal control problem

The optimal control system which involves the system of state variables in (18), costate variables in (21), and optimal characterization in (22) is analyzed using the Runge–Kutta forward–backward iterative numerical approximation method [21], [23], [41]. The idea of this method as follows. First, give an initial guess for control variables and for all time . Using this value, and the initial condition for state variables, solve system (18) forward in time to find values of state variables in all-time , and calculate the cost function (19). Next, solve the costate system (21) backward in time using the transversality condition of costate variables, an initial guess of and solution of state variables from the previous step. Next, update the control variables using Eq. (22). Repeat this scheme until the convergence criteria achieved. In this article, the terminate condition is until the error of the optimal solution in iteration- is less than small constant , or in this case:

For the base-case, we choose the weight cost: Note that and are larger than because it is quite difficult for the policymaker to handle the undetected crises. Furthermore, we have that since the rapid test for the asymptomatic individual is easier to implement if the candidate had already shown the symptoms. We run our simulation for .

For the base case, we use parameter value as shown in Table 1, and the numerical results are given in Fig. 13, and the final number of infected individuals at , total cost function (19), and a number of an averted infected individual given in Table 2. Please note that we use two different y-axes to identify the number of the infected individual, with and without control, since the scale of an infected individual without control is almost ten times larger than with controls. It can be seen that the time-dependent control succeeds in suppressing the number of infected individuals in all classes. The Control profile for and is relatively different. It can be seen that rapid test for the symptomatic individual should be given at a maximum rate since the early time of simulation, and start to decrease to it minimum value after day 70, where the total of infected individuals are already decreasing. On the other hand, rapid-test for asymptomatic individual start to be implemented at day 35 when a total of individual infected start to increase, and it starts to decrease when the number of asymptomatic individuals also shows a decreasing trend. It is interesting to see that the dynamic of the symptomatic individual has three outbreaks. The first outbreak occurs as a consequence of start to increase in day 35, while the second outbreak occurs when starts to decrease to its minimum value in day 62. When starts to decrease in day 70, the infected population will start to increase, and the symptomatic class will reach the third outbreak on day 90. Using the control profile, which depends on time, as shown in Fig. 13, we can avoid new cases as much as 138 thousand cases with the cost of intervention more than 50% cheaper than without the intervention.

Fig. 13

Numerical simulation of the optimal control problem of the system (18) for the base-case.

Table 2

Numerical result for the base-case in Fig. 13.

Case	J	New infection averted
No control	2.43×10^6	–
Time-dependent control	4.87×10^5	138000

Our second simulation conducted to see the impact of a cheaper rapid-test cost on the dynamic of the infected population. To do this, we redefine the weight parameter of the base case as , while the other parameters remain the same. It can be seen that and for this scenario are ten times smaller than in the base-case. The dynamic of the infected population and the control profile can be seen in Fig. 14, and the numerical results are shown in Table 3. Similar to the base-case scenario, the time-dependent control could avoid more than one hundred thousand new cases. As a consequence of a cheaper rapid-test cost, it can be seen that and can remain at the maximum rate for a longer period than in the base-case. Furthermore, we can see that instead of having three outbreak as in the base-case, the symptomatic undetected case only have two outbreak in this scenario. Therefore, we can conclude that a massive implementation of rapid-test as a consequence of a cheaper cost for the implementation could prevent a future outbreak of COVID-19.

Fig. 14

Numerical simulation of the optimal control problem of system (18) when and is smaller ten times than in the base-case.

Table 3

Numerical result for Fig. 14.

Case	J	New infection averted
No control	2.45×10^6	–
Time-dependent control	5.5×10^5	157740

Application of the proposed model to forecast some COVID-19 incidence data

In this section, we present some examples of how our model in (17) can fit COVID-19 incidence data. The incidence data used in this article can be accessed in [42]. We fit the daily infected data to compartment  in the proposed model (17) for the early outbreak period. We use parameter value as shown in Table 1, except , and which we estimated from the real data. We use COVID-19 incidence data from China (Date of 22 January 2020–27 May 2020), Italy (Date of 6 March 2020–27 May 2020) and Pakistan (6 March 2020–27 May 2020). To fit these real incidence data, we use the software MATLAB. The results are shown in Fig. 15, and the parameters that been estimated can be seen in Table 4. It can be seen that our model could fit the Incidence data in China. For data of Italy, our model suggests that the disease In Italy shows a decreasing trend and will tend to zero cases approximately after day 150 of simulation (Late of July 2020). On the other hand, our numerical results for Pakistan data show that our proposed model can capture the dynamics of COVID-19 in the early period of the outbreak.

Fig. 15

Confirmed COVID-19 data (dotted) vs. simulation of model (17) in (a) China, (b) Italy and (c) Pakistan.

Table 4

Numerical result for Fig. 15.

Country	β_1	β_2	β_3
China	4.015×10^−8	7.4827×10^−8	8.04×10^−9
Italy	3.768×10^−8	5.4767×10^−8	7.536×10^−9
Italy	1.817×10^−7	3.634×10^−7	3.5×10^−8

Conclusions

Let us summarized the discussion started with the mathematical modeling and concluded with the prime results, while all the data for studies are obtained from the WHO situation reports NHCRC.

Here we have modeled the dynamics of all possible cases of human to human transmission, i.e., susceptible, asymptomatic infected, reported symptomatic infected, unreported symptomatic infected, and recovered individuals to analyze accurate transmission dynamics of the COVID-19 outbreak.

The solutions of the model equations for different parameters and initial populations have been numerically approximated using System Biology Toolbox (SBedit) for MATLAB.

The modeling and simulation based on the suggested sensitivity analysis indicate that almost all model parameters may have a role in spreading this virus among susceptible, infected, recorded symptomatic, unrecorded symptomatic, and recovered individuals.

The effect of the control strategies on the model is analyzed graphically and analytically. From the analysis of the basic reproduction number, we found that rapid-test intervention, which aims to detect infected individuals among the human population, is promising to suppress the spread of the COVID-19. The success level of rapid-test also depending on how the government follows up the cases that have been detected with the help of rapid-test, for instance, with self-quarantine monitored for the infected individual if the case is not yet serious. Increasing the quality of the hospital by upgrading the capacity of the hospital or increase the number of medical staff can also increase the chance of COVID-19 eradication program.

The prospect of rapid test intervention to help the eradication program of COVID-19 analyzed using the optimal control theory. We find that a time-dependent rapid test intervention could reduce the number of new COVID-19 infection at a lower cost. We also find that rapid tests could reduce the size of future outbreak, delay the time of it appearance, furthermore, it also could eliminate the possibility of outbreak occurrence if the implementation of rapid-test is set to be adapted to the increasing number of infections.

The coronavirus, which is the cause of COVID-19, is very easy to mutate. Based on , until 2022, the virus variants for COVID-19 are divided into three, namely variants of concern (4 serotypes), Variants of Interest (3 serotypes), Variants under monitoring (7 serotypes), and De-escalated variants (27 serotypes). These serotype differences are shown in the virus’s physiology and its consequences on people infected with COVID-19, such as the speed it spreads, the response to the environment, and how dangerous the serotype is to humans. Therefore, further analysis is needed regarding these serotype differences. Differences in vaccine efficacy against different serotypes make it difficult to predict the dynamics of COVID-19 with simple models. Multiple variants model can be considered to understand this issue better.

Many types of vaccines have been introduced in various parts of the world. Not all communities in various countries accept this vaccination policy. In addition, several countries have not yet achieved the high vaccination coverage recommended by WHO. Based on this, the mathematical model in this paper needs to be developed by considering the complexity of this vaccination problem, such as differences in vaccine efficacy, multiple phases and vaccine booster, and others.

Results in this study provide a good step forward in predicting the model dynamics in the future for development programs, interventions, and health care strategies. For further development of the model, the reader could consider the existence of the maximum capacity of the hospital since the low capacity of the hospital will prolong the infection period. An application of optimal control problems to model the rate of rapid-test intervention as the time-dependent variable could be considered to handle the budget limitation on the COVID-19 eradication program.

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.