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

Bishal Chhetri^1,a, D. K. K. Vamsi^∗,a, Carani B Sanjeevi^b,c
^aDepartment of Mathematics and Computer Science, Sri Sathya Sai Institute of Higher Learning, Prasanthi Nilayam,
Puttaparthi, Anantapur District - 515134, Andhra Pradesh, India
^b Vice-Chancellor, Sri Sathya Sai Institute of Higher Learning - SSSIHL, India
^c Department of Medicine, Karolinska Institute, Stockholm, Sweden
bishalchhetri@sssihl.edu.in, dkkvamsi@sssihl.edu.in^∗,
sanjeevi.carani@sssihl.edu.in, sanjeevi.carani@ki.se
¹ First Authors, ^∗ Corresponding Author

Abstract

COVID pandemic has catalysed the development of novel coronavirus vaccines across pharmaceutical companies and research organizations. In a situation where the vaccine is still unavailable and the disease is spreading exponentially, an effective control measures seems to be 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. This disease has caused the most severe health issues for adults over the age of 60 with particularly fatal results for those 80 years and old. In this situation an age structured modelling could be helpful in understanding the spread of the disease across different age groups. In this study initially, we propose an age structured model and calculate the equilibrium points and basic reproduction number. Later we propose an optimal control problem to understand the roles of treatment in controlling the epidemic. From the Stability analysis we see that the infection free equilibrium remains asymptotically stable whenever $R_{0}<1$ and as $R_{0}$ crosses unity we have the infected equilibrium to be stable. From the sensitivity analysis parameters $u_{11}$ , $b_{1}$ , $\beta_{1}$ , $d_{1}$ and $\mu$ were found to be sensitive. Findings from the Optimal Control studies suggests that the infection among the adult population(age $\geq 30)$ is least considering the second control $u_{12}$ whereas, when both the controls $u_{11}$ and $u_{12}$ are considered together the infectives is minimum in case of young populations(age $\leq 30$ ). The cumulative infected population reduced the maximum when the second control was considered followed by considering both the controls together. The control $u_{12}$ was effective for mild epidemic $(R_{0}\in(1,2))$ whereas control $u_{11}$ was found to be highly effective when epidemic was severe $(R_{0}\in(2,7))$ for the population of age group $(\leq 30)$ . Whereas for age group $(\geq 30)$ the control $u_{12}$ was highly effective for the entire range of basic reproduction number. The effect of saturation level in treatment is also explored numerically.

Introduction

Coronavirus disease 2019 (COVID-19) is a contagious respiratory and vascular disease caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). The Coronavirus caused disease COVID-19 has been declared a pandemic by WHO. As on 20 November 2020, around 55 million have been affected and around 1.3 million have lost their worldwide [3]. The social and economic life of people is drastically affected by this COVID-19 Pandemic.

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 [1, 2, 4, 18]. 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 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 populations [5]. A study with 425 patients indicated that the median age was 59 years (range, 15 to 89) and reported different case distributions in four age groups groups: 0 – 14, 15 – 44 , 45 73 – 64 , and ¿ 65 years [6]. 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 ([16], [11]).

Motivated by the above, in this study, we have proposed a non-linear age structured compartmental model in which the population is divided into two age groups. The first group with age $\leq$ 30(we call this group by young population) and the other group with age $\geq$ 30(we call this group by adult population) and in each of these age groups we have three compartments, namely susceptible denoted by $S_{1}$ (young) and $S_{2}$ (adult), infected denoted by $I_{1}$ (young) and $I_{2}$ (adult) and recovered denoted by $R_{1}$ (young) and $R_{2}$ (adult). Later, we frame an optimal control problem to understand the role of treatments by considering them as controls in the two age groups considered.

We have used Holling type III recovery rate function of infected adult individuals 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.

The section-wise split-up of the paper is as follows: In the next section we talk about the formulation of age structured model followed by the Positivity and Boundedness of the model. In section 3 the equilibrium points and the basic reproduction number is calculated following which the stability of equilibrium is analysed. In section 5 we perform the sensitivity analysis of the parameters of the model. We formulate an optimal control problem and do the simulations to see the role of the controls in reducing the infection in section 6. Finally, we present the discussions and conclusions.

1 Model Formulation

We have divided the entire population into two age groups, one with age $\leq 30$ (we call this group as young population) and the other ones with age $\geq 30$ (we call this group as adult population). In each of these groups we have three compartments namely susceptible( $S_{1},S_{2}$ ), Infected( $I_{1},I_{2}$ ) and the Recovered( $R_{1},R_{2}$ ). At any time an individual can be in one of the states namely susceptible, infected, or removed. Susceptibles are the ones who are not infected but can get infected when it comes in contact with he 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 $\beta_{1},\beta_{2},\beta_{3},\beta_{4}$ are the infection rates corresponding to two age groups considered. We assume that there is a constant birth rate of the young susceptible population denoted by parameter $b_{1}$ and also the recovered can again become susceptible again with rate $\delta_{1}$ and $\delta_{2}$ . Transition rate or maturation rate describing transition from age group to another is used and is denoted by parameters $m$ . We have assumed that the natural mortality rate is common across all the compartments and is denoted by parameter $\mu$ , whereas the disease induced death rate differs among age groups denoted by $d_{1}$ and $d_{2}$ . As the recovery from a disease differs among different age groups we denote the recovery rate for the first age group by $u_{11}$ and for the second group by $u_{12}$ . This model leads to the system of ordinary differential equations as follows:

$\displaystyle\frac{dS_{1}}{dt}$	$\displaystyle=$	$\displaystyle b_{1}+\delta_{1}R_{1}\ -\beta_{1}S_{1}I_{1}-\beta_{2}S_{1}I_{2}-\mu S_{1}-mS_{1}$	(1)
$\displaystyle\frac{dI_{1}}{dt}$	$\displaystyle=$	$\displaystyle\beta_{1}S_{1}I_{1}+\beta_{2}S_{1}I_{2}\ -d_{1}I_{1}\ -\mu I_{1}-\mu_{11}I_{1}$	(2)
$\displaystyle\frac{dR_{1}}{dt}$	$\displaystyle=$	$\displaystyle\mu_{11}I_{1}\ -\mu R_{1}-\delta_{1}R_{1}-mR_{1}$	(3)
$\displaystyle\frac{dS_{2}}{dt}$	$\displaystyle=$	$\displaystyle mS_{1}+\delta_{2}R_{2}\ -\beta_{3}S_{2}I_{1}-\beta_{4}S_{2}I_{2}-\mu S_{2}$	(4)
$\displaystyle\frac{dI_{2}}{dt}$	$\displaystyle=$	$\displaystyle\beta_{3}S_{2}I_{1}+\beta_{4}S_{2}I_{2}\ -d_{2}I_{2}\ -\mu I_{2}-\frac{\mu_{12}I_{2}^{2}}{1+\alpha I_{2}^{2}}$	(5)
$\displaystyle\frac{dR_{2}}{dt}$	$\displaystyle=$	$\displaystyle mR_{1}+\frac{\mu_{12}I_{2}^{2}}{1+\alpha I_{2}^{2}}\ -\mu R_{2}-\delta_{2}R_{2}$	(6)

Table 1: Parameters and their Meanings

Parameters	Biological Meaning
$S_{1}$	Suceptible young population
$S_{2}$	Suceptible 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
$\delta_{1}$	rate of recovered young becoming suceptible
$\beta_{1},\beta_{2}$	Rate at which Suceptible young population are infected
	because of infected young ones and infected adults
$\mu$	natural death rate
$d_{1},d_{2}$	disease induced death rate of young and adult popultaion
$\mu_{11}$	recovery rate of young ones because of treatment
$\delta_{2}$	rate of recovered adult becoming suceptible again
$\beta_{3},\beta_{4}$	Rate at which Suceptible adult population are infected
	because of infected young ones and infected adults
$\mu_{12},\alpha$	$\mu_{12}$ is the treatment rate of infected adult individuals
	and $\alpha$ is the non-negative constant related with saturation in treatment.
$m$	maturation rate

2 Positivity and Boundedness

We now show that if the initial conditions of the system (1)-(6) are positive, then the solution remain positive for any future time. Using the equations (1)-(6), we get,

$\displaystyle\frac{dS_{1}}{dt}\bigg{\|}_{S_{1}=0}$	$\displaystyle=\bigg{(}b_{1}S_{2}+\delta_{1}R_{1}\bigg{)}\geq 0,$	$\displaystyle\frac{dI_{1}}{dt}\bigg{\|}_{I_{1}=0}$	$\displaystyle=\beta_{2}S_{1}I_{2}\geq 0,$
$\displaystyle\frac{dR_{1}}{dt}\bigg{\|}_{R_{1}=0}$	$\displaystyle=\mu_{11}I_{1}\geq 0,$	$\displaystyle\frac{dS_{2}}{dt}\bigg{\|}_{S_{2}=0}$	$\displaystyle=\bigg{(}mS_{1}+\delta_{2}R_{2}\bigg{)}\geq 0$
$\displaystyle\frac{dI_{2}}{dt}\bigg{\|}_{I_{2}=0}$	$\displaystyle=\beta_{3}S_{2}I_{1}\geq 0,$	$\displaystyle\frac{dR_{2}}{dt}\bigg{\|}_{R_{2}=0}$	$\displaystyle=mR_{1}+\frac{\mu_{12}I_{2}^{2}}{1+\alpha I_{2}^{2}}\geq 0$

Thus all the above rates are non-negative on the bounding planes (given by $S_{1}=0$ , $I_{1}=0$ , $R_{1}=0$ , $S_{2}=0$ , $I_{2}=0$ , and $R_{2}=0$ ) 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 (1)-(6) remain positive for any time $t>0$ provided that the initial conditions are positive. This establishes the positivity of the solutions of the system (1)-(6). Next we will show that the solution is bounded.

Boundedness: Let $N(t)=S_{1}(t)+I_{1}(t)+R_{1}(t)+S_{2}(t)+I_{2}(t)+R_{2}(t)$

Now,

\begin{split}\frac{dN}{dt}&=\frac{dS_{1}}{dt}+\frac{dS_{2}}{dI_{1}}{dt}+\frac{dR_{1}}{dt}+\frac{dS_{2}(t)}{dt}+\frac{dI_{2}(t)}{dt}+\frac{dR_{2}(t)}{dt}\\[4.0pt] &\leq b_{1}-\mu(S_{1}+I_{1}+R_{1}+S_{2}+I_{2}+R_{2})\\ \end{split}

Here the integrating factor is $e^{\mu t}.$ Therefore after integration we get,

N(t)\leq\frac{b_{1}}{\mu}

Thus we have shown that the system (1)-(6) is positive and bounded. Therefore the biologically feasible region is given by the following set,

\Omega=\bigg{\{}\bigg{(}S_{1}(t),I_{1}(t),R_{1}(t),S_{2}(t),I_{2}(t),R_{2}(t)\bigg{)}\in\mathbb{R}^{6}_{+}:S_{1}(t)+I_{1}(t)+R_{1}(t)+S_{2}(t)+I_{2}(t)+R_{2}(t)\leq\frac{b_{1}}{\mu},\ t\geq 0\bigg{\}}

3 Equilibrium points and Reproduction Number

The system (1.1)-(1.6) admits two equilibrium namely the disease free equilibrium and the infected equilibrium. The disease free equilibrium denoted by $E_{0}$ is calculated to be,
$E_{0}=(S_{1}^{*},S_{2}^{*},0,0,0,0)$ where,

S_{1}^{*}=\frac{b_{1}}{(\mu+m)}

S_{2}^{*}=\frac{b_{1}m}{\mu(\mu+m)}

The infected equilibrium is denoted by $E_{1}$ and is given by,
$E_{1}=(S_{1}^{*},I_{1}^{*},R_{1}^{*},S_{2}^{*},I_{2}^{*},R_{2}^{*})$ where,

