Optimal Control Studies on Age Structured Modeling of COVID-19 in Presence of Saturated Medical Treatment of Holling Type III.

COVID-19 pandemic has caused the most severe health problems to adults over 60 years of age, with particularly fatal consequences for those over 80. In this case, age-structured mathematical modeling could be useful to determine the spread of the disease and to develop a better control strategy for different age groups. In this study, we first propose an age-structured model considering two different age groups, the first group with population age below 30 years and the second with population age above 30 years, and discuss the stability of the equilibrium points and the sensitivity of the model parameters. In the second part of the study, we propose an optimal control problem to understand the age-specific role of treatment in controlling the spread of COVID -19 infection. From the stability analysis of the equilibrium points, it was found that the infection-free equilibrium point remains locally asymptotically stable when R 0 < 1 , and when R 0 is greater than one, the infected equilibrium point remains locally asymptotically stable. The results of the optimal control study show that infection decreases with the implementation of an optimal treatment strategy, and that a combined treatment strategy considering treatment for both age groups is effective in keeping cumulative infection low in severe epidemics. Cumulative infection was found to increase with increasing saturation in medical treatment. © Foundation for Scientific Research and Technological Innovation 2022.

Introduction

Mathematical modeling of infectious diseases such as COVID-19, influenza, dengue, HIV/AIDS etc. are one of the most important research areas today. Mathematical epidemiology has contributed to a better understanding of the dynamical behavior of these infectious diseases, its impacts, and possible future predictions about its spreading. Mathematical models are used in comparing, planning, implementing, evaluating, and optimizing various detection, prevention, therapy, and control programs. COVID-19 is one such contagious respiratory and vascular disease that has shaken the world today. It is caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). On 30 january it was declared as a Public Health Emergency of International Concern. As of latest world statistics of COVID-19, around 219 million cases has been reported and around 4.55 million have died worldwide. Several mathematical models has been developed to understand the dynamics of the disease. A SAIU (susceptible–asymptomatic–reported symptomatic infectious–unreported symptomatic infectious) compartmental model that explains the transmission dynamics of COVID-19 is developed in [47]. Recent works on COVID-19 includes [1–4, 48]. The role of some of the control policies such as treatment, quarantine, isolation, screening, etc. are also applied to control the spread of infectious diseases [7, 26, 38]. COVID-19 has caused the most severe health issues for adults over the age of 60 with particularly fatal results for those 80 years and older. This is due to the number of underlying health conditions present in older population [11]. A mathematical model for estimating the age-specific transmissibility of a novel coronavirus is developed in [56]. In this study, the age specific SEIARW model was fitted with the reported data by dividing the population into four age groups and the results from this study suggested that the highest transmissibility occurred from age group (1–14) to (15–44).

This pandemic has catalysed the development of novel coronavirus vaccines across pharmaceutical companies and research organizations. Many potential vaccines for COVID-19 are being studied and some are under clinical trials. Drugs such as remdesivir, favipiravir, ivermectin, lopinavir/ritonavir, mRNA-1273, phase I trial (NCT04280224) and AVT technology are being used as therapeutic agents by different countries for treating COVID-19 [5, 30, 44, 46, 51, 53]. Therefore, in situation where the vaccine is still unavailable and the disease is spreading exponentially, an effective control measures has become the need of an hour to at least contain the size of the disease.

In this context age related transmissibility of COVID-19 has become a public health concern. It has caused the most severe health issues to adults over 60 years with particularly fatal results for those above 80 years. This is due to the number of underlying health conditions present in older populations [11]. A study with 425 patients indicated that the median age was 59 years (range 15–89) and reported different case distributions in four age groups: 0–14, 15–44 , 45–64 , and > 65 years [8]. However, there is no enough epidemiological evidence to classify the age groups in transmission. Age structure modelling studies for diseases such as influenza and Dengue can be found in [33, 45].

It has been noted by many authors that the saturation in medical treatment leads to complex behavior of the diseases such as existence of multiple endemic equilibria, backward bifurcation, Hopf bifurcation, etc. [12, 29, 41, 55]. In [55] backward bifurcation and global dynamics of an epidemic model with saturated incidence rate and saturated treatment function is studied. Because of these rich dynamics shown by models that consider saturation in medical treatment, we consider type III recovery rate due to treatment to incorporate saturation in medical treatment. Motivated by the above, in this study we propose a nonlinear age-structured compartmental model and first study the stability of the equilibrium points and the sensitivity of the model parameters. Later, we formulate an optimal control problem to understand the role of treatments by considering them as controls in the age-structured model.

