SARS-CoV-2 infection with lytic and non-lytic immune responses: A fractional order optimal control theoretical study.

In this research article, we establish a fractional-order mathematical model to explore the infections of the coronavirus disease (COVID-19) caused by the novel SARS-CoV-2 virus. We introduce a set of fractional differential equations taking uninfected epithelial cells, infected epithelial cells, SARS-CoV-2 virus, and CTL response cell accounting for the lytic and non-lytic effects of immune responses. We also include the effect of a commonly used antiviral drug in COVID-19 treatment in an optimal control-theoretic approach. The stability of the equilibria of the fractional ordered system using qualitative theory. Numerical simulations are presented using an iterative scheme in Matlab in support of the analytical results.

Introduction

The novel coronavirus disease (COVID-19) caused by severe acute respiratory syndrome coronavirus (SARS-CoV-2) has rapidly spread across the globe. This disease causes more than 9.9 lakhs cases and lakhs death as of 27th June 2020. The death rate of COVID-19 is approximately 5%. Most of the cases of COVID-19 are asymptomatic and can be recovered before testing. 20% of the patients [1] are found to be symptomatic. Common symptoms for COVID-19 include fever, cough, shortness of breath, sore throat, chills, fatigue, etc. But in the case of the patient suffering from diabetes, high blood pressure, chronic pneumonia, COVID-19 causes various organ collapses and also mortality.

The angiotensin-converting enzyme (ACE2) receptors on the surface of the epithelial cells [1] are the major target area of the SARS-CoV-2 virus. The density of epithelial cells is the highest in the lungs, followed by the nose, and lastly, the trachea/bronchi tissues. Thus pneumonia is quite common for COVID-19 affected patients.

The immune system acts against the virus during viral infections in presence of inherent and antigen-specific immune responses [2]. These forms of immune responses classified into two categories (lytic and nonlytic). The lytic mechanism represents the killing process of infected cells and the non-lytic effector mechanism inhibit the viral replication by soluble mediators. Cytotoxic T Lymphocyte cells (CTLs) play a pivotal role to kill infected cells, though antibodies deactivate free virus which results from the inhibition of infected cells. In this article, a COVID-19 model is constructed to study the cell-to-cell transmission of the virus. We have considered uninfected, infected epithelial cells, SARS-CoV-2 population, and CTL population. We have also studied the lytic and non-lytic immune responses effect on the model dynamics.

It has been observed that for COVID-19 infection, the average maturation period is days to days and the median time to discharge from the hospital of symptomatic patients is days [1], however the average time to death is days [3]. From clinical findings, it is observed among patients who were undergone through antiviral drug lopinavir/ritonavir and were released. The average time from infection start to commencement of antiviral treatment was days [4]. Also, the intermediate period of viral coming off was days. Considering the effect of commonly used drugs, we have studied the effect of commonly used drugs to explore the best possible treatment regimen.

Fractional order differential equation is used as an important tool in nonlinear mathematical modeling [5], [6], [7], [8], [9], [10], [11], [12].

Fractional differential equation (FDE) have recently used to model real-life phenomena and it stands as a strong tool to study the dynamics of the system. In recent studies, FDEs are being used as new and effective tools to model real-life phenomena instead of Integer order differential equations.

The fractional calculus is initially proposed by Leibniz [13]. The most frequently used definitions for fractional derivatives are the Riemann–Liouville and Caputo fractional derivatives. There are some other definitions of fractional derivatives [14] but they are hard to implement in studying real-life problems. In order to overcome the difficulties, the Caputo type derivative is introduced which is closely related to the Riemann–Liouville derivative.

FDE provides a good tool for describing the memory feature for the biological system in mathematical epidemiology. Actually, a differential operator with an integer-order is a local operator, on the other hand, a differential operator with a fractional-order is nonlocal. It takes into account all the history of its previous states along with the present state. A large number of advanced techniques have been constructed for solving numerically the integer-order differential equations numerically. Whereas, for FDEs, the numerical techniques are not well developed. One of the effective numerical methods, so far, to solve FDEs, is a generalized Adams–Bashforth–Moulton algorithm [?].

In the present study, we consider Caputo type fractional order differential equations (FODE) of order . In a mathematical model, the idea of FODE is used in order to explore the concept of memory, one of the genetic behavior in a host-cell system [6]. The impact of memory between an integer-order derivative and a power of time was first explored by Caputo and Riemann through Liouville derivatives [15]. Cell biological fractals structure is closely related to the fractional-order differential equations [16].

