# Optimal Allocation of Vaccine and Antiviral Drugs for Influenza Containment over Delayed Multiscale Epidemic Model considering Time-Dependent Transmission Rate

Authors: Zohreh Abbasi; Iman Zamani; Amir Hossein Amiri Mehra; Asier Ibeas; Mohsen Shafieirad
Journal: Computational and Mathematical Methods in Medicine
 (2021)

---

## Abstract

In this study, two types of epidemiological models called “within host” and “between hosts” have been studied. The within-host model represents the innate immune response, and the between-hosts model signifies the SEIR (susceptible, exposed, infected, and recovered) epidemic model. The major contribution of this paper is to break the chain of infectious disease transmission by reducing the number of susceptible and infected people via transferring them to the recovered people group with vaccination and antiviral treatment, respectively. Both transfers are considered with time delay. In the first step, optimal control theory is applied to calculate the optimal final time to control the disease within a host’s body with a cost function. To this end, the vaccination that represents the effort that converts healthy cells into resistant-to-infection cells in the susceptible individual’s body is used as the first control input to vaccinate the susceptible individual against the disease. Moreover, the next control input (antiviral treatment) is applied to eradicate the concentrations of the virus and convert healthy cells into resistant-to-infection cells simultaneously in the infected person’s body to treat the infected individual. The calculated optimal time in the first step is considered as the delay of vaccination and antiviral treatment in the SEIR dynamic model. Using Pontryagin’s maximum principle in the second step, an optimal control strategy is also applied to an SEIR mathematical model with a nonlinear transmission rate and time delay, which is computed as optimal time in the first step. Numerical results are consistent with the analytical ones and corroborate our theoretical results.


---

## Body

## 1. Introduction

Seasonal and pandemic influenza A virus (IAV) infection causes severe morbidity and mortality worldwide [1]. There are many drugs to treat influenza, but the immune system is naturally the first defensive line against the disease. Therefore, it is essential to understand the mechanisms that trigger the immune system response and how this response affects the overall spread of the disease. The innate immune response plays an important role in the control and clearance of pathogens following viral infection inside an individual’s body. However, in the majority of virus-infected individuals, the innate immune response is insufficient because viruses use different evasion strategies to escape the immune response. Innate immunity plays a critical role in the control of viral infection because most infectious pathogens are eliminated through innate immune response without necessarily requiring the activation of adaptive immunity [2].Since innate immunity is an important part of the body for dealing with viral infections, some methods have been determined to model different aspects of the human immune system. Differential equation models are the most popular methods to simulate critical diseases’ immunological and epidemiological dynamics [3]. In addition, the use of mathematical modeling of the innate immune response interpreting experimental results has made a significant contribution to this field. Many mathematical models of the immune response to IAV infection have been introduced that considered both different aspects of the immune system and different detail levels. The first mathematical model to illustrate the within-host dynamic of IAV in infected mice was developed in 1976 by Larson et al. [4]. The proposed compartmental model includes seven compartments with five associated rate parameters. Also, the virus population dynamic was introduced in mathematical terms. Another model introduced by Bocharov and Romanyukha [5] studied immune responses to viral infections in human infections, such as the influenza A virus, and considered 12 immune populations. Also, Baccam et al. [6] described the viral kinetics of influenza A during infection in the human body, including target cell and the innate interferon response with delayed virus production through infected cells. However, as opposed to previous papers, which did not consider a complete model of the innate immune system, we used a more comprehensive one adapted from [7].However, the disease’s impact within the community and the relationship between people in society also need to be examined. Mathematical models are instrumental in explaining and quantifying epidemic dynamics. Hence, various mathematical models have been proposed in the following for between-hosts epidemic models to provide a general overview of strategies to be employed in the period of influenza virus infection. Kermack and McKendrick created one of the simplest types to explain epidemic dynamics that is called SIR (susceptible, infected, and recovered individuals, respectively). The proposed model was described in three papers in 1927 [8], 1932 [9], and 1933 [10]. Moreover, in [11], another type of epidemic model (SIS) is considered with a nonlinear incidence rate and time delay. In addition, four different control strategies are investigated and compared, including optimal control to decrease the density of infected individuals, increase that of susceptible ones, and reduce relevant costs. Besides these, in [12], a model of the relationship between two types of systems is formulated and analyzed. The first dynamic is the within-host model dynamic and focuses on cellular interactions. The second is the between-hosts model dynamic that focuses on transmission and infection statuses that are governed by the susceptible-exposed-infected-recovered (SEIR) model. Likewise, in this study, motivated by the previous paper, a combined model with different relationships is presented that will be discussed in detail in the following sections. The authors in [13] have used the sliding mode control on the SEIAR dynamic to overcome the uncertainty of parameters. However, they have not considered the optimal time to eradicate the infection in society, and they only refer to the between-hosts dynamic. For this reason, since time is one of the priorities in order to control disease, we use an optimal control strategy to get the optimal time to cure the infected person and recover the susceptible person.Optimal control is a mathematical method derived from the calculus of variations. There are different methods for computing the optimal control for mathematical models [14]. Pontryagin’s maximum principle, first formulated in 1956 by L. S. Pontryagin [15], is to calculate the optimal control for an ordinary differential equation model system with a given constraint. In [16], an optimal control method is applied to a generic model of a pathogenic attack on the innate immune system. The minimization of a quadratic cost function is generated by numerical optimization. Terminal optimal time control is also investigated. For further studies, the readers are referred to [2], and the references therein introduce an optimal control model based on an ordinary differential equation system. Some studies have investigated the optimal control of SIR and SEIR models, but we mention just some of them in the following. The authors in [17] used an optimal control approach to an SIR epidemic model with a time delay to minimize the spread of infected individuals and maximize the number of susceptible and recovered individuals. Likewise, in [18, 19], optimal control is considered to minimize the level of infection at the terminal time for an SIR and SEIR model, respectively, where the incidence rate is an unspecified nonlinear function, while their models only consider one aspect of the impact of the disease that affects the community and the relationship between susceptible, infected, and recovered people. They do not take into account the effects of the disease inside the body, and the relationship between cells is not considered.Epidemic mathematical models can be used for several epidemic diseases like COVID-19 by modifying the parameter and adding or eliminating some states. In this regard, the authors in [20] considered a SEIAR epidemic model for the COVID-19 pandemic. Also, some new SIR-type models are introduced in [21–24]. The authors in [21] have considered “isolated,” “vaccinated,” and “quarantined” people groups due to several different conditions (after and before vaccine development). Moreover, “hospitalized” and “ICU-admitted” people groups are introduced in [22] to express the spread of coronavirus between the healthcare workers. Also, in [23], a new epidemic model (called SIDARTHE) is studied that considers the distinction between diagnosed and nondiagnosed infected individuals with and without symptoms detected acutely symptomatic infected ones. Besides, the discrete-time SEIR epidemic models can be used in such studies instead of continuous ones [25, 26], to name a few. To further study about optimal control with time delay, readers refer to [27–31] and references therein.Optimal control theory can be used to control other diseases like HIV [32] and especially COVID-19. In this regard, the authors in [33] used optimal control theory based on the SQEIAR epidemic model (Susceptible, Quarantined, Exposed, Infected, Asymptomatic, and Recovered) for the eradication of infection considering quarantine and treatment policies in China and Spain for different types of disease (COVID-19, Ebola, and influenza). Also, according to travel between cities that cause an increase in the outbreak, the authors in [33] investigated the impact of travel or immigration and impulsive change of population. These changes lead to generate impulsive epidemic models and impulsive control strategies [34]. Also, the authors in [35] have considered an optimal control strategy in the age-structured SEIRQ model for the COVID-19 epidemic in Brazil, where the quarantine entrance parameters are the control efforts of the system. In the same way, to explore the dynamic of the COVID-19 pandemic in Pakistan, the necessary optimality conditions are introduced in [36] using the optimal control theory in a mathematical model. Some optimal control strategies with application to COVID-19 are also introduced in [37], for further study.This paper investigates the correlation between two modeling aspects: the effect of the disease on society and inside the individual’s body, which is an important issue. The community includes four groups (SEIR), all of which relate to each other in a community. All communications occur only in the community and among individuals. In contrast, there are some interactions among the cells inside an individual’s body. Recovering the susceptible by vaccination takes time that is considered as delay caused by the length of converting the healthy cells into resistant-to-infectious cells. Also, the infected people need to be treated with antiviral treatment to recover, which is also considered as a delay (the length of time to kill the infected and partially infected cells and viruses and convert the healthy cells into resistant-to-infection cells simultaneously).The goal of the present paper will be presented clearly and thoroughly regarding the reviewed papers. This paper introduces a multiscale model in which both innate immune responses (within-host model) and the SEIR dynamic (between-hosts model) are used simultaneously to define a single system. A multiscale-type model has already been presented in [12], with the difference that the proposed model in the present study is more comprehensive and concisely describes the relationships between the virus, cells, and the immune system. The model is also more realistic and closer to the real world. Moreover, a nonlinear sinusoidal transmission rate (contact rate) is considered, which varies in time according to changes in contacts between people throughout day and night. In this study, the optimal control is used to eradicate the infection at the optimal time in the community. As discussed in Introduction, different papers have been published on the subject of optimal control. However, to the best of the authors’ knowledge, an optimal control approach that simultaneously tackles both types of dynamics (within-host and between-hosts) for epidemic models has never been considered in the literature before.1This paper is organized as follows. Section2 describes the preliminaries. After that, we propose an innate immune response dynamic in Section 3. The optimal control of innate immune response is introduced in Section 4, followed by the SEIR dynamic model with delay and nonlinear transmission rate in Section 5. In Section 6, we apply optimal control of the SEIR epidemic model. The results of simulations of different cases and sensitivity analysis are presented in Section 7. Finally, Section 8 summarizes the conclusions.


## 2. Preliminaries

In this section, we first look at the relationship between the two models considered in this paper. According to Figures1 and 2, our multiscale model considers two scales: the first is the interaction of the virus and innate immune response within a host (inside the body), and the second is the contact between susceptible, exposed, infected, and recovered individuals in society. The dynamic of the innate immune response and relationship between cells and other components are now described in detail. The within-host model dynamic is the innate immune response to the influenza virus using a mathematical model. This model is based on interferon-induced resistance to infection of respiratory epithelial cells and the clearance of infected cells by natural killers with seven state variables (the numbers of healthy (UH), partially infected (UE), infected (UI), and resistant-to-infection cells (UR) and IFN-I molecules ([IFN]), natural killer cells K, and virus particles V). Also, a nonlinear SEIR epidemiological model for influenza is given where the nonnegative state variables St,Et,It, and Rt are the susceptible, exposed, infected (symptomatic), and recovered compartments, respectively. Here, St represents the number of people that are susceptible but not infected with influenza, Et denotes the number of people exposed to influenza (infected but not yet infectious), It denotes the population of infected humans with infectious influenza symptoms, Rt is the number of recovered people from influenza, and N is the total population size. The relationship between the two epidemic models is shown in Figures 1 and 2. The dynamics of the innate immune response and SEIR epidemiological model is described further in the following sections.Figure 1
Block diagram of the relationship between innate immune response and SEIR model.Figure 2
The disease modeling with two scales: within a host and between hosts.The general strategy for this study is described by the following Algorithm1.Algorithm 1: General strategy algorithm.
Step1. Optimal control design (vaccination strategy) for within-host model to convert healthy cells to resistant-to-infection cells in the optimal timet1.Step2. Optimal control design (antiviral treatment strategy) for within-host model to convert healthy cells to resistant-to-infection cells at the same time as destroying the virus particles in the optimal timet2.Step3. Optimal control design for between-hosts model to eradicate the infection, considering delayst1,t2 in the community.


## 3. The Innate Immune Response Dynamic Model