The paper is organised as follows: In the next section, we discuss the formulation of the age structured model, followed by the positivity and boundedness of the solutions. In the next section, we compute the equilibrium points and the basic reproduction number, after which the stability of the equilibrium points are analysed. In the following section, we perform the sensitivity analysis of the parameters of the model. In the next section, we formulate an optimal control problem and perform numerical simulations to see what role the controls play in reducing infection. Finally, we present the conclusions.

Model Formulation

We divide the entire population into two age groups, first with age below 30 years (we call this group as young population) and second with age years (we call this group as adult population). In each of these groups we have three compartments namely susceptible denoted by , Infected denoted by and Recovered denoted by for . At any time an individual can be in one of the states. Susceptible are the ones who are not infected but can get infected when it come in contact with the infected ones. Infected are those who spread the infection to others and ones an infected individual recovers from the disease then he develops immunity against the disease and is said to be recovered. The parameters and are the rates at which susceptible in first group are subjected to infection. Here, susceptible get infected by infected population within the same group or from the other group. Similarly, denote the rates at which susceptible in second group are infected. We assume that there is a constant birth rate of the young susceptible population and denote it by . The recovered population in each group can become susceptible again with rate and . Transition rate or maturation rate describing transition from first age group to the second is denoted by parameters m. We assume that the natural mortality rate is common across all the compartments and denote it by parameter , whereas the disease induced death rate and recovery rate differs among age groups. The disease induced death rate and recovery rate for the first age group is denoted by and and similarly, and denote the same for the second age group. To incorporate the limitations in medical resources, we consider type III recovery rate function due to treatment, wherein, the treatment provided initially is less, owing to less infected individuals, and as the epidemic progresses, the treatment increases accordingly and non-linearly saturates due to limitations on medical facilities. Considering the fact that the infection rate and hospitalization rate for COVID-19 is comparatively higher in case of adult population [56], we consider type III recovery rate function only for the second age group. For the first age group, recovery due to treatment is considered as a linear function. Considering all the above the model describing COVID-19 dynamics among two age groups is given by the following system of ordinary differential equations:With and for . The initial population is considered to be positive. The rate describes the type III recovery rate function. All the parameters considered in the model are positive and their biological meaning is given in Table 1.

Table 1

Parameters and their meanings

Symbols	Biological meaning
S_1	Susceptible young population
S_2	Susceptible adult population
I_1	Infected young population
I_2	Infected adult population
R_1	Recovered young population
R_2	Recovered adult population
b_1	Constant birth rate of young population
δ_1	Rate of recovered young becoming susceptible
β_1,β_2	Rate at which susceptible young population are infected because of infected young ones and infected adults
μ	Natural death rate
d_1,d_2	Disease induced death rate of young and adult population
μ_11	Recovery rate due to treatment of young population
δ_2	Rate of recovered adult becoming susceptible again
β_3,β_4	Rate at which susceptible adult population are infected because of infected young ones and infected adults
μ_12,α	μ_12 is the recovery rate due to treatment of infected adult population and α is the non-negative constant related with saturation in treatment
m	Maturation rate

Positivity, Boundedness and Existence of Solutions

Positivity

We now show that if the initial conditions of the system (2.1)–(2.6) are positive, then the solution remains positive for any future time. Positivity of the solutions are established in similar lines to the method discussed in [34, 41]. Using the Eqs. (2.1)–(2.6), we get,Thus all the above rates are non-negative on the bounding planes (given by , , , , , and ) of the non-negative region of the real space. So, if a solution begins in the interior of this region, it will remain inside it throughout time t. This happens because the direction of the vector field is always in the inward direction on the bounding planes as indicated by the above inequalities. Hence, we conclude that all the solutions of the the system (2.1)–(2.6) remain positive for any time provided that the initial conditions are positive. This establishes the positivity of the solutions of the system (2.1)–(2.6). Next we will show that the solution is bounded.