S_{1}^{*}=\frac{(b_{1}-d_{1}-\mu-u_{11})+\delta_{1}R_{1}^{*}}{(\mu+m)}

I_{1}^{*}=\frac{\beta_{1}(A+BR_{1}^{*})+\beta_{2}(A+BR_{1}^{*})}{(d_{1}+\mu+u_{11})}

R_{1}^{*}=\frac{u_{11}}{\mu+\delta_{1}+m}

S_{2}^{*}=\frac{m(A+BR_{1}^{*})+\delta_{2}R_{2}^{*}}{\beta_{3}+\beta_{4}I_{2}^{*}+\mu}

I_{2}^{*}=\frac{b_{1}+\delta_{1}R_{1}^{*}-\beta_{1}(A+BR_{1}^{*})-(\mu+m)(A+BR_{1}^{*})}{\beta_{2}(A+BR_{1}^{*})}

R_{2}^{*}=\frac{mR_{1}^{*}+\frac{u_{12}I_{2}^{*2}}{1+\alpha I_{2}^{*2}}}{(\mu+\delta_{2})}

Here,

A=\frac{(b_{1}-d_{1}-\mu-u_{11})}{(\mu+m)}

B=\frac{\delta_{1}}{(\mu+m)}

3.1 Calculation of $R_{0}$

The basic reproduction number which is the average number of secondary cases produced per primary case is calculated using the next generation matrix method described in [8]. Our system (1.1)-(1.6) has two infected states and four uninfected states. Calculating the jacobian matrix at infection free equilibrium $E_{0}$ considering infected states we have,

J(E_{0})=\begin{bmatrix}\beta_{1}S_{1}^{*}-d_{1}-\mu-u_{11}&\beta_{2}S_{2}^{*}\\[6.0pt] \beta_{3}S_{2}^{*}&\beta_{4}S_{2}^{*}-d_{2}-\mu\\[6.0pt] \end{bmatrix}

or,

J(E_{0})=F+V

where,

F=\begin{bmatrix}\beta_{1}S_{1}^{*}&\beta_{2}S_{2}^{*}\\[6.0pt] \beta_{3}S_{2}^{*}&\beta_{4}S_{2}^{*}\\[6.0pt] \end{bmatrix}

V=\begin{bmatrix}-(d_{1}+\mu+u_{11})&0\\[6.0pt] 0&-d_{2}-\mu\\[6.0pt] \end{bmatrix}

Calculating the inverse of $V$ we get,

V^{-1}=\begin{bmatrix}\frac{-1}{(d_{1}+\mu+u_{11})}&0\\[6.0pt] 0&\frac{-1}{d_{2}+\mu}\\[6.0pt] \end{bmatrix}

Now

-FV^{-1}=\begin{bmatrix}\beta_{1}S_{1}^{*}p&\beta_{2}S_{1}^{*}q\\[6.0pt] \beta_{3}S_{2}^{*}p&\beta_{4}S_{2}^{*}q\\[6.0pt] \end{bmatrix}

where,

p=\frac{1}{d_{1}+\mu+u_{11}}

q=\frac{1}{d_{2}+\mu}

The characterstics equation of $-FV^{-1}$ is calculated as,

\lambda^{2}-\bigg{(}\beta_{1}S_{1}^{*}p+\beta_{4}S_{2}^{*}q\bigg{)}\lambda+S_{1}^{*}S_{2}^{*}pq\bigg{(}\beta_{1}\beta_{4}-\beta_{2}\beta_{3}\bigg{)}

The spectral radius of $-FV^{-1}$ is given by,

\ \mathbf{\frac{(\beta_{1}S_{1}^{*}p+\beta_{4}S_{2}^{*}q)+\sqrt{M}}{2}}

where,

M=(\beta_{1}S_{1}^{*}p-\beta_{4}S_{2}^{*}q)^{2}+4S_{1^{*}}S_{2}^{*}\beta_{2}\beta_{3}pq

Therefore the basic reproduction number is given by,

\mathbf{R_{0}}=\mathbf{\frac{(\beta_{1}S_{1}^{*}p+\beta_{4}S_{2}^{*}q)+\sqrt{M}}{2}}

4 Numerical Simulation

Parameters Values In the following table we provide the parameter values and the source from which the parameter values are taken. Since the Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) mimics the influenza virus regarding clinical presentation, transmission mechanism, and seasonal coincidence [7] and also due to the lack of data on the age structure modelling of COVID-19, we have taken some of the parameter values from the work on influenza disease [11]. In similar lines to ([11],[12]) for numerical simulation we take the values of $\delta_{1}$ and $\delta_{2}$ and $u_{11}$ and $u_{12}$ to be same. Taking parameter values from table 2, we will numerically show the asymptotic stability of $E_{0}$ and $E_{1}$ based on the values of basic reproduction number in similar lines to [11].

Table 2: Parameters values and their Source

Parameters	Values	Source
$b_{1}$	7.192	[17]
$\delta_{1}$	.0714	[11]
$\beta_{1},\beta_{2}$	4/3, 2	[11]
$\mu$	.062	[17]
$d_{1},d_{2}$	.000073, .0000913	[11]
$\mu_{11}$	.1	-
$\delta_{2}$	.0714	[11]
$\beta_{3},\beta_{4}$	4, .00000008	[11]
$\mu_{12}$	.1	[12]
$\alpha$	0.4	[12]
$m$	0.000182	[11]

Stability of infection free equilibrium
We will numerically show that whenever we have $R_{0}<1$ the infection free equilibrium $E_{0}$ is locally asymptotically stable. For the parameter values from table 3, the value of $R_{0}$ was calculated as 0.98 and $E_{0}=(0.1157,.00039,0,0,0,0)$ . From figure 1 we see that the disease free equilibrium $E_{0}$ is locally asymptotically stable. The initial values were taken as $(S1,S2,I1,I2,R1,R2)=(100,150,5,70,10,30)$