Recently many mathematicians and biologists formulate the mathematical model and study the COVID-19 dynamics by using FDE, for example, see [17], [18], [19] and the references therein. Atangana [17] used Italy’s COVID-19 data to formulate the mathematical model and they have considered FDE to include the non-locality into the modeling. Tuan et al. [19] investigated the transmission of COVID-19 with the help of fractional–order Caputo derivative. With the help of fixed point theory, the prove the existence and uniqueness of the solution. Baleanu et al. [20] investigate a model of the COVID-19 transmission in different groups of people using the Caputo–Fabrizio fractional derivative. They have solved the system by using the homotopy analysis transform method (HATM) and obtain approximate solutions in convergent series. Khan and Atangana [18] formulate the fractional-order mathematical model on considering cases for January 21, 2020, to January 28, 2020, and estimated the parameters of the model. They compute the basic reproduction number for the data is and numerically illustrate their findings.

The article is organized in the following order. We formulate a fractional-order differential equation system in “Formulation of mathematical model’. In ‘Mathematical analysis’ we study the non-negativity and local and global existence of the model. The stability of the system is studied in ‘Equilibrium and stability analysis’. In ‘The fractional optimal control problem’ the fractional order optimal control therapy problem is analyzed. In ‘Numerical findings’, the numerical simulations are shown and the discussions are given in ‘Discussion’.

Formulation of mathematical model

We first present the definitions of fractional-order derivatives and fractional-order integration [21]. According to the concept of fractional order derivative, we have to adopt Caputo’s definition which is the modification of the Riemann–Liouville definition [22]. This definition has the more advantage of dealing properly with the initial value problems.

The left-sided Caputo fractional order derivative is defined as:

The right-sided Caputo fractional order derivative is defined as: where is the order of the derivative and . Also is the gamma function, whereas is an integer.

The left-sided Riemann–Liouville fractional order derivative is defined as: Also, the right-sided Riemann–Liouville fractional derivative is defined as: where and are positive constants.

In general, most of the dynamic systems have been analyzed by using integer-order differential equations. Recently fractional order equations have been used widely. The dynamical model depends on the chronological situation along with the present state [23]. Thus we study the system by using fractional order differential equation. Recently, Chatterjee and Basir [24] has introduced a basic SARS-CoV-2 model with lytic immune response and its antiviral treatment. Here we have studied both the non-lytic and lytic effect of immune responses in the presences of antiviral drug therapy.

We begin with the basic dynamics between uninfected epithelial cells , the infected cells , free virus population , and CTL population . Fig. 1 represent the model graphically. The uninfected epithelial cells are generated at a rate of and die at a rate of . The disease transmission rate is and the rate of inhibition of viral replication due to non-lytic antiviral activity is . The death rate of the infected cells is and the killing rate of the infected cells due to the lytic factor of CTL is denoted by .

Fig. 1

Transmission dynamics of COVID-19 within the host of the mathematical model (5).

The CTL responses are expected to proliferate at a rate proportionate to the infected cells and CTL responses population. Here we consider that CTL is proliferated at a rate of and is the maximum CTL responses at which proliferation shuts off. The removal rate of CTLs is . The strength of the non-lytic component is represented by and represents the strength of lytic components. We have considered that the antiviral drug blocks the infection at a rate and inhibits the viral production at a rate of .

The changes of uninfected cells, infected cells, virus population and CTL responses in presence of antiviral drug effect over the time is described by the following set of differential equations: with non-negative initial condition

Here and are the effectiveness of the antiviral drug with . represents no antiviral drug effect while represents that the drug is 100% effective.

To study the nonlocal property of the model (5), we use the fractional-order differential equations. Since the present state of a cell population not only depends on its current state but also depends on all of its historical states, hence fractional-order differential equation is applicable in the study of the dynamics of a disease.

The model (5) in form of fractional order differential equations is as follows with non-negative initial condition are the initial conditions and the Left-Caputo fractional derivative is indicated as . Here we consider all the parameter values are non-negative. Also, all the state variables in our model (7) are non-negative for

antiviral drug block infection at a rate and inhibit the viral production at a rate of

Mathematical analysis

The study of existence and uniqueness of the mathematical model (7) with the initial condition (8) is discussed in this section.