In this section, we consider innate immune system dynamics for dealing with the influenza virus, which is derived from [7]. There are seven components to the dynamic states, and the following equations govern the dynamics. For brevity, we suppress the notation “t” for variables at the current time.
(1)U̇H=SH−kIUHV−kRUHIFN−δHUH,(2)U̇E=kIUHV−kEUE−qKUEK,(3)U̇I=kEUE−δIUI−qKUIK,(4)U̇R=kRUHIFN−δRUR,(5)V̇=ρVUI−δVV,(6)IFṄ=aIUI−δIFNIFN,(7)K̇=SK+ΦKUI−δKK,where UH is the number of healthy cells that increase with cell growth rate SH. These cells become partially infected cells UE via the virus particles V with rate kI. They also convert to resistant-to-infection cells UR by the rate kR with IFN-I molecules IFN that increase by the rate of aI. Partially infected cells are in the eclipse phase. This means that they cannot spread the virus in the body and they infect the other healthy cells by the term kE. When the partially infected cells pass through the eclipse phase, they will become infected cells UI and can spread the virus and cause other healthy cells to become infected. Partially infected UE and infected cells UI are destroyed by natural killers K at a rate of qK. Once the cells are infected UI, they broadcast thousands of similar cells a day, which is denoted by the parameter ρV. The natural death rates of healthy, infected, infection-resistant cells, and IFN-I molecules are denoted byδH,δI,δR,andδIFN, respectively, and the clear rate of virus particles V is δV. The number of NKs (natural killers) increases with a constant cell growth rate SK. In the same way, NKs die at the rate of δK. After inflammatory stimulation caused by infected cells, natural killers are recruited from the blood indicated by ΦKUI. Figure 3 shows a schematic diagram of the innate immune system response to the influenza virus, which shows the cells’ relation very well.Figure 3
Conceptual flow diagram of innate immune dynamics with vaccination.After introducing the model, we now introduce the controller and investigate how to influence the vaccine and antiviral treatment on the cells. As mentioned earlier, the vaccination is used to make healthy cells resistant to infection in a susceptible person’s body and by using the antiviral treatment, the virus particles, infected, and partially infected cells are removed inside the infected body as discussed in detail in the next section.


### 3.1. Equilibrium Points and the Basic Reproduction Number

To determine the disease spread and the number of secondary infections, an infected individual would produce in his/her infectious period the most important threshold criterion named “the basic reproduction number”R0 can be calculated as the follows.According to the definition ofR0=ρFV−1, where
(8)F=00kISHδH000000,V=kE+qKSKδK00−kEδI+qKSKδK00−ρVδV,R0 is calculated as R0=δHSHkIkEρVδK2/δVkEδK+qKSKδIδK+qKSK.The disease-free equilibrium pointDFE is calculated by equating Equations (1)–(7) to zero and solving for UE∗=UI∗=V∗=IFN∗=0 as SH/δH,0,0,0,0,0,0,SK/δK. To find the endemic equilibrium point EE, Equations (1)–(7) are equated to zero, UH=SH/kIV+kRIFN+δH, UE=kIUHV/kE+qKK, UI=kEUE/δI+qKK, UR=kRUHIFN/δR, V=ρVUI/δV, and IFN=aIUI/δIFN. Therefore, by replacing the states, the following equations can be calculated: U¯H=SH/a1U¯I+δH, U¯E=a2/a3δH+a5U¯I−a1a4U¯I2U¯I, U¯R=kRSHδVaI/a7U¯I−δHδVδRδIFNU¯I, V¯=ρV/δVU¯I, IFN¯=aI/δIFNU¯I, K¯=SK+ΦKU¯I/δK, and α1U¯I4+α2U¯I3+α3U¯I2+α4U¯I=0. The last equation has four solutions that one of them is DFE and can be achieved by U¯I=0. The other three solutions are obtained using the parameter values in Table 1, one of which is 1.84×105, and the other two solutions are negative and unacceptable, because the number of cells is a positive value. Therefore, the other equilibrium points are given as EE=U¯H,U¯E,U¯I,U¯R,V¯,IFN¯,K¯=SH/a1U¯I+δH,a2/a3δH+a5U¯I−a1a4U¯I2U¯I,U¯I,kRSHδVaI/a7U¯I−δHδVδRδIFNU¯I,ρV/δVU¯I,aI/δIFNU¯I,SK+ΦKU¯I/δK, in which, α1=qK2ΦK2/δK2δIFNδIFNkIρV+kRδIδV, α2=qKΦK/δK2δIFN2δIFNkIρVqKSK+δKδIFNkIρVkE+δI+2qKkRSKδVaI+kRδVδKaIkE+δI+δHδVδIFNqKΦK, α3=1/δK2δIFNkIρVδK2δIFNkEδI+δKδIFNqKSKkE+δI+qK2δIFNSK2+kRaIδK2δVkEδI+δKδVqKSKkE+δI+qK2δVSK2+δHδVqKΦK2δIFNqKSK+δKδIFNkE+δI, and α4=1/δK2kEδIδHδVδK2+δKδHδVqKSKkE+δI+δHqK2δVSK2−kEkIρVSH.Table 1
Equilibrium values of the model.
State variableDFEEEU¯HSHδH2.44×106U¯E08.94×105U¯I01.84×105U¯R02.37×108V¯0354.65IFN¯00.138K¯SKδK804.610×103Then, the Jacobian of Equations (1)–(7) is given by
(9)J=−kIV−kRIFN−δH000−kIUH−kRUH0kIV−kE−qKK00−kIUH0−qKUE0kE−δI−qKK000−qKUIkRIFN00−δR0kRUH000ρV0−δV0000aI00−δIFN000ΦK000−δK.Theorem 1.
The DFE is locally stable ifR0<1 while it is unstable if R0>1.Proof.
By calculating the seventh-order characteristic equation ofJ about its DFE and according to the Routh-Hurwitz criterion, the following inequality must be satisfied:
(10)δVkEδK+qKSKδIδK+qKSK−δHSHkIkEρVδK2>0.
Therefore, it can be concluded thatR0<1. For R0>1, there exists a positive eigenvalue for the characteristic equation and the equilibrium point is unstable.Theorem 2.
The EE is locally stable ifR0>1 and unstable if R0<1.Proof.
In the same way and by calculating the characteristic equation ofJ about its EE and according to the negative eigenvalues obtained from Jacobian, the EE is stable if R0>1, otherwise it is unstable.


### 3.2. Sensitivity Analysis

In this subsection, the sensitivity ofR0 is investigated. Therefore, to check the R0 sensitivity, SαR0=∂R0/∂αα/R0, we have
(11)SδHR0=1,SSHR0=1,SkIR0=1,SρVR0=1,SkER0=qKSKkEδK+qKSK,SδKR0=qKSKkEδK+δIδK+2SKqKkEδK+qKSKδKδI+qKSK,SδVR0=−δVkEδK+qKSKδIδK+qKSK,SqKR0=−qKkEδKSKδV+SKδIδKδV+2δVSK2qKδVkEδK+qKSKδIδK+qKSK,SSKR0=−SKkEδKqKδV+qKδIδKδV+2δVSKqK2δVkEδK+qKSKδIδK+qKSK.It can be concluded thatR0 is directly related to the change of parameters δH, SH, kI,ρV, kE, and δK. That is, if the parameters δH, SH, kI,andρVincrease by 10%, R0 also changes accordingly (10%) and is inversely related to the parameters δV,qK, and SK.


### 3.3. Positivity Analysis

In this subsection, the positivity of solutions is investigated with Theorem3.Theorem 3.
IfUH0≥0,UE0≥0,UI0≥0,UR0≥0,V0≥0,IFN0≥0, and K0≥0, then solutions of the model are positive for all t≥0.Proof.
According to the positivity of all the parameters of the model, Equation (4) can be written as U̇Rt≥−δRURt. Multiplying e∫0tδRdr into two both side gives U̇Rte∫0tδRdr+δRURte∫0tδRdr≥0; then, d/dte∫0tδRdrURt≥0; integrating from 0 to t gives URt≥UR0e−∫0tδRdr≥0. Therefore URt is positive URt≥0. In the same way, the other states Vt,IFNt,andKt (Equations (5)–(7)) are also positive states and can be written as Vt≥V0e−∫0tδVdr≥0,IFNt≥IFN0e−∫0tδIFNdr≥0, and Kt≥K0e−∫0tδKdr≥0. According to the positivity of URt,Vt,IFNt and Kt and according to Equation (3), UIt≥UI0e−∫0tδI+qKKdr≥0; therefore, UIt≥0. Similarly, UEt≥UE0e−∫0tkE+qKKdr≥0 and UHt≥UH0e−∫0t−kIV−kRIFN−δHdr≥0. Therefore, all solutions are positive.


### 3.4. Existence of Solutions

In this subsection, the existence of solutions is investigated with Theorem4.Theorem 4.
The proposed model with the initial conditionsUH0≥0,UE0≥0,UI0≥0,UR0≥0,V0≥0,IFN0≥0,K0≥0, and E0≥0 has a unique solution.Proof.
The model can be rewritten as(12)fX=AX+f1Xwhere X=UHtUEtUItURtVtIFNtKtT and A=diag −δH,−kE,−δI,−δR,−δV,−δIFN,−δKand f1X=SH−kIUHV−kRUHIFNkIUHV−qKUEKkEUE−qKUIKkRUHIFNρVUIaIUISK+ΦKUIT.
The functionf1X satisfies
(13)f1X1−f1X2=kIUH2V2−UH1V1+kIUH1V1−UH2V2+kRUH2IFN2−UH1IFN1+qKUE2K2−UE1K1+kE−UE2+UE1+qKUI2K2−UI1K1+kRUH1IFN1−UH2IFN2+ρV+aI+ΦKUI1−UI2
Assumed that all parameters have the maximum valueM which is a bounded value:
(14)f1X1−f1X2≤MUH1V1−UH2V2+MUH1IFN1−UH2IFN2+MUE1K1−UE2K2+MUI1K1−UI2K2+MUE1−UE2+MUI1−UI2.
ReplacingUH1V1−UH2V2 by UH1V1−UH2V1+UH2V1−UH2V2 and similarly for others,
(15)f1X1−f1X2≤MUH1V1−UH2V1+UH2V1−UH2V2+MUH1IFN1−UH2IFN1+UH2IFN1−UH2IFN2+MUE1K1−UE2K1+UE2K1−UE2K2+MUI1K1−UI2K1+UI2K1−UI2K2+MUE1−UE2+MUI1−UI2
Therefore,(16)f1X1−f1X2≤MV1UH1−UH2+MUH2V1−V2+MIFN1UH1−UH2+MUH2IFN1−IFN2+MK1UE1−UE2+MUE2K1−K2+MK1UI1−UI2+MUI2K1−K2+MUE1−UE2+MUI1−UI2.
The maximum value of states is equal toR (bounded value); therefore,
(17)f1X1−f1X2≤2RMUH1−UH2+RMV1−V2+RMIFN1−IFN2+2RMUI1−UI2+2RMK1−K2+2RMUE1−UE2≤2RMX1−X2.
It can be concluded thatf1X1−f1X2≤2RMX1−X2. Therefore, fX1−fX2≤LX1−X2, where L=max2RM,A. Thus, f is uniformly Lipschitz continuous, and the solutions exist.


## 3.1. Equilibrium Points and the Basic Reproduction Number