Refer to caption — Figure 1: Stability of $E_{0}$

Table 3: Parameters values for

R_{0}<1

Parameters	Values
$b_{1}$	.007192
$\delta_{1}$	.0714
$\beta_{1},\beta_{2}$	4/3, 2
$\mu$	.062
$d_{1},d_{2}$	.000073, .0000913
$\mu_{11}$	.1
$\delta_{2}$	.0714
$\beta_{3},\beta_{4}$	4, .00000008
$\mu_{12}$	.1
$\alpha$	0.4
$m$	0.000182

Stability of infected equilibrium
Here we shall numerically show that whenever the value of basic reproduction number denoted by $R_{0}$ crosses unity we have the infected equilibrium $E_{1}$ to be locally asymptotically stable. The value of $R_{0}$ was calculated to be 2.7615 and $E_{1}=(10.42,1,0.144,0.0041,0.03,0.0005)$ for the parameter values from table 4 . From figure 2 we see that the infected equilibrium $E_{1}$ is locally asymptotically stable provided $R_{0}>1$ . The initial values were taken as $(S_{1},S_{2},I_{1},I_{2},R_{1},R_{2})=(100,200,1,2,1,2)$ .

Table 4: Parameters values for

R_{0}>1

Parameters	Values
$b_{1}$	7.192
$\delta_{1}$	.0714
$\beta_{1},\beta_{2}$	.0133, 2
$\mu$	.62
$d_{1},d_{2}$	.000073, .0000913
$\mu_{11}$	.1
$\delta_{2}$	.0714
$\beta_{3},\beta_{4}$	4, .00000008
$\mu_{12},\alpha$	.1, .5
$m$	0.00182

Corona virus has caused the most severe health issues for adults over the age of 60 — with particularly fatal results for those 80 years and older[5]. With the parameter values from table 2, now we will numerically illustrate this fact. From figure 3 we see that the number of infected cases in case of adult population is comparatively higher than that in case of young.

5 Sensitivity Analysis

In this section we check the sensitivity of all the parameters of the model (1.1)-(1.6). As each parameter is varied in different intervals, the total infected cell population, mean infected cell population and the mean square error are plotted with respect to time. These plots are used to determine the sensitivity of the parameter. The different intervals chosen are given in the following table 5. The fixed parameter values are taken from table 2.

Table 5: Sensitivity Analysis

Parameter	Interval	Step Size
$u_{11}$	0 to 0.5	0.01
	1.5 to 2
$b_{1}$	6.5 to 7.1920	0.01
	7.1920 to 8
	0.1 to 0.5
$m$	0 to 0.00182	0.0001
	.00182 to 1	0.01
$u_{12}$	0 to .5	0.05
	0.5 to 2
$\beta_{1}$	0 to 1.33	0.01
	1.33 to 2
$\beta_{2}$	0 to 2	0.01
	2 to 3
$\beta_{3}$	0 to 2.5	0.01
	2.5 to 5
$\beta_{4}$	0 to 0.5	0.01
	0.5 to 1
$\alpha$	0 to 0.5	0.01
	0.5 to 2
$d_{1}$	0 to 0.000073	0.00001
	0.000073 to 1	0.01
$d_{2}$	0 to 0.0000913	0.00001
	0.0000913 to 2
$\mu$	0 to 0.5	0.01
	0.5 to 2
$\delta_{1}$	0 to 0.0714	0.001
	0.0714 to 1
$\delta_{2}$	0 to 0.0714	0.001
	0.0714 to 1

5.0.1 Parameter ${b_{1}}$

The results related to sensitivity of $b_{1}$ , varied in three intervals as mentioned in table 5, are given in figure (4-6). The plots of infected population for each varied value of the parameter $b_{1}$ 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 $b_{1}$ is not sensitive in interval III and sensitive in I and II.

5.0.2 Parameter $\boldsymbol{\mu}$

The results related to sensitivity of $\mu$ , varied in two intervals as mentioned in table 5, are given in figure (7-8). The plots of infected population for each varied value of the parameter $\mu$ 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 $\mu$ is sensitive in interval I and insensitive in II.

5.0.3 Parameter $\boldsymbol{\beta_{1}}$

The results related to sensitivity of $\beta_{1}$ , varied in two intervals as mentioned in table 5, are given in figure (9-10). The plots of infected population for each varied value of the parameter $\beta_{1}$ 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 $\beta_{1}$ 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.

Table 6: Summary of Sensitivity Analysis

Parameter	Interval	Step Size
$u_{11}$	0 to 0.5	✓
	1.5 to 2	$\times$
$b_{1}$	6.5 to 7.1920	✓
	7.1920 to 8	✓
	0.1 to 0.5	$\times$
$m$	0 to 0.00182	$\times$
	.00182 to 1	$\times$
$u_{12}$	0 to .5	$\times$
	0.5 to 2	$\times$
$\beta_{1}$	0 to 1.33	✓
	1.33 to 2	$\times$
$\beta_{2}$	0 to 2	$\times$
	2 to 3	$\times$
$\beta_{3}$	0 to 2.5	$\times$
	2.5 to 5	$\times$
$\beta_{4}$	0 to 0.5	$\times$
	0.5 to 1	$\times$
$\alpha$	0 to 0.5	$\times$
	0.5 to 2	$\times$
$d_{1}$	0 to 0.000073	$\times$
	0.000073 to 1	✓
$d_{2}$	0 to 0.0000913	$\times$
	0.0000913 to 2	$\times$
$\mu$	0 to 0.5	✓
	0.5 to 2	$\times$
$\delta_{1}$	0 to 0.0714	$\times$
	0.0714 to 1	$\times$
$\delta_{2}$	0 to 0.0714	$\times$
	0.0714 to 1	$\times$

6 OPTIMAL CONTROL PROBLEM

In this section, we will formulate an optimal control problem to see the role of treatment in reducing the number of infection. The controls that we considered are:

1.