Boundedness

Let .

Now,Here the integrating factor is Therefore, after integration we get, , where, c is a constant. Now as we get,Thus we have shown that the system (2.1)–(2.6) is positive and bounded. Therefore the biologically feasible region is given by the following set,

Existence and Uniqueness of Solution

For the general first order ODE of the formOne would have interest in knowing the answers to the following questions: We use the following theorem to established the existence and uniqueness of solution of the system (2.1)–(2.6).

Under what conditions solution exists for the (3.1)?

Under what conditions unique solution exists for the system (3.1)?

Theorem 3.1

Let D denote the domain:and suppose that f(t, x) satisfies the Lipschitz condition:and whenever the pairs and belong to the domain D , where k is used to represent a positive constant. Then, there exist a constant such that a unique (exactly one) continuous vector solution x(t) exists for the system (3.1) in the interval . It is important to note that condition (3.2) is satisfied by requirement that:be continuous and bounded in the domain D.

We use boundedness of the solutions proved above and show that a unique solution exists for system (2.1)–(2.6) by showing partial derivative of right hand side of Eqs. (2.1)–(2.6) are continuous and bounded with respect to each of the variables and .

LetFrom Eq. (3.3) we haveSimilarly, from Eq. (3.4) we have,From (3.5) we have,From (3.6) we have,From (3.7) we have,Finally, from (3.8) we have,Hence, we have shown that the partial derivatives of are continuous and bounded. Hence, by Theorem 2.1, there exists a unique solution of system (2.1)–(2.6).

Equilibrium points and Basic Reproduction Number

System (2.1)–(2.6) admits two equilibrium namely the disease free equilibrium and the infected equilibrium. The disease free equilibrium denoted by is calculated to be, where,The infected equilibrium is denoted by and is given by, where,Here,

Calculation of

The basic reproduction number is one of the most important quantities in disease modeling. It is defined as the average number of secondary cases produced per primary case. In epidemic models, the basic reproduction number is calculated using the next generation matrix method [25]. The system (2.1)–(2.6) has two infected states and four non-infected states. Computing the Jacobian matrix of the infected states of the system (2.1)–(2.6) at infection-free equilibrium we have,or,where,Calculating the inverse of V we get,Nowwhere,The characteristics equation of is calculated as,The spectral radius of is given by,where,Therefore the basic reproduction number is given by,In the following we numerically show that the infection free equilibrium point remains locally asymptotically stable when and the infected equilibrium point becomes locally asymptotically stable whenever .

Numerical Simulation

Parameters Values

In Table 2 we give the parameter values and the source from which the parameter values are taken. With these parameter values, we will numerically show the asymptotic stability of and based on the values of basic reproduction number in similar lines to [33].

Table 2

Parameters values and their source

Parameters	Values	Source
b_1	μN(0)	[47]
δ_1	0.0714	[33]
β_1,β_2	0.0028	[50]
μ	0.062	[47]
d_1,d_2	0.013, 0.014	[50]
μ_11	0.071	[50]
δ_2	0.0714	[33]
β_3,β_4	0.0086	[50]
μ_12	0.071	[50]
α	0.4	[34]
m	0.000182	[33]

Stability of Infection Free Equilibrium

We show numerically that whenever , the infection-free equilibrium is locally asymptotically stable. To make the value of less than one, we change some of the parameter values of Table 2. When the value of and were taken to be 0.00028 and 0.62 respectively, the value of was calculated as 0.38 and . Taking parameter values from Table 3, the system of Eqs. (2.1)–(2.6) was numerically solved in math-lab. The solutions of the system (2.1)–(2.6) are depicted in Fig. 1 with the initial values . From Fig. 1, we see that the solutions eventually reaches to the infection free state, . Therefore, from this numerical study, we conclude that the infection-free equilibrium of the system (2.1)–(2.6) remains locally asymptotically stable when the value of is less than unity. The parameter values are given in Table 3.

Table 3

Parameters values for

Parameters	Values
b_1	310
δ_1	0.0714
β_1,β_2	0.00028
μ	0.62
d_1,d_2	0.013, 0.014
μ_11	0.071
δ_2	0.0714
β_3,β_4	0.0086
μ_12	0.071
α	0.4
m	0.0182