To determine the disease spread and the number of secondary infections, an infected individual would produce in his/her infectious period the most important threshold criterion named “the basic reproduction number”R0 can be calculated as the follows.According to the definition ofR0=ρFV−1, where
(8)F=00kISHδH000000,V=kE+qKSKδK00−kEδI+qKSKδK00−ρVδV,R0 is calculated as R0=δHSHkIkEρVδK2/δVkEδK+qKSKδIδK+qKSK.The disease-free equilibrium pointDFE is calculated by equating Equations (1)–(7) to zero and solving for UE∗=UI∗=V∗=IFN∗=0 as SH/δH,0,0,0,0,0,0,SK/δK. To find the endemic equilibrium point EE, Equations (1)–(7) are equated to zero, UH=SH/kIV+kRIFN+δH, UE=kIUHV/kE+qKK, UI=kEUE/δI+qKK, UR=kRUHIFN/δR, V=ρVUI/δV, and IFN=aIUI/δIFN. Therefore, by replacing the states, the following equations can be calculated: U¯H=SH/a1U¯I+δH, U¯E=a2/a3δH+a5U¯I−a1a4U¯I2U¯I, U¯R=kRSHδVaI/a7U¯I−δHδVδRδIFNU¯I, V¯=ρV/δVU¯I, IFN¯=aI/δIFNU¯I, K¯=SK+ΦKU¯I/δK, and α1U¯I4+α2U¯I3+α3U¯I2+α4U¯I=0. The last equation has four solutions that one of them is DFE and can be achieved by U¯I=0. The other three solutions are obtained using the parameter values in Table 1, one of which is 1.84×105, and the other two solutions are negative and unacceptable, because the number of cells is a positive value. Therefore, the other equilibrium points are given as EE=U¯H,U¯E,U¯I,U¯R,V¯,IFN¯,K¯=SH/a1U¯I+δH,a2/a3δH+a5U¯I−a1a4U¯I2U¯I,U¯I,kRSHδVaI/a7U¯I−δHδVδRδIFNU¯I,ρV/δVU¯I,aI/δIFNU¯I,SK+ΦKU¯I/δK, in which, α1=qK2ΦK2/δK2δIFNδIFNkIρV+kRδIδV, α2=qKΦK/δK2δIFN2δIFNkIρVqKSK+δKδIFNkIρVkE+δI+2qKkRSKδVaI+kRδVδKaIkE+δI+δHδVδIFNqKΦK, α3=1/δK2δIFNkIρVδK2δIFNkEδI+δKδIFNqKSKkE+δI+qK2δIFNSK2+kRaIδK2δVkEδI+δKδVqKSKkE+δI+qK2δVSK2+δHδVqKΦK2δIFNqKSK+δKδIFNkE+δI, and α4=1/δK2kEδIδHδVδK2+δKδHδVqKSKkE+δI+δHqK2δVSK2−kEkIρVSH.Table 1
Equilibrium values of the model.
State variableDFEEEU¯HSHδH2.44×106U¯E08.94×105U¯I01.84×105U¯R02.37×108V¯0354.65IFN¯00.138K¯SKδK804.610×103Then, the Jacobian of Equations (1)–(7) is given by
(9)J=−kIV−kRIFN−δH000−kIUH−kRUH0kIV−kE−qKK00−kIUH0−qKUE0kE−δI−qKK000−qKUIkRIFN00−δR0kRUH000ρV0−δV0000aI00−δIFN000ΦK000−δK.Theorem 1.
The DFE is locally stable ifR0<1 while it is unstable if R0>1.Proof.
By calculating the seventh-order characteristic equation ofJ about its DFE and according to the Routh-Hurwitz criterion, the following inequality must be satisfied:
(10)δVkEδK+qKSKδIδK+qKSK−δHSHkIkEρVδK2>0.
Therefore, it can be concluded thatR0<1. For R0>1, there exists a positive eigenvalue for the characteristic equation and the equilibrium point is unstable.Theorem 2.
The EE is locally stable ifR0>1 and unstable if R0<1.Proof.
In the same way and by calculating the characteristic equation ofJ about its EE and according to the negative eigenvalues obtained from Jacobian, the EE is stable if R0>1, otherwise it is unstable.


## 3.2. Sensitivity Analysis

In this subsection, the sensitivity ofR0 is investigated. Therefore, to check the R0 sensitivity, SαR0=∂R0/∂αα/R0, we have
(11)SδHR0=1,SSHR0=1,SkIR0=1,SρVR0=1,SkER0=qKSKkEδK+qKSK,SδKR0=qKSKkEδK+δIδK+2SKqKkEδK+qKSKδKδI+qKSK,SδVR0=−δVkEδK+qKSKδIδK+qKSK,SqKR0=−qKkEδKSKδV+SKδIδKδV+2δVSK2qKδVkEδK+qKSKδIδK+qKSK,SSKR0=−SKkEδKqKδV+qKδIδKδV+2δVSKqK2δVkEδK+qKSKδIδK+qKSK.It can be concluded thatR0 is directly related to the change of parameters δH, SH, kI,ρV, kE, and δK. That is, if the parameters δH, SH, kI,andρVincrease by 10%, R0 also changes accordingly (10%) and is inversely related to the parameters δV,qK, and SK.


## 3.3. Positivity Analysis

In this subsection, the positivity of solutions is investigated with Theorem3.Theorem 3.
IfUH0≥0,UE0≥0,UI0≥0,UR0≥0,V0≥0,IFN0≥0, and K0≥0, then solutions of the model are positive for all t≥0.Proof.
According to the positivity of all the parameters of the model, Equation (4) can be written as U̇Rt≥−δRURt. Multiplying e∫0tδRdr into two both side gives U̇Rte∫0tδRdr+δRURte∫0tδRdr≥0; then, d/dte∫0tδRdrURt≥0; integrating from 0 to t gives URt≥UR0e−∫0tδRdr≥0. Therefore URt is positive URt≥0. In the same way, the other states Vt,IFNt,andKt (Equations (5)–(7)) are also positive states and can be written as Vt≥V0e−∫0tδVdr≥0,IFNt≥IFN0e−∫0tδIFNdr≥0, and Kt≥K0e−∫0tδKdr≥0. According to the positivity of URt,Vt,IFNt and Kt and according to Equation (3), UIt≥UI0e−∫0tδI+qKKdr≥0; therefore, UIt≥0. Similarly, UEt≥UE0e−∫0tkE+qKKdr≥0 and UHt≥UH0e−∫0t−kIV−kRIFN−δHdr≥0. Therefore, all solutions are positive.


## 3.4. Existence of Solutions

In this subsection, the existence of solutions is investigated with Theorem4.Theorem 4.
The proposed model with the initial conditionsUH0≥0,UE0≥0,UI0≥0,UR0≥0,V0≥0,IFN0≥0,K0≥0, and E0≥0 has a unique solution.Proof.
The model can be rewritten as(12)fX=AX+f1Xwhere X=UHtUEtUItURtVtIFNtKtT and A=diag −δH,−kE,−δI,−δR,−δV,−δIFN,−δKand f1X=SH−kIUHV−kRUHIFNkIUHV−qKUEKkEUE−qKUIKkRUHIFNρVUIaIUISK+ΦKUIT.
The functionf1X satisfies
(13)f1X1−f1X2=kIUH2V2−UH1V1+kIUH1V1−UH2V2+kRUH2IFN2−UH1IFN1+qKUE2K2−UE1K1+kE−UE2+UE1+qKUI2K2−UI1K1+kRUH1IFN1−UH2IFN2+ρV+aI+ΦKUI1−UI2
Assumed that all parameters have the maximum valueM which is a bounded value:
(14)f1X1−f1X2≤MUH1V1−UH2V2+MUH1IFN1−UH2IFN2+MUE1K1−UE2K2+MUI1K1−UI2K2+MUE1−UE2+MUI1−UI2.
ReplacingUH1V1−UH2V2 by UH1V1−UH2V1+UH2V1−UH2V2 and similarly for others,
(15)f1X1−f1X2≤MUH1V1−UH2V1+UH2V1−UH2V2+MUH1IFN1−UH2IFN1+UH2IFN1−UH2IFN2+MUE1K1−UE2K1+UE2K1−UE2K2+MUI1K1−UI2K1+UI2K1−UI2K2+MUE1−UE2+MUI1−UI2
Therefore,(16)f1X1−f1X2≤MV1UH1−UH2+MUH2V1−V2+MIFN1UH1−UH2+MUH2IFN1−IFN2+MK1UE1−UE2+MUE2K1−K2+MK1UI1−UI2+MUI2K1−K2+MUE1−UE2+MUI1−UI2.
The maximum value of states is equal toR (bounded value); therefore,
(17)f1X1−f1X2≤2RMUH1−UH2+RMV1−V2+RMIFN1−IFN2+2RMUI1−UI2+2RMK1−K2+2RMUE1−UE2≤2RMX1−X2.
It can be concluded thatf1X1−f1X2≤2RMX1−X2. Therefore, fX1−fX2≤LX1−X2, where L=max2RM,A. Thus, f is uniformly Lipschitz continuous, and the solutions exist.


## 4. Optimal Control of Within-Host Epidemic Model

Optimal control techniques are used to improve optimal strategies to control numerous kinds of diseases. The flexibility and relative simplicity of optimal control techniques can lead to the development of the optimal strategies to control the monitoring and treating of various kinds of diseases [14]. Therefore, the optimal control theory is applied to the innate immune dynamic in this part of the paper. In the first subsection, the optimal control is used to vaccinate susceptible people and cure infected people.


### 4.1. Vaccination Strategy