The first control that we consider is the treatment based recovery of young infected individuals as a result of which an infected young individual recovers. This treatment could be any common drugs 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 [1, 2, 4, 18]. We denote this intervention by the control variable $u_{11}$ .
2.

The second control that we consider is the treatment based recovery of adult Population. We consider a limited treatment as $\frac{\mu_{2}I_{2}^{2}}{1+\alpha I_{2}^{2}}$ where $\mu_{2}(t)$ is a constant treatment rate. Here $\frac{\mu_{2}I_{2}^{2}}{1+\alpha I_{2}^{2}}$ represents Holling type III response where the initial treatment is less because of less infectives and as the disease spreads, the treatment also increases and non-linearly saturates due to the limitation in the medical facilities.

The set of all admissible controls is given by

$U=\left\{(u_{11}(t),u_{12}(t)):u_{11}(t)\in[0,u_{11}max],u_{12}(t)\in[0,u_{12}max],t\in[0,T]\right\}$

Without medical interventions, $u_{11}$ and $\;u_{12},\;$ are just constant parameters.

In order to reduce the complexity of the problem here we choose to model the control efforts via a linear combination of the quadratic terms. Also when the objective function is quadratic with respect to the control, diﬀerential equations arising from optimization have a known solution. Other functional forms sometimes lead to systems of diﬀerential equations that are diﬃcult to solve ([9], [13]). Based on these we now propose and define the optimal control problem with the goal to reduce the cost functional defines as follows,

J(u_{11}(t),u_{12}(t))=\int_{0}^{T}(A_{1}u_{11}(t)^{2}+A_{2}u_{12}(t)^{2}+I_{1}(t)+I_{2}(t))dt

(7)

where $u=(u_{11}(t),u_{12}(t))\in U$

subject to the system

\displaystyle\begin{aligned} \frac{dS_{1}}{dt}&=&b_{1}+\delta_{1}R_{1}\ -\beta_{1}S_{1}I_{1}-\beta_{2}S_{1}I_{2}-\mu S_{1}-mS_{1}\\ \frac{dI_{1}}{dt}&=&\beta_{1}S_{1}I_{1}+\beta_{2}S_{1}I_{2}\ -d_{1}I_{1}\ -\mu I_{1}-\mu_{11}(t)I_{1}\\ \frac{dR_{1}}{dt}&=&\mu_{11}(t)I_{1}\ -\mu R_{1}-\delta_{1}R_{1}-mR_{1}\\ \frac{dS_{2}}{dt}&=&mS_{1}+\delta_{2}R_{2}\ -\beta_{3}S_{2}I_{1}-\beta_{4}S_{2}I_{2}-\mu S_{2}\\ \frac{dI_{2}}{dt}&=&\beta_{3}S_{2}I_{1}+\beta_{4}S_{2}I_{2}\ -d_{2}I_{2}\ -\mu I_{2}-\frac{\mu_{12}(t)I_{2}^{2}}{1+\alpha I_{2}^{2}}\\ \frac{dR_{2}}{dt}&=&mR_{1}+\frac{\mu_{12}(t)I_{2}^{2}}{1+\alpha I_{2}^{2}}\ -\mu R_{2}-\delta_{2}R_{2}\end{aligned}

(8)

The integrand of the cost function (6.1), denoted by

L(S,I,V,u_{11},u_{12})=(I_{1}(t)+I_{2}(t)+A_{1}u_{11}(t)^{2}+A_{2}u_{12}(t)^{2})

is called the Lagrangian or the running cost.

Here, the cost function represents the number of infected population throughout the observation period, and the overall cost of implementation of the treatments. Effectively, we want to minimize the infected population and the cost. Here, $A_{1}$ and $A_{2}$ 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

$\Omega=\left\{(S_{1},I_{1},R_{1},S_{2},I_{2},R_{2},u_{11},u_{12})\;|\;S_{1},I_{1},S_{2},I_{2},R_{1},R_{2}\text{that satisfy }(8)\right\}$

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 [10].

Theorem 1.

There exists a 2-tuple of optimal controls $(u_{11}^{*}(t),u_{12}^{*}(t))$ in the set of admissible controls U such that the cost functional is minimized i.e.,

J[u_{11}^{*}(t),u_{12}^{*}(t)]=\min_{(u_{11}^{*},u_{12}^{*})\in U}\bigg{\{}J[u_{11},u_{12}]\bigg{\}}

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 :

1.

The solution set for the system (6.2) along with bounded controls must be non-empty, $i.e.$ , $\Omega\neq\phi$ .
2.

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.
3.

The Lagrangian L should be convex on U and $L(S,I,V,u_{11},u_{12})\geq g(u_{11},u_{12})$ , where $g(u_{11},u_{12})$ is a continuous function of control variables such that $|(u_{11},u_{12})|^{-1}g(u_{11},u_{12})\to\infty$ whenever $|(u_{11},u_{12})|\to\infty$ , where $|.|$ is an $l^{2}(0,T)$ norm.

Now we will show that each of the conditions are satisfied :

1. 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-Lindelof Theorem[15], we have satisfied condition 1.

2. $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.

3. Choosing $g(u_{11},u_{12})=c(u_{11}^{2}+u_{12}^{2})$ such that $c=min\left\{A_{1},A_{2}\right\}$ , we can satisfy the condition 3.

Hence there exists a control 2-tuple $(u_{11}^{*},u_{12}^{*})\in U$ 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 [14] and also obtain the characteristics of the optimal controls.

The Hamiltonian for this problem is given by

H(S,I,V,u_{11},u_{12},u_{2},\lambda):=L(S_{1},I_{1},R_{1},S_{2},I_{2},R_{2},u_{11},u_{12})+\lambda_{1}\frac{\mathrm{d}S_{1}}{\mathrm{d}t}+\lambda_{2}\frac{\mathrm{d}S_{2}}{\mathrm{d}t}+\lambda_{3}\frac{\mathrm{d}I_{1}}{\mathrm{d}t}+\lambda_{4}\frac{\mathrm{d}I_{2}}{\mathrm{d}t}+\lambda_{5}\frac{\mathrm{d}R_{1}}{\mathrm{d}t}+\lambda_{6}\frac{\mathrm{d}R_{2}}{\mathrm{d}t}