Fig. 1

Figure depicting local asymptotic stability of whenever

Stability of Infected Equilibrium

Here we numerically show that whenever the value of basic reproduction number denoted by crosses unity we have the infected equilibrium to be locally asymptotically stable. The value of was calculated to be 1.1514 and for the parameter values from Table 4 . From Fig. 2 we see that the infected equilibrium is locally asymptotically stable provided . The initial values were taken as .

Table 4

Parameters values for

Parameters	Values
b_1	210.84
δ_1	0.0714
β_1,β_2	0.00028
μ	0.62
d_1,d_2	0.013, 0.014
μ_11	0.071
δ_2	0.0714
β_3,β_4	0.0086
μ_12	0.071
α	0.4
m	0.182

Fig. 2

Figure depicting local asymptotic stability of whenever

If the parameters are chosen very small in magnitude then some of the matlab based numerical solution may not converge to the exact solution. For alternate numerical analysis and convergence one can look into the methods discussed in [18, 20, 21, 35, 49].

Sensitivity Analysis

From the previous sections, it is evident that the infection dies out in the population when . Therefore, it becomes very important to control the model parameters so that becomes less than one. Therefore, determining the intervals at which the model parameters are sensitive is crucial. In similar lines to the sensitivity analysis done in [43], in this section we carry out sensitivity analysis of the model parameters. As each parameter is varied, we plot the infected population, the mean infected population, and the mean squared error as a function of time. These plots can be used to determine whether or not the parameter is sensitive in the given interval. The different intervals chosen are given in Table 5. The fixed parameter values are taken from Table 2.

Table 5

Sensitivity analysis

Parameter	Interval	Step size
u_11	0 to 0.05	0.001
	0.05 to 0.2
b_1	300 to 305	0.01
	305 to 310
	280 to 300
m	0 to 0.00182	0.0001
	0.00182 to 0.1	0.0001
u_12	0 to .05	0.001
	0.05 to 1.5
β_1	0 to 0.0028	0.0001
	0.0028 to 0.1
β_2	0 to 0.0028	0.0001
	0.0028 to 0.1
β_3	0 to 0.0.0086	0.0001
	0.0086 to 0.1
β_4	0 to 0.0086	0.0001
	0.0086 to 0.1
α	0 to 0.5	0.01
	0.5 to 2
d_1	0 to 0.013	0.001
	0.013 to 0.5	0.001
d_2	0 to 0.014	0.001
	0.014 to 0.5
μ	0 to 0.5	0.01
	0.5 to 2
δ_1	0 to 0.0714	0.001
	0.0714 to 1
δ_2	0 to 0.0714	0.001
	0.0714 to 1

Parameter

The results related to sensitivity of , varied in two intervals as mentioned in Table 5, are given in Fig. 3. The plots of infected population for each value of the parameter , the mean infected population and the mean square error are used to determine the sensitivity. Figure 3, we see that total infected population remains same for all the values of varied in both the interval given in Table 5. The mean infection decreases rapidly and gets cleared within few days. It is clearly seen that the mean square error of the total infected population increases in the beginning. However, this variation around the mean exists only for a short period and then the mean square error converges to zero. Since there is only one predictable fluctuation and standard deviation goes to negligible levels, we infer from this observation that is insensitive in both the interval I and II.

Fig. 3

Figure depicting the sensitivity Analysis of varied in three intervals in Table 5. The plots depict the infected population for each varied value of the parameter per interval along with the mean infected population and the mean square error in the same interval

The results related to sensitivity of , varied in three different intervals as mentioned in Table 5, are given in Fig. 4. The plots of infected population for each varied value of the parameter per interval, the mean infected population and the mean square error are used to determine the sensitivity. Figure 4, we see that total infected population varies with different values of and the mean square error also increases with time in the first two intervals. But in the third interval we see that the total infection do not vary with variations in the values of and the mean square error also decreases. It is inferred from this observation that insensitive in interval III and sensitive in I and II.

Fig. 4

Figure depicting the sensitivity analysis of varied in three intervals in Table 5. The plots depict the infected population for each varied value of the parameter per interval along with the mean infected population and the mean square error in the same interval