In this subsection, the optimal control is applied to the innate immune response to convert healthy cells to resistant-to-infection cells. The vaccine works with the innate immune to develop protection (immunity) to the disease. This subsection is aimed at minimizing the number of healthy cellsUH by converting them to resistant-to-infection cells UR with the control input (vaccination) at a rate of ν in the optimal time. It is obvious that injecting a vaccine into the susceptible person’s bodies takes time to vaccinate them (recovered) against the disease. Unlike [13, 19], in the current paper, this time is considered as a delay in moving susceptible people to the recovered people group in the society, which is investigated in Section 5 in detail. As long as the healthy cells have not become resistant-to-infection cells in the susceptible body, the person has not been vaccinated against the disease, and there is still a risk of infection. The above description is given by the following model:
(18)U̇H=SH−kIUHV−kRUHIFN−δHUH−νUH,(19)U̇E=kIUHV−kEUE−qKUEK,(20)U̇I=kEUE−δIUI−qKUIK,(21)U̇R=kRUHIFN−δRUR+νUH,(22)V̇=ρVUI−δVV,(23)IFṄ=aIUI−δIFNIFN,(24)K̇=SK+ΦKUI−δKK.Figure3 shows the diagram of the within-host dynamic model with vaccination. The susceptible individual’s body does not have the virus, infected, or partially infected cells. Therefore, the goal is to convert healthy cells into resistant-to-infection cells using the vaccination strategy in the optimal time.Therefore, the optimal control problem is to minimize the cost (objective) function applied to Equations (18)–(24) given by
(25)Jν,tf=∫0tfA1UH−A2UR+A32ν2dt+∅tf,in which a free terminal time is investigated, which gives a minimum duration of a vaccination. tf represent the duration of the vaccination. ∅ is a positive increasing function such that limt⟶∞∅t=+∞. In other words, the vaccination strategy is applied using the optimal control ν∗ in optimal time tf∗ such that
(26)Jν∗,tf∗=minJν,tf∣ν∈U1,tf∈ℝ+,where ν is control and U1 is the set of admissible controls defined by
(27)U1=U1ismeasurable,0≤ν≤νmax=1,t∈0,tf.The necessary conditions for optimality are expressed by(28)∇Xf∅−ptfTδXf+HXtf,utf,ptf,tf+∂∅∂tfδtf=0,where Xtf or Xf is vector of constant values (the final value of all states is known). The gradient ∇Xf∅ is the derivative of ∅tf. tf is unknown but Xtf is defined; therefore, δXf=0. Also, a stationarity condition on the control is as follows:
(29)∇νH=0,which are expressed using the Hamiltonian of the system. The Hamiltonian for the control problem is defined as follows:
(30)H=g+pTf,where f=ẋ=U̇H,U̇E,U̇I,U̇R,V̇,IFṄ,K̇T and g=A1UH−A2UR+A3/2ν2.Therefore,(31)H=A1UH−A2UR+A32ν2+p1U̇H+p2U̇E+p3U̇I+p4U̇R+p5V̇+p6IFṄ+p7K,̇and p1,p2,p3,p4,p5,p6, and p7 are the adjoint functions to be determined suitably. Then, the necessary conditions are computed as follows:
(32)ẋ=∂H∂p=U̇H,U̇E,U̇I,U̇R,V̇,IFṄ,K̇T.Then, the necessary conditions for optimality are computed, which are expressed by the mentioned Hamiltonian equations [38]. These equations include a linear, ordinary differential equation that are shown by
(33)ṗ=−∂H∂x=−∂H∂UH,−∂H∂UE,−∂H∂UI,−∂H∂UR,−∂H∂V,−∂H∂IFN,−∂H∂KT.In order to derive the necessary conditions for this optimal control, the Pontryagin maximum principle [39] is applied to characterize the optimal control problem and the optimal final time, which is given as the following theorem.Theorem 5.
Given an optimal controlν∗, solutions UH∗,UE∗,UI∗,UR∗,V∗,IFN∗, and K∗ for the optimal control problem Jν∗,tf∗=minJν,tf∣ν∈U1,tf∈ℝ+. Then, there are adjoint variables p1,p2,p3,p4,p5,p6, and p7 that satisfy
(34)ṗ1=−A1+p1δH−kIVp2−p1−kRIFN+νp4−p1,(35)ṗ2=p2kE+qKK−p3kE,(36)ṗ3=p3δI+qKK−p5ρV−p6aI−p7ΦK,(37)ṗ4=A2+p4δR,(38)ṗ5=kIUHp2+p1+p5ρV,(39)ṗ6=kRUHp1−p4+p6δIFN,(40)ṗ7=p2qKUE+p3qKUI+p7δK,with transversality conditions as
(41)pitf=0,i=1,2,3,4,5,6,7.
As a result, an optimal controlν is given by
(42)ν∗=maxminUH∗p1−p4A3,1,0,and the optimal final time is given by
(43)∂ϕ∂tftf∗=−A1UH∗+A2UR∗−A32ν∗2.Proof.
Using Pontryagin’s principle, the adjoint equations and transversality conditions are obtained such that(44)ṗ1=−∂H∂UH,p1tf=0,(45)ṗ2=−∂H∂UE,p2tf=0,(46)ṗ3=−∂H∂UI,p3tf=0,(47)ṗ4=−∂H∂UR,p4tf=0,(48)ṗ5=−∂H∂V,p5tf=0,(49)ṗ6=−∂H∂IFN,p6tf=0,(50)ṗ7=−∂H∂K,p7tf=0,
and using the optimality conditions, the optimal control is solved as∇νH=0 in the following form:
(51)A3ν+p1∂U̇H∂ν+p2∂U̇E∂ν+p3∂U̇I∂ν+p4∂U̇R∂ν+p5∂V̇∂ν+p6∂IFṄ∂ν+p7∂K̇∂ν=A3ν+p1−UH+p4UH=0.
Thus, the control law is written as(52)ν∗=UH∗p1−p4A3.Using the property of the control space, the following equation is obtained:
(53)ν∗=1UH∗p1−p4A3≥1,UH∗p1−p4A30<UH∗p1−p4A3<1,0UH∗p1−p4A3≤0,so, the optimal control is characterized as
(54)ν∗=maxminUH∗p1−p4A3,1,0.
Now, according to Equation (28) and knowing δXf=0, the optimal final time tf∗ is computed using the following equation:
(55)∂ϕ∂tftf∗=−HUH∗,UE∗,UI∗,UR∗,V∗,IFN∗,K∗,ν∗,tf∗.
Therefore, using the characterization of the optimal control, the following optimality system is presented to obtain the optimal control and the state. This optimal control problem consists of the state system Equations (18)–(24) with initial conditions at t=0, the adjoint system Equations (34)–(40) with the final conditions Equation (41), and the characterization of the optimal control Equation (42).
(56)U̇H∗=SH−kIUH∗V∗−kRUH∗IFN∗−δHUH∗−maxminUH∗p1−p4A3,1,0UH∗,U̇E∗=kIUH∗V∗−kEUE∗−qKUE∗K∗,U̇I∗=kEUE∗−δIUI∗−qKUI∗K∗,U̇R∗=kRUHIFN∗−δRUR∗+maxminUH∗p1−p4A3,1,0UH∗,V̇∗=ρVUI∗−δVV∗,IFṄ∗=aIUI∗−δIFNIFN∗,K̇∗=SK+ΦKUI∗−δKK∗,and the adjoint system is as follows:
(57)ṗ1=−A1+p1δH−kIV∗p2−p1−kRIFN∗+maxminUH∗p1−p4A3,1,0p4−p1,ṗ2=p2kE+qKK∗−p3kE,ṗ3=p3δI+qKK∗−p5ρV−p6aI−p7ΦK,ṗ4=A2+p4δR,ṗ5=kIUH∗p2+p1+p5ρV,ṗ6=kRUH∗p1−p4+p6δIFN,ṗ7=p2qKUE∗+p3qKUI∗+p7δK,with p1tf=0, p2tf=0, p3tf=0, p4tf=0, p5tf=0, p6tf=0,p7tf=0, UH0=UH0, UE0=UE0, UI0=UI0, UR0=UR0, V0=V0, IFN0=IFN0, K0=K0, and T∗which can be rewritten as in Equation (78). So, the proof is completed.After examining the effect of the vaccine on the person’s body and the control strategy applied to the model for the best effect, it is necessary to investigate the efficiency of the optimal antiviral treatment strategy in an optimal time, which is considered in the following subsection.


### 4.2. Antiviral Treatment Strategy

This subsection is aimed at destroying the virus particles and converting healthy cells into resistant-to-infection cells simultaneously using the optimal antiviral treatment strategy. There are infected and partially infected cells within the infected individual’s body. Therefore, the control input (antiviral treatment) is used to eradicate the disease by converting healthy cells to resistant-to-infection cells at the rate ofT1 and destroying virus particles at a rate of T2. Figure 4 shows the conceptual flow diagram of the innate immune dynamics with the proposed control action (antiviral treatment). Moreover, the dynamic of the innate immune with antiviral treatment is considered as
(58)U̇H=SH−kIUHV−kRUHIFN−δHUH−T1UH,(59)U̇E=kIUHV−kEUE−qKUEK,(60)U̇I=kEUE−δIUI−qKUIK,(61)U̇R=kRUHIFN−δRUR+T1UH,(62)V̇=ρVUI−δVV−T2V,(63)IFṄ=aIUI−δIFNIFN,(64)K̇=SK+ΦKUI−δKK.Figure 4
The diagram of the innate immune dynamics with the proposed control action (antiviral treatment).To reach the proposed goal and minimize the cost of the optimal antiviral treatment strategy within an optimal duration, the following objective function with free terminal time has to be minimized:(65)JT1,T2,tf=∫0tfA1UH−A2UR+A3V+A42T12+A52T22dt+∅tf,in which tf represent the duration of the antiviral treatment. ∅ is a positive increasing function such that limt⟶∞∅t=+∞. Also, to reach the optimal controls T1∗,T2∗ and an optimal terminal time tf∗, we considered
(66)JT1∗,T2∗,tf∗=minJT1,T2,tf∣T1,T2∈U2,tf∈ℝ+,where U2 is the set of admissible controls defined by
(67)U2=T1,T2∣T1,T2aremeasurable,0≤T1≤1,0≤T2≤1,t∈0,tf.Then, the necessary conditions for optimality are written as follows:(68)ṗ=−∂H∂x,(69)∇Xf∅−ptfTδXf+HXtf,utf,ptf,tf+∂∅∂tfδtf=0,in which tf is unknown while Xtf is fixed; therefore, δXf=0. Then, the Hamiltonian for the control problem is defined as
(70)H=A1UH−A2UR+A3V+A42T12+A52T22+p1U̇H+p2U̇E+p3U̇I+p4U̇R+p5V̇+p6IFṄ+p7K̇,where p1,p2,p3,p4,p5,p6, and p7 are also the adjoint functions to be obtained. As discussed earlier, the following theorem is obtained to apply Pontryagin’s principle to the Hamiltonian.Theorem 6.
Given optimal control pairT1∗,T2∗ and solutions UH∗,UE∗,UI∗,UR∗,V∗,IFN∗, and K∗ for the optimal control problem JT1∗,T2∗,tf∗=minJT1,T2,tf∣T1,T2∈U,tf∈ℝ+, let UH∗,UE∗,UI∗,UR∗,V∗,IFN∗, and K∗ be optimal state solutions with an associated optimal control pair T1,T2. There are adjoint variables that satisfy the following equations:
(71)ṗ1=−A1+p1δH+kIVp1−p2+kRIFN+νp1−p4,ṗ2=p2kE+qKK−p3kE,ṗ3=p3δI+qKK−p5ρV−p6aI−p7ΦK,ṗ4=A2+p4δR,ṗ5=−A3+kIUHp1−p2+p4δv+T2,ṗ6=kRUHp1−p4+p6δIFN,ṗ7=p2qKUE+p3qKUI+p7δK,with transversality conditions:
(72)pitf=0,i=1,2,3,4,5,6,7.
Furthermore, the optimal controlsT1,T2 are given by
(73)T1∗=maxminUH∗p1−p4A4,1,0,T2∗=maxminp5V∗A5,1,0,and the optimal final time is given by
(74)HUH∗,UE∗,UI∗,UR∗,V∗,IFN∗,K∗,T1∗,T2∗,tf∗+∂ϕ∂tftf∗=0.Proof.
The adjoint equation and transversality conditions can be obtained using Pontryagin’s principle as similar to Equations (44)–(50). Then, by using the optimality conditions, the optimal controls are obtained using ∇T1H=0 and ∇T2H=0. As a result, the control inputs are computed as
(75)T1∗=UH∗p1−p4A3,T2∗=p5V∗A5.
Using the property of the control space, we obtain(76)T1∗=1UH∗p1−p4A4≥1,UH∗p1−p4A30<UH∗p1−p4A4<1,0UH∗p1−p4A4≤0,T2∗=1p5V∗A5≥1,p5V∗A50<p5V∗A5<1,0p5V∗A5≤0.
So, the optimal control is characterized as(77)T1∗=maxminUH∗p1−p4A4,1,0,T2∗=maxminp5V∗A5,1,0,and according to Equation (74) and knowing that δXf=0, the optimal final time tf∗ is computed using the following equation:
(78)∂ϕ∂tftf∗=−A1UH∗−A2UR∗+A3V∗+A42T1∗2+A52T2∗2.
Therefore, using the characterization of the optimal control, we have the following optimality system to obtain the optimal states:(79)U̇H∗=SH−kIUH∗V∗−kRUH∗IFN∗−δHUH∗−maxminUH∗p1−p4A4,1,0UH∗,U̇E∗=kIUH∗V∗−kEUE∗−qKUE∗K∗,U̇I∗=kEUE∗−δIUI∗−qKUI∗K∗,U̇R∗=kRUHIFN∗−δRUR∗+maxminUH∗p1−p4A4,1,0UH∗,V̇∗=ρVUI∗−δVV∗−maxminp5V∗A5,1,0V∗,IFṄ∗=aIUI∗−δIFNIFN∗,K̇∗=SK+ΦKUI∗−δKK∗,and the adjoint equations are as follows:
(80)ṗ1=−A1+p1δH−kIV∗p2−p1−kRIFN∗+maxminUH∗p1−p4A4,1,0p4−p1,ṗ2=p2kE+qKK∗−p3kE,ṗ3=p3δI+qKK∗−p5ρV−p6aI−p7ΦK,ṗ4=A2+p4δR,ṗ5=−A3+kIUH∗p1−p2+p4δv+maxminp5V∗A5,1,0,ṗ6=kRUH∗p1−p4+p6δIFN,ṗ7=p2qKUE∗+p3qKUI∗+p7δK,where p1tf=0, p2tf=0, p3tf=0, p4tf=0, p5tf=0, p6tf=0,p7tf=0, UH0=UH0, UE0=UE0, UI0=UI0, UR0=UR0, V0=V0, IFN0=IFN0, and K0=K0. Then, the proof is completed.Remark 7.
In the vaccination strategy, the healthy cells convert to resistant-to-infection cells using the proposed controller and their number converges to zero whent⟶∞UH⟶0 and the rest remain at their initial zero value, while in the antiviral treatment strategy, in addition to converting healthy cells to resistant-to-infection cells, the virus particles decrease and converge to zero using the proposed controller. Meanwhile, the number of natural killers converges to the initial constant value. Then, after a while, Equation (59) becomes U̇E=−kE+qKKUE in which kE+qKK is constant and is considered as first-order linear differential equation that yields UE=UE0e−kE+qKKt. As a result, UE converges exponentially to zero.In this section, the effect of the optimal control on the within-host model is examined. At the same time, it is necessary to explore the controller’s influence on the between-hosts model (SEIR epidemic model). Thus, the impact of the delay computed in this section is discussed in the following section in detail.


