Hindawi Computational and Mathematical Methods in Medicine Volume 2021, Article ID 4348910, 27 pages https://doi.org/10.1155/2021/4348910 Research Article Optimal Allocation of Vaccine and Antiviral Drugs for Influenza Containment over Delayed Multiscale Epidemic Model considering Time-Dependent Transmission Rate Zohreh Abbasi ,1 Iman Zamani ,2 Amir Hossein Amiri Mehra ,1,2 Asier Ibeas,3,4 and Mohsen Shafieirad 1 1 Department of Electrical and Computer Engineering, University of Kashan, Iran 2 Electrical and Electronic Engineering Department, Shahed University, Tehran, Iran 3 Departament de Telecomunicació i Enginyeria de Sistemes, Escolad’Enginyeria, Universitat Autònoma de Barcelona, Barcelona, Spain 4 Universidad de Bogotá Jorge Tadeo Lozano, Bogotá D.C., Colombia Correspondence should be addressed to Iman Zamani; zamaniiman@shahed.ac.ir Received 26 June 2021; Revised 19 July 2021; Accepted 26 August 2021; Published 18 October 2021 Academic Editor: Andrei Korobeinikov Copyright © 2021 Zohreh Abbasi et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. 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. 1. Introduction ual’s body. However, in the majority of virus-infected individuals, the innate immune response is insufficient Seasonal and pandemic influenza A virus (IAV) infection because viruses use different evasion strategies to escape causes severe morbidity and mortality worldwide [1]. There the immune response. Innate immunity plays a critical role are many drugs to treat influenza, but the immune system is in the control of viral infection because most infectious naturally the first defensive line against the disease. There- pathogens are eliminated through innate immune response fore, it is essential to understand the mechanisms that trigger without necessarily requiring the activation of adaptive the immune system response and how this response affects immunity [2]. the overall spread of the disease. The innate immune Since innate immunity is an important part of the body response plays an important role in the control and clear- for dealing with viral infections, some methods have been ance of pathogens following viral infection inside an individ- determined to model different aspects of the human 2 Computational and Mathematical Methods in Medicine immune system. Differential equation models are the most Optimal control is a mathematical method derived from popular methods to simulate critical diseases’ immunologi- the calculus of variations. There are different methods for cal and epidemiological dynamics [3]. In addition, the use computing the optimal control for mathematical models of mathematical modeling of the innate immune response [14]. Pontryagin’s maximum principle, first formulated in interpreting experimental results has made a significant con- 1956 by L. S. Pontryagin [15], is to calculate the optimal con- tribution to this field. Many mathematical models of the trol for an ordinary differential equation model system with immune response to IAV infection have been introduced a given constraint. In [16], an optimal control method is that considered both different aspects of the immune system applied to a generic model of a pathogenic attack on the and different detail levels. The first mathematical model to innate immune system. The minimization of a quadratic illustrate the within-host dynamic of IAV in infected mice cost function is generated by numerical optimization. Ter- was developed in 1976 by Larson et al. [4]. The proposed minal optimal time control is also investigated. For further compartmental model includes seven compartments with studies, the readers are referred to [2], and the references five associated rate parameters. Also, the virus population therein introduce an optimal control model based on an dynamic was introduced in mathematical terms. Another ordinary differential equation system. Some studies have model introduced by Bocharov and Romanyukha [5] studied investigated the optimal control of SIR and SEIR models, immune responses to viral infections in human infections, but we mention just some of them in the following. The such as the influenza A virus, and considered 12 immune authors in [17] used an optimal control approach to an populations. Also, Baccam et al. [6] described the viral kinet- SIR epidemic model with a time delay to minimize the ics of influenza A during infection in the human body, spread of infected individuals and maximize the number of including target cell and the innate interferon response with susceptible and recovered individuals. Likewise, in [18, 19], delayed virus production through infected cells. However, as optimal control is considered to minimize the level of infec- opposed to previous papers, which did not consider a com- tion at the terminal time for an SIR and SEIR model, respec- plete model of the innate immune system, we used a more tively, where the incidence rate is an unspecified nonlinear comprehensive one adapted from [7]. function, while their models only consider one aspect of However, the disease’s impact within the community the impact of the disease that affects the community and and the relationship between people in society also need to the relationship between susceptible, infected, and recovered be examined. Mathematical models are instrumental in people. They do not take into account the effects of the dis- explaining and quantifying epidemic dynamics. Hence, var- ease inside the body, and the relationship between cells is ious mathematical models have been proposed in the follow- not considered. ing for between-hosts epidemic models to provide a general Epidemic mathematical models can be used for several overview of strategies to be employed in the period of influ- epidemic diseases like COVID-19 by modifying the param- enza virus infection. Kermack and McKendrick created one eter and adding or eliminating some states. In this regard, of the simplest types to explain epidemic dynamics that is the authors in [20] considered a SEIAR epidemic model called SIR (susceptible, infected, and recovered individuals, for the COVID-19 pandemic. Also, some new SIR-type respectively). The proposed model was described in three models are introduced in [21–24]. The authors in [21] have papers in 1927 [8], 1932 [9], and 1933 [10]. Moreover, in considered “isolated,” “vaccinated,” and “quarantined” people [11], another type of epidemic model (SIS) is considered groups due to several different conditions (after and before with a nonlinear incidence rate and time delay. In addition, vaccine development). Moreover, “hospitalized” and “ICU- four different control strategies are investigated and com- admitted” people groups are introduced in [22] to express pared, including optimal control to decrease the density of the spread of coronavirus between the healthcare workers. infected individuals, increase that of susceptible ones, and Also, in [23], a new epidemic model (called SIDARTHE) is reduce relevant costs. Besides these, in [12], a model of the studied that considers the distinction between diagnosed and relationship between two types of systems is formulated nondiagnosed infected individuals with and without symp- and analyzed. The first dynamic is the within-host model toms detected acutely symptomatic infected ones. Besides, dynamic and focuses on cellular interactions. The second is the discrete-time SEIR epidemic models can be used in such the between-hosts model dynamic that focuses on transmis- studies instead of continuous ones [25, 26], to name a few. sion and infection statuses that are governed by the suscep- To further study about optimal control with time delay, tible-exposed-infected-recovered (SEIR) model. Likewise, in readers refer to [27–31] and references therein. this study, motivated by the previous paper, a combined Optimal control theory can be used to control other model with different relationships is presented that will be diseases like HIV [32] and especially COVID-19. In this discussed in detail in the following sections. The authors in regard, the authors in [33] used optimal control theory based [13] have used the sliding mode control on the SEIAR on the SQEIAR epidemic model (Susceptible, Quarantined, dynamic to overcome the uncertainty of parameters. How- Exposed, Infected, Asymptomatic, and Recovered) for the ever, they have not considered the optimal time to eradicate eradication of infection considering quarantine and treat- the infection in society, and they only refer to the between- ment policies in China and Spain for different types of dis- hosts dynamic. For this reason, since time is one of the pri- ease (COVID-19, Ebola, and influenza). Also, according to orities in order to control disease, we use an optimal control travel between cities that cause an increase in the outbreak, strategy to get the optimal time to cure the infected person the authors in [33] investigated the impact of travel or immi- and recover the susceptible person. gration and impulsive change of population. These changes Computational and Mathematical Methods in Medicine 3 lead to generate impulsive epidemic models and impulsive 2. Preliminaries control strategies [34]. Also, the authors in [35] have consid- ered an optimal control strategy in the age-structured SEIRQ In this section, we first look at the relationship between the model for the COVID-19 epidemic in Brazil, where the two models considered in this paper. According to quarantine entrance parameters are the control efforts of Figures 1 and 2, our multiscale model considers two scales: the system. In the same way, to explore the dynamic of the the first is the interaction of the virus and innate immune COVID-19 pandemic in Pakistan, the necessary optimality response within a host (inside the body), and the second is conditions are introduced in [36] using the optimal control the contact between susceptible, exposed, infected, and theory in a mathematical model. Some optimal control strat- recovered individuals in society. The dynamic of the innate egies with application to COVID-19 are also introduced in immune response and relationship between cells and other [37], for further study. components are now described in detail. The within-host This paper investigates the correlation between two model dynamic is the innate immune response to the influ- modeling aspects: the effect of the disease on society and enza virus using a mathematical model. This model is based inside the individual’s body, which is an important issue. on interferon-induced resistance to infection of respiratory The community includes four groups (SEIR), all of which epithelial cells and the clearance of infected cells by natural relate to each other in a community. All communications killers with seven state variables (the numbers of healthy occur only in the community and among individuals. In (U H ), partially infected (U E ), infected (U I ), and resistant- contrast, there are some interactions among the cells to-infection cells (U R ) and IFN-I molecules ([IFN]), natural inside an individual’s body. Recovering the susceptible by killer cells ðKÞ, and virus particles ðVÞ). Also, a nonlinear vaccination takes time that is considered as delay caused SEIR epidemiological model for influenza is given where by the length of converting the healthy cells into the nonnegative state variables SðtÞ, EðtÞ, IðtÞ, and RðtÞ are resistant-to-infectious cells. Also, the infected people need the susceptible, exposed, infected (symptomatic), and recov- to be treated with antiviral treatment to recover, which is ered compartments, respectively. Here, SðtÞ represents the also considered as a delay (the length of time to kill the number of people that are susceptible but not infected with infected and partially infected cells and viruses and con- influenza, EðtÞ denotes the number of people exposed to vert the healthy cells into resistant-to-infection cells influenza (infected but not yet infectious), IðtÞ denotes the simultaneously). population of infected humans with infectious influenza The goal of the present paper will be presented clearly symptoms, RðtÞ is the number of recovered people from and thoroughly regarding the reviewed papers. This paper influenza, and N is the total population size. The relation- introduces a multiscale model in which both innate ship between the two epidemic models is shown in immune responses (within-host model) and the SEIR Figures 1 and 2. The dynamics of the innate immune dynamic (between-hosts model) are used simultaneously response and SEIR epidemiological model is described fur- to define a single system. A multiscale-type model has ther in the following sections. already been presented in [12], with the difference that The general strategy for this study is described by the fol- the proposed model in the present study is more compre- lowing Algorithm 1. hensive and concisely describes the relationships between the virus, cells, and the immune system. The model is also 3. The Innate Immune Response more realistic and closer to the real world. Moreover, a Dynamic Model nonlinear sinusoidal transmission rate (contact rate) is considered, which varies in time according to changes in In this section, we consider innate immune system dynamics contacts between people throughout day and night. In this for dealing with the influenza virus, which is derived from study, the optimal control is used to eradicate the infection [7]. There are seven components to the dynamic states, at the optimal time in the community. As discussed in and the following equations govern the dynamics. For brev- Introduction, different papers have been published on the ity, we suppress the notation “ðtÞ” for variables at the current subject of optimal control. However, to the best of the time. authors’ knowledge, an optimal control approach that simultaneously tackles both types of dynamics (within-host U_ H = SH − kI U H V − kR U H ½I FN  − δH U H , ð1Þ and between-hosts) for epidemic models has never been considered in the literature before. U_ E = kI U H V − kE U E − qK U E K, ð2Þ 1This paper is organized as follows. Section 2 describes the preliminaries. After that, we propose an innate U_ I = kE U E − δI U I − qK U I K, ð3Þ immune response dynamic in Section 3. The optimal con- trol of innate immune response is introduced in Section 4, U_ R = kR U H ½I FN  − δR U R , ð4Þ followed by the SEIR dynamic model with delay and non- linear transmission rate in Section 5. In Section 6, we V_ = ρV U I − δV V, ð5Þ apply optimal control of the SEIR epidemic model. The results of simulations of different cases and sensitivity _  = a U − δ ½IFN , ½I FN ð6Þ I I I FN analysis are presented in Section 7. Finally, Section 8 sum- marizes the conclusions. K_ = SK + ΦK U I − δK K, ð7Þ 4 Computational and Mathematical Methods in Medicine Susceptible people UH UR Recovered people Infected people UH UR V 0 Figure 1: Block diagram of the relationship between innate immune response and SEIR model. Healthy Virus Infected cells particles cells Within-host (Exposed) Within-host (Infected) Susceptible Exposed Infected Recovered Between-hosts Between-hosts Antiviral tratment for infected Vaccination for susceptible Figure 2: The disease modeling with two scales: within a host and between hosts. Step1. Optimal control design (vaccination strategy) for within-host model to convert healthy cells to resistant-to-infection cells in the optimal time ðt 1 Þ. 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 time ðt 2 Þ. Step3. Optimal control design for between-hosts model to eradicate the infection, considering delays ðt 1 , t 2 Þ in the community. Algorithm 1: General strategy algorithm. where U H is the number of healthy cells that increase with NKs (natural killers) increases with a constant cell growth cell growth rate SH . These cells become partially infected rate SK . In the same way, NKs die at the rate of δK . After cells ðU E Þ via the virus particles ðVÞ with rate kI . They also inflammatory stimulation caused by infected cells, natural convert to resistant-to-infection cells ðU R Þ by the rate kR killers are recruited from the blood indicated by ΦK U I . with IFN-I molecules ð½IFNÞ that increase by the rate of aI Figure 3 shows a schematic diagram of the innate immune . Partially infected cells are in the eclipse phase. This means system response to the influenza virus, which shows the that they cannot spread the virus in the body and they infect cells’ relation very well. the other healthy cells by the term kE . When the partially After introducing the model, we now introduce the con- infected cells pass through the eclipse phase, they will troller and investigate how to influence the vaccine and anti- become infected cells ðU I Þ and can spread the virus and viral treatment on the cells. As mentioned earlier, the cause other healthy cells to become infected. Partially vaccination is used to make healthy cells resistant to infec- infected ðU E Þ and infected cells ðU I Þ are destroyed by natu- tion in a susceptible person’s body and by using the antiviral ral killers ðKÞ at a rate of qK . Once the cells are infected ðU I Þ, treatment, the virus particles, infected, and partially infected they broadcast thousands of similar cells a day, which is cells are removed inside the infected body as discussed in denoted by the parameter ρV . The natural death rates of detail in the next section. healthy, infected, infection-resistant cells, and IFN-I mole- cules are denoted by δH , δI , δR ,and δIFN , respectively, and 3.1. Equilibrium Points and the Basic Reproduction Number. the clear rate of virus particles ðVÞ is δV . The number of To determine the disease spread and the number of Computational and Mathematical Methods in Medicine 5 R0 is calculated as R0 = δH SH kI kE ρV δK 2 /ðδV ðkE δK + qK SK Þð (IFN) δI δK + qK SK ÞÞ. SH The disease-free equilibrium point ðDFEÞ is calculated δIFN(IFN) aIUI by equating Equations (1)–(7) to zero and solving for kR kI kEUE U E ∗ = U I ∗ = V ∗ = ½IFN∗ = 0 as ðSH /δH , 0, 0, 0, 0, 0, 0, SK / UR UH UE UI δK Þ. To find the endemic equilibrium point ðEEÞ, Equa- qk tions (1)–(7) are equated to zero, U H = SH /ðkI V + kR ½IFN qk UHv δIUI + δH Þ, U E = kI U H V/ðkE + qK KÞ, U I = kE U E /ðδI + qK KÞ, δRUR δHUH U R = ðkR U H ½IFNÞ/δR , V = ρV U I /δV , and ½IFN = aI U I / δIFN . Therefore, by replacing the states, the following 𝜌VUI SK 𝛷KUI equations can be calculated: U  H = SH /ða1 U  I + δH Þ, U E = K  I − a1 a4 U 2  I ÞU I, U  R = kR SH δV aI /ða7 U  I − δH a2 /ða3 δH + a5 U δVV δKK  I, V = ðρ /δV ÞU  I , ½IFN = ðaI /δIFN ÞU  I, K = ðSK δV δR δIFN ÞU V   4  3  2  Figure 3: Conceptual flow diagram of innate immune dynamics + ΦK U I Þ/δK , and α1 U I + α2 U I + α3 U I + α4 U I = 0. The last with vaccination. equation has four solutions that one of them is DFE and can be achieved by U  I = 0. The other three solutions are secondary infections, an infected individual would produce obtained using the parameter values in Table 1, one of in his/her infectious period the most important threshold which is 1:84 × 105 , and the other two solutions are nega- criterion named “the basic reproduction number” ðR0 Þ can tive and unacceptable, because the number of cells is a be calculated as the follows. positive value. Therefore, the other equilibrium points are According to the definition of R0 = ρðFV −1 Þ, where given as EE = ðU  H, U  E, U I, U R , V,  ½IFN, KÞ  = ðSH /ða1 U I +  I − a1 a4 U 2  I ÞU I, U  I , kR SH δV aI /ða7 UI − δH Þ, a2 /ða3 δH + a5 U 2 SH 3 δH δV δR δIFN ÞU  I , ðρ /δV ÞU  I , ðaI /δIFN ÞU  I , ðSK + ΦK U  I Þ/δK Þ, 0 0 kI V 6 δH 7 in which, α1 = qK 2 ΦK 2 /δK 2 δIFN ðδIFN kI ρV + kR δI δV Þ, α2 = F=6 60 7, qK ΦK /δK 2 δIFN ð2δIFN kI ρV qK SK + δK δIFN kI ρV ðkE + δI Þ + 2qK 4 0 0 75 kR SK δV aI + kR δV δK aI ðkE + δI Þ + δH δV δIFN qK ΦK Þ, α3 = 1/ 0 0 0 2 3 δK 2 δIFN ðkI ρV ðδK 2 δIFN kE δI + δK δIFN qK SK ðkE + δI Þ + qK 2 δIFN SK ð8Þ SK 2 Þ + kR aI ðδK 2 δV kE δI + δK δV qK SK ðkE + δI Þ + qK 2 δV SK 2 Þ 6 kE + qK δ 0 0 7 + δH δV qK ΦK ð2δIFN qK SK + δK δIFN ðkE + δI ÞÞÞ, and α4 = ð1/ 6 K 7 6 7 6 V =6 SK 7, δK 2 ÞðkE δI δH δV δK 2 + δK δH δV qK SK ðkE + δI Þ + δH qK 2 δV SK 2 6 −kE δI + q K 0 7 7 − kE kI ρV SH Þ. 4 δK 5 Then, the Jacobian of Equations (1)–(7) is given by 0 −ρV δV 2 3 −kI V − kR ½I FN  − δH 0 0 0 −kI U H −kR U H 0 6 7 6 6 kI V −kE − qK K 0 0 −kI U H 0 −qK U E 7 7 6 7 6 6 0 kE −δI − qK K 0 0 0 −qK U I 7 7 6 7 J =6 6 kR ½I FN  0 0 −δR 0 kR U H 0 7 7: ð9Þ 6 7 6 6 0 0 ρV 0 −δV 0 0 7 7 6 7 6 4 0 0 aI 0 0 −δI FN 0 7 5 0 0 ΦK 0 0 0 −δK Theorem 1. The DFE is locally stable if R0 < 1 while it is Hurwitz criterion, the following inequality must be satisfied: unstable if R0 > 1. δV ðkE δK + qK SK ÞðδI δK + qK SK Þ − δH SH kI kE ρV δK 2 > 0: Proof. By calculating the seventh-order characteristic equa- ð10Þ tion of J about its DFE and according to the Routh- 6 Computational and Mathematical Methods in Medicine Table 1: Equilibrium values of the model. 3.3. Positivity Analysis. In this subsection, the positivity of solutions is investigated with Theorem 3. State variable DFE EE H SH Theorem 3. If U H 0 ≥ 0, U E0 ≥ 0, U I 0 ≥ 0, U R0 ≥ 0, V 0 ≥ 0, U 2:44 × 106 δH ½IFN0 ≥ 0, and K 0 ≥ 0, then solutions of the model are posi- E U 0 8:94 × 105 tive for all t ≥ 0: I U 0 1:84 × 105 Proof. According to the positivity of all the parameters of the R U 0 2:37 × 108 _  Ð (4) can be written as U R ðtÞ ≥ −δÐR U R ðtÞ. model, Equation V 0 354.65 δR dr t δ dr t   MultiplyingÐt e 0 into two bothÐ side gives U_ R ðtÞe 0 R + IFN 0 0.138 δ dr t δ dr δR U R ðtÞe 0 R ≥ 0; then, ðd/dtÞðe 0 R UÐ R ðtÞÞ ≥ 0 ; integrat- SK  t K 804:610 × 103 − ing from 0 to t gives U R ðtÞ ≥ U R ð0Þe 0 R ≥ 0. Therefore δ dr δK U R ðtÞ is positive ðU R ðtÞ ≥ 0Þ. In the same way, the other states VðtÞ, ½IFNðtÞ,and KðtÞ (Equations (5)–(7)) Ðare also t Therefore, it can be concluded that R0 < 1: For R0 > 1, positive states and can Ðbe written as VðtÞ ≥ Vð0Þe 0 ÐV ≥ − δ dr there exists a positive eigenvalue for the characteristic equa- t − δ dr t − δ dr tion and the equilibrium point is unstable. 0, ½IFNðtÞ ≥ ½IFNð0Þe 0 IFN ≥ 0, and KðtÞ ≥ Kð0Þe 0 K ≥ 0: According to the positivity of U R ðtÞ, VðtÞ, ½IFNðtÞ Theorem 2. The EE is locally stable if R0 > 1 and unstable if R0 < 1. Ð t KðtÞ and according to Equation (3), U I ðtÞ ≥ U I ð0Þ and − ðδI +qK KÞdr e 0 Ð ≥ 0; therefore, U I ðtÞ ≥ 0: Similarly, U E ðtÞ ≥ t Proof. In the same way and by calculating the characteristic − ðkE +qK KÞdr U EÐ ð0Þe 0 ≥0 and U H ðtÞ ≥ U H ð0Þ equation of J about its EE and according to the negative t − ð−kI V−kR ½IFN−δH Þdr eigenvalues obtained from Jacobian, the EE is stable if R0 e 0 ≥ 0. Therefore, all solutions are posi- > 1, otherwise it is unstable. tive. 3.2. Sensitivity Analysis. In this subsection, the sensitivity of 3.4. Existence of Solutions. In this subsection, the existence of R0 is investigated. Therefore, to check the R0 sensitivity, solutions is investigated with Theorem 4. SRα 0 = ð∂R0 /∂αÞðα/R0 Þ, we have Theorem 4. The proposed model with the initial conditions U H ð0Þ ≥ 0, U E ð0Þ ≥ 0, U I ð0Þ ≥ 0, U R ð0Þ ≥ 0, Vð0Þ ≥ 0, ½IFNð R SδH0 = 1, 0Þ ≥ 0, Kð0Þ ≥ 0, and Eð0Þ ≥ 0 has a unique solution. R SSH0 = 1, Proof. The model can be rewritten as R SkI0 = 1, f ðX Þ = AX + f 1 ðX Þ ð12Þ SRρV0 = 1, q K SK where X = ½U H ðtÞ U E ðtÞ U I ðtÞ U R ðtÞ VðtÞ½IFNðtÞ KðtÞT R SkE0 = , and A =diag ½ −δH , −kE , −δI , −δR , −δV , −δIFN ,−δK  ðkE δK + qK SK Þ ð11Þ and f 1 ðXÞ = ½ SH − kI U H V − kR U H ½IFNkI U H V − qK U E K R qK SK ðkE δK + δI δK + 2SK qK Þ kE U E − qK U I K kR U H ½IFN ρV U I aI U I SK + ΦK U I T . SδK0 = , ðkE δK + qK SK ÞðδK δI + qK SK Þ The function f 1 ðXÞ satisfies R SδV0 = −δV ðkE δK + qK SK ÞðδI δK + qK SK Þ,   j f 1 ðX 1 Þ − f 1 ðX 2 Þ j = j k I ðU H 2 V 2 − U H 1 V 1 Þj −qK kE δK SK δV + SK δI δK δV + 2δV S2K qK + jkI ðU H 1 V 1 − U H 2 V 2 Þj SqK = R0 ,    δV ðkE δK + qK SK ÞðδI δK + qK SK Þ + kR U H 2 ½IFN − U H 1 ½IFN    2 1 R0 −SK kE δK qK δV + qK δI δK δV + 2δV SK q2K + jqK ðU E2 K 2 − U E1 K 1 Þj SSK = : ð13Þ δV ðkE δK + qK SK ÞðδI δK + qK SK Þ + jkE ð−U E2 + U E1 Þj + jqK ðU I 2 K 2 − U I 1 K 1 Þj    It can be concluded that R0 is directly related to the + kR U H1 ½IFN1 − U H 2 ½IFN2  change of parameters δH , SH , kI , ρV , kE , and δK . That is, if the parameters δH , SH , kI ,and ρV increase by 10%, R0 also + jðρV + aI + ΦK ÞðU I 1 − U I 2 Þj changes accordingly (10%) and is inversely related to the parameters δV ,qK , and SK . Assumed that all parameters have the maximum Computational and Mathematical Methods in Medicine 7 value M which is a bounded value: subsection, the optimal control is used to vaccinate suscepti-   ble people and cure infected people. j f 1 ðX 1 Þ − f 1 ðX 2 Þj ≤ jMðU H 1 V 1 − U H 2 V 2 Þj + M U H 1 ½IFN1  4.1. Vaccination Strategy. In this subsection, the optimal − U H 2 ½IFN2 j + jMðU E1 K 1 − U E2 K 2 Þj control is applied to the innate immune response to convert + jMðU I 1 K 1 − U I 2 K 2 Þj + jMðU E1 − U E2 Þj healthy cells to resistant-to-infection cells. The vaccine + jMðU I 1 − U I 2 Þj: works with the innate immune to develop protection (immunity) to the disease. This subsection is aimed at min- ð14Þ imizing the number of healthy cells ðU H Þ by converting them to resistant-to-infection cells ðU R Þ with the control Replacing ðU H 1 V 1 − U H 2 V 2 Þ by ðU H 1 V 1 − U H 2 V 1 + input (vaccination) at a rate of ν in the optimal time. It is U H 2 V 1 − U H 2 V 2 Þ and similarly for others, obvious that injecting a vaccine into the susceptible person’s bodies takes time to vaccinate them (recovered) against the j f 1 ðX 1 Þ − f 1 ðX 2 Þj ≤ jMðU H 1 V 1 − U H 2 V 1 + U H 2 V 1 − U H 2 V 2 Þj disease. Unlike [13, 19], in the current paper, this time is   + M U H 1 ½IFN1 − U H 2 ½IFN1 + U H 2 ½IFN1 considered as a delay in moving susceptible people to the  recovered people group in the society, which is investigated − U H 2 ½IFN2 j + jMðU E1 K 1 − U E2 K 1 in Section 5 in detail. As long as the healthy cells have not + U E2 K 1 − U E2 K 2 Þj + jMðU I 1 K 1 − U I 2 K 1 become resistant-to-infection cells in the susceptible body, + U I 2 K 1 − U I 2 K 2 Þj + jMðU E1 − U E2 Þj the person has not been vaccinated against the disease, and + jMðU I 1 − U I 2 Þj there is still a risk of infection. The above description is given ð15Þ by the following model: Therefore, U_ H = SH − kI U H V − kR U H ½IFN − δH U H − νU H , ð18Þ j f 1 ðX 1 Þ − f 1 ðX 2 Þj ≤ jMV 1 ðU H 1 − U H 2 Þ + MU H 2 ðV 1 − V 2 Þj U_ E = kI U H V − kE U E − qK U E K, ð19Þ  + M½IFN1 ðU H 1 − U H 2 Þ   U_ I = kE U E − δI U I − qK U I K, ð20Þ + MU H 2 ½IFN − ½IFN  1 2 + jMK 1 ðU E1 − U E2 Þ + MU E2 ðK 1 − K 2 Þj U_ R = kR U H ½IFN − δR U R + νU H , ð21Þ + jMK 1 ðU I 1 − U I 2 Þ + MU I 2 ðK 1 − K 2 Þj V_ = ρV U I − δV V, ð22Þ + jMðU E1 − U E2 Þj + jMðU I 1 − U I 2 Þj: _  = a U − δ ½IFN, ½IFN ð23Þ ð16Þ I I IFN K_ = SK + ΦK U I − δK K: ð24Þ The maximum value of states is equal to R (bounded value); therefore, Figure 3 shows the diagram of the within-host dynamic model with vaccination. The susceptible individual’s body j f 1 ðX 1 Þ − f 1 ðX 2 Þj ≤ 2RMjU H 1 − U H 2 j + RMjV 1 − V 2 j does not have the virus, infected, or partially infected cells.   + RM½IFN1 − ½IFN2  + 2RMjU I 1 − U I 2 j Therefore, the goal is to convert healthy cells into resistant-to-infection cells using the vaccination strategy in + 2RMjK 1 − K 2 j + 2RMjU E1 − U E2 j the optimal time. ≤ 2RMjX 1 − X 2 j: Therefore, the optimal control problem is to minimize ð17Þ the cost (objective) function applied to Equations (18)–(24) given by It can be concluded that k f 1 ðX 1 Þ − f 1 ðX 2 Þk ≤ 2RMkX 1 ðt f   − X 2 k. Therefore, k f ðX 1 Þ − f ðX 2 Þk ≤ LkX 1 − X 2 k, where L   A3 2   J ν, t f = A1 U H − A2 U R + ν dt+∅ t f , ð25Þ = max ð2RM, kAkÞ. Thus, f is uniformly Lipschitz contin- 0 2 uous, and the solutions exist. in which a free terminal time is investigated, which gives a 4. Optimal Control of Within-Host minimum duration of a vaccination. t f represent the dura- Epidemic Model tion of the vaccination. ∅ is a positive increasing function such that lim ∅ðtÞ = +∞. In other words, the vaccination Optimal control techniques are used to improve optimal t⟶∞ strategies to control numerous kinds of diseases. The flexi- strategy is applied using the optimal control ðν∗ Þ in optimal bility and relative simplicity of optimal control techniques time ðt f ∗ Þ such that can lead to the development of the optimal strategies to con-     trol the monitoring and treating of various kinds of diseases J ν∗ , t f ∗ = min J ν, t f ∣ ν ∈ U 1 , t f ∈ ℝ+ , ð26Þ [14]. Therefore, the optimal control theory is applied to the innate immune dynamic in this part of the paper. In the first where ν is control and U 1 is the set of admissible controls 8 Computational and Mathematical Methods in Medicine defined by Theorem 5. Given an optimal control ν∗ , solutions U H ∗ ,   U E ∗ , U I ∗ , U R ∗ , V ∗ , ½IFN∗ , and K ∗ for the optimal control U 1 = U 1 is measurable, 0 ≤ ν ≤ νmax = 1, t ∈ 0, t f : problem Jðν∗ , t f ∗ Þ = min fJðν, t f Þ ∣ ν ∈ U 1 , t f ∈ ℝ+ g. Then, ð27Þ there are adjoint variables p1 , p2 , p3 , p4 , p5 , p6 , and p7 that satisfy The necessary conditions for optimality are expressed by " # p_ 1 = −A1 + p1 δH − kI V ½p2 − p1  − ½kR ½IFN  + ν½p4 − p1 , h   iT         ∂∅ ∇∅ − p tf δX f + H X t f , u t f , p t f , t f + δt f = 0, ð34Þ Xf ∂t f ð28Þ p_ 2 = p2 ½kE + qK K  − p3 kE , ð35Þ where Xðt f Þ or X f is vector of constant values (the final p_ 3 = p3 ½δI + qK K  − p5 ρV − p6 aI − p7 ΦK , ð36Þ value of all states is known). The gradient ð∇∅ X f Þ is the deriv- p_ 4 = A2 + p4 δR , ð37Þ ative of ∅ðt f Þ. t f is unknown but Xðt f Þ is defined; therefore, δX f = 0. Also, a stationarity condition on the control is as p_ 5 = kI U H ½p2 + p1  + p5 ρV , ð38Þ follows: p_ 6 = kR U H ½p1 − p4  + p6 δIFN , ð39Þ ∇H ν = 0, ð29Þ p_ 7 = p2 qK U E + p3 qK U I + p7 δK , ð40Þ which are expressed using the Hamiltonian of the system. with transversality conditions as The Hamiltonian for the control problem is defined as fol-   lows: pi t f = 0, i = 1, 2, 3, 4, 5, 6, 7: ð41Þ H = g + pT f , ð30Þ As a result, an optimal control ðνÞ is given by where f = x_ = ½U_ H , U_ E , U_ I , U_ R , V, _ K _ ½IFN, _ T and g = ½A1 U H ∗ ½p 1 − p 4  ν∗ = max min ,1 ,0 , ð42Þ U H − A2 U R + ðA3 /2Þν2 . A3 Therefore, and the optimal final time is given by A H = A1 U H − A2 U R + 3 ν2 + p1 U_ H + p2 U_ E + p3 U_ I + p4 U_ R 2 ∂ϕ  ∗  A _ _ _ + p5 V + p6 ½IFN + p7 K, t = −A1 U H ∗ + A2 U R ∗ − 3 ν∗ 2 : ð43Þ ∂t f f 2 ð31Þ Proof. Using Pontryagin’s principle, the adjoint equations and p1 , p2 , p3 , p4 , p5 , p6 , and p7 are the adjoint functions to and transversality conditions are obtained such that be determined suitably. Then, the necessary conditions are computed as follows: ∂H p_ 1 = − , ∂U H ∂H h _ iT _ , K_ :   ð44Þ x_ = = U H , U_ E , U_ I , U_ R , V, _ ½IFN ð32Þ p1 t f = 0, ∂p Then, the necessary conditions for optimality are com- ∂H puted, which are expressed by the mentioned Hamiltonian p_ 2 = − , ∂U E equations [38]. These equations include a linear, ordinary   ð45Þ differential equation that are shown by p2 t f = 0,   ∂H ∂H ∂H ∂H ∂H ∂H ∂H ∂H T ∂H p_ = − = − ,− ,− ,− ,− ,− ,− : p_ 3 = − , ∂x ∂U H ∂U E ∂U I ∂U R ∂V ∂½IFN ∂K ∂U I   ð46Þ ð33Þ p3 t f = 0, In order to derive the necessary conditions for this opti- ∂H mal control, the Pontryagin maximum principle [39] is p_ 4 = − , ∂U R applied to characterize the optimal control problem and   ð47Þ the optimal final time, which is given as the following p4 t f = 0, theorem. Computational and Mathematical Methods in Medicine 9 ∂H the optimal control and the state. This optimal control prob- p_ 5 = − , lem consists of the state system Equations (18)–(24) with ∂V ð48Þ   initial conditions at t = 0, the adjoint system Equations p5 t f = 0, (34)–(40) with the final conditions Equation (41), and the characterization of the optimal control Equation (42). ∂H p_ 6 = − , ∂½IFN ∗ U_ H = SH − kI U H ∗ V ∗ − kR U H ∗ ½IFN∗ − δH U H ∗   ð49Þ  U H ∗ ½p1 − p4  p6 t f = 0, − max min , 1 , 0 U H ∗, A3 ∗ ∂H U_ E = kI U H ∗ V ∗ − kE U E ∗ − qK U E ∗ K ∗ , p_ 7 = − , ∗ ∂K ð50Þ U_ I = kE U E ∗ − δI U I ∗ − qK U I ∗ K ∗ ,    p7 t f = 0, ∗ U_ R = kR U H ½IFN∗ − δR U R ∗ + max min U H ∗ ½p1 − p4  ,1 ,0 U H ∗, A3 and using the optimality conditions, the optimal control ∗ V_ = ρV U I ∗ − δV V ∗ , is solved as ∇H ν = 0 in the following form: _ ∗ = a U ∗ − δ ½IFN∗ , ½IFN IFN ! ! ! ! I I ∗ ∂U_ H ∂U_ E ∂U_ I ∂U_ R K_ = SK + ΦK U I ∗ − δK K ∗ , A3 ν + p 1 + p2 + p3 + p4 ∂ν ∂ν ∂ν ∂ν ð56Þ ! ! ! ∂V_ _  ∂½IFN ∂K_ and the adjoint system is as follows: + p5 + p6 + p7 ∂ν ∂ν ∂ν p_ 1 = −A1 + p1 δH − kI V ∗ ½p2 − p1  = A3 ν + p1 ð−U H Þ + p4 ðU H Þ = 0:   U H ∗ ½p1 − p4  − kR ½IFN∗ + max min , 1 , 0 ½p4 − p1 , ð51Þ A3 p_ 2 = p2 ½kE + qK K ∗  − p3 kE , Thus, the control law is written as p_ 3 = p3 ½δI + qK K ∗  − p5 ρV − p6 aI − p7 ΦK , p_ 4 = A2 + p4 δR , U H ∗ ½p1 − p4  ν∗ = : ð52Þ p_ 5 = kI U H ∗ ½p2 + p1  + p5 ρV , A3 p_ 6 = kR U H ∗ ½p1 − p4  + p6 δIFN , Using the property of the control space, the following equa- p_ 7 = p2 qK U E ∗ + p3 qK U I ∗ + p7 δK , tion is obtained: ð57Þ 8 > U H ∗ ½p 1 − p 4  > > 1 ≥ 1, with p1 ðt f Þ = 0, p2 ðt f Þ = 0, p3 ðt f Þ = 0, p4 ðt f Þ = 0, p5 ðt f Þ = 0, > > A3 > > p6 ðt f Þ = 0, p7 ðt f Þ = 0, U H ð0Þ = U H 0 , U E ð0Þ = U E0 , U I ð0Þ = > < U ∗ ½p − p  ∗ 1 4 U ∗ ½p − p4  U I 0 , U R ð0Þ = U R0 , Vð0Þ = V 0 , ½IFNð0Þ = ½IFN0 , Kð0Þ = K 0 , ν = H 0< H 1 < 1, ð53Þ > A3 > A3 and T ∗ which can be rewritten as in Equation (78). So, the > > > > U H ∗ ½p 1 − p 4  proof is completed. > > :0 ≤ 0, A3 After examining the effect of the vaccine on the person’s body and the control strategy applied to the model for the so, the optimal control is characterized as best effect, it is necessary to investigate the efficiency of the optimal antiviral treatment strategy in an optimal time, U H ∗ ½p1 − p4  ν∗ = max min ,1 ,0 : ð54Þ which is considered in the following subsection. A3 4.2. Antiviral Treatment Strategy. This subsection is aimed at Now, according to Equation (28) and knowing δX f = 0, destroying the virus particles and converting healthy cells the optimal final time ðt f ∗ Þ is computed using the following into resistant-to-infection cells simultaneously using the equation: 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 t = −H U H ∗ , U E ∗ , U I ∗ , U R ∗ , V ∗ , ½IFN∗ , K ∗ , ν∗ , t f ∗ : eradicate the disease by converting healthy cells to ∂t f f resistant-to-infection cells at the rate of T 1 and destroying ð55Þ virus particles at a rate of T 2 . Figure 4 shows the conceptual flow diagram of the innate immune dynamics with the pro- Therefore, using the characterization of the optimal con- posed control action (antiviral treatment). Moreover, the trol, the following optimality system is presented to obtain dynamic of the innate immune with antiviral treatment is 10 Computational and Mathematical Methods in Medicine ∂H (IFN) p_ = − , ð68Þ ∂x SH aIUI δIFN(IFN) " # h   iT         ∂∅ kR kI kEUE ∇∅ − p tf δX f + H X t f , u t f , p t f , t f + δt f = 0, UR UH UE UI Xf ∂t f qk ð69Þ qk δIUI δRUR UH 1 δHUH in which t f is unknown while Xðt f Þ is fixed; therefore, δX f 𝜌VUI SK 𝛷KUI = 0. Then, the Hamiltonian for the control problem is V K defined as V δVV δKK 2 A4 A H = A 1 U H − A2 U R + A3 V + T 2 + 5 T 2 2 + p1 U_ H Figure 4: The diagram of the innate immune dynamics with the 2 1 2 proposed control action (antiviral treatment). _  + p K, + p2 U_ E + p3 U_ I + p4 U_ R + p5 V_ + p6 ½IFN _ 7 considered as ð70Þ U_ H = SH − kI U H V − kR U H ½IFN − δH U H − T 1 U H , ð58Þ where p1 , p2 , p3 , p4 , p5 , p6 , and p7 are also the adjoint func- tions to be obtained. As discussed earlier, the following the- U_ E = kI U H V − kE U E − qK U E K, ð59Þ orem is obtained to apply Pontryagin’s principle to the Hamiltonian. U_ I = kE U E − δI U I − qK U I K, ð60Þ Theorem 6. Given optimal control pair T 1 ∗ , T 2 ∗ and U_ R = kR U H ½IFN − δR U R + T 1 U H , ð61Þ solutions U H ∗ , U E ∗ , U I ∗ , U R ∗ , V ∗ , ½IFN∗ , and K ∗ for the optimal control problem JðT 1 ∗ , T 2 ∗ , t f ∗ Þ = min fJðT 1 , T 2 , V_ = ρV U I − δV V − T 2 V, ð62Þ t f Þ ∣ T 1 , T 2 ∈ U, t f ∈ ℝ+ g, let U H ∗ , U E ∗ , U I ∗ , U R ∗ , V ∗ , _  = a U − δ ½IFN, ½IFN ð63Þ ½IFN∗ , and K ∗ be optimal state solutions with an associated I I IFN optimal control pair ðT 1 , T 2 Þ. There are adjoint variables that K_ = SK + ΦK U I − δK K: ð64Þ satisfy the following equations: To reach the proposed goal and minimize the cost of the p_ 1 = −A1 + p1 δH + kI V ½p1 − p2  + ½kR ½IFN  + ν½p1 − p4 , optimal antiviral treatment strategy within an optimal dura- tion, the following objective function with free terminal time p_ 2 = p2 ½kE + qK K  − p3 kE , has to be minimized: p_ 3 = p3 ½δI + qK K  − p5 ρV − p6 aI − p7 ΦK , ðt f     A A   J T 1, T 2, t f = A1 U H − A2 U R + A3 V + 4 T 1 2 + 5 T 2 2 dt+∅ t f , p_ 4 = A2 + p4 δR , 0 2 2 ð65Þ p_ 5 = −A3 + kI U H ½p1 − p2  + p4 ½δv + T 2 , p_ 6 = kR U H ½p1 − p4  + p6 δIFN , in which t f represent the duration of the antiviral treatment. ∅ is a positive increasing function such that lim ∅ðtÞ = + p_ 7 = p2 qK U E + p3 qK U I + p7 δK , t⟶∞ ∞. Also, to reach the optimal controls ðT 1 ∗ , T 2 ∗ Þ and an ð71Þ optimal terminal time ðt f ∗ Þ, we considered with transversality conditions:     J T 1 ∗ , T 2 ∗ , t f ∗ = min J T 1 , T 2 , t f ∣ ðT 1 , T 2 Þ ∈ U 2 , t f ∈ ℝ+ ,   ð66Þ pi t f = 0, i = 1, 2, 3, 4, 5, 6, 7: ð72Þ where U 2 is the set of admissible controls defined by Furthermore, the optimal controls ðT 1 , T 2 Þ are given by U 2 = fðT 1 , T 2 Þ ∣ T 1 , T 2 are measurable,   U H ∗ ½p 1 − p 4  0 ≤ T 1 ≤ 1, 0 ≤ T 2 ≤ 1, t ∈ 0, t f : ð67Þ T 1 ∗ = max min ,1 ,0 , A4 ð73Þ ∗ p5 V ∗ Then, the necessary conditions for optimality are written T 2 = max min ,1 ,0 , as follows: A5 Computational and Mathematical Methods in Medicine 11 ∗ and the optimal final time is given by U_ H = SH − kI U H ∗ V ∗ − kR U H ∗ ½IFN∗ − δH U H ∗    U H ∗ ½p 1 − p 4  H U H ∗ , U E ∗ , U I ∗ , U R ∗ , V ∗ , ½IFN ∗ , K ∗ , T 1 ∗ , T 2 ∗ , t f ∗ − max min , 1 , 0 U H ∗, A4 ∂ϕ  ∗  + t = 0: ∂t f f ∗ U_ E = kI U H ∗ V ∗ − kE U E ∗ − qK U E ∗ K ∗ , ð74Þ ∗ U_ I = kE U E ∗ − δI U I ∗ − qK U I ∗ K ∗ , Proof. The adjoint equation and transversality conditions can ∗ be obtained using Pontryagin’s principle as similar to Equa- U_ R = kR U H ½IFN∗ − δR U R ∗ tions (44)–(50). Then, by using the optimality conditions,  U H ∗ ½p 1 − p 4  the optimal controls are obtained using ∇H T 1 = 0 and ∇T 2 = 0 H + max min , 1 , 0 U H ∗, A4 . As a result, the control inputs are computed as U ∗ ½p − p 4  ∗ p5 V ∗ T1 ∗ = H 1 , V_ = ρV U I ∗ − δV V ∗ − max min , 1 , 0 V ∗, A3 A5 ð75Þ ∗ p V∗ _ ∗ = a U ∗ − δ ½IFN∗ , ½IFN T2 = 5 : I I IFN A5 ∗ K_ = SK + ΦK U I ∗ − δK K ∗ , ð79Þ Using the property of the control space, we obtain and the adjoint equations are as follows: 8 > U H ∗ ½p1 − p4  > >1 ≥ 1, > > A4 p_ 1 = −A1 + p1 δH − kI V ∗ ½p2 − p1  > >   > < U ∗ ½p − p  U H ∗ ½p 1 − p 4  ∗ 1 4 U ∗ ½p − p 4  − kR ½IFN∗ + max min ,1 ,0 T1 = H 0< H 1 < 1, A4 > A3 > A4 > > > >  ½p4 − p 1 , > > U H ∗ ½p1 − p4  :0 ≤ 0, A4 8 ð76Þ p_ 2 = p2 ½kE + qK K ∗  − p3 kE , > p5 V ∗ > > 1 ≥ 1, > > A5 p_ 3 = p3 ½δI + qK K ∗  − p5 ρV − p6 aI − p7 ΦK , > > > < p V∗ ∗ 5 p V∗ T2 = 0< 5 < 1, p_ 4 = A2 + p4 δR , > > A5 A5 > > > > p5 V ∗ > > p_ 5 = −A3 + kI U H ∗ ½p1 − p2  :0 ≤ 0:   A5 p5 V ∗ + p4 δv + max min ,1 ,0 , A5 So, the optimal control is characterized as U H ∗ ½p 1 − p 4  p_ 6 = kR U H ∗ ½p1 − p4  + p6 δIFN , ∗ T 1 = max min ,1 ,0 , A4 p_ 7 = p2 qK U E ∗ + p3 qK U I ∗ + p7 δK , ð80Þ ∗ ð77Þ p5 V T 2 ∗ = max min ,1 ,0 , where p1 ðt f Þ = 0, p2 ðt f Þ = 0, p3 ðt f Þ = 0, p4 ðt f Þ = 0, p5 ðt f Þ = 0, A5 p6 ðt f Þ = 0, p7 ðt f Þ = 0, U H ð0Þ = U H 0 , U E ð0Þ = U E0 , U I ð0Þ = and according to Equation (74) and knowing that δX f = 0, the U I 0 , U R ð0Þ = U R0 , Vð0Þ = V 0 , ½IFNð0Þ = ½IFN0 , and Kð0Þ optimal final time ðt f ∗ Þ is computed using the following = K 0 . Then, the proof is completed. equation: Remark 7. In the vaccination strategy, the healthy cells con- ∂ϕ  ∗  A A vert to resistant-to-infection cells using the proposed con- t = −A1 U H ∗ − A2 U R ∗ + A3 V ∗ + 4 T 1 ∗2 + 5 T 2 ∗2 : troller and their number converges to zero when t ⟶ ∞ ∂t f f 2 2 ðU H ⟶ 0Þ and the rest remain at their initial zero value, ð78Þ while in the antiviral treatment strategy, in addition to converting healthy cells to resistant-to-infection cells, the Therefore, using the characterization of the optimal con- virus particles decrease and converge to zero using the pro- trol, we have the following optimality system to obtain the posed controller. Meanwhile, the number of natural killers optimal states: converges to the initial constant value. Then, after a while, 12 Computational and Mathematical Methods in Medicine Equation (59) becomes U_ E = −ðkE + qK KÞU E in which ðkE uv + qK KÞ is constant and is considered as first-order linear 𝜇bE+𝜇bS+𝜇bR+𝜇b(1-𝜌)I 𝜇b𝜌I differential equation that yields U E = U E0 e−ðkE +qK KÞt . As a result, U E converges exponentially to zero. 𝛽SI kE gI S E I R In this section, the effect of the optimal control on the within-host model is examined. At the same time, it is nec- essary to explore the controller’s influence on the between- 𝜇dS 𝜇dE 𝜇dI u 𝜇dR hosts model (SEIR epidemic model). Thus, the impact of the delay computed in this section is discussed in the follow- Figure 5: Conceptual flow diagram of the influenza SEIAR ing section in detail. dynamics with the proposed controllers. Table 2: Initial values of the SEIR epidemic model. 5. SEIR Dynamic Model State variable Initial value In this section, a SEIR epidemic model is introduced with two additional properties. First is two time-delay that aroused S0 9375 from this fact that vaccination and antiviral treatment take E0 125 time to impact a person’s body. Second is a nonlinear contact I0 313 rate (transmission rate) function that is followed by this fact R0 0 that the contact rate is varying throughout day and night. In this regard, the susceptible and infected individuals move to N0 10,000 the recovered people group with time delays t 1 and t 2 (optimal times of vaccination and antiviral treatment within the body, Table 3: Initial values of innate immune response. which is computed in the previous section), respectively. The proposed SEIR epidemic model is presented as State Susceptible Infected Exposed variable person person person S_ ðt Þ = −βðt ÞSðt ÞI ðt Þ + μd ð3ð1 − ρI ðt ÞÞ − Sðt ÞÞ, UH0 5 × 108 5 × 108 5 × 108 E_ ðt Þ = βðt ÞSðt ÞI ðt Þ − ðμd + κÞEðt Þ, UE0 0 3 × 103 0 ð81Þ UI0 0 1․5 × 10 3 0 I_ ðt Þ = κEðt Þ − αI ðt Þ, UR0 0 0 0 R_ ðt Þ = gI ðt Þ − μd Rðt Þ, V0 0 10−3 10−3 ½IFN0 0 0 0 where α = ðg − 3μd ρ + μd Þ. As shown in Figure 5, SðtÞ repre- 5 5 sents the number of susceptible individuals who are vulnerable K0 8 × 10 8 × 10 8 × 105 and not infected, EðtÞ or exposed individual denotes the num- ber of people that are infected but not infectious yet, IðtÞ sig- nifies the population of infectious people with symptoms that viduals in the population per unit time” [40], the contact can transmit the disease to others, and RðtÞ is the number of rate must be more than one during the day due to the completely recovered (immune) people. N is the total popula- more relations between people in the community. More- tion size ðSðtÞ + EðtÞ + IðtÞ + RðtÞ = NÞ. μd is the rate of deaths over, late at night and early in the morning, the transmis- from causes unrelated to the infection, and μb is the birth rate sion rate is less than one. Therefore, with trial and error and taking into account the assumptions below, the contact in each group, which is triple the mortality rate ðμb = 3μd Þ. t 1 is the vaccination time of susceptible people, and t 2 is the rate is obtained as βðtÞ = 0:5ð3 − sin ð2πtÞ2 − sin ð2πtÞ − cos recovery time of infected people. In the case of births, in the ð2πtÞÞ. susceptible, exposed, and recovered groups, newborns are sus- ceptible to the disease and add to the susceptible group at the (a) 00:00 midnight: βð00 : 00Þ = 1 rate of μb . Fraction ρ of the infected individuals’ infants goes (b) 06:00 morning: βð06 : 00Þ = 0:5 to the infected compartment, and the others (fraction ð1 − ρÞ of infected individuals’ infants) move to the susceptible group. (c) 12:00 midday: βð12 : 00Þ = 2 The exposed people get infected at a rate of κ, and infected peo- (d) 18:00 evening: βð18 : 00Þ = 1:5 ple leave their compartment at the rate of α. The contact rate is βðtÞ that is purposefully variable Thus, the contact rate is less than one late at night and at throughout day and night. It is also obvious that this rate daybreak and is greater than one in the middle of the day at night is higher than during the day, which makes it because there are more social gatherings and offices and closer to the real world and is more effective in practice schools are also open. [40]. According to the basic definition of βðtÞ that is In the next section, an optimal control is applied to the defined as “the average number of contacts between indi- proposed model to overcome the disease outbreak. Computational and Mathematical Methods in Medicine 13 10000 9000 8000 7000 Number of humans 6000 5000 4000 3000 2000 1000 0 0 5 10 15 20 25 Time (days) Susceptible with 𝛽(t) Susceptible with 1.4×𝛽(t) Susceptible with 0.6×𝛽(t) (a) 4500 4000 3500 3000 Number of humans 2500 2000 1500 1000 500 0 0 5 10 15 20 25 Time (days) Exposed with 𝛽(t) Exposed with 1.4×𝛽(t) Exposed with 0.6×𝛽(t) (b) Figure 6: Continued. 14 Computational and Mathematical Methods in Medicine 4000 3500 3000 Number of humans 2500 2000 1500 1000 500 0 0 5 10 15 20 25 Time (days) Infected with 𝛽(t) Infected with 1.4×𝛽(t) Infected with 0.6×𝛽(t) (c) 14000 12000 10000 Number of humans 8000 6000 4000 2000 0 0 5 10 15 20 25 Time (days) Recovered with 𝛽(t) Recovered with 1.4×𝛽(t) Recovered with 0.6×𝛽(t) (d) Figure 6: The impact of the contact rate changes on the (a) susceptible, (b) exposed, (c) infected, and (d) recovered people. 6. Optimal Control of the Between-Hosts ral treatment in 25 days, as shown in Figure 5. The SEIR Epidemic Model dynamic model in the presence of the controller is pre- sented as 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 of uv ð0 ≤ uv ≤ 1 S_ ðt Þ = −βðt ÞSðt ÞI ðt Þ + μd ð3ð1 − ρI ðt ÞÞ − Sðt ÞÞ − uv ðt − t 1 ÞSðt Þ, Þ in 25 days and second, moving the infected people to the recovered group at a rate of uτ ð0 ≤ uτ ≤ 1Þ by antivi- ð82Þ Computational and Mathematical Methods in Medicine 15 10000 9000 8000 7000 Number of humans 6000 5000 4000 3000 2000 1000 0 0 5 10 15 20 25 Time (days) Susceptible with g Susceptible with 2×g Susceptible with 0.5×g (a) 4500 4000 3500 3000 Number of humans 2500 2000 1500 1000 500 0 0 5 10 15 20 25 Time (days) Exposed with g Exposed with 2×g Exposed with 0.5×g (b) Figure 7: Continued. 16 Computational and Mathematical Methods in Medicine 6000 5500 5000 4500 4000 Number of humans 3500 3000 2500 2000 1500 1000 500 0 0 5 10 15 20 25 Time (days) Infected with g Infected with 2×g Infected with 0.5×g (c) 14000 12000 10000 Number of humans 8000 6000 4000 2000 0 0 5 10 15 20 25 Time (days) Recovered with g Recovered with 2×g Recovered with 0.5×g (d) Figure 7: The impact of the recovery rate changes on the (a) susceptible, (b) exposed, (c) infected, and (d) recovered people. E_ ðt Þ = βðt ÞSðt ÞI ðt Þ − ðμd + κÞEðt Þ, ð83Þ R_ ðt Þ = gI ðt Þ − μd Rðt Þ + uv ðt − t 1 ÞSðt Þ + uτ ðt − t 2 ÞI ðt Þ, ð85Þ I_ ðt Þ = κEðt Þ − αI ðt Þ − uτ ðt − t 2 ÞI ðt Þ, ð84Þ where uv ðtÞ = 0 for t ∈ ½−t 1 , 0 and uτ ðtÞ = 0 for t ∈ ½−t 2 , 0. Computational and Mathematical Methods in Medicine 17 9000 8000 7000 6000 Number of humans 5000 4000 3000 2000 1000 0 0 5 10 15 20 25 Time (days) Susceptible with K Susceptible with 2×K Susceptible with 0.5×K (a) 6000 5500 5000 4500 4000 Number of humans 3500 3000 2500 2000 1500 1000 500 0 0 5 10 15 20 25 Time (days) Exposed with K Exposed with 2×K Exposed with 0.5×K (b) Figure 8: Continued. 18 Computational and Mathematical Methods in Medicine 4500 4000 3500 3000 Number of humans 2500 2000 1500 1000 500 0 0 5 10 15 20 25 Time (days) Infected with K Infected with 2×K Infected with 0.5×K (c) 14000 12000 10000 Number of humans 8000 6000 4000 2000 0 0 5 10 15 20 25 Time (days) Recovered with K Recovered with 2×K Recovered with 0.5×K (d) Figure 8: The impact of the κ rate changes on the (a) susceptible, (b) exposed, (c) infected, and (d) recovered people. Computational and Mathematical Methods in Medicine 19 To this end, the optimal control strategy is used to Notice that the control uv ðt − t 1 Þ can only influence achieve these goals. The objective function is considered the state variable on the time interval ½t f − t 1 , t f . Thus, as it suffices to compute the optimal control on the interval ½0, t f − t 1 . Therefore, we have ðt f     B4 B J uv ðt Þ, uτ ðt Þ, t f = B1 I ðt Þ + B2 Sðt Þ − B3 Rðt Þ + uv ðt Þ2 + 5 uτ ðt Þ2 dt, 0 2 2 ð86Þ  ∂H ∂H  +  = 0, ð89Þ where B1 , B2 , and B3 are gains of the infected, susceptible, ∂uv ðt Þ ∂uv ðt − t 1 Þt+t 1 and recovered individuals. Also, B4 and B5 are the gains of control inputs uv , and uτ , respectively. Moreover, t f and Xðt f Þ are fixed. Tacking the Hamiltonian yields which yields 8 S∗ t + t 1 Þ½p1 ðt + t 1 Þ − p4 ðt + t 1 Þ   < uv ∗ ð t Þ = ð B4 B > H = B 1 I ðt Þ + B 2 S ð t Þ − B 3 R ð t Þ + u ðt Þ2 + 5 uτ ðt Þ2 , t ∈ 0, t f − t 1 , 2 v 2 ð87Þ B4 _ _ _ _ + p1 Sðt Þ + p2 Eðt Þ + p3 I ðt Þ + p4 Rðt Þ, > : ∗   uv ðt Þ = 0, t ∈ t f − t1 , t f : where p1 :p2 ,p3 , and p4 are adjoint functions that satisfy the ð90Þ adjoint equations. Pontryagin’s maximum principle is applied to the Hamiltonian H to derive the necessary con- ditions for the optimal control. The adjoint equation and Similarly, transversality conditions are obtained by using Pontrya- gin’s maximum principle such that  ∂H ∂H  +  =0 ð91Þ ∂H ∂uτ ðt Þ ∂uτ ðt − t 2 Þt+t p_ 1 = − = −B2 + βðt ÞI ðt Þ½p1 − p2  + uv ðt − t 1 Þ½p1 − p4  2 ∂Sðt Þ   + p1 μd , p1 t f = 0, that can be written as ∂H   p_ 2 = − = κ ðp2 − p3 Þ + p2 μd , p2 t f = 0, ∂Eðt Þ 8 ∂H I ∗ t + t 2 Þ½p3 ðt + t 2 Þ − p4 ðt + t 2 Þ   p_ 3 = − = −B1 + βðt ÞSðt Þ½p1 − p2  + uτ ðt − t 2 Þ½p3 − p4  < uτ ∗ ð t Þ = ð > , t ∈ 0, t f − t 2 , ∂I ðt Þ B5   > : ∗   + p3 α − p4 g, p3 t f = 0, uτ ðt Þ = 0, t ∈ t f − t2, t f , ∂H   ð92Þ p_ 4 = − = B3 + p4 μd , p4 t f = 0: ∂Rðt Þ ð88Þ which yields 8 > S∗ ðt + t 1 Þ½p1 ðt + t 1 Þ − p4 ðt + t 1 Þ > > 1 ≥ 1, > > B4 > > > < S∗ ðt + t Þ½p ðt + t Þ − p ðt + t Þ S∗ ðt + t 1 Þ½p1 ðt + t 1 Þ − p4 ðt + t 1 Þ uv ∗ ð t Þ = 1 1 1 4 1 0< < 1, > > B4 B4 > > > > > > S∗ ðt + t 1 Þ½p1 ðt + t 1 Þ − p4 ðt + t 1 Þ :0 ≤ 0, B4 8 ð93Þ > I ∗ ðt + t 2 Þ½p3 ðt + t 2 Þ − p4 ðt + t 2 Þ > >1 ≥ 1, > > B5 > > > < I ∗ ðt + t Þ½p ðt + t Þ − p ðt + t Þ ∗ 2 3 2 4 2 I ∗ ðt + t 2 Þ½p3 ðt + t 2 Þ − p4 ðt + t 2 Þ uτ ð t Þ = 0< < 1, > > B5 B5 > > > > I ∗ ðt + t 2 Þ½p3 ðt + t 2 Þ − p4 ðt + t 2 Þ > > :0 ≤ 0: B5 20 Computational and Mathematical Methods in Medicine So, the optimal controls are given as 50 S∗ ðt + t 1 Þ½p1 ðt + t 1 Þ − p4 ðt + t 1 Þ   uv ∗ ðt Þ = max min ,1 ,0 , t ∈ 0, t f − t 1 , B4 0 Y–0.58881 I ∗ ðt + t 2 Þ½p3 ðt + t 2 Þ − p4 ðt + t 2 Þ   uτ ∗ ðt Þ = max min ,1 ,0 , t ∈ 0, t f − t 2 : B5 ð94Þ –50 Therefore, using the characterization of the optimal con- trol and considering Equations (90) and (92), the following –100 optimality system is obtained as ∗ S_ ðt Þ = −βðt ÞS∗ ðt ÞI ∗ ðt Þ + μd ð3ð1 − ρI ∗ ðt ÞÞ − S∗ ðt ÞÞ − uv ∗ ðt ÞS∗ ðt Þ, ∗ –150 E_ ðt Þ = βðt ÞS∗ ðt ÞI ∗ ðt Þ − ðμd + κÞE∗ ðt Þ, 0 2 4 6 8 10 12 ∗ I_ ðt Þ = κE∗ ðt Þ − gI ∗ ðt Þ − uτ ∗ ðt ÞI ∗ ðt Þ, Time (days) ∗ R_ ðt Þ = gI ∗ ðt Þ − μd R∗ ðt Þ + uv ∗ ðt ÞS∗ ðt Þ + uτ ∗ ðt ÞI ∗ ðt Þ: t1 ð95Þ Figure 9: The optimal final time. The adjoint equations are also obtained as values of the SEIR epidemic model states and innate p_ 1 = −B2 + βðt ÞI ∗ ðt Þ½p1 − p2  + u∗v ðt − t 1 Þ½p1 − p4  + p1 μd , immune response model are shown in Tables 2 and 3, p_ 2 = κ ðp2 − p3 Þ + p2 μd , respectively. This section is categorized as follows: p_ 3 = −B1 + βðt ÞS ðt Þ½p1 − p2  + u∗τ ðt − t 2 Þ½p3 − p4  + p3 α − p4 g, ∗ p_ 4 = B3 + p4 μd , (a) The influence of model parameters ð96Þ (b) The effect of the optimal controller on the within- host model. This section is divided into the following where p1 ðt f Þ = 0, p2 ðt f Þ = 0, p3 ðt f Þ = 0, p4 ðt f Þ = 0, p5 ðt f Þ = 0, two subsections: p6 ðt f Þ = 0, p7 ðt f Þ = 0, Sð0Þ = S0 , Eð0Þ = E0 , Ið0Þ = I 0 , and Rð (i) Within-host vaccination 0Þ = R0 . (ii) Within-host therapy Remark 8. According to the proposed controller, susceptible (c) The effect of the optimal controller on the between- and infected persons converge to zero as ⟶∞ hosts model. This section is divided into the follow- ðSðtÞ, IðtÞ ⟶ 0Þ; then, Equation (83) becomes a first-order ing two subsections: _ = −ðμ + κÞEðtÞ in which linear differential equation as EðtÞ d ðμd + κÞ is a positive constant and yields EðtÞ = E0 e−ðμd +κÞt . (i) Between-hosts vaccination Thereupon, EðtÞ converges exponentially to zero indirectly (ii) Between-hosts therapy 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 7.1. Influence of Model Parameters. This section is aimed at be illustrated. Therefore, the simulation results of the pro- conducting the sensitivity analysis to identify the factors that posed work are outlined in the next section. noticeably affect the number of people engaged in influenza. Therefore, the changes in the value of parameters βðtÞ, g, 7. Simulation Results and k are investigated, respectively. The tremendous influ- ence of contact rate cannot be neglected. The susceptible In this section, a numerical simulation is given to prove the people get infected if they have more contact with infected accuracy of the results obtained in the previous sections. people; as a result, the epidemic widely spreads. Thus, the First, the effect of the optimal controllers on the innate contact rate has to be considered as a critical parameter in immune response and then on the SEIR epidemic model is the epidemic containment process. To this end, the value assessed. The value of innate immune response parameters of the contact rate is changed from 0:6 βðtÞ to 1:4 βðtÞ to ðkI , kR , kE , qK , ΦK , aI , ρV , δH , δI , δR , δIFN , δK , δV , SH , SK Þ is observe its impact on the population size in each engaged based on experimental data from Hernandez-Vargas and group. As shown in Figures 6(a)–6(d), by increasing in the Meyer-Hermann [7], and the parameter values of the SEIR value of contact rate to 1:4 βðtÞ, as a result, the more contact epidemic model are κ = g = ρ = 0:3 and μd = 0:02. The initial between people, the number of susceptible people reduces Computational and Mathematical Methods in Medicine 21 5.5 ×10 5 ×10 8 8 5 4.5 4.5 4 4 3.5 3.5 Number of cells 3 Number of cells 3 2.5 2.5 2 2 1.5 1.5 1 1 0.5 0.5 0 0 0 1 2 3 4 5 6 7 8 0 1 2 3 4 5 6 7 8 Time (days) Time (days) UH without controller UR without controller UH with controller UR with controller (a) (b) Figure 10: Evolution of healthy cells and resistant-to-infection cells with and without the controller. regarding increment in the number of infected people, while 40 the number of susceptible people does not decrease more because of the less contact between infected individuals ð 30 0:6 βðtÞÞ. Obviously, the number of exposed and infected people decreases followed by nearly half decrement in the 20 contact rate with susceptible people. Thus, few people get 10 infected; as a result, few people have recovered. In Figure 7, the impact of changes in the value of g, Y 0.93807 0 which is the recovery rate of infected people, is presented. As shown in Figure 7(a), when the rate of g is doubled, a –10 greater number of susceptible people have recovered; as a result, fewer people get infected and the number of suscepti- –20 ble people is increased. Moreover, the number of infected –30 people increases when the recovery rate decreases and decreases with increment in the value of recovery rate –40 because more infected people are recovered. As a result, it 0 1 2 3 4 5 6 7 8 is obvious that the number of recovered people is grown up followed by the rise in the recovery rate. Time (days) According to Figure 8, the exposed people are infected at t2 a rate of κ; therefore, increment in this rate leads to growth in the infected people number; as a result, there is growth Figure 11: The optimal final time. in the number of recovered people. It is obvious that the number of susceptible people is decreased when κ increases. As observed, for different values of βðtÞ, g,and k, the dis- ease does not disappear in about 25 days; therefore, it is nec- become immune to the infection following vaccination essary to use vaccination and antiviral treatment strategies to because it takes time for all healthy cells (within the suscep- overcome the epidemic. Thus, in the next sections, the tible individual’s body) to become resistant-to-infection model in the presence of vaccination and antiviral treatment cells. Moreover, the destruction of the virus (and at the same is simulated compared to the model without the controller. time the transformation of healthy cells into resistant-to- infection cells) within the infected person’s body is not 7.2. Optimal Control Results in the Within-Host Model. In immediate with antiviral therapy; therefore, the infected per- the real world, susceptible people will not immediately son is treated after a while. To this end, two-time delay (the 22 Computational and Mathematical Methods in Medicine ×108 5 3.5 ×10 8 4.5 4 3 3.5 2.5 Number of cells Number of cells 3 2 2.5 2 1.5 1.5 1 1 0.5 0.5 0 0 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.5 1 1.5 2 2.5 3 3.5 4 Time (days) Time (days) UH without controller UR without controller UH with controller UR with controller (a) (b) Figure 12: Evolution of healthy cells and resistant-to-infection cells with and without the controller. length of vaccination and antiviral treatment) is considered 7.2.2. Within-Host Therapy. There are virus particles (as a in the SEIR model. 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 7.2.1. Within-Host Vaccination. Vaccine injection into the process, the infected and partially infected cells and virus susceptible individual’s body makes the healthy cells particles are killed, and the healthy cells become resistant- resistant-to-infection. Consequently, the susceptible individ- to-infection simultaneously, which also takes time. ual, which includes healthy cells, will be recovered against In Figure 11, the final optimal time to eradicate disease the infection. Converting the healthy cells into the computed in Section 4 is shown. As shown in resistant-to-infection cells takes time, which is considered Figures 12(a) and 12(b), using antiviral drugs, healthy cells as delay in recovering process of susceptible individuals in become resistant-to-infection cells. The drug’s effect on the SEIR model. In Section 4, initially, the optimal time of infected and partially infected cells and viruses can also be healthy cells converting into resistant-to-infection cells in seen in Figures 13 and 14, respectively. It is clear that the the susceptible individual’s body has been obtained. Numer- infected and partially infected cells are destroyed within six ical simulations suggest eight days as the optimal time to days, and the number of virus particles is also converged to convert the healthy cells into the resistant ones in a suscep- zero, which indicates that the infected person is cured. tible person’s body (see Figure 9). The proposed optimal Moreover, the evolution of IFN-I molecules and natural time will be considered as a delay of moving susceptible peo- killer cells with and without the controller is shown in ple to the recovered people group in the SEIR model in the Figures 15(a) and 15(b). next section. Figures 10(a) and 10(b) show the effect of control input (vaccination) on the number of healthy and resistant-to- 7.3. Optimal Control Results between Hosts. In this section, infection cells. There is a surge in the number of resistant- we examine the effect of the optimal controller on the SEIR to-infection cells during vaccination in eight days. After vac- epidemic model, which was discussed in Section 7, in detail. cination, it can be deduced that rising in the resistant-to- First, the changes in the number of susceptible, exposed, infection cell numbers is associated with the reduction in infected, and recovered individuals in the absence of appro- the healthy cell numbers, indicating the controller efficacy. priate control are investigated. Then, the effect of the vacci- As mentioned, as the number of resistant-to-infection cells nation strategy and antiviral treatment on individuals in a increases moment by moment, the susceptible person even- community is examined. For the numerical simulations, tually becomes immune to the disease after eight days the delay differential equation solver (dde23) is used. This because the healthy cells in his/her body have all become solver does not solve the advance equations ðt + t i Þ directly; resistant. therefore, it is necessary to convert the advance equations Computational and Mathematical Methods in Medicine 23 ×108 ×107 2 2 1.8 1.8 1.6 1.6 1.4 1.4 Number of cells Number of cells 1.2 1.2 1.0 1.0 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Time (days) Time (days) UE without controller UI without controller UE with controller UI with controller (a) (b) Figure 13: Evolution of partially infected cells and infected cells with and without the controller. ×104 3.5 3 2.5 Concentration per mg 2 1.5 1 0.5 0 0 1 2 3 4 5 6 Time (days) V without controller V with controller Figure 14: Evolution of virus particle with and without the controller. into delay differential equations. To this end, a new descends and converts to infected people. Eventually, it reversed-time variable ς = t f − t is defined and new adjoint remains stationary at a constant value (less than 2,000 peo- variables λi ðςÞ = pi ðt f − ςÞ are introduced. Thus, this new ple). Besides, the number of recovered people without any system is solved numerically using dde23. 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, 7.3.1. Between-Hosts Vaccination. As shown in Figures 16(a) more people become infected, and as a result, more people and 16(b), the number of susceptible individuals in the will be recovered. By injecting the vaccine into the suscepti- absence of a controller (vaccination has not been applied) ble person’s body, during the conversion of healthy cells to 24 Computational and Mathematical Methods in Medicine 14 12 10 Number of molecules 8 6 4 2 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Time (days) IFN without controller IFN with controller (a) ×105 8.2 8.18 8.16 8.14 Number of cell 8.12 8.1 8.08 8.06 8.04 8.02 8 0 10 20 30 40 50 60 70 80 90 100 Time (days) K without controller K with controller (b) Figure 15: Evolution of IFN-I molecules and natural killer cells with and without the controller. resistant cells, the susceptible person also vaccinates (direct 7.3.2. Between-Hosts Therapy. When the infection spreads transfer of individuals from the susceptible group to the in the community without control over its prevalence, recovered group). As a result, the number of susceptible the number of infected people is increasing every day, individuals decreases and moves to the recovered individual as shown in Figure 17(a). It rose by over 3,000 people. group at the same time. This means the susceptible people Finally, their number decreases until the 25th day and vaccinate against the disease in almost 25 days ðt 3 Þ (the reaches a little under 2,000 and remains unchanged, healthy cells of the susceptible individual convert to indicating the community’s disease persistence. Also, the resistant-to-infection cells in eight days as the delay time exposed people number reaches a high number of 4,000 with vaccination), while in [13], the final time of the vaccina- in only twelve days and ultimately converges to a stable tion campaign of nearly 40 days is suggested without consid- population and remains at almost 1,500 (see ering the delay in the vaccination of susceptible people. Figure 17(b)). Computational and Mathematical Methods in Medicine 25 15000 10000 9000 8000 Number of humans 7000 10000 Number of humans 6000 5000 4000 5000 3000 2000 1000 0 0 0 5 10 15 20 25 0 5 10 15 20 25 Time (days) Time (days) Susceptible without controller Recovered without controller Susceptible withcontroller Recovered with controller (a) (b) Figure 16: The population of susceptible and recovered individuals with and without the controller. 3500 4000 3000 3500 2500 3000 Number of humans Number of humans 2500 2000 2000 1500 1500 1000 1000 500 500 0 0 0 5 10 15 20 25 0 5 10 15 20 25 Time (days) Time (days) Infected without controller Exposed without controller Infected with controller Exposed with controller (a) (b) Figure 17: The population of infected and exposed individuals with and without the controller. In contrast, the number of infected people exceeds 8. Conclusion almost 500 with antiviral therapy, which is used to treat infected people (transfer them to a recovered group). During In this paper, the optimal control theory was applied to two the computed optimal time (the same delay as six days), the types of dynamic epidemic models: the innate immune healthy cells convert to resistant-to-infection cells; at the response dynamic and the SEIR epidemic model. A nonlin- same time, virus particles are eradicated by antiviral treat- ear transmission rate and time delay were considered in the ment in the infected person’s body. Therefore, the infected SEIR epidemic model. The aim of this paper was twofold: people move to the recovered people group in almost 25 firstly, preventing the virus from spreading within the body days. According to Remark 8, the number of exposed people of a host and secondly, curb the epidemic in the society also converges to zero as a result of decrement of infected and between hosts. To this end, two control strategies were people in 25 days. introduced to satisfy the first goal in the susceptible and 26 Computational and Mathematical Methods in Medicine infected person’s bodies. In this regard, vaccination was used [9] W. O. Kermack and A. G. McKendrick, “Contributions to the to convert healthy cells to resistant-to-infection cells inside mathematical theory of epidemics. II. The problem of ende- the susceptible individual’s body that recovered them against micity,” Proceedings of the Royal Society of London. Series A, the disease in society. In the following, antiviral treatment Containing Papers of a Mathematical and Physical Character, was used to reduce the concentration of viruses and convert vol. 138, no. 834, pp. 55–83, 1932. healthy cells to resistant-to-infection cells at the same time [10] W. O. Kermack and A. G. McKendrick, “Contributions to the inside the infected individual’s body to recover them. More- mathematical theory of epidemics. III. Further studies of the over, to control the epidemic in society, the optimal control problem of endemicity,” Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and was used to increase the number of recovered people by Physical Character, vol. 141, no. 843, pp. 94–122, 1933. reducing the number of susceptible and infected people in [11] L. Zhu, G. Guan, and Y. Li, “Nonlinear dynamical analysis and the fixed 25 days. Transferring occurred with a delay com- control strategies of a network-based SIS epidemic model with puted as the optimal time of disease eradication inside the time delay,” Applied Mathematical Modelling, vol. 70, pp. 512– infected individual’s body and the recovery of the susceptible 531, 2019. individual’s body. Finally, numerical simulation was used to [12] A. E. S. Almocera, V. K. Nguyen, and E. A. Hernandez-Vargas, illustrate the usefulness of the obtained results. “Multiscale model within-host and between-host for viral infectious diseases,” Journal of Mathematical Biology, vol. 77, Data Availability pp. 1035–1057, 2018. [13] A. H. Amiri Mehra, I. Zamani, Z. Abbasi, and A. Ibeas, The data used are included within the paper and cited “Observer-based adaptive PI sliding mode control of devel- accordingly. oped uncertain SEIAR influenza epidemic model considering dynamic population,” Journal of Theoretical Biology, vol. 482, Conflicts of Interest article 109984, 2019. [14] M. Lhous, M. Rachik, and A. Larrache, “Free optimal time The authors declare that they have no conflicts of interest. control problem for a SEIR-epidemic model with immigration of infective,” International Journal of Computer Applications, Acknowledgments vol. 159, no. 3, pp. 1–5, 2017. [15] A. B. Kurzhanski, Control Systems, Robotics and Automation– This work was supported in part by the Basque Government Volume VIII, Advanced Control Systems-II, 2009. through project IT1207-19. [16] R. F. Stengel, R. Ghigliazza, N. Kulkarni, and O. Laplace, “Optimal control of innate immune response,” Optimal Con- References trol Applications and Methods, vol. 23, no. 2, pp. 91–104, 2002. [1] A. Boianelli, V. Nguyen, T. Ebensen et al., “Modeling influenza [17] G. Zaman, Y. H. Kang, and I. H. Jung, “Optimal treatment of virus infection: a roadmap for influenza research,” MDPI, an SIR epidemic model with time delay,” Biosystems, vol. 98, vol. 7, no. 10, pp. 5274–5304, 2015. no. 1, pp. 43–50, 2009. [2] J. Tan and X. Zou, “Optimal control strategy for abnormal [18] E. Grigorieva, E. Khailov, and A. Korobeinikov, “Optimal con- innate immune response,” Computational and Mathematical trol for a SIR epidemic model with nonlinear incidence rate,” Methods in Medicine, vol. 2015, Article ID 386235, 16 pages, Mathematical Modelling of Natural Phenomena, vol. 11, 2015. no. 4, pp. 89–104, 2016. [3] C. Beauchemin, J. Samuel, and J. Tuszynski, “A simple cellular [19] E. Grigorieva, E. Khailov, and A. Korobeinikov, “Optimal con- automaton model for influenza A viral infections,” Journal of trol for an SEIR epidemic model with nonlinear incidence Theoretical Biology, vol. 232, no. 2, pp. 223–234, 2005. rate,” Studies in Applied Mathematics, vol. 141, no. 3, [4] E. W. Larson, J. W. Dominik, A. H. Rowberg, and G. A. Hig- pp. 353–398, 2018. bee, “Influenza virus population dynamics in the respiratory [20] M. Chen, M. Li, Y. Hao, Z. Liu, L. Hu, and L. Wang, “The tract of experimentally infected mice,” Infection and Immu- introduction of population migration to SEIAR for COVID- nity, vol. 13, no. 2, pp. 438–447, 1976. 19 epidemic modeling with an efficient intervention strategy,” [5] G. Bocharov and A. Romanyukha, “Mathematical model of Information Fusion, vol. 64, pp. 252–258, 2020. antiviral immune response III. Influenza A virus infection,” [21] A. H. Amiri Mehra, M. Shafieirad, Z. Abbasi, and I. Zamani, Journal of Theoretical Biology, vol. 167, no. 4, pp. 323–360, “Parameter estimation and prediction of COVID-19 epi- 1994. demic turning point and ending time of a case study on [6] P. Baccam, C. Beauchemin, C. A. Macken, F. G. Hayden, and SIR/SQAIR epidemic models,” Computational and Mathe- A. Perelson, “Kinetics of influenza A virus infection in humans,” matical Methods in Medicine, vol. 2020, Article ID Journal of Virology, vol. 80, no. 15, pp. 7590–7599, 2006. 1465923, 13 pages, 2020. [7] A. E. Hernandez-Vargas and M. Meyer-Hermann, “Innate [22] N. Gozalpour, E. Badfar, and A. Nikoofard, “Transmission immune system dynamics to influenza virus,” IFAC Proceed- dynamics of novel coronavirus SARS-Cov-2 among healthcare ings Volumes, vol. 45, no. 18, pp. 260–265, 2012. workers, a case study in Iran,” Nonlinear Dynamics, vol. 105, [8] W. O. Kermack and A. G. McKendrick, “A contribution to the no. 4, pp. 3749–3761, 2021. mathematical theory of epidemics,” Proceedings of the Royal [23] G. Giordano, F. Blanchini, R. Bruno et al., “Modelling the Society of London. Series A, Containing Papers of a Mathemat- COVID-19 epidemic and implementation of population- ical and Physical Character, vol. 115, no. 772, pp. 700–721, wide interventions in Italy,” Nature Medicine, vol. 26, no. 6, 1927. pp. 855–860, 2020. Computational and Mathematical Methods in Medicine 27 [24] J. L. Gevertz, J. M. Greene, C. H. Sanchez-Tapia, and E. D. Son- tag, “A novel COVID-19 epidemiological model with explicit susceptible and asymptomatic isolation compartments reveals unexpected consequences of timing social distancing,” Journal of Theoretical Biology, vol. 510, article 110539, 2021. [25] A. Ibeas, M. de la Sen, S. Alonso-Quesada, and I. Zamani, “Sta- bility analysis and observer design for discrete-time SEIR epi- demic models,” Advances in Difference Equations, vol. 2015, no. 1, 2015. [26] H. Boutayeb, S. Bidah, O. Zakary, and M. Rachik, “A new sim- ple epidemic discrete-time model describing the dissemination of information with optimal control strategy,” Discrete Dynamics in Nature and Society, vol. 2020, Article ID 7465761, 11 pages, 2020. [27] L. Göllmann and H. Maurer, “Optimal control problems with time delays: two case studies in biomedicine,” Mathematical Biosciences & Engineering, vol. 15, no. 5, pp. 1137–1154, 2018. [28] M. Elhia, M. Rachik, and E. Benlahmar, “Optimal control of an SIR model with delay in state and control variables,” Interna- tional Scholarly Research Notices, vol. 2013, Article ID 403549, 7 pages, 2013. [29] S. Sharma, A. Mondal, A. K. Pal, and G. P. Samanta, “Stability analysis and optimal control of avian influenza virus a with time delays,” International Journal of Dynamics and Control, vol. 6, no. 3, pp. 1351–1366, 2018. [30] L. Göllmann and H. Maurer, “Theory and applications of opti- mal control problems with multiple time-delays,” Journal of Industrial & Management Optimization, vol. 10, no. 2, pp. 413–441, 2014. [31] J. E. Forde, S. M. Ciupe, A. Cintron-Arias, and S. Lenhart, “Optimal control of drug therapy in a hepatitis B model,” Applied Sciences, vol. 6, no. 8, p. 219, 2016. [32] P. A. Naik, J. Zu, and K. M. Owolabi, “Global dynamics of a fractional order model for the transmission of HIV epidemic with optimal control,” Solitons & Fractals, vol. 138, article 109826, 2020. [33] Z. Abbasi, I. Zamani, A. H. Amiri Mehra, M. Shafieirad, and A. Ibeas, “Optimal control design of impulsive SQEIAR epi- demic models with application to COVID-19,” Chaos, Solitons & Fractals, vol. 139, article 110054, 2020. [34] J. Li and Y. Yang, “SIR-SVS epidemic models with continuous and impulsive vaccination strategies,” Journal of Theoretical Biology, vol. 280, no. 1, pp. 108–116, 2011. [35] J. A. Gondim and L. Machado, “Optimal quarantine strategies for the COVID-19 pandemic in a population with a discrete age structure,” 2020, https://arxiv.org/abs/2005.09786. [36] S. Ullah and M. A. Khan, “Modeling the impact of non- pharmaceutical interventions on the dynamics of novel coro- navirus with optimal control analysis with a case study,” Chaos, Solitons & Fractals, vol. 139, article 110075, 2020. [37] A. T. Azar and A. E. Hassanien, Modeling, control and drug development for COVID-19 outbreak prevention, Springer, 2021. [38] R. F. Stengel, Optimal Control and Estimation, Courier Corpo- ration, 1994. [39] M. I. Kamien and N. L. Schwartz, Dynamic Optimization: the Calculus of Variations and Optimal Control in Economics and Management, Courier Corporation, 2012. [40] X. Liu and P. Stechlinski, Infectious Disease Modeling: A Hybrid System Approach, Springer, 2017.