The results related to sensitivity of , varied in two intervals as mentioned in Table 5, are given in Fig. 5. The plots of infected population for each varied value of the parameter per interval, the mean infected population and the mean square error are used to determine the sensitivity. We conclude from these plots that the parameter is sensitive in interval I and insensitive in II.

Fig. 5

Figure depicting the sensitivity analysis of varied in three intervals in Table 5. The plots depict the infected population for each varied value of the parameter per interval along with the mean infected population and the mean square error in the same interval

The results related to sensitivity of , varied in two intervals as mentioned in Table 5, are given in Fig. 6. The plots of infected population for each varied value of the parameter per interval, the mean infected population and the mean square error are used to determine the sensitivity. We conclude from these plots that the parameter is sensitive in interval I and insensitive in II. In similar lines, the sensitivity analysis is done for other parameters. The results are summarized in Table 6. The corresponding plots are given in Appendix A owing to the brevity of the manuscript.

Fig. 6

Figure depicting the sensitivity analysis of varied in two intervals in Table 5. The plots depict the infected population for each varied value of the parameter per interval along with the mean infected population and the mean square error in the same interval

Table 6

Summary of sensitivity analysis

Parameter	Interval	Step size
u_11	0 to 0.05	✓
	0.05 to 0.2	×
b_1	300 to 305	✓
	305 to 310	✓
	280 to 300	×
m	0 to 0.00182	×
	.00182 to 1	×
u_12	0 to 0.05	×
	0.05 to 1.5	×
β_1	0 to 0.0028	✓
	0.0028 to 0.1	×
β_2	0 to 0.0028	×
	0.0028 to 0.1	×
β_3	0 to 0.0086	×
	0.0086 to 0.1	×
β_4	0 to 0.0086	×
	0.0086 to 0.1	×
α	0 to 0.5	×
	0.5 to 2	×
d_1	0 to 0.013	×
	0.013 to 0.5	✓
d_2	0 to 0.014	×
	0.014 to 0.5	×
μ	0 to 0.5	✓
	0.5 to 2	×
δ_1	0 to 0.0714	×
	0.0714 to 1	×
δ_2	0 to 0.0714	×
	0.0714 to 1	×

Optimal Control Problem

Optimal utilization of control interventions is one of the fundamental concern to the policy makers, health agencies and government. As there is always a trade-off in the implantation of suitable and available control interventions and also they consume huge amount of money, therefore, it is not only important to control the spread of the disease but also to minimize the total cost. In this section, we will formulate an optimal control problem to obtain an optimal treatment rate that minimizes COVID-19 infection as well as the overall cost incurred in the implementation of treatment strategy. The controls that we considered are:

The controls that we consider in our study is the treatment based recovery of young and adult infected population as a result of which an infected individual recovers. This treatment could be any common drug interventions that are recommended in the treatment of COVID-19. Many potential vaccines for COVID-19 are being studied and some are under clinical trials. Drugs such as remdesivir, favipiravir, ivermectin, lopinavir/ritonavir, mRNA-1273, phase I trial (NCT04280224) and AVT technology are being used as therapeutic agents by different countries for treating Covid-19 [44, 46, 51]. As resources and funds are limited, the efforts should be optimally used. Thus in optimal control study, the constant recovery rate due to treatment and , are taken as a functions of time, , (t). The aim is to investigate an optimal way of treatment that minimizes new infection as well as the total cost

The set of all admissible controls is given byWithout medical interventions, and are just constant parameters. We assume that the control variables are measurable bounded functions and the upper bounds and are based on the resource limitations. Based on these we now propose and define the optimal control problem with the goal to reduce the cost functional defined as follows,where subject to the systemThe drugs used for treating COVID-19 infected patients could have multiple effects, this justifies the quadratic terms in the definition of the cost functional [31]. Also when the objective function is defined as a linear combination of the quadratic terms of control variables the complexity of the problem reduces. Studies where objective function is considered as a linear combination of the quadratic terms of control variables can be found in [32, 39]. The integrand of the cost function (6.1), denoted byis called the Lagrangian or the running cost.

Here, the cost functional (6.1) represents the cumulative number of infected population and the overall cost incurred in implementation of the treatment strategy, throughout the observation period. Effectively, we want to minimize the infected population and the cost. Here, and are positive weight constants which not only balance units of integrand but also related cost.