## 4.1. Vaccination Strategy

In this subsection, the optimal control is applied to the innate immune response to convert healthy cells to resistant-to-infection cells. The vaccine works with the innate immune to develop protection (immunity) to the disease. This subsection is aimed at minimizing the number of healthy cellsUH by converting them to resistant-to-infection cells UR with the control input (vaccination) at a rate of ν in the optimal time. It is obvious that injecting a vaccine into the susceptible person’s bodies takes time to vaccinate them (recovered) against the disease. Unlike [13, 19], in the current paper, this time is considered as a delay in moving susceptible people to the recovered people group in the society, which is investigated in Section 5 in detail. As long as the healthy cells have not become resistant-to-infection cells in the susceptible body, the person has not been vaccinated against the disease, and there is still a risk of infection. The above description is given by the following model:
(18)U̇H=SH−kIUHV−kRUHIFN−δHUH−νUH,(19)U̇E=kIUHV−kEUE−qKUEK,(20)U̇I=kEUE−δIUI−qKUIK,(21)U̇R=kRUHIFN−δRUR+νUH,(22)V̇=ρVUI−δVV,(23)IFṄ=aIUI−δIFNIFN,(24)K̇=SK+ΦKUI−δKK.Figure3 shows the diagram of the within-host dynamic model with vaccination. The susceptible individual’s body does not have the virus, infected, or partially infected cells. Therefore, the goal is to convert healthy cells into resistant-to-infection cells using the vaccination strategy in the optimal time.Therefore, the optimal control problem is to minimize the cost (objective) function applied to Equations (18)–(24) given by
(25)Jν,tf=∫0tfA1UH−A2UR+A32ν2dt+∅tf,in which a free terminal time is investigated, which gives a minimum duration of a vaccination. tf represent the duration of the vaccination. ∅ is a positive increasing function such that limt⟶∞∅t=+∞. In other words, the vaccination strategy is applied using the optimal control ν∗ in optimal time tf∗ such that
(26)Jν∗,tf∗=minJν,tf∣ν∈U1,tf∈ℝ+,where ν is control and U1 is the set of admissible controls defined by
(27)U1=U1ismeasurable,0≤ν≤νmax=1,t∈0,tf.The necessary conditions for optimality are expressed by(28)∇Xf∅−ptfTδXf+HXtf,utf,ptf,tf+∂∅∂tfδtf=0,where Xtf or Xf is vector of constant values (the final value of all states is known). The gradient ∇Xf∅ is the derivative of ∅tf. tf is unknown but Xtf is defined; therefore, δXf=0. Also, a stationarity condition on the control is as follows:
(29)∇νH=0,which are expressed using the Hamiltonian of the system. The Hamiltonian for the control problem is defined as follows:
(30)H=g+pTf,where f=ẋ=U̇H,U̇E,U̇I,U̇R,V̇,IFṄ,K̇T and g=A1UH−A2UR+A3/2ν2.Therefore,(31)H=A1UH−A2UR+A32ν2+p1U̇H+p2U̇E+p3U̇I+p4U̇R+p5V̇+p6IFṄ+p7K,̇and p1,p2,p3,p4,p5,p6, and p7 are the adjoint functions to be determined suitably. Then, the necessary conditions are computed as follows:
(32)ẋ=∂H∂p=U̇H,U̇E,U̇I,U̇R,V̇,IFṄ,K̇T.Then, the necessary conditions for optimality are computed, which are expressed by the mentioned Hamiltonian equations [38]. These equations include a linear, ordinary differential equation that are shown by
(33)ṗ=−∂H∂x=−∂H∂UH,−∂H∂UE,−∂H∂UI,−∂H∂UR,−∂H∂V,−∂H∂IFN,−∂H∂KT.In order to derive the necessary conditions for this optimal control, the Pontryagin maximum principle [39] is applied to characterize the optimal control problem and the optimal final time, which is given as the following theorem.Theorem 5.
Given an optimal controlν∗, solutions UH∗,UE∗,UI∗,UR∗,V∗,IFN∗, and K∗ for the optimal control problem Jν∗,tf∗=minJν,tf∣ν∈U1,tf∈ℝ+. Then, there are adjoint variables p1,p2,p3,p4,p5,p6, and p7 that satisfy
(34)ṗ1=−A1+p1δH−kIVp2−p1−kRIFN+νp4−p1,(35)ṗ2=p2kE+qKK−p3kE,(36)ṗ3=p3δI+qKK−p5ρV−p6aI−p7ΦK,(37)ṗ4=A2+p4δR,(38)ṗ5=kIUHp2+p1+p5ρV,(39)ṗ6=kRUHp1−p4+p6δIFN,(40)ṗ7=p2qKUE+p3qKUI+p7δK,with transversality conditions as
(41)pitf=0,i=1,2,3,4,5,6,7.
As a result, an optimal controlν is given by
(42)ν∗=maxminUH∗p1−p4A3,1,0,and the optimal final time is given by
(43)∂ϕ∂tftf∗=−A1UH∗+A2UR∗−A32ν∗2.Proof.
Using Pontryagin’s principle, the adjoint equations and transversality conditions are obtained such that(44)ṗ1=−∂H∂UH,p1tf=0,(45)ṗ2=−∂H∂UE,p2tf=0,(46)ṗ3=−∂H∂UI,p3tf=0,(47)ṗ4=−∂H∂UR,p4tf=0,(48)ṗ5=−∂H∂V,p5tf=0,(49)ṗ6=−∂H∂IFN,p6tf=0,(50)ṗ7=−∂H∂K,p7tf=0,
and using the optimality conditions, the optimal control is solved as∇νH=0 in the following form:
(51)A3ν+p1∂U̇H∂ν+p2∂U̇E∂ν+p3∂U̇I∂ν+p4∂U̇R∂ν+p5∂V̇∂ν+p6∂IFṄ∂ν+p7∂K̇∂ν=A3ν+p1−UH+p4UH=0.
Thus, the control law is written as(52)ν∗=UH∗p1−p4A3.Using the property of the control space, the following equation is obtained:
(53)ν∗=1UH∗p1−p4A3≥1,UH∗p1−p4A30<UH∗p1−p4A3<1,0UH∗p1−p4A3≤0,so, the optimal control is characterized as
(54)ν∗=maxminUH∗p1−p4A3,1,0.
Now, according to Equation (28) and knowing δXf=0, the optimal final time tf∗ is computed using the following equation:
(55)∂ϕ∂tftf∗=−HUH∗,UE∗,UI∗,UR∗,V∗,IFN∗,K∗,ν∗,tf∗.
Therefore, using the characterization of the optimal control, the following optimality system is presented to obtain the optimal control and the state. This optimal control problem consists of the state system Equations (18)–(24) with initial conditions at t=0, the adjoint system Equations (34)–(40) with the final conditions Equation (41), and the characterization of the optimal control Equation (42).
(56)U̇H∗=SH−kIUH∗V∗−kRUH∗IFN∗−δHUH∗−maxminUH∗p1−p4A3,1,0UH∗,U̇E∗=kIUH∗V∗−kEUE∗−qKUE∗K∗,U̇I∗=kEUE∗−δIUI∗−qKUI∗K∗,U̇R∗=kRUHIFN∗−δRUR∗+maxminUH∗p1−p4A3,1,0UH∗,V̇∗=ρVUI∗−δVV∗,IFṄ∗=aIUI∗−δIFNIFN∗,K̇∗=SK+ΦKUI∗−δKK∗,and the adjoint system is as follows:
(57)ṗ1=−A1+p1δH−kIV∗p2−p1−kRIFN∗+maxminUH∗p1−p4A3,1,0p4−p1,ṗ2=p2kE+qKK∗−p3kE,ṗ3=p3δI+qKK∗−p5ρV−p6aI−p7ΦK,ṗ4=A2+p4δR,ṗ5=kIUH∗p2+p1+p5ρV,ṗ6=kRUH∗p1−p4+p6δIFN,ṗ7=p2qKUE∗+p3qKUI∗+p7δK,with p1tf=0, p2tf=0, p3tf=0, p4tf=0, p5tf=0, p6tf=0,p7tf=0, UH0=UH0, UE0=UE0, UI0=UI0, UR0=UR0, V0=V0, IFN0=IFN0, K0=K0, and T∗which can be rewritten as in Equation (78). So, the proof is completed.After examining the effect of the vaccine on the person’s body and the control strategy applied to the model for the best effect, it is necessary to investigate the efficiency of the optimal antiviral treatment strategy in an optimal time, which is considered in the following subsection.


## 4.2. Antiviral Treatment Strategy