Here $\lambda$ = ( $\lambda_{1}$ , $\lambda_{2}$ , $\lambda_{3}$ , $\lambda_{4}$ , $\lambda_{5}$ , $\lambda_{6}$ ) is called co-state vector or adjoint vector.

Now the Canonical equations that relate the state variables to the co-state variables are given by

$\displaystyle\frac{\mathrm{d}\lambda_{1}}{\mathrm{d}t}$	$\displaystyle=-\frac{\partial H}{\partial S_{1}}$	(9)
$\displaystyle\frac{\mathrm{d}\lambda_{2}}{\mathrm{d}t}$	$\displaystyle=-\frac{\partial H}{\partial S_{2}}$
$\displaystyle\frac{\mathrm{d}\lambda_{3}}{\mathrm{d}t}$	$\displaystyle=-\frac{\partial H}{\partial I_{1}}$
$\displaystyle\frac{\mathrm{d}\lambda_{4}}{\mathrm{d}t}$	$\displaystyle=-\frac{\partial H}{\partial I_{2}}$
$\displaystyle\frac{\mathrm{d}\lambda_{5}}{\mathrm{d}t}$	$\displaystyle=-\frac{\partial H}{\partial R_{1}}$
$\displaystyle\frac{\mathrm{d}\lambda_{6}}{\mathrm{d}t}$	$\displaystyle=-\frac{\partial H}{\partial R_{2}}$

Substituting the Hamiltonian value gives the canonical system

$\displaystyle\frac{\mathrm{d}\lambda_{1}}{\mathrm{d}t}$	$\displaystyle=\lambda_{1}(\beta_{1}I_{1}+\beta_{2}I_{2}+m+\mu)+\lambda_{2}m+\lambda_{3}(\beta_{1}I_{1}+\beta_{2}I_{2})$	(10)
$\displaystyle\frac{\mathrm{d}\lambda_{2}}{\mathrm{d}t}$	$\displaystyle=\lambda_{2}(\beta_{3}I_{1}+\beta_{4}I_{2}+\mu)-\lambda_{4}(\beta_{3}I_{1}+\beta_{4}I_{2})$
$\displaystyle\frac{\mathrm{d}\lambda_{3}}{\mathrm{d}t}$	$\displaystyle=-1+\lambda_{1}\beta_{1}S_{1}+\lambda_{2}\beta_{3}S_{2}-\lambda_{3}(\beta_{1}S_{1}-d_{1}-\mu-\mu_{11}(t))-\lambda_{4}\beta_{3}S_{2}-\lambda_{5}\mu_{11}(t)$
$\displaystyle\frac{\mathrm{d}\lambda_{4}}{\mathrm{d}t}$	$\displaystyle=\lambda_{1}\beta_{2}S_{1}+\lambda_{2}\beta_{4}S_{2}-\lambda_{3}\beta_{2}S_{1}-\lambda_{4}\bigg{(}\beta_{4}S_{2}-d_{2}-\mu-\frac{2\mu_{12}(t)I_{2}}{1+\alpha I_{2}^{2}}\bigg{)}-\lambda_{6}\bigg{(}\frac{2\mu_{12}(t)I_{2}}{1+\alpha I_{2}^{2}}\bigg{)}$
$\displaystyle\frac{\mathrm{d}\lambda_{5}}{\mathrm{d}t}$	$\displaystyle=-\lambda_{1}\delta_{1}+\lambda_{5}(\mu+\delta_{1}+m)-\lambda_{6}m$
$\displaystyle\frac{\mathrm{d}\lambda_{6}}{\mathrm{d}t}$	$\displaystyle=-\lambda_{2}\delta_{2}+\lambda_{6}(\mu+\delta_{2})$

along with transversality conditions $\lambda_{1}(T)=0,\ \lambda_{2}(T)=0,\ \lambda_{3}(T)=0.\lambda_{4}(T)=0,\ \lambda_{5}(T)=0,\ \lambda_{6}(T)=0.$

Now, to obtain the optimal controls, we will use the Hamiltonian minimization condition $\frac{\partial H}{\partial u_{i}}$ = 0 , at $u_{i}=u_{i}^{*}$ for i = 11, 12 .

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

	$\displaystyle u_{11}^{*}$	$\displaystyle=$	$\displaystyle\min\bigg{\{}\max\bigg{\{}\frac{(\lambda_{3}-\lambda_{5})I_{1}}{2A_{1}},0\bigg{\}},u_{11}max\bigg{\}}$
	$\displaystyle u_{12}^{*}$	$\displaystyle=$	$\displaystyle\min\bigg{\{}\max\bigg{\{}\frac{(\lambda_{4}-\lambda_{6})I_{2}^{2}}{2A_{2}(1+\alpha I_{2}^{2})},0\bigg{\}},u_{12}max\bigg{\}}$

7 Simulations for Optimal Control Problem

In this section, we perform numerical simulations to understand the role of the treatments in reducing the infection. Entire simulations is done using MATLAB software.

The various combinations of controls considered are:

1. Implementation of treatment only for young population

2. Implementation of treatment only for adult population.

3. Implementation of treatment for both young and the adult population.

For our simulations, we have taken the total number of days as $T=100$ 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 $(S_{1},S_{2},I_{1},I_{2},R_{1},R_{2})=(100,100,10,10,5,5)$

Now, to simulate the system with controls, we use the Forward-Backward Sweep method stating 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 [14]. In similar lines to [19] the positive weights chosen for objective coefficients are $A_{1}$ = .0001, $A_{2}$ = .005. 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 $(u_{11}=0,u_{12}=0)$ along with same initial population size. The corresponding count of the infectives $(I_{1},I_{2})$ are shown in figures 11 and 12 with blue color. From figure 11 we see that the number of infected adult cases considering the implementation of controls are less than the case considering no controls and best result is obtained when the control $u_{12}$ is considered alone as the number of infected adult population is least compared to the other two cases.