In this section we write the model (7) in the form as follows; where, and with the initial conditions for the derivative of Left-sided Caputo. Here, , , , and are the right hand side of the model (7), that is as we get:

The function with dimension is a vector field.

Non-negativity solutions

In this subsection we have analyzed non-negativity solution of our model (7).

We define and . For later use, we now recall the following theorem (known as Generalized Mean Value Theorem) [25] with a corollary in order to prove our main results.

Suppose thatand, then for, the following result is truewhere, for all(Generalized Mean Value Theorem, [25] ).

Suppose that and , for . If , for all , then is non-decreasing for each . If , , then is non-increasing for each .

Local existence and uniqueness

In this subsection first we define the function , as follows Putting the values of , , we get

Writing, , and , we can derive from Eq. (15) the following:

Therefore, the existence condition for our fractional order model (7) can be proved with the help of the theorem [26] stated below.

Assuming that and the function satisfies the following three conditions:

is Lebesgue measurable with ,

is continuous with ,

there is a real-valued function such that for almost every and for all .

Then for , there exists at least one solution of the initial value problem (9) in the interval , .

Here of model (7) are continuous on , bounded and measurable on for all . From the relation (11) it is clear that satisfies the conditions of Theorem 2 where . Therefore, we can conclude that the model (7) has a solution for . Also by using the theorem [26] we can prove that the model has a unique solution.

All the assumption of the Theorem 2 hold. Assume that there exists a real-valued function such that for almost every and all . Then there exists a unique solution of the initial value problem (12) on with some positive number .

Using (11) for the model, we get with, . Since the model (7) fulfill the condition of the Theorem 3, thus, we can conclude that the system has a unique solution.

Global existence