This subsection is aimed at destroying the virus particles and converting healthy cells into resistant-to-infection cells simultaneously using the optimal antiviral treatment strategy. There are infected and partially infected cells within the infected individual’s body. Therefore, the control input (antiviral treatment) is used to eradicate the disease by converting healthy cells to resistant-to-infection cells at the rate ofT1 and destroying virus particles at a rate of T2. Figure 4 shows the conceptual flow diagram of the innate immune dynamics with the proposed control action (antiviral treatment). Moreover, the dynamic of the innate immune with antiviral treatment is considered as
(58)U̇H=SH−kIUHV−kRUHIFN−δHUH−T1UH,(59)U̇E=kIUHV−kEUE−qKUEK,(60)U̇I=kEUE−δIUI−qKUIK,(61)U̇R=kRUHIFN−δRUR+T1UH,(62)V̇=ρVUI−δVV−T2V,(63)IFṄ=aIUI−δIFNIFN,(64)K̇=SK+ΦKUI−δKK.Figure 4
The diagram of the innate immune dynamics with the proposed control action (antiviral treatment).To reach the proposed goal and minimize the cost of the optimal antiviral treatment strategy within an optimal duration, the following objective function with free terminal time has to be minimized:(65)JT1,T2,tf=∫0tfA1UH−A2UR+A3V+A42T12+A52T22dt+∅tf,in which tf represent the duration of the antiviral treatment. ∅ is a positive increasing function such that limt⟶∞∅t=+∞. Also, to reach the optimal controls T1∗,T2∗ and an optimal terminal time tf∗, we considered
(66)JT1∗,T2∗,tf∗=minJT1,T2,tf∣T1,T2∈U2,tf∈ℝ+,where U2 is the set of admissible controls defined by
(67)U2=T1,T2∣T1,T2aremeasurable,0≤T1≤1,0≤T2≤1,t∈0,tf.Then, the necessary conditions for optimality are written as follows:(68)ṗ=−∂H∂x,(69)∇Xf∅−ptfTδXf+HXtf,utf,ptf,tf+∂∅∂tfδtf=0,in which tf is unknown while Xtf is fixed; therefore, δXf=0. Then, the Hamiltonian for the control problem is defined as
(70)H=A1UH−A2UR+A3V+A42T12+A52T22+p1U̇H+p2U̇E+p3U̇I+p4U̇R+p5V̇+p6IFṄ+p7K̇,where p1,p2,p3,p4,p5,p6, and p7 are also the adjoint functions to be obtained. As discussed earlier, the following theorem is obtained to apply Pontryagin’s principle to the Hamiltonian.Theorem 6.
Given optimal control pairT1∗,T2∗ and solutions UH∗,UE∗,UI∗,UR∗,V∗,IFN∗, and K∗ for the optimal control problem JT1∗,T2∗,tf∗=minJT1,T2,tf∣T1,T2∈U,tf∈ℝ+, let UH∗,UE∗,UI∗,UR∗,V∗,IFN∗, and K∗ be optimal state solutions with an associated optimal control pair T1,T2. There are adjoint variables that satisfy the following equations:
(71)ṗ1=−A1+p1δH+kIVp1−p2+kRIFN+νp1−p4,ṗ2=p2kE+qKK−p3kE,ṗ3=p3δI+qKK−p5ρV−p6aI−p7ΦK,ṗ4=A2+p4δR,ṗ5=−A3+kIUHp1−p2+p4δv+T2,ṗ6=kRUHp1−p4+p6δIFN,ṗ7=p2qKUE+p3qKUI+p7δK,with transversality conditions:
(72)pitf=0,i=1,2,3,4,5,6,7.
Furthermore, the optimal controlsT1,T2 are given by
(73)T1∗=maxminUH∗p1−p4A4,1,0,T2∗=maxminp5V∗A5,1,0,and the optimal final time is given by
(74)HUH∗,UE∗,UI∗,UR∗,V∗,IFN∗,K∗,T1∗,T2∗,tf∗+∂ϕ∂tftf∗=0.Proof.
The adjoint equation and transversality conditions can be obtained using Pontryagin’s principle as similar to Equations (44)–(50). Then, by using the optimality conditions, the optimal controls are obtained using ∇T1H=0 and ∇T2H=0. As a result, the control inputs are computed as
(75)T1∗=UH∗p1−p4A3,T2∗=p5V∗A5.
Using the property of the control space, we obtain(76)T1∗=1UH∗p1−p4A4≥1,UH∗p1−p4A30<UH∗p1−p4A4<1,0UH∗p1−p4A4≤0,T2∗=1p5V∗A5≥1,p5V∗A50<p5V∗A5<1,0p5V∗A5≤0.
So, the optimal control is characterized as(77)T1∗=maxminUH∗p1−p4A4,1,0,T2∗=maxminp5V∗A5,1,0,and according to Equation (74) and knowing that δXf=0, the optimal final time tf∗ is computed using the following equation:
(78)∂ϕ∂tftf∗=−A1UH∗−A2UR∗+A3V∗+A42T1∗2+A52T2∗2.
Therefore, using the characterization of the optimal control, we have the following optimality system to obtain the optimal states:(79)U̇H∗=SH−kIUH∗V∗−kRUH∗IFN∗−δHUH∗−maxminUH∗p1−p4A4,1,0UH∗,U̇E∗=kIUH∗V∗−kEUE∗−qKUE∗K∗,U̇I∗=kEUE∗−δIUI∗−qKUI∗K∗,U̇R∗=kRUHIFN∗−δRUR∗+maxminUH∗p1−p4A4,1,0UH∗,V̇∗=ρVUI∗−δVV∗−maxminp5V∗A5,1,0V∗,IFṄ∗=aIUI∗−δIFNIFN∗,K̇∗=SK+ΦKUI∗−δKK∗,and the adjoint equations are as follows:
(80)ṗ1=−A1+p1δH−kIV∗p2−p1−kRIFN∗+maxminUH∗p1−p4A4,1,0p4−p1,ṗ2=p2kE+qKK∗−p3kE,ṗ3=p3δI+qKK∗−p5ρV−p6aI−p7ΦK,ṗ4=A2+p4δR,ṗ5=−A3+kIUH∗p1−p2+p4δv+maxminp5V∗A5,1,0,ṗ6=kRUH∗p1−p4+p6δIFN,ṗ7=p2qKUE∗+p3qKUI∗+p7δK,where p1tf=0, p2tf=0, p3tf=0, p4tf=0, p5tf=0, p6tf=0,p7tf=0, UH0=UH0, UE0=UE0, UI0=UI0, UR0=UR0, V0=V0, IFN0=IFN0, and K0=K0. Then, the proof is completed.Remark 7.
In the vaccination strategy, the healthy cells convert to resistant-to-infection cells using the proposed controller and their number converges to zero whent⟶∞UH⟶0 and the rest remain at their initial zero value, while in the antiviral treatment strategy, in addition to converting healthy cells to resistant-to-infection cells, the virus particles decrease and converge to zero using the proposed controller. Meanwhile, the number of natural killers converges to the initial constant value. Then, after a while, Equation (59) becomes U̇E=−kE+qKKUE in which kE+qKK is constant and is considered as first-order linear differential equation that yields UE=UE0e−kE+qKKt. As a result, UE converges exponentially to zero.In this section, the effect of the optimal control on the within-host model is examined. At the same time, it is necessary to explore the controller’s influence on the between-hosts model (SEIR epidemic model). Thus, the impact of the delay computed in this section is discussed in the following section in detail.


## 5. SEIR Dynamic Model

In this section, a SEIR epidemic model is introduced with two additional properties. First is two time-delay that aroused from this fact that vaccination and antiviral treatment take time to impact a person’s body. Second is a nonlinear contact rate (transmission rate) function that is followed by this fact that the contact rate is varying throughout day and night. In this regard, the susceptible and infected individuals move to the recovered people group with time delayst1 and t2 (optimal times of vaccination and antiviral treatment within the body, which is computed in the previous section), respectively. The proposed SEIR epidemic model is presented as
(81)Ṡt=−βtStIt+μd31−ρIt−St,Ėt=βtStIt−μd+κEt,İt=κEt−αIt,Ṙt=gIt−μdRt,where α=g−3μdρ+μd. As shown in Figure 5, St represents the number of susceptible individuals who are vulnerable and not infected, Et or exposed individual denotes the number of people that are infected but not infectious yet, It signifies the population of infectious people with symptoms that can transmit the disease to others, and Rt is the number of completely recovered (immune) people. N is the total population size St+Et+It+Rt=N. μdis the rate of deaths from causes unrelated to the infection, and μb is the birth rate in each group, which is triple the mortality rate μb=3μd. t1 is the vaccination time of susceptible people, and t2 is the recovery time of infected people. In the case of births, in the susceptible, exposed, and recovered groups, newborns are susceptible to the disease and add to the susceptible group at the rate of μb. Fraction ρ of the infected individuals’ infants goes to the infected compartment, and the others (fraction 1−ρ of infected individuals’ infants) move to the susceptible group. The exposed people get infected at a rate of κ, and infected people leave their compartment at the rate of α.Figure 5
Conceptual flow diagram of the influenza SEIAR dynamics with the proposed controllers.The contact rate isβt that is purposefully variable throughout day and night. It is also obvious that this rate at night is higher than during the day, which makes it closer to the real world and is more effective in practice [40]. According to the basic definition of βt that is defined as “the average number of contacts between individuals in the population per unit time” [40], the contact rate must be more than one during the day due to the more relations between people in the community. Moreover, late at night and early in the morning, the transmission rate is less than one. Therefore, with trial and error and taking into account the assumptions below, the contact rate is obtained as βt=0.53−sin2πt2−sin2πt−cos2πt.
(a)
00:00 midnight:β00:00=1(b)
06:00 morning:β06:00=0.5(c)
12:00 midday:β12:00=2(d)
18:00 evening:β18:00=1.5Thus, the contact rate is less than one late at night and at daybreak and is greater than one in the middle of the day because there are more social gatherings and offices and schools are also open.In the next section, an optimal control is applied to the proposed model to overcome the disease outbreak.


## 6. Optimal Control of the Between-Hosts Epidemic Model

The aim of this section is twofold: first, minimizing the number of susceptible individuals and moving them to the recovered group by vaccination at a rate ofuv0≤uv≤1 in 25 days and second, moving the infected people to the recovered group at a rate of uτ0≤uτ≤1 by antiviral treatment in 25 days, as shown in Figure 5. The SEIR dynamic model in the presence of the controller is presented as
(82)Ṡt=−βtStIt+μd31−ρIt−St−uvt−t1St,(83)Ėt=βtStIt−μd+κEt,(84)İt=κEt−αIt−uτt−t2It,(85)Ṙt=gIt−μdRt+uvt−t1St+uτt−t2It,where uvt=0fort∈−t1,0 and uτt=0for t∈−t2,0.To this end, the optimal control strategy is used to achieve these goals. The objective function is considered as(86)Juvt,uτt,tf=∫0tfB1It+B2St−B3Rt+B42uvt2+B52uτt2dt,where B1, B2, and B3are gains of the infected, susceptible, and recovered individuals. Also, B4 and B5 are the gains of control inputs uv, and uτ, respectively. Moreover, tf and Xtf are fixed. Tacking the Hamiltonian yields
(87)H=B1It+B2St−B3Rt+B42uvt2+B52uτt2+p1Ṡt+p2Ėt+p3İt+p4Ṙt,where p1.p2,p3, and p4 are adjoint functions that satisfy the adjoint equations. Pontryagin’s maximum principle is applied to the Hamiltonian H to derive the necessary conditions for the optimal control. The adjoint equation and transversality conditions are obtained by using Pontryagin’s maximum principle such that
(88)ṗ1=−∂H∂St=−B2+βtItp1−p2+uvt−t1p1−p4+p1μd,p1tf=0,ṗ2=−∂H∂Et=κp2−p3+p2μd,p2tf=0,ṗ3=−∂H∂It=−B1+βtStp1−p2+uτt−t2p3−p4+p3α−p4g,p3tf=0,ṗ4=−∂H∂Rt=B3+p4μd,p4tf=0.Notice that the controluvt−t1 can only influence the state variable on the time interval tf−t1,tf. Thus, it suffices to compute the optimal control on the interval 0,tf−t1. Therefore, we have
(89)∂H∂uvt+∂H∂uvt−t1t+t1=0,which yields
(90)uv∗t=S∗t+t1p1t+t1−p4t+t1B4,t∈0,tf−t1,uv∗t=0,t∈tf−t1,tf.Similarly,(91)∂H∂uτt+∂H∂uτt−t2t+t2=0that can be written as
(92)uτ∗t=I∗t+t2p3t+t2−p4t+t2B5,t∈0,tf−t2,uτ∗t=0,t∈tf−t2,tf,which yields
(93)uv∗t=1S∗t+t1p1t+t1−p4t+t1B4≥1,S∗t+t1p1t+t1−p4t+t1B40<S∗t+t1p1t+t1−p4t+t1B4<1,0S∗t+t1p1t+t1−p4t+t1B4≤0,uτ∗t=1I∗t+t2p3t+t2−p4t+t2B5≥1,I∗t+t2p3t+t2−p4t+t2B50<I∗t+t2p3t+t2−p4t+t2B5<1,0I∗t+t2p3t+t2−p4t+t2B5≤0.So, the optimal controls are given as(94)uv∗t=maxminS∗t+t1p1t+t1−p4t+t1B4,1,0,t∈0,tf−t1,uτ∗t=maxminI∗t+t2p3t+t2−p4t+t2B5,1,0,t∈0,tf−t2.Therefore, using the characterization of the optimal control and considering Equations (90) and (92), the following optimality system is obtained as
(95)Ṡ∗t=−βtS∗tI∗t+μd31−ρI∗t−S∗t−uv∗tS∗t,Ė∗t=βtS∗tI∗t−μd+κE∗t,İ∗t=κE∗t−gI∗t−uτ∗tI∗t,Ṙ∗t=gI∗t−μdR∗t+uv∗tS∗t+uτ∗tI∗t.The adjoint equations are also obtained as(96)ṗ1=−B2+βtI∗tp1−p2+uv∗t−t1p1−p4+p1μd,ṗ2=κp2−p3+p2μd,ṗ3=−B1+βtS∗tp1−p2+uτ∗t−t2p3−p4+p3α−p4g,ṗ4=B3+p4μd,where p1tf=0, p2tf=0, p3tf=0, p4tf=0, p5tf=0, p6tf=0, p7tf=0, S0=S0, E0=E0, I0=I0, and R0=R0.Remark 8.
According to the proposed controller, susceptible and infected persons converge to zero as⟶∞St,It⟶0; then, Equation (83) becomes a first-order linear differential equation as Ėt=−μd+κEt in which μd+κ is a positive constant and yields Et=E0e−μd+κt. Thereupon, Et converges exponentially to zero indirectly and as a result of zeroing infected and susceptible people.After investigating the effect of the proposed controller on both models, the results of the controller actions should be illustrated. Therefore, the simulation results of the proposed work are outlined in the next section.