The admissible solution set for the optimal control problem (6.1)–(6.2) is given by

Existence of Optimal Control

We will show the existence of optimal control functions that minimize the cost functions within a finite time span [0, T] showing that we satisfy the conditions stated in Theorem 4.1 of [27].

Theorem 6.1

There exists a 2-tuple of optimal controls in the set of admissible controls U such that the cost functional is minimized i.e.,corresponding to the optimal control problem (6.1)–(6.2).

Proof

In order to show the existence of optimal control functions, we will show that the following conditions are satisfied : Now we will show that each of the conditions are satisfied:

The solution set for the system (6.2) along with bounded controls must be non-empty, i.e., .

U is closed and convex and system should be expressed linearly in terms of the control variables with coefficients that are functions of time and state variables.

The Lagrangian should be convex on U and , where is a continuous function of control variables such that whenever , where |.| is an norm.

From positivity and boundedness of solutions of the system (6.2), all solutions are bounded for each bounded control variable in U. Also, the right hand side of the system (6.2) satisfies Lipschitz condition with respect to state variables. Hence, using the positivity and boundedness condition and the existence of solution from Picard–Lindel of Theorem [40], we have satisfied condition 1.

U is closed and convex by definition. Also, the system (6.2) is clearly linear with respect to controls such that coefficients are only state variables or functions dependent on time. Hence condition 2 is satisfied.

By definition Lagrangian is convex on U. Choosing such that , we can satisfy the condition 3. Hence there exists a control 2-tuple that minimizes the cost function (6.1).

Characterization of Optimal Control

We will obtain the necessary conditions for optimal control functions using the Pontryagin’s maximum principle [37] and also obtain the characteristics of the optimal controls.

The Hamiltonian for this problem is given byHere = (, , , , , ) is called co-state vector or adjoint vector.

Now the Canonical equations that relate the state variables to the co-state variables are given bySubstituting the Hamiltonian value gives the canonical systemalong with transversality conditions

Now, to obtain the optimal controls, we will use the Hamiltonian minimization condition = 0 , at for i = 11, 12 .

Differentiating the Hamiltonian and solving the equations, we obtain the optimal controls as

Simulations for Optimal Control Problem

In this section, we perform numerical simulations to understand the dynamics of the infection in different age groups along with treatments. Entire simulations are done using MATLAB software.

The various combinations of controls considered are: For our simulations, we have taken the total number of days as and the other fixed parameters from Table 2. We first solve the state system numerically using fourth order Runge–Kutta method in MATLAB without any interventions along with the initial values as

Implementation of treatment only for young population

Implementation of treatment only for adult population.

Implementation of treatment for both young and the adult population.

Now, to simulate the system with controls, we use the forward-backward sweep method with the initial values of controls and solve the state system forward in time. Following this we solve the adjoint state system backward in time due to the transversality conditions, using the optimal state variables and initial values of optimal control.

Now, using the values of adjoint state variables, the values of optimal control are updated and with these updated control variables, we go through this process again. We continue this till the convergence criterion is met [37]. In similar lines to [54] the positive weights chosen for objective coefficients are = 0.0001, = 0.001. We have chosen the weights related to the treatment of adults 10 times more than that of young ones because the cost related to the treatment of adult population is generally higher compared to the cost of treatment of young.

First, we solve system (6.2) in the absence of the controls and then simulate the system (6.2) with controls. The count of the infected population for both the age groups without the controls are shown in Figs. 7 and 8 with blue color. From Fig. 7 we see that the number of infected adult population decreases with the implementation of treatment strategy and the decrease in the infection is more when the second control is considered alone followed by considering both optimal controls and . The decrease in the infected adult population is found to be more considering alone compared to the combined use of optimal controls and This can be attributed to the type III recovery rate function considered for the second age group.

Fig. 7

under optimal controls ,

Fig. 8

under optimal controls ,

We note that when parameter values of the model (6.2) are taken arbitrary small, the standard numerical methods like Runge-Kutta method may not always converge. Therefore, in this cases adaptive numerical methods become necessary which works even for arbitrary small parameters values. Some of the studies in these lines that discusses numerical methods with convergence can be found in [9, 10, 15, 17, 19].