We also calculate the average values of infected adult population in table 7 over the time period $t=100$ to support the above fact.

Table 7: Table depicting the average values of the infected adult population.

Control Combinations	Avg Infected Population(I2)
$u_{11}^{}=0,u_{12}^{}$	0.7185
$u_{11}^{},u_{12}^{}$	24.5626
$u_{11}^{},u_{12}^{}=0$	30.1420
$u_{11}^{}=u_{12}^{}=0$	32.6219

In figure 12 we simulate the infected young population over the time considering different control combinations. We see from the figure that the number of infectives is higher when no control is implemented and as treatment is implemented the number of infectives reduces and reduces the maximum when both $u_{11}$ and $u_{12}$ is considered followed by considering only $u_{11}$ .

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 in case when both the controls are implemented.

Table 8: Table depicting the average values of the infected young population.

Control Combinations	Avg Infected Population(I1)
$u_{11}^{},u_{12}^{}$	15.1129
$u_{11}^{},u_{12}^{}=0$	15.2773
$u_{11}^{}=0,u_{12}^{}$	22.9808
$u_{11}^{}=u_{12}^{}=0$	31.9349

We also simulate the recovered population over the time with and without controls in figure 13 and figure 14. From Figure 13 we see that the recovered young population is highest when the control $u_{11}$ is considered. Whereas the recovered young population remains almost constant at constant value 5 in case of second control $u_{12}$ consideration. From figure 13 we also see that the recovered population curve considering $u_{12}$ alone and the curve without controls coincide with each other. The reason for this could be that the second control $u_{12}$ features only in the second age group(adult) population considered in the model.

From figure 14 we see that the recovered adult population is highest considering the control $u_{12}$ whereas it is least in case of no control consideration and implementation of only $u_{11}$ . Similar to the previous case we see that the recovered population curve considering $u_{11}$ alone and the curve without controls coincide with each other and the reason for this could be that the control $u_{11}$ features only in the first age group(young) population considered in the model. In figure 15 we plot the cumulative infected population over the time considering different controls and we see that the cumulative count of the infectives is minimum when the second control $u_{12}$ is considered followed by considering both the controls. Average values of the recovered population are calculated for young and adult population considering different combinations of the controls in table 9 and table 10.

Table 9: Table depicting the average values of the Recovered young population.

Control Combinations	Avg Recovered Population(R1)
$u_{11}^{}=0,u_{12}^{}$	34.9995
$u_{11}^{},u_{12}^{}$	10.4040
$u_{11}^{},u_{12}^{}=0$	4.9834
$u_{11}^{}=u_{12}^{}=0$	4.9834

Table 10: Table depicting the average values of the Recovered adult population.

Control Combinations	Avg Recovered Population(R2)
$u_{11}^{},u_{12}^{}=0$	20.4898
$u_{11}^{},u_{12}^{}$	19.9415
$u_{11}^{}=0,u_{12}^{}$	4.9833
$u_{11}^{}=u_{12}^{}=0$	4.9833

8 The Effect of $R_{0}$ 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 $(R_{0})$ varies. Since the severity of the epidemic characterized by the high epidemic peaks which is measured by the higher values of $R_{0}$ , 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 $R_{0}$ in absence of any controls i.e. $u_{11}=0$ and $u_{12}=0$ as follows for the system (1.1)-(1.6):

\mathbf{R_{0}}=\mathbf{\frac{(\beta_{1}S_{1}^{*}p+\beta_{4}S_{2}^{*}q)+\sqrt{M}}{2}}

where,

p=\frac{1}{d_{1}+\mu}

q=\frac{1}{d_{2}+\mu}

We choose the parameters same as taken above for the comparative study to explore the impact of $R_{0}$ 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 $R_{0}$ for different control strategies.

In figure (16,17,18) we plot the effect of $R_{0}$ on the infected young population, infected adult and the cumulative infected population considering different controls. Our findings suggest that when the epidemic is mild $(R_{0}\in(1;1.5))$ , the second control $u_{12}$ treatment works better than the control $u_{11}$ as given in yellow and violet colored curves in figure 16. Whereas when severity of the epidemic increases $(R_{0}\in(1.5;7))$ , the effect of the controls $u_{11}$ and $u_{11}$ and $u_{12}$ together is found to be better than the effect of $u_{12}$ alone. From figure 17 and figure 18 we see that the treatment $u_{12}$ works better in keeping the severity very low throughout the range of the basic reproduction number $(R_{0}\in(1;7))$ as shown in yellow colored curve in the figure 17 followed by the combined effect of both the controls. In figure 19 we plot the cumulative recovered population and see that with the control $u_{12}$ , we have the highest recovered population throughout the range of the basic reproduction number $(R_{0}\in(1;7))$ .

9 The Effect of $\alpha$ and $R_{0}$ 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 $\alpha$ on the disease burden for mild as well as severe epidemic for each designed control strategy. In figure 20 we plot the effect of $\alpha$ on the cumulative count of infected population considering both the controls $u_{11}$ and $u_{12}$ . We observed from figure 20 that the cumulative count of the disease grows significantly as the level of saturation increases for entire range of $R_{0}$ . In figure 21 we consider the control $u_{12}$ alone and study the effect of $\alpha$ by varying the value of $\alpha$ on the disease burden. Here again we see that the cumulative count of the disease grows significantly as the level of saturation increases but as epidemic becomes severe $(R_{0}>5)$ the effect of $\alpha$ remains the same as the cumulative count of infection saturates.

10 Discussion and Conclusion

In this study, initially we have proposed a non-linear age structured compartmental model in which the population is divided into two age groups the young (age $\leq$ 30) and the adult ones ( $\geq$ 30) and in each of to age groups we have three compartments, namely susceptible denoted by $S_{1}$ (young) and $S_{2}$ (adult), infected denoted by $I_{1}$ (young) and $I_{2}$ (adult) and recovered denoted by $R_{1}$ (young) and $R_{2}$ (adult). We have used Holling type III recovery rate function of infected adult individuals 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. Later we proposed an optimal control problem to investigate the role of the control strategies namely the role of treatments in reducing the infection.