## 7. Simulation Results

In this section, a numerical simulation is given to prove the accuracy of the results obtained in the previous sections. First, the effect of the optimal controllers on the innate immune response and then on the SEIR epidemic model is assessed. The value of innate immune response parameterskI,kR,kE,qK,ΦK,aI,ρV,δH,δI,δR,δIFN,δK,δV,SH,SK is based on experimental data from Hernandez-Vargas and Meyer-Hermann [7], and the parameter values of the SEIR epidemic model are κ=g=ρ=0.3 and μd=0.02. The initial values of the SEIR epidemic model states and innate immune response model are shown in Tables 2 and 3, respectively.Table 2
Initial values of the SEIR epidemic model.
State variableInitial valueS09375E0125I0313R00N010,000Table 3
Initial values of innate immune response.
State variableSusceptible personInfected personExposed personUH05×1085×1085×108UE003×1030UI001․5×1030UR0000V0010−310−3IFN0000K08×1058×1058×105This section is categorized as follows:(a)
The influence of model parameters(b)
The effect of the optimal controller on the within-host model. This section is divided into the following two subsections:(i)
Within-host vaccination(ii)
Within-host therapy(c)
The effect of the optimal controller on the between-hosts model. This section is divided into the following two subsections:(i)
Between-hosts vaccination(ii)
Between-hosts therapy


### 7.1. Influence of Model Parameters

This section is aimed at conducting the sensitivity analysis to identify the factors that noticeably affect the number of people engaged in influenza. Therefore, the changes in the value of parametersβt,g,andkare investigated, respectively. The tremendous influence of contact rate cannot be neglected. The susceptible people get infected if they have more contact with infected people; as a result, the epidemic widely spreads. Thus, the contact rate has to be considered as a critical parameter in the epidemic containment process. To this end, the value of the contact rate is changed from 0.6βtto1.4βt to observe its impact on the population size in each engaged group. As shown in Figures 6(a)–6(d), by increasing in the value of contact rate to 1.4βt, as a result, the more contact between people, the number of susceptible people reduces regarding increment in the number of infected people, while the number of susceptible people does not decrease more because of the less contact between infected individuals 0.6βt. Obviously, the number of exposed and infected people decreases followed by nearly half decrement in the contact rate with susceptible people. Thus, few people get infected; as a result, few people have recovered.Figure 6
The impact of the contact rate changes on the (a) susceptible, (b) exposed, (c) infected, and (d) recovered people.
(a)(b)(c)(d)In Figure7, the impact of changes in the value of g, which is the recovery rate of infected people, is presented. As shown in Figure 7(a), when the rate of g is doubled, a greater number of susceptible people have recovered; as a result, fewer people get infected and the number of susceptible people is increased. Moreover, the number of infected people increases when the recovery rate decreases and decreases with increment in the value of recovery rate because more infected people are recovered. As a result, it is obvious that the number of recovered people is grown up followed by the rise in the recovery rate.Figure 7
The impact of the recovery rate changes on the (a) susceptible, (b) exposed, (c) infected, and (d) recovered people.
(a)(b)(c)(d)According to Figure8, the exposed people are infected at a rate of κ; therefore, increment in this rate leads to growth in the infected people number; as a result, there is growth in the number of recovered people. It is obvious that the number of susceptible people is decreased when κ increases.Figure 8
The impact of theκ rate changes on the (a) susceptible, (b) exposed, (c) infected, and (d) recovered people.
(a)(b)(c)(d)As observed, for different values ofβt,g,andk, the disease does not disappear in about 25 days; therefore, it is necessary to use vaccination and antiviral treatment strategies to overcome the epidemic. Thus, in the next sections, the model in the presence of vaccination and antiviral treatment is simulated compared to the model without the controller.


### 7.2. Optimal Control Results in the Within-Host Model

In the real world, susceptible people will not immediately become immune to the infection following vaccination because it takes time for all healthy cells (within the susceptible individual’s body) to become resistant-to-infection cells. Moreover, the destruction of the virus (and at the same time the transformation of healthy cells into resistant-to-infection cells) within the infected person’s body is not immediate with antiviral therapy; therefore, the infected person is treated after a while. To this end, two-time delay (the length of vaccination and antiviral treatment) is considered in the SEIR model.


#### 7.2.1. Within-Host Vaccination

Vaccine injection into the susceptible individual’s body makes the healthy cells resistant-to-infection. Consequently, the susceptible individual, which includes healthy cells, will be recovered against the infection. Converting the healthy cells into the resistant-to-infection cells takes time, which is considered as delay in recovering process of susceptible individuals in the SEIR model. In Section4, initially, the optimal time of healthy cells converting into resistant-to-infection cells in the susceptible individual’s body has been obtained. Numerical simulations suggest eight days as the optimal time to convert the healthy cells into the resistant ones in a susceptible person’s body (see Figure 9). The proposed optimal time will be considered as a delay of moving susceptible people to the recovered people group in the SEIR model in the next section.Figure 9
The optimal final time.Figures10(a) and 10(b) show the effect of control input (vaccination) on the number of healthy and resistant-to-infection cells. There is a surge in the number of resistant-to-infection cells during vaccination in eight days. After vaccination, it can be deduced that rising in the resistant-to-infection cell numbers is associated with the reduction in the healthy cell numbers, indicating the controller efficacy. As mentioned, as the number of resistant-to-infection cells increases moment by moment, the susceptible person eventually becomes immune to the disease after eight days because the healthy cells in his/her body have all become resistant.Figure 10
Evolution of healthy cells and resistant-to-infection cells with and without the controller.
(a)(b)


#### 7.2.2. Within-Host Therapy

There are virus particles (as a result, infected and partially infected cells) in the infected person’s body. Therefore, the infected person needs to be treated with antiviral treatment. Thus, over the treatment process, the infected and partially infected cells and virus particles are killed, and the healthy cells become resistant-to-infection simultaneously, which also takes time.In Figure11, the final optimal time to eradicate disease computed in Section 4 is shown. As shown in Figures 12(a) and 12(b), using antiviral drugs, healthy cells become resistant-to-infection cells. The drug’s effect on infected and partially infected cells and viruses can also be seen in Figures 13 and 14, respectively. It is clear that the infected and partially infected cells are destroyed within six days, and the number of virus particles is also converged to zero, which indicates that the infected person is cured. Moreover, the evolution of IFN-I molecules and natural killer cells with and without the controller is shown in Figures 15(a) and 15(b).Figure 11
The optimal final time.Figure 12
Evolution of healthy cells and resistant-to-infection cells with and without the controller.
(a)(b)Figure 13
Evolution of partially infected cells and infected cells with and without the controller.
(a)(b)Figure 14
Evolution of virus particle with and without the controller.Figure 15
Evolution of IFN-I molecules and natural killer cells with and without the controller.
(a)(b)


### 7.3. Optimal Control Results between Hosts

In this section, we examine the effect of the optimal controller on the SEIR epidemic model, which was discussed in Section7, in detail. First, the changes in the number of susceptible, exposed, infected, and recovered individuals in the absence of appropriate control are investigated. Then, the effect of the vaccination strategy and antiviral treatment on individuals in a community is examined. For the numerical simulations, the delay differential equation solver (dde23) is used. This solver does not solve the advance equations t+ti directly; therefore, it is necessary to convert the advance equations into delay differential equations. To this end, a new reversed-time variable ς=tf−t is defined and new adjoint variables λiς=pitf−ς are introduced. Thus, this new system is solved numerically using dde23.


#### 7.3.1. Between-Hosts Vaccination

As shown in Figures16(a) and 16(b), the number of susceptible individuals in the absence of a controller (vaccination has not been applied) descends and converts to infected people. Eventually, it remains stationary at a constant value (less than 2,000 people). Besides, the number of recovered people without any controller is greater than their number in the control-based strategy. The number of recovered people has reached more than 10,000 because in the absence of effective treatment, more people become infected, and as a result, more people will be recovered. By injecting the vaccine into the susceptible person’s body, during the conversion of healthy cells to resistant cells, the susceptible person also vaccinates (direct transfer of individuals from the susceptible group to the recovered group). As a result, the number of susceptible individuals decreases and moves to the recovered individual group at the same time. This means the susceptible people vaccinate against the disease in almost 25 days t3 (the healthy cells of the susceptible individual convert to resistant-to-infection cells in eight days as the delay time with vaccination), while in [13], the final time of the vaccination campaign of nearly 40 days is suggested without considering the delay in the vaccination of susceptible people.Figure 16
The population of susceptible and recovered individuals with and without the controller.
(a)(b)


#### 7.3.2. Between-Hosts Therapy

When the infection spreads in the community without control over its prevalence, the number of infected people is increasing every day, as shown in Figure17(a). It rose by over 3,000 people. Finally, their number decreases until the 25th day and reaches a little under 2,000 and remains unchanged, indicating the community’s disease persistence. Also, the exposed people number reaches a high number of 4,000 in only twelve days and ultimately converges to a stable population and remains at almost 1,500 (see Figure 17(b)).Figure 17
The population of infected and exposed individuals with and without the controller.
(a)(b)In contrast, the number of infected people exceeds almost 500 with antiviral therapy, which is used to treat infected people (transfer them to a recovered group). During the computed optimal time (the same delay as six days), the healthy cells convert to resistant-to-infection cells; at the same time, virus particles are eradicated by antiviral treatment in the infected person’s body. Therefore, the infected people move to the recovered people group in almost 25 days. According to Remark8, the number of exposed people also converges to zero as a result of decrement of infected people in 25 days.


## 7.1. Influence of Model Parameters

This section is aimed at conducting the sensitivity analysis to identify the factors that noticeably affect the number of people engaged in influenza. Therefore, the changes in the value of parametersβt,g,andkare investigated, respectively. The tremendous influence of contact rate cannot be neglected. The susceptible people get infected if they have more contact with infected people; as a result, the epidemic widely spreads. Thus, the contact rate has to be considered as a critical parameter in the epidemic containment process. To this end, the value of the contact rate is changed from 0.6βtto1.4βt to observe its impact on the population size in each engaged group. As shown in Figures 6(a)–6(d), by increasing in the value of contact rate to 1.4βt, as a result, the more contact between people, the number of susceptible people reduces regarding increment in the number of infected people, while the number of susceptible people does not decrease more because of the less contact between infected individuals 0.6βt. Obviously, the number of exposed and infected people decreases followed by nearly half decrement in the contact rate with susceptible people. Thus, few people get infected; as a result, few people have recovered.Figure 6
The impact of the contact rate changes on the (a) susceptible, (b) exposed, (c) infected, and (d) recovered people.
(a)(b)(c)(d)In Figure7, the impact of changes in the value of g, which is the recovery rate of infected people, is presented. As shown in Figure 7(a), when the rate of g is doubled, a greater number of susceptible people have recovered; as a result, fewer people get infected and the number of susceptible people is increased. Moreover, the number of infected people increases when the recovery rate decreases and decreases with increment in the value of recovery rate because more infected people are recovered. As a result, it is obvious that the number of recovered people is grown up followed by the rise in the recovery rate.Figure 7
The impact of the recovery rate changes on the (a) susceptible, (b) exposed, (c) infected, and (d) recovered people.
(a)(b)(c)(d)According to Figure8, the exposed people are infected at a rate of κ; therefore, increment in this rate leads to growth in the infected people number; as a result, there is growth in the number of recovered people. It is obvious that the number of susceptible people is decreased when κ increases.Figure 8
The impact of theκ rate changes on the (a) susceptible, (b) exposed, (c) infected, and (d) recovered people.
(a)(b)(c)(d)As observed, for different values ofβt,g,andk, the disease does not disappear in about 25 days; therefore, it is necessary to use vaccination and antiviral treatment strategies to overcome the epidemic. Thus, in the next sections, the model in the presence of vaccination and antiviral treatment is simulated compared to the model without the controller.