We also calculate the average values of infected adult population in Table 7 over the time period to support the above fact.

Table 7

Table depicting the average values of the infected adult population

Control combinations	Avg infected population (I_2)
μ11∗=0,μ12∗	0.2778
μ11∗,μ12∗	5.7477
μ11∗,μ12∗=0	10.2981
μ11∗=μ12∗=0	10.4165

In Fig. 8, we simulate the infected young population over the time considering different optimal control combinations. From Fig. 8 we see that the number of infection is higher when no treatment strategy is implemented and as optimal treatment strategy is implemented the number of infection reduces. In case of the young population we see that the reduction in infection is maximum when both and are considered followed by considering alone.

We also calculate the average values of the infected young population over the time period considered. From Table 8, we see that the average value of infected young population is least when both optimal treatment strategies are considered.

Table 8

Table depicting the average values of the infected young population

Control combinations	Avg infected population (I_1)
μ11∗,μ12∗	3.6180
μ11∗,μ12∗=0	3.6348
μ11∗=0,μ12∗	10.0538
μ11∗=μ12∗=0	10.1237

We also simulate the recovered population over the time with and without optimal controls in Figs. 9 and 10. From Fig. 9, we see that the recovered young population is higher when the control alone or combination of optimal controls and are considered. The recovered young population is found to remain almost constant at constant value in case of second control consideration. We also see that the recovered population curve considering alone and the curve without controls coincide with each other. The reason for this could be that the second control features only in the second age group (adult) population considered in the model.

Fig. 9

under optimal controls ,

Fig. 10

under optimal controls ,

From Fig. 10, we see that the recovery of adult population is higher considering the optimal control whereas, it is least in case of no control case and implementation of alone. Similar to the previous case, here too we observe that the recovered population curve considering alone and without controls coincide with each other and the reason for this could be that the control features only in the first age group (young) population considered in the model. Average values of the recovered population are calculated for young and adult population considering different combinations of the controls in Tables 9 and 10.

Table 9

Table depicting the average values of the recovered young population

Control combinations	Avg recovered population (R_1)
μ11∗=0,μ12∗	4.9833
μ11∗,μ12∗	11.4254
μ11∗,μ12∗=0	11.4195
μ11∗=μ12∗=0	4.9833

Table 10

Table depicting the average values of the recovered adult population

Control combinations	Avg recovered population (R_2)
μ11∗,μ12∗=0	4.9836
μ11∗,μ12∗	9.4601
μ11∗=0,μ12∗	14.8945
μ11∗=μ12∗=0	4.9834

In Fig. 11 we plot the cumulative infected population considering different optimal controls and we see that initially the cumulative count of the infection is minimum when second control is considered but over the time combination of optimal controls and is found to be effective.

Fig. 11

Cumulative infection () under different controls

The Effect of on Optimal Controls and on Disease Burden

This section devotes to study the dynamics of disease and the effect of optimal controls on the disease when the basic reproduction number varies. Since the severity of the epidemic characterized by the high epidemic peaks which is measured by the higher values of , therefore we will observe the prevalence of the cumulative count of the disease by varying the basic reproduction number. For this purpose, we will consider the basic reproduction number in absence of any controls i.e. and as follows for the system (2.1)–(2.6):where,We choose the parameters same as taken above for the comparative study to explore the impact of on optimal controls and on the disease. For this purpose, the total infected population is plotted to depict the severity of the disease by varying the degree of transmissibility as reflected by the basic reproduction number for different control strategies.

In Figs. 12, 13, 14, 15 we numerically study the effect of on the infected young population, infected adult and cumulative infected population considering different optimal controls. The findings from these simulations suggest the optimal controls alone works better than the combined control and for the entire range of for the second age groups, whereas, the combined treatment strategy consisting of and together is found to work better in minimizing the infection for the first age group (Figs. 12, 13).

Fig. 12

Effect of on under different controls

Fig. 13

Effect of on under different controls

Fig. 14

Effect of on cumulative infected population under different controls

Fig. 15

Effect of on cumulative recovered population under different controls

From Figs. 14 and 15 we see that when the epidemic is mild , the optimal control works better than the combined control and in minimizing the infection and maximizing the recovery of infected population. But, as severity of the epidemic increases , the combined treatment strategy consisting of and is found to be effective in keeping the infection low.