From the Stability analysis we conclude that the infection free equilibrium remains asymptotically stable whenever $R_{0}<1$ and as $R_{0}$ crosses unity we have the infected equilibrium to be stable. From the sensitivity analysis of the parameters of the model we conclude that the parameters $u_{11},b_{1},\beta_{1},d_{1}$ and $\mu$ are the only sensitive parameters in some intervals as discussed in table 6. Findings from the Optimal Control studies suggests that the infection among the adult population(age $\geq 30)$ is least considering the second control $u_{12}$ whereas, when both the controls $u_{11}$ and $u_{12}$ are considered together the infectives is minimum in case of young populations(age $\leq 30$ ). The cumulative infected population reduced the maximum when the second control was considered followed by considering both the controls together.

We have also studied the effect of $R_{0}$ and $\alpha$ on the disease burden considering different control strategies. Our findings suggests that when the epidemic is mild $(R_{0}\in(1;1.5))$ , the second control $u_{12}$ treatment works better than the control $u_{11}$ as given in yellow and violet colored curves in figure 16. Whereas when severity of the epidemic increased $(R_{0}\in(1.5;7))$ , the effect of the controls $u_{11}$ and $u_{11}$ and $u_{12}$ together was found to be better than the effect of $u_{12}$ alone. From figure 17 and figure 18 we see that the treatment $u_{12}$ works better in keeping the severity very low throughout the range of the basic reproduction number $(R_{0}\in(1;7))$ as shown in yellow colored curve in the figure 17 followed by the combined effect of both the controls.

From the study on the effect of $\alpha$ on disease burden we saw that the cumulative count of the disease grew significantly as the level of saturation increased for entire range of $R_{0}$ considering both the controls together. Whereas when we considered the second control alone we saw that the cumulative count of the disease grew significantly as the level of saturation increased but as epidemic became severe $(R_{0}>5)$ the effect of $\alpha$ remained the same.

References

[1] https://pib.gov.in/newsite/printrelease.aspx?relid=201174.
[2] https://www.bloomberg.com/news/articles/2020-04-10/two-thirds-of-severe-covid-19-improved-on-gilead-s-remdesivir.
[3] https://www.google.com/search?client=firefox-b-dq=covid-19+world+statistics.
[4] https://www.nih.gov/news-events/news-releases/nih-clinical-trial-investigational-vaccine-covid-19-begins.
[5] https://www.nwhn.org/how-does-covid-19-affect-different-age-groups/.
[6] Yan Bai, Lingsheng Yao, Tao Wei, Fei Tian, Dong-Yan Jin, Lijuan Chen, and Meiyun Wang, Presumed asymptomatic carrier transmission of covid-19, Jama 323 (2020), no. 14, 1406–1407.
[7] Elena Cuadrado-Payán, Enrique Montagud-Marrahi, Manuel Torres-Elorza, Marta Bodro, Miquel Blasco, Esteban Poch, Alex Soriano, and Gaston J Piñeiro, Sars-cov-2 and influenza virus co-infection, Lancet (London, England) 395 (2020), no. 10236, e84.
[8] Odo Diekmann, JAP Heesterbeek, and Michael G Roberts, The construction of next-generation matrices for compartmental epidemic models, Journal of the Royal Society Interface 7 (2010), no. 47, 873–885.
[9] Ramses Djidjou-Demasse, Yannis Michalakis, Marc Choisy, Micea T Sofonea, and Samuel Alizon, Optimal covid-19 epidemic control until vaccine deployment, medRxiv (2020).
[10] Wendell H Fleming and Raymond W Rishel, Deterministic and stochastic optimal control, vol. 1, Springer Science & Business Media, 2012.
[11] Etienne Kouokam, Jean-Daniel Zucker, Franklin Fondjo, and Marc Choisy, Disease control in age structure population, International Scholarly Research Notices 2013 (2013).
[12] Anuj Kumar and Prashant K Srivastava, Role of optimal screening and treatment on infectious diseases dynamics in presence of self-protection of susceptible, Differential Equations and Dynamical Systems (2019), 1–29.
[13] Sunmi Lee, Gerardo Chowell, and Carlos Castillo-Chávez, Optimal control for pandemic influenza: the role of limited antiviral treatment and isolation, Journal of Theoretical Biology 265 (2010), no. 2, 136–150.
[14] Daniel Liberzon, Calculus of variations and optimal control theory: a concise introduction, Princeton University Press, 2011.
[15] Evgeny Makarov and Bas Spitters, The picard algorithm for ordinary differential equations in coq, International Conference on Interactive Theorem Proving, Springer, 2013, pp. 463–468.
[16] P Pongsumpun and IM Tang, Transmission of dengue hemorrhagic fever in an age structured population, Mathematical and Computer Modelling 37 (2003), no. 9-10, 949–961.
[17] Piu Samui, Jayanta Mondal, and Subhas Khajanchi, A mathematical model for covid-19 transmission dynamics with a case study of india, Chaos, Solitons & Fractals 140 (2020), 110173.
[18] Yung-Fang Tu, Chian-Shiu Chien, Aliaksandr A Yarmishyn, Yi-Ying Lin, Yung-Hung Luo, Yi-Tsung Lin, Wei-Yi Lai, De-Ming Yang, Shih-Jie Chou, Yi-Ping Yang, et al., A review of sars-cov-2 and the ongoing clinical trials, International Journal of Molecular Sciences 21 (2020), no. 7, 2657.
[19] Muhmmad Zamir, Zahir Shah, Fawad Nadeem, Arif Memood, Hussam Alrabaiah, and Poom Kumam, Non pharmaceutical interventions for optimal control of covid-19, Computer methods and programs in biomedicine 196 (2020), 105642.