Here, we study the global existence of the solution of (7) using the theorem (Theorem 3.1 of [26] stated below.

Suppose that the first two conditions of the Theorem 3 are true for the vector field and for two constants and , for all most every and each , then there exists a function on which satisfies the initial value problem (9) .

Using Eq. (11) along with norm property of norm, we get

Hence, the model (7) satisfies the global existence theorem with , . Hence the global existence of the model (7) is satisfied.

Equilibrium and stability analysis

In this section, we only deal with the non negative equilibrium’s of the model with there stability. The three steady state of the model (7) are given below  . : the Disease free equilibrium .

. : the CTL responses free equilibrium when with

and .

. The endemic equilibrium point which we get, and comes from the equation where,

Note that, if and then Eq. (15) has at least one positive real root. Again and then equation (15) can have two positive roots. For a feasible endemic equilibrium we also get

Stability analysis

For the stability of , the following theorem is obtained.

The disease-free equilibrium is locally asymptotically stable if and unstable if

To determine the local stability of , we calculate the Jacobian matrix of the model (7) at where

Now we consider be the eigenvalue of the matrix . From the characteristics equation corresponding to , there exist two eigenvalues which are real, negative that is and and the other two eigenvalues satisfy the equation where, We observe that implies . Hence, the quadratic equation (17) has two strictly negative real roots if , that is, if . Hence, we can conclude that if disease-free equilibrium is locally asymptotically stable and when the system becomes unstable . □

CTL responses free equilibrium point is stable locally asymptotically iff

The Jacobian matrix of the system (7) around , is given by

Now we consider be the eigenvalue of the matrix . From the characteristics equation corresponding to , we have one eigenvalue is real, negative that is and the other three eigenvalues can be expressed in the form as follows: where, Following the Lienard–Chipard test [27], [28], we claim that each root is negative or having negative real part if the following conditions are verified:

,

.

We observe that is always positive, thus the Lienard–Chipard conditions are satisfied and local stability of CTL responses free equilibrium point is established. □

The endemic equilibrium is asymptotically stable if each eigenvalue of satisfies

Denoting we have the following result.

If the discriminant defined above is positive i.e.  , then the infected equilibrium is asymptotically stable iff Routh–Hurwitz conditions written below are satisfied

If the discriminant defined above is negative i.e.  , then the infected equilibrium is asymptotically stable and

If and then the infected steady state is unstable.

Endemic

To determine the local stability of , we compute the Jacobian matrix of the system (7) around the DFE is given by where The characteristic equation of which is given by where, We now have the following proposition:

The endemic equilibrium is asymptotically stable if for each eigenvalue of the following condition holds

Let us define the following discriminant: , Using this, we have the following proposition.

If holds, then the endemic steady state is asymptotically stable if the following conditions hold

If and then for , is asymptotically stable.

If , , , , and then is unstable.

The fractional optimal control problem

The Optimal control approach is a convenient mathematical technique used to control a mathematical model of biological or medical progressions. Several mathematicians have established control based mathematical models to attain optimal drug dose.

Here the aim is to analyze the effect of control measures. For this we introduced two control functions and into the system (7). is used to block the infection and is applied to inhibits the viral production. Biologically, represents the drug influence in the model dynamics. Thus our control based model (7) becomes with initial conditions: , , , and represents the Caputo fractional derivative of order . The system (20) can be rewritten as the following compact form: Since the aim is to minimize the infected epithelial cells, we define our objective function as below where represents the cost of the drug to block infection and are cost of the drug to block infection the weight constants on benefit of the cost. Here A and B represents the positive balancing coefficients (i.e. weights) that regularize the optimal controls. The quadratic expressions for the controls included in (20) indicate nonlinear costs potentially arising at high intervention levels [29].

We find the optimal controls to minimize the objective function (21) over the set subject to the state system (20) such that

Optimal conditions derivation

Fractional order control problem (FOCC) as in (20) describes the fractional dynamic system (FDS). General formulation of the optimality conditions of the FOCP is derived below using the methodology proposed by Agarwal [30].

We rewrite our fractional order control system as below: assuming as the state vector and as the control vector. We further write the functional as Thus in a compact form, the optimal control problem is the following: subject to the state system The co-state vector satisfy the following relation: and the optimal control pair is determined using the following relation:, The FOCP with Caputo fractional derivatives gives the Euler–Lagrange optimality conditions (25), (26), and (27). Note that when , the FOCP becomes the classical optimal control problem.

Using the above relations we solve our optimal control problem (21). In order to find a solution of the FOCP, the Hamiltonian function is considered: where , , , . We can minimized the objective function (21) by using the Euler–Lagrange conditions satisfying the optimality conditions (25), (26), (27).

(20) represents the state system. By using Eqs. (22), (24) the costate system can be obtained as: with the boundary conditions: .

By using (27), (28), we can obtain the optimal function stated below: The compact form of the optimal control can be written as The solution of the FOCP can be obtained by replacing and respectively by and in then system (20) and system (29).

Numerical findings

Here, we provide the numerical simulations on the basis of our analytical findings of the previous sections. To solve our model (7) we use the iterative scheme [31].

The above equation, last term of right hand side represents the memory. Here is defined as and , and , , and are the positive initial conditions.

The solution of the system (7) for the different values of represents by Fig. 2. Here we have observed that in case of integer order system (i.e., ) the system shows its periodic nature. But in case of fractional order () the system converges to its steady state in shorter time.

Fig. 2

The system trajectories (7) varies with . Other parameters are taken from Table 1.

Table 1

Parameter values of the system (7) and (20).

Parameter	Definition	Value
Π	Production rate of uninfected epithelial cell	5
β	Disease transmission rate	0.0001
δ_T	Death rate of epithelial cells	0.1
δ_I	Death rate of infected epithelial cells	0.1
δ_V	Removal rate of virus	0.1
δ_C	Decay rate of CTL	0.1
p	Efficacy of nonlytic component	0.01
q	Efficacy of lytic component
κ	Number of new virus produced	10–100
γ	CTL proliferation rate	0.22
C_max	Maximum proliferation of CTL	100
ε_1	Drug effectiveness of blocking infection	0–1
ε_2	Drug effectiveness of inhibiting infection	0–1

The solution of the system (7) for the different values of non-lytic factor () represents by Fig. 3, Fig. 4 shows the solution trajectories for different values of lytic factor . These numerical findings show two basic results.

Fig. 3

The system trajectories (7) for different value of with the set of parameters as in Table 1.

Fig. 4

The system trajectories (7) varies with with the parameters values from Table 1.

Fig. 3 shows that decreasing the value of viral transmission rate by nonlytic effector mechanism is useful to the epithelial cells in case of SARS-CoV-2 infection. Also if the killing rate by a lytic mechanism of infected cell increases, then lytic effector can be beneficial to control the virus replication. Fig. 4 shows that for low lytic effectors () the system attains a periodic nature. But if the lytic effector increases to , the system moves to its stable region within a short time. Thus we can apply these findings to fractional order optimal control study.

To study the fractional-order control problem (20), we numerically solve the state system (20) as an initial value problem and the costate system (29) as a boundary value problem. Thus forward iteration method has been used in solving the state system (20) and a backward iteration solves the costate system (29) with the help of MATLAB source code (see [6] for more detail).

Fig. 5, Fig. 6 is plotted with memory in presence of antiviral drug effects. Fig. 5 shows the path of the model variables in the presence of control. Here we observe that the uninfected epithelial cells attain their maximum value after days and infected and the virus population reduces with the effect of antiviral drugs. Also proper choice of drug weightage results in the extinction of infection. Fig. 6 shows the optimal control profile of the antiviral drugs and .

Fig. 5

Numerical solution of FOCP.

Fig. 6

Optimal control profile of control variables.

Discussion

CTL responses play a crucial role both in the death rate of infected cells due to viral infection and on the infection rate via the lytic and nonlytic mechanisms. In order to explore the effect of CTL responses in presence of the antiviral drugs, we have proposed a fractional-order mathematical model considering human host cells as model variables and analyzed the biological dynamics of the system analytically and numerically.

In the case of a fractional differential equation, the memory kernel declines as a law of power. However, the memory kernel converts the Dirac Delta function in the case of the ordinary differential equation. Thus the mathematical model enables us to predict the optimal drug dosage. Thus, by selecting the appropriate fractional index following the actual clinical records best treatment regimen may be obtained.

In this article, we have proved that the solution of the model is nonnegative and bounded. We have found that there exist three equilibriums, namely, the disease-free equilibrium , the CTL responses free equilibrium , and the endemic equilibrium . Here we have found the basic reproduction ratio of . It is clearly observed that if the infection rate reaches above a certain value, then becomes greater than the unity and the system attains its endemic state. On the other hand, if the death rate of infected cells attains above a certain value, then reduces than unity, and the system moves to its disease-free state which results in wiped out of infected cells. Analytically, we have derived that when , the disease-free state is asymptotically stable and for it is unstable, and the endemic equilibrium becomes feasible.

Using numerical simulation by changing the values of the fractional-order (Fig. 2), we study a better understanding of the biological process considered here. Reducing virus reproduction rate by nonlytic mechanism and increasing the death rate of infected cells by lytic mechanism is always beneficial to control the disease (Fig. 3, Fig. 4). Our numerical findings reveal that in the presence of antiviral drug therapy the lytic and nonlytic mechanism of CTL responses acts more successfully. We also study the control-theoretic methodology in view of optimal antiviral drug dosing in our proposed system. From the fractional-order optimal control-theoretic study along with numerical findings (see Fig. 5, Fig. 6) it is clearly observed that optimal treatment procedures can reduce the viral load and eradicate the infection. From our theoretical and numerical findings, we can predict that the lytic and nonlytic effector of the CTL responses plays a pivotal role in governing and eliminate the disease in presence of antiviral drug therapy.

In this article, we have revealed that the memory of cells plays a pivotal role to sustain the overall stability for long-term dynamics. It also gives us a more realistic way to reveal the viral dynamics. Thus regulating the memory of cells helps the system to explore more effectively. We have analyzed the non-negativity and stability of the system according to the value of . Numerical simulations confirmed our analytical findings. At the same time, the optimal drug dosing into the fractional-order system plays a crucial role to reduce the viral as well as infected cell load. Hence, the FDE model is more suitable for understanding the host cell dynamics and CTL responses against viruses in comparison to the ODE system. Further, this study helps us to develop effective treatments against COVID-19 infections in order to reduce the patients’ basic reproduction number to less than one. This study also leads for activation of immune response represents the main edge for reducing and eliminating viral load and infected cells from the human body.

CRediT authorship contribution statement

Amar Nath Chatterjee: Conceptualization, Funding acquisition, Methodology, Resources, Supervision, Visualization, Writing - review & editing. Fahad Al Basir: Formal analysis, Investigation, Project administration, Software, Validation, Writing - original draft. Muqrin A. Almuqrin: Formal analysis, Validation, Writing - review & editing. Jayanta Mondal: Conceptualization, Funding acquisition, Methodology, Resources, Supervision, Visualization, Writing - review & editing. Ilyas Khan: Formal analysis, Investigation, Project administration, Software, Validation, Writing - original draft.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.