Eﬀect of an antiviral drug control in host COVID-19 kinetics

The article proposes and analyzes a system of diﬀerential equations modeling the interaction of the SARS-CoV-2 virus and the epithelial cells of the human lungs. Optimal control strategies representing antiviral drug treatment eﬀects of this model are explored here. The Pontryagin’s maximum principle is used to clarify the optimal control strategies. The existence of optimal control is proved and eﬀective strategies are illustrated. Numerical simulations, eﬃciency analysis, and cost-eﬀectiveness analysis reveal that time-dependent antiviral drug with other control mechanisms, would reduce the viral load and control the infection process at low cost.


Introduction
COVID-19 is the disease that affects the whole world and this disease is caused by Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) [1]. Coronavirus can be spread and infects humans through droplets via coughing and sneezing from infected humans. A total of 13 million people worldwide living with COVID-19 today and more than half a million people died to date. Scientists and researchers of the whole world are searching for the anti-COVID-19 vaccine. But to date, no such fruitful results have shown yet. In view of that, the use of available antiviral drugs and immunosuppressive drugs combination is the only way to fight against COVID-19 infection [2]. Therefore non-pharmaceutical interventions like self-isolation, using masks, hand washing, sanitation are suggested by the WHO (World Health Organisation) and ICMR (Indian Council for Medical Research). Also, the Government of states has implemented lockdown, travel restrictions, quarantine measures and testing to control the disease. Several epidemiological mathematical models [3][4][5][6][7][8][9] have been developed to make the right decisions in these measures. These have highlighted that social distancing intervention to mitigate the epidemic is a key aspect.
However, a few mathematical models have been developed to study the COVID-19 dynamics within host levels based [2,[10][11][12]. The epithelial cells in the upper respiratory tract and upper divisions of bronchi are the major target area of the SARS-CoV-2 virus [2,13]. If the infections spread to the lower lungs, the patients' conditions become severe [2,14,15]. The cell infection process by the SARS CoV-2 virus has mainly two stages. At first, the SARS-CoV-2 virus binds to the angiotensin-converting enzyme 2 (ACE2) receptors on the surface epithelial cells [16] and viral genome penetrate the cell and cells become infected. The infected cells turned into virus-producing cells and the new virions reproduced through the eclipse and burst phase [2].
In this research article, we proposed mathematical modeling to investigate the SARS-CoV-2 replication within host, specifically the interaction between susceptible epithelial and SARS-CoV-2 virus. We also study the current pharmaceutical intervention in our model. Here we study the effect of commonly used antiviral and immune-suppressing drugs. It has been observed that antiviral drugs and immunosuppressive drugs play a pivotal role in the treatment of high -risk patients [2]. Immunosuppressive drugs reduce the adaptive immune responses (AIR) to a level low enough not to interfere with the innate immune response during the early phase of infection [2]. Due to the effect of the immunosuppressive drugs, the infected cells become inadequate to produce new virions, the drug should be removed so that the AIRs can clear the remaining virus. This new insight looks at the drug regimen which can lead to fast and complete recovery. In this perspective, we have formulated the basic mathematical model of COVID-19 based on drug effectiveness.
The article is designed as follows: first we provide the formulation of our mathematical model centered with drug effectiveness in the Section 2 and studied the model analytically. In Section 3, We have framed the optimal control problem in the light of cost function for optimal drug dosage. In Section 4, the numerical studies of both the models (time independent and time dependent control) are depicted. Also the drug efficiency along with cost effectiveness have been verified. In Section 5, we have drawn discussions regarding our analytical and numerical findings.

Analysis of basic model with constant control
In terms of biological view we have considered an ordinary differential equation model (taken from [2]) which describes the interaction between epithelial cells and SARS-CoV-2. We have explored the optimal control strategy of this ODE model. Here we consider E S , E I and V representing the concentration of the susceptible epithelial cells, infected epithelial cells and free virus particles respectively.
with this non negative initial conditions: The model considers that t 0 ∈ [0, ∞) is the starting day of COVID-19 infection.

Positivity and boundedness
Theorem 1. All the solutions of the system (1) subject to initial condition (2) are positive and bounded in R 3 + , for all t > 0.
Proof. The first equation of the system (1) can be written as; where a = [r − (1 − 1 )βV ] and b = r/K. Then on integration we get, In the similar way, it can be observed that Consequently, all the solution trajectories (E S , E I , V ) of system (1) with initial conditions (2) are positive. Next we show the boundedness of the solutions. According to [17], we can construct from the first and second equations of system (1) lim sup Accordingly, E S ≤ µ and E I ≤ µ. Similarly, from the third equation of the system (1) we can write Therefore, we get the positively invariant bounded region The region Φ is well posed and attracting since all the solutions of the system (1) with initial conditions (2) will enter this region and will never leave it.

Equilibria
In order to study, we determine the biologically feasible equilibrium points. The system (1) has the following equilibria-1. The model (1) executes two disease-free equilibrium points, one is P 0 (0, 0, 0) and second one is P 0 (K, 0, 0). Here we consider only second one disease-free equilibrium P 0 (K, 0, 0) because this is very relevant for our system (1) with initial condition (2) to be biologically feasible. Disease-free equilibrium P 0 (K, 0, 0) always exists.
2. P * (E * S , E * I , V ) is the endemic equilibrium where COVID-19 infection persists in the system.

Basic reproduction number
In this subsection we now determine the basic reproduction number (R 0 ) of the system (1) using the next generation matrix method anticipated by Diekmann et. al. [18] and van den Driessche and Watmough [19]. It is notable that the infected compartments of the system (1) are E I and V . Let M represents the appearance of new infections while N stands for transfer of the infections between the compartments and x = [E I , V ] T . Then at the disease free steady state P 0 (K, 0, 0) the linearized infectious subsystem of the system (1) is given Accordingly we can calculate that Therefore next generation matrix is M N −1 and The basic reproduction number (R 0 ) is the spectral radius of the next generation matrix that is

Existence of the endemic steady state
To acquire feasible solutions of the system (1), the equilibria of the system must be positive. The positive components of the unique endemic equilibrium Therefore it is observed that the unique endemic equilibrium P * (E * S , E * I , V * ) would exist only if R 0 > 1.

Local stability analysis
Theorem 2. The disease free equilibrium point P 0 (K, 0, 0) of the system (1) is locally asymptomatic steady state point whenever R 0 < 1 and would be unstable Proof. The Jacobian matrix calculated at the disease free equilibrium point P 0 (K, 0, 0) is The characteristic equation corresponding of the Jacobian matrix J P0 corresponding to the eigenvalue ρ is given by Therefore we get ρ = −r and the real parts of the other two eigenvalues must be negative for the locally asymptomatic stability of the equilibrium P 0 and this will be possible only when If one of the eigenvalue or its real part would be positive, then R 0 > 1 and consequently P 0 would be unstable.
Proof. The Jacobian matrix around the equilibrium point P * (E * S , E * I , V * ) is The characteristic equation of the Jacobian matrix J P * corresponding to the eigen value γ is given by where, It is notable that 3 > 0, 2 > 0, 1 > 0 and 1 2 > 3 only for R 0 > 1 that is the Routh -Hurwitz criterion is satisfied for the cubic equation (6) and this ensures that all the three roots of the characteristic equation (6) have negative real parts when R 0 > 1. Therefore the endemic equilibrium point P * (E * S , E * I , V * ) is locally asymptotic stable provided R 0 > 1.

Persistence of the system
Proof. Previously we already proved that, lim sup t→∞ E S ≤ rK Let us consider that U = max(E S , E I , V ). Thus there exists a U > 0 such that The third equation of the system (1) indicates, Since is arbitrary small, From the first equation of the system (1), we can obtain Thus we get, Again from the second equation of the system (1), we get Therefore, Accordingly, E I > 0 if the condition (7) holds. Let us consider χ = min(E S , E I , V ). Therefore there exists a χ > 0 such that Therefore the system (1) is uniformly persistent.

Global stability analysis of disease free equilibrium
Theorem 5. The disease free equilibrium point P 0 (K, 0, 0) is a globally asymptomatic steady state point whenever R 0 < 1 and would be unstable if R 0 > 1.
Proof. Let us define a Lyapunov function L : Φ → R for the system (1) as The Lyapunov function L is positive in the non-negative cone Φ and attains zero at P 0 (K, 0, 0). Taking orbital derivative of L along positive solutions of the system (1), we get It is notable that in Φ, E S ≤ K and consequently we get, The expression (8) reveals that the dL/dt would be zero or negative only if R 0 ≤ 1. Let us define the subset S of Φ where dL/dt = 0 as t → ∞ and is given by the equations From the biological point of view, it is observed that the maximum invariant set Q contained in the set S is the plane E I = 0 and V = 0. Therefore S = P 0 (K, 0, 0) and solutions started in Q approaches to P 0 (K, 0, 0) as t → ∞. Using Lyapunov-Lasalle theorem and LaSalle's Invariant Principle [20], [21], we can conclude that the disease-free steady state is globally asymptomatic stable and all the solution trajectories started in Φ approaches to P 0 (K, 0, 0) as t → ∞. Hence the proof is completed.

Sensitivity
In this subsection, sensitivity analysis has been executed to elucidate the relative weightage and the robustness of predictions to the baseline parameters associated to the basic reproduction number (R 0 ) [22]. Sensitivity analysis advantages to emphasize those parameters which are highly sensitive comparative to others in COVID-19 infection and we should pay attention to them to yield proper intervention schemes. The normalized forward sensitivity indices of R 0 with respect to its all the six associated parameters are given by The normalized sensitivity indices can be independent or dependent on the baseline parameters. In the sensitivity analysis of R 0 , the most sensitive parameters are carrying capacity of epithelial cell (K), rate of infection (β) and number of free virus produced (η). Since sensitivity index of these parameters are all equal to 1, hence we can conclude that 10% increases (or decreases) of if the value of carrying capacity of epithelial cell (K), infection rate or number of free virus increases (or decreases) the basic reproduction number by 10%. Other important parameter is virus removal rate (δ V ). The decreases (or increases) of virus removal rate by 10% increases (or decreases) the basic reproduction number by 6.67%. Increases (or decreases) of the effectiveness of immunosupressive drug by 10% decreases (or increases) the basic reproduction number by 9.23%. Similarly increasing (or decreasing) of antiviral drug effectiveness by 10% decreases (or increases) the basic reproduction number by 7.24%. Hence immunosuppressive drug plays a pivotal role to control the disease progression in comparison with antiviral drug.
Next, we consider (time dependent) optimal control strategies associated with case holding and case finding based on host SARS-CoV-2 / COVID-19 model.  Table 1. Table 2 shows the Partial Rank Correlation Coefficients between R 0 and each parameters of model

Model with time dependent control
Optimal control technique; this is one of the most powerful mathematical tools which make the result involving infectious disease dynamical systems [23], [24], [25]. By using this mathematical tool we are trying to reduce the spread of corona virus infection in our body. Our main objective is to minimize the number infected cell and to maximize the number of uninfected cells [26]. In this optimal control problem, we introduce two control variables.
(i). The first control u 1 (t) represents the antiviral drug strategy (Strategy I) which block infection.
(ii). The second control u 2 (t) represents the drug strategies (Strategy II) which block the production of virus particles.
Our main aim is to maximize the level of susceptible epithelial cells. Also we want to keep cost -as measured in terms of antiviral drug therapy as low as possible. Our control is a function of u 1 (t) and u 2 (t) with values normalized between 0 and 1, where u i (t) = 1, i = 1, 2 represents totally effective chemother-apy and u i (t) = 0, i = 1, 2 represents no treatment. We choose our control class as follows: Our goal is to minimize the objective functional: subject to the state system with this positive initial conditions: All the coefficients A, B 1 , and B 2 are non-negative and they represent weights on the different terms of objective functional.

Existence of an optimal control
In this subsection, the existence of optimal control for the system (11) is analyzed. In order to use Pontryagin Maximum Principle [26] to characterize an optimal control, first we have to prove the existence of an optimal control [23,27]. Theorem 6. There exists an optimal control vector (u * 1 (t), u * 2 (t)) ∈ M with corresponding state solutions (E * S , E * I , V * ) that minimizes the objective functional J (u 1 (t), u 2 (t)) defined by (9) [23].
Proof. All the solutions of the system (11) are non-negative and uniformly bounded for given initial conditions (12). Non-negativity of objective functional values is obvious and non-negativity implies the boundedness of the objective functional. Hence a minimizing sequence of controls (u v 1 (t), u v 2 (t)) ∈ M exists such that; Therefore the controls belonging to M are uniformly bounded in L ∞ and consequently they are uniform bounded in L 2 ([0, T ]). Since the space L 2 ([0, T ]) is reflexive [28], there exists (u * 1 (t), u * 2 (t)) ∈ M such that on a sub sequence, Thus obviously the state sequence (E v S , E v I , V v ) is uniformly bounded corresponding to (u v 1 (t), u v 2 (t)). It can be seen that the uniform boundedness of the right-hand sides of the system (11) implies the uniform boundedness of the derivatives for (E v S , E v I , V v ) and equicontinuity of the corresponding state sequence (E v S , E v I , V v ). According to the Arzelï-Ascoli Theorem, there exists (E * S , E * I , V * ) such that on a subsequence, Next, we choose the proper subsequence passing the limit to system (11) and corresponding to the controls u * 1 (t), u * 2 (t), we can achieve the state solution (E * S , E * I , V * ). Then the lower semi-continuity of the L 2 -norm with respect to L 2 weak convergence implies that, Therefore, (E * S , E * I , V * ) is an optimal control. Theorem 7. Define an optimal control vector (u * 1 (t), u * 2 (t)) ∈ M and the corresponding state solutions (E * S , E * I , V * ) in the model (11), there exist adjoint variables ϑ 1 (t), ϑ 2 (t), and ϑ 3 (t), satisfying, with the transversality conditions: ϑ 1 (T ) = 0, ϑ 2 (T ) = 0, ϑ 3 (T ) = 0. Furthermore, the optimal control vector is given by (u * 1 (t), u * 2 (t)), where u * 1 (t) = min max 0, Proof. Following the Pontryagin's Minimum Principle, we can obtain the Hamiltonian as: Now, adjoint variables ϑ 1 (t) = 0, ϑ 2 (t) = 0, ϑ 3 (t) = 0 by: with the transversality conditions ϑ 1 (T ) = 0, ϑ 2 (T ) = 0, ϑ 3 (T ) = 0. We get the characterization of optimal controls by saying ∂H ∂u 1 (t) = 0 and ∂H ∂u 2 (t) = 0.
From ∂H/∂u 1 (t) = 0 and ∂H/∂u 2 (t) = 0, we get By holding the upper and lower bounds for u 1 (t) and u 2 (t) into account, the characterization of optimal controls are as follows: Thus the proof is completed.

Numerical simulation
In this present section, we develop graphical displays of the fixed control system (1) and optimal control system (11), (12) and present the results. For fixed control problem the drug efficacy 1 and 2 are updated. Also for time dependent control, the controls u 1 (t) and u 2 (t) are updated and used to sole the state system and adjoint system together. Optimal control is a two-points boundary value problem with separate boundary conditions at times t = 0 and t = T . Here our aim is to solve this problem for the value of T = 100. This is the value (in days) at which treatment is stopped. We use a finite difference approach to solve the optimality system (11) -(12) by applying the package BVP4C of MATLAB. To solve the problem easily we use T = 1 and t = 0. The solution of the problem for T = 1 can be used as an initial guess to the solution nearby such that T = 1 + ∆T for ∆T sufficiently small. This process can be continued until the required problem is solved and this referred to as a homotopy path. By selecting different weight factors we can generate several treatment schedules for various time periods. Here we illustrate the 100 -days treatment schedule. The optimal solution trajectories are displayed below in Figures (6, 8, 9) with different weight constants keeping all other parameters constant.
Here we study the following control strategies To study the efficiency of the drug therapy, we consider three strategies as follows:   Strategy I: For this strategy, we observe that for fixed control (in figure (5)) the total number of uninfected cells increases and infected cell number decreases with increasing value of the antiviral drug efficacy ( 1 ) . Figure (6) represents the effect of time dependent antiviral drug which block infection. In this figure, we observe that the number of infected cell along virus decreases as we keep the drug level at its maximum value throughout the treatment period.
Strategy II: For this strategy, we observe that for fixed control (in figure (7)) the total number of uninfected cells increases and infected cell number decreases with increasing value of the antiviral drug efficacy ( 2 ). Figure (8) represents the effect of time dependent antiviral drug which block new virus production. In this figure, we observe that the number of virus decreases as we keep the drug level at its maximum value throughout the treatment period. But in this strategy no such significant effect observed in case of infected cells.
Strategy III: Figure (9) represents the effect of time dependent combination of antiviral drug effect. In this figure, we observe that the number of infected cell along virus decreases as we keep the drug level at its maximum value throughout the treatment period. The outcome of this strategy is almost similar to the Strategy I.  (11) showing the effect of the optimal strategies (u 1 = 0, u 2 = 0) for A = 100, B 1 = 300.

Efficiency analysis
For Strategy I we consider u 1 (t) = 0, u 2 (t) = 0, for Strategy II, u 1 (t) = 0, u 2 (t) = 0 and for Strategy III we consider u 1 (t) = 0, u 2 (t) = 0. Here we shall calculate the efficiency index E which is defined as where A C represents the area under the infected epithelial cell concentration as a function of time when the control is used and A S is the area under the infected total population curve in absence of control input. The cumulative number of infected epithelial cells during the time interval [0, 1] is defined by By calculating the efficacy index we can adopt the best strategy whose efficiency index will be the biggest [29,30]. The values of A C and the efficacy index for two strategies are given in Table 4.2.
From the Table 4.2, we can conclude that Strategy I is more effective than Strategy II. But Strategy III (which is the combination of two drug therapy) is the best strategy.

Cost effectiveness analysis
To compare the differences between the costs and health outcomes of two alternative strategies, Incremental Cost-Effectiveness Ratio (ICER) plays a crucial role to select the best strategies. In the case of ICER, when we compare two competing strategies, the best strategy will be selected from less effective ICER values. In order to quantify the cost effectiveness the Incremental Cost-Effectiveness Ratio (ICER) measures is defined as ICER = The differences in intervention costs The differences in averted infection .
Based on the control model simulation, we calculate the ICER for the strategies given below The ICER is calculated as follows:    (11) showing the effect of the optimal strategies (u 1 = 0, u 2 = 0) for A = 100, B 1 = 300, B 2 = 70000 .

Discussion and conclusion
We proposed the time-independent and time-dependent control problem to keep the susceptible epithelial cells level high and infected epithelial cells and virus level low. Also, our aim is to minimize drug costs. Here we use two control strategies. The first control input is described in terms of both susceptible epithelial cells and virus and their corresponding adjoint variables. It mainly blocks the infection or disease transmission. The second control input is described in terms of free virus generation and it acts to block new virus production.
We have derived the expression of the basic reproduction number by using the next-generation approach. By calculating the sensitivity index we have obtained that the most sensitive parameters are the carrying capacity of epithelial cell (K), rate of infection (β), and the number of free viruses produced (η). PRCC results showed the negative correlation with the parameters δ V , 1 , and strategies. Thus WHO suggests usage of the commonly used antiviral drugs to fight against the COVID-19 infection. As a result, the death rate would be restricted below 5% worldwide. Also, we have investigated the cost-effectiveness of the controls to determine the best strategies with minimum costs. By calculating ICER, we found that Strategy III i.e. combination of drug therapy is the most effective therapy. Also, it has been observed that the Strategy I is effective in eliminating the disease. This strategy has the same effect as a combination of drug therapy. But this single drug strategy is not cost-effective as like as combination therapy. In a nutshell, we can conclude that: (i) The model has a disease free equilibrium that is locally asymptotically stable if R 0 < 1.
(ii) When R 0 > 1 the system attain its endemic equilibrium.