The Effect of and on Optimal Controls and Disease Burden

In this section we focus on studying the effect of saturation level in treatment by varying the value of on the disease burden for mild as well as severe epidemic for designed control strategy. In Fig. 16 we plot the effect of on the cumulative count of infected population considering both the controls and . From Fig. 16 we observe that the cumulative count of the disease grows as the level of saturation increases for entire range of . In Fig. 17 we plot the cumulative recovered population consider optimal controls and . From this figure we observe that as the saturation level in medical treatment increases, the recovery of infected population decreases for the entire range of basic reproduction number .

Fig. 16

Effect of on the cumulative disease burden with and

Fig. 17

Effect of on the cumulative disease burden with

Discussion and Conclusion

In this study, we first proposed a nonlinear age-structured compartmental model in which the total population was divided into two age groups, first group with age below 30 years ( we call this as young population) and the second with age above 30 years ( we call this group as adult population), and in each of these groups there were three compartments, namely the susceptible population, denoted by (young) and (adult), the infected population, denoted by (young) and (adult), and the recovered population, denoted by (young) and (adult). For the second age group, we consider the type II recovery rate function due to treatment to account for limited medical resources. Later, we proposed an optimal control problem to study the role of control strategies, namely the role of treatments in reducing infection. Pontryagin’s minimum principle was used to determine the optimal treatment rates that minimize both the number of infections and the cumulative costs. Numerical simulations were used to investigate the individual and combined roles of the optimal treatment strategies.

From the stability analysis of the equilibrium points it was concluded that the infection free equilibrium remains locally asymptotically stable whenever and as crosses unity the infected equilibrium is locally asymptotically stable. From the sensitivity analysis of the parameters of the model it was concluded that the parameters and are the sensitive in some intervals as discussed in Table 6.

The results of the optimal control studies indicate that infection is lowest in the adult population (age when optimal control is considered, whereas infection is lowest in the young population (age ) when a combination of optimal controls and is considered. Cumulative infection decreased most when the combination of optimal controls and was considered, followed by consideration of the second control alone.

The effects of and on disease burden considering different control strategies were also numerically investigated. The results suggest that in the case of a mild epidemic , the second control works better to minimize cumulative infection. On the other hand, when the severity of the epidemic increases , the combination of optimal controls and works better to keep infection low and maximize recovery. The study on the effect of on disease burden found that the cumulative number of infections increases with increasing saturation level for the whole range of considering the combination of optimal controls. In contrast, the cumulative recovery decreases with the increase of the saturation level for medical treatment.

Recovery of the infected population by treatment minimizes the severity of the outbreak. However, since the implementation of these control strategies involves costs, optimal use is the main concern. Optimal control theory has been used to obtain optimal control strategies that minimize the disease burden and also the total cost of implementing these strategies [6, 7, 24, 26, 28, 38, 42]. The optimal treatment strategy has been shown to be effective in keeping COVID-19 infection low [32, 36]. With the effective implementation of various control strategies like vaccination and treatment, almost all the countries of the world including India are successful in minimizing the COVID-19 cases [52]. The results obtained in this study are consistent with the real world data and can be useful to the physicians in treating COVID -19 infected individuals.

In conclusion we state the results of our work are in line and qualitatively consistent with the real data as in [52]. For instance from Figs. 7 and 8 we see the decrease of the infected are in similar pattern to decrease of infected in the Indian scenario [52].

We make a note here that the standard numerical methods on a uniform mesh may sometimes fail to provide an approximate solutions due to the presence of boundary layers. One effective approach in this case is to produce the suitable layer adapted meshes to capture the layer phenomena. Some of the studies in this line can be found in [13, 14, 16, 23]. In [22] the approximation of solutions of a type of fractional-order Volterra–Fredholm integro-differential equations is discussed along with the convergence analysis of the method. The present work is focused on understanding the dynamics of COVID-19 and the age-specific role of optimal treatment strategies. For this purpose, the standard Runge-Kutta method is used for performing the numerical studies. For further analysis, the methods discussed in [13, 14, 16, 22, 23] can be used and the results can be compared with the results of the present work.