## 7.2. Optimal Control Results in the Within-Host Model

In the real world, susceptible people will not immediately become immune to the infection following vaccination because it takes time for all healthy cells (within the susceptible individual’s body) to become resistant-to-infection cells. Moreover, the destruction of the virus (and at the same time the transformation of healthy cells into resistant-to-infection cells) within the infected person’s body is not immediate with antiviral therapy; therefore, the infected person is treated after a while. To this end, two-time delay (the length of vaccination and antiviral treatment) is considered in the SEIR model.


### 7.2.1. Within-Host Vaccination

Vaccine injection into the susceptible individual’s body makes the healthy cells resistant-to-infection. Consequently, the susceptible individual, which includes healthy cells, will be recovered against the infection. Converting the healthy cells into the resistant-to-infection cells takes time, which is considered as delay in recovering process of susceptible individuals in the SEIR model. In Section4, initially, the optimal time of healthy cells converting into resistant-to-infection cells in the susceptible individual’s body has been obtained. Numerical simulations suggest eight days as the optimal time to convert the healthy cells into the resistant ones in a susceptible person’s body (see Figure 9). The proposed optimal time will be considered as a delay of moving susceptible people to the recovered people group in the SEIR model in the next section.Figure 9
The optimal final time.Figures10(a) and 10(b) show the effect of control input (vaccination) on the number of healthy and resistant-to-infection cells. There is a surge in the number of resistant-to-infection cells during vaccination in eight days. After vaccination, it can be deduced that rising in the resistant-to-infection cell numbers is associated with the reduction in the healthy cell numbers, indicating the controller efficacy. As mentioned, as the number of resistant-to-infection cells increases moment by moment, the susceptible person eventually becomes immune to the disease after eight days because the healthy cells in his/her body have all become resistant.Figure 10
Evolution of healthy cells and resistant-to-infection cells with and without the controller.
(a)(b)


### 7.2.2. Within-Host Therapy

There are virus particles (as a result, infected and partially infected cells) in the infected person’s body. Therefore, the infected person needs to be treated with antiviral treatment. Thus, over the treatment process, the infected and partially infected cells and virus particles are killed, and the healthy cells become resistant-to-infection simultaneously, which also takes time.In Figure11, the final optimal time to eradicate disease computed in Section 4 is shown. As shown in Figures 12(a) and 12(b), using antiviral drugs, healthy cells become resistant-to-infection cells. The drug’s effect on infected and partially infected cells and viruses can also be seen in Figures 13 and 14, respectively. It is clear that the infected and partially infected cells are destroyed within six days, and the number of virus particles is also converged to zero, which indicates that the infected person is cured. Moreover, the evolution of IFN-I molecules and natural killer cells with and without the controller is shown in Figures 15(a) and 15(b).Figure 11
The optimal final time.Figure 12
Evolution of healthy cells and resistant-to-infection cells with and without the controller.
(a)(b)Figure 13
Evolution of partially infected cells and infected cells with and without the controller.
(a)(b)Figure 14
Evolution of virus particle with and without the controller.Figure 15
Evolution of IFN-I molecules and natural killer cells with and without the controller.
(a)(b)


## 7.2.1. Within-Host Vaccination

Vaccine injection into the susceptible individual’s body makes the healthy cells resistant-to-infection. Consequently, the susceptible individual, which includes healthy cells, will be recovered against the infection. Converting the healthy cells into the resistant-to-infection cells takes time, which is considered as delay in recovering process of susceptible individuals in the SEIR model. In Section4, initially, the optimal time of healthy cells converting into resistant-to-infection cells in the susceptible individual’s body has been obtained. Numerical simulations suggest eight days as the optimal time to convert the healthy cells into the resistant ones in a susceptible person’s body (see Figure 9). The proposed optimal time will be considered as a delay of moving susceptible people to the recovered people group in the SEIR model in the next section.Figure 9
The optimal final time.Figures10(a) and 10(b) show the effect of control input (vaccination) on the number of healthy and resistant-to-infection cells. There is a surge in the number of resistant-to-infection cells during vaccination in eight days. After vaccination, it can be deduced that rising in the resistant-to-infection cell numbers is associated with the reduction in the healthy cell numbers, indicating the controller efficacy. As mentioned, as the number of resistant-to-infection cells increases moment by moment, the susceptible person eventually becomes immune to the disease after eight days because the healthy cells in his/her body have all become resistant.Figure 10
Evolution of healthy cells and resistant-to-infection cells with and without the controller.
(a)(b)


## 7.2.2. Within-Host Therapy

There are virus particles (as a result, infected and partially infected cells) in the infected person’s body. Therefore, the infected person needs to be treated with antiviral treatment. Thus, over the treatment process, the infected and partially infected cells and virus particles are killed, and the healthy cells become resistant-to-infection simultaneously, which also takes time.In Figure11, the final optimal time to eradicate disease computed in Section 4 is shown. As shown in Figures 12(a) and 12(b), using antiviral drugs, healthy cells become resistant-to-infection cells. The drug’s effect on infected and partially infected cells and viruses can also be seen in Figures 13 and 14, respectively. It is clear that the infected and partially infected cells are destroyed within six days, and the number of virus particles is also converged to zero, which indicates that the infected person is cured. Moreover, the evolution of IFN-I molecules and natural killer cells with and without the controller is shown in Figures 15(a) and 15(b).Figure 11
The optimal final time.Figure 12
Evolution of healthy cells and resistant-to-infection cells with and without the controller.
(a)(b)Figure 13
Evolution of partially infected cells and infected cells with and without the controller.
(a)(b)Figure 14
Evolution of virus particle with and without the controller.Figure 15
Evolution of IFN-I molecules and natural killer cells with and without the controller.
(a)(b)


## 7.3. Optimal Control Results between Hosts

In this section, we examine the effect of the optimal controller on the SEIR epidemic model, which was discussed in Section7, in detail. First, the changes in the number of susceptible, exposed, infected, and recovered individuals in the absence of appropriate control are investigated. Then, the effect of the vaccination strategy and antiviral treatment on individuals in a community is examined. For the numerical simulations, the delay differential equation solver (dde23) is used. This solver does not solve the advance equations t+ti directly; therefore, it is necessary to convert the advance equations into delay differential equations. To this end, a new reversed-time variable ς=tf−t is defined and new adjoint variables λiς=pitf−ς are introduced. Thus, this new system is solved numerically using dde23.


### 7.3.1. Between-Hosts Vaccination

As shown in Figures16(a) and 16(b), the number of susceptible individuals in the absence of a controller (vaccination has not been applied) descends and converts to infected people. Eventually, it remains stationary at a constant value (less than 2,000 people). Besides, the number of recovered people without any controller is greater than their number in the control-based strategy. The number of recovered people has reached more than 10,000 because in the absence of effective treatment, more people become infected, and as a result, more people will be recovered. By injecting the vaccine into the susceptible person’s body, during the conversion of healthy cells to resistant cells, the susceptible person also vaccinates (direct transfer of individuals from the susceptible group to the recovered group). As a result, the number of susceptible individuals decreases and moves to the recovered individual group at the same time. This means the susceptible people vaccinate against the disease in almost 25 days t3 (the healthy cells of the susceptible individual convert to resistant-to-infection cells in eight days as the delay time with vaccination), while in [13], the final time of the vaccination campaign of nearly 40 days is suggested without considering the delay in the vaccination of susceptible people.Figure 16
The population of susceptible and recovered individuals with and without the controller.
(a)(b)


### 7.3.2. Between-Hosts Therapy

When the infection spreads in the community without control over its prevalence, the number of infected people is increasing every day, as shown in Figure17(a). It rose by over 3,000 people. Finally, their number decreases until the 25th day and reaches a little under 2,000 and remains unchanged, indicating the community’s disease persistence. Also, the exposed people number reaches a high number of 4,000 in only twelve days and ultimately converges to a stable population and remains at almost 1,500 (see Figure 17(b)).Figure 17
The population of infected and exposed individuals with and without the controller.
(a)(b)In contrast, the number of infected people exceeds almost 500 with antiviral therapy, which is used to treat infected people (transfer them to a recovered group). During the computed optimal time (the same delay as six days), the healthy cells convert to resistant-to-infection cells; at the same time, virus particles are eradicated by antiviral treatment in the infected person’s body. Therefore, the infected people move to the recovered people group in almost 25 days. According to Remark8, the number of exposed people also converges to zero as a result of decrement of infected people in 25 days.


## 7.3.1. Between-Hosts Vaccination

As shown in Figures16(a) and 16(b), the number of susceptible individuals in the absence of a controller (vaccination has not been applied) descends and converts to infected people. Eventually, it remains stationary at a constant value (less than 2,000 people). Besides, the number of recovered people without any controller is greater than their number in the control-based strategy. The number of recovered people has reached more than 10,000 because in the absence of effective treatment, more people become infected, and as a result, more people will be recovered. By injecting the vaccine into the susceptible person’s body, during the conversion of healthy cells to resistant cells, the susceptible person also vaccinates (direct transfer of individuals from the susceptible group to the recovered group). As a result, the number of susceptible individuals decreases and moves to the recovered individual group at the same time. This means the susceptible people vaccinate against the disease in almost 25 days t3 (the healthy cells of the susceptible individual convert to resistant-to-infection cells in eight days as the delay time with vaccination), while in [13], the final time of the vaccination campaign of nearly 40 days is suggested without considering the delay in the vaccination of susceptible people.Figure 16
The population of susceptible and recovered individuals with and without the controller.
(a)(b)


## 7.3.2. Between-Hosts Therapy

When the infection spreads in the community without control over its prevalence, the number of infected people is increasing every day, as shown in Figure17(a). It rose by over 3,000 people. Finally, their number decreases until the 25th day and reaches a little under 2,000 and remains unchanged, indicating the community’s disease persistence. Also, the exposed people number reaches a high number of 4,000 in only twelve days and ultimately converges to a stable population and remains at almost 1,500 (see Figure 17(b)).Figure 17
The population of infected and exposed individuals with and without the controller.
(a)(b)In contrast, the number of infected people exceeds almost 500 with antiviral therapy, which is used to treat infected people (transfer them to a recovered group). During the computed optimal time (the same delay as six days), the healthy cells convert to resistant-to-infection cells; at the same time, virus particles are eradicated by antiviral treatment in the infected person’s body. Therefore, the infected people move to the recovered people group in almost 25 days. According to Remark8, the number of exposed people also converges to zero as a result of decrement of infected people in 25 days.


## 8. Conclusion

In this paper, the optimal control theory was applied to two types of dynamic epidemic models: the innate immune response dynamic and the SEIR epidemic model. A nonlinear transmission rate and time delay were considered in the SEIR epidemic model. The aim of this paper was twofold: firstly, preventing the virus from spreading within the body of a host and secondly, curb the epidemic in the society and between hosts. To this end, two control strategies were introduced to satisfy the first goal in the susceptible and infected person’s bodies. In this regard, vaccination was used to convert healthy cells to resistant-to-infection cells inside the susceptible individual’s body that recovered them against the disease in society. In the following, antiviral treatment was used to reduce the concentration of viruses and convert healthy cells to resistant-to-infection cells at the same time inside the infected individual’s body to recover them. Moreover, to control the epidemic in society, the optimal control was used to increase the number of recovered people by reducing the number of susceptible and infected people in the fixed 25 days. Transferring occurred with a delay computed as the optimal time of disease eradication inside the infected individual’s body and the recovery of the susceptible individual’s body. Finally, numerical simulation was used to illustrate the usefulness of the obtained results.



---

*Source: 4348910-2021-10-18.xml*