key: cord-0790158-tvukkd7v
authors: Chhetri, Bishal; Bhagat, Vijay M.; Vamsi, D.K.K.; Ananth, V. S.; Prakash D, Bhanu; Mandale, Roshan; Muthusamy, Swapna; Sanjeevi, Carani B
title: Within-host mathematical modeling on crucial inflammatory mediators and drug interventions in COVID-19 identifies combination therapy to be most effective and optimal
date: 2021-04-30
journal: Alexandria Engineering Journal
DOI: 10.1016/j.aej.2020.12.011
sha: d166e3064c0c0b033a9b346fe0abee3b393f4f32
doc_id: 790158
cord_uid: tvukkd7v

The unprecedented Covid-19 pandemic has resulted in more than 14.75 million infections and 6, 10, 839 deaths in 212 countries. Appropriate interventions can decrease the rate of Covid-19 related mortality. Fast track clinical trials around the world are addressing the efficacy of individual pharmaceutical agent acting at various stages of pathogenesis. However, lessons learnt while dealing with past viral epidemics mandates, simultaneous use of such drugs in combination amongst different populations. Mathematical modelling studies can be extremely helpful in understanding the efficacy of drugs both individually and in combination. The current within-host mathematical model studies the natural history of Covid-19 in terms of complex interplay of virus replication and behaviour of host immune response. Additionally it studies the role of various drugs at various stages of pathogenesis. The model was validated by generating two-parameter heat plots, representing the characteristics of Covid-19, the sensitivity analysis identified the crucial parameters. The efficacy of interventions was assessed by optimal control problem. The model dynamics exhibited disease-free equilibrium and the infected equilibrium with their stability, based on the reproduction number R 0 , the transcritical bifurcation observed at R 0 = 1 . The burst rate and the natural death rate of the virus were observed as most significant parameters in the life-threatening Covid-19 pneumonia. The antiviral drugs affecting viral replication and those modulating the immune response, reduce the infected cells and viral load significantly. However, the yield was optimal and most effective when the combination therapy involving one or more antiviral and one or more immunomodulating drugs were administered together. These findings may help physicians with early decision making in treatment of life-threatening Covid-19 infection.

The Novel Corona Virus and its associated symptoms have involved more than 212 countries and been declared a pandemic by WHO. This unprecedented pandemic resulted into more than 14.75 million infections and 6, 10, 839 deaths as of July 20th, 2020, and continues to worsen [2] . The causative agent of Covid-19 was identified as Severe Acute Respiratory Syndrome Coronavirus-2 (SARS-CoV-2). Out of four coronavirus genera (a; b; c; d), b CoV strain showed 88% identity of genetic sequence with two bat derived SARS corona viruses (bat-SL-CoVZXC45,bat-SL-CoVZXC21) and 50% with the virus causing Middle East Respiratory Syndrome(MERSCoV) [28] . Therefore the functional mechanism in pathogenesis closely resembles SARS CoV and MERSCoV.

Human to human transmission of SARS-CoV-2 occurs either through droplet infection or direct contact from an infected host. Transmission from asymptomatic carriers and through faeco-oral route have also been reported. It is reported that Angiotensin Converting Enzyme-2 (ACE-2) plays crucial role in providing binding site to the viral structural protein(Spike protein-S) on the cell surface and subsequent entry into the host cell [35] . The SARS-CoV-2 have much higher affinity towards ACE-2 receptors as compared with its previous counterparts(SARS-CoV and MERS-CoV). ACE-2 is expressed in alveolar Type-2 cells(AT2), cholangiocytes, colonocytes, keratinocytes, endothelial cells of the esophagus, stomach, ileum, rectum and proximal tubules of kidneys. AT2 secret surfactant, which reduces surface tension preventing the collapse of alveoli and plays crucial role in oxygen diffusion across lungs and blood vessels. Viral antigens are presented by antigen presenting cells (APC) in Human Leucocyte Antigen (HLA) to T cells, resulting in generation of cytotoxic T cell (also known as T-killer cell, cytotoxic T-lymphocyte, CD8+) [23] . The virus enters the cell by fusion of its membranes with host cell and begins transcription with ssRNA acting as a template. Synthesis of viral proteins takes place in the cytoplasm of pneumocytes, and new copies are released by budding. These copies are ready to infect new cells which is confirmed by presence of abundant viral antigens in the cytoplasm of the pneumocytes(AT2) in case of SARS-CoV [45] . Viremia was associated with very high levels of IL-6 especially among severe cases of Covid-19, leading to increased vascular permeability and organs impairment [11] . Acute Respiratory Distress Syndrome(ARDS) is the common immunopathological event for all the aforesaid corona viral diseases. One of the main mechanisms in causation of ARDS is a deadly uncontrolled systemic inflammatory response due to the release of large quantity of pro-inflammatory cytokines and chemokines by immune effector cells. These cytokines are identified as IFN-a, IFN-c, IL-1b, IL-6, IL-2, IL-18, IL-33, TNF-a, TGF-b and chemokines as CCL-2, CCL-3, CCL-5, CXCL-8, CXCL-9, CXCL-10. Severe infections correlated high levels of IL-6, IFN a, CCL-5, CXCL-8 and CXCL-10 [11, 23] . Furthermore, the viral load was noted to be crucial in determining severity of the disease and strongly correlated with the lung injury Murray score [27] . The cytokine storm is a violent attack by the immune system causing ARDS, multiorgan failure and eventually death [23] .

Within-host mathematical models have gained considerable focus on decision making for therapeutic regimen and prophylactic interventions in several infectious and non communicable diseases. Works such as [21, 32] used within-host mathematical models to simulate the combination chemotherapy including reverse transcriptase inhibitors and protease inhibitors for HIV infection. These studies highlighted that the survival time is proportional to the CD4 + cells. Now, the combination of these interventions is mainstay of the treatment of HIV infection and survival outcome focused on the counts of CD4 + cells. Another example of successful application of within-host mathematical modeling in medical understanding is of within-host modeling studies done in [37] . The authors in this work described the determinants for severity of dengue disease (Dengue shock syndrome) due to cytokine storm and their role to predict the onset of severe life threatening disease.

Therefore, in order to identify possible interventions to break the natural history of disease, it is crucial to identify the interplay of viral multiplication and the response of host immune system in the form of inflammatory mediators in the dynamics of pneumonia inflicted by Covid-19. Carefully designed antiviral interventions acting at specific stages of pathogenesis and preventing complications, will be a vital breakthrough in reducing Covid-19 related mortality.

The clinical trials testing pharmacological interventions have several issues such as inadequate patient numbers to meet the sample size, internal and external variation of the behaviour of the virus among different subgroups like gender, race, nutritional status etc. These issues not only make the interventions time consuming but also do not adequately address potential confounders and biases.

In this context, within-host mathematical modelling can be extremely helpful in understanding the pathogenesis of the disease and the efficacy of interventions. The current study was aimed to develop mathematical models to understand the natural history of Covid-19 in terms of crucial inflammatory mediators and the behaviour of host immune response further to establish the role of various pharmaceutical interventions in reducing viral load and infected cells. Some of the mathematical modelling studies that deal with transmission and spread of COVID-19 at the population level can be found in [10, 22, 26, 44] . In the recent work [16] , an in-host modelling study deals with the qualitative characteristics and estimation of standard parameters of corona viral infections. Modelling the withinhost dynamics of Covid-19 involving the crucial biomarkers, which is being attempted here is the first of its kind for Covid-19.

2. Objectives of the proposed study 1. To study the natural history of Covid-19 with reference to the pathogenesis and changes in the dynamics of viral and infected populations and its correlation to the clinical condition. 2. To investigate the role of pharmaceutical interventions by incorporating them as controls at specific sites in the pathogenesis. 3. To study and compare the dynamics of infected cells and viral load with and without these control interventions.

Evidence before the study: After a thorough literature review in the relevant databases, it was noticed that current evidences were specific to the available intervention for Covid-19, but there were only a few instances for use of combination therapy. Also, there were limited scientific description for pathogenesis and the natural history of Covid-19.

Added value of this study: This is the first study to describe pathogenesis, with focus on the interplay between virus replication and immune response with and without interventions. The study revealed that immune modulation interventions required to be used in combination with chemotherapeutic agents which controls viral replications for optimal desired results.

Implications of all the available evidence: The study will help treating physicians with decision making and can be referred to by researchers, and epidemiologists in formulating appropriate policies and guidelines.

A mathematical model dealing with natural history was formulated and the stability theory of dynamical systems was used to understand the asymptotic behaviour of the models. The technique of two parameter heat plot was used to validate the model and sensitivity analysis was done by varying the model parameters in a specific range and for each varied value of the parameter, the infected population is plotted with respect to time. Next an within-host model incorporating antiviral drug interventions was studied as an optimal control problem. The Fillipov Existence Theorem was used to obtain an optimal solution. Numerical simulations involving with and without control interventions were performed.

From the above discussed pathogenesis, it can be understood that the study of both the pneumocytes (healthy and infected) and virus population levels and their changes over the time due to inflammatory mediators play a crucial in understanding the dynamics of Covid-19 associated pneumonia.

Motivated by these we first consider the following mathematical model.

Model 1: Model without Interventions/Medication The purpose of drug interventions can be two fold, the first to target the virus replication cycle and the second based on immunotherapy approaches either aimed to boost innate antiviral immune responses or alleviate damage induced by dysregulated inflammatory responses. Based on this the therapeutic agents for viral infections can be divided into two categories each serving the designated purpose [41] .

Drugs such as remdesivir, favipiravir inhibit RNAdependent RNA polymerase and drugs ivermectin, lopinavir/ritonavir inhibit the viral protease there by reducing viral replication [41] . On the other hand clinical trials such as phase I trial (NCT04280224) in China aim to enhance the innate immune system [3] . Interventions for improving the acquired immune response mediated via activation of T and B lymphocytes that help in reducing viral load are discussed in [4, 38] .

Motivated by the above clinical findings we consider a control problem with the following drug interventions as controls.

Here the controls u 11 and u 12 incorporate the effect of drug interventions which enhance the innate and acquired immune response which in turn leads to decrease in the infected population and viral load. The control u 2 incorporates the effect of drug interventions which prevent viral replication thereby reducing the virus birth rate.

Model 2: Model with Interventions/Medication as Controls

In this section we consider model 1 which deals with the natural history and course of infection. We study its dynamics.

The positivity and boundedness of the solutions of the model 1 is fundamental for further analysis. We establish these now. Positivity: We now show that if the initial conditions of the system (1)-(3) are positive, then the solution remain positive for any future time. Using the (1)-(3), we get,

dV dt V¼0 ¼ aI P 0: Thus all the above rates are non-negative on the bounding planes (given by S ¼ 0; I ¼ 0, and V ¼ 0) of the non-negative region of the real space. So, if a solution begins in the interior of this region, it will remain inside it throughout time t. This happens because the direction of the vector field is always in the inward direction on the bounding planes as indicated by the above inequalities. Hence, we conclude that all the solutions of the system (1)-(3) remain positive for any time t > 0 provided that the initial conditions are positive. This establishes the positivity of the solutions of the system (1)-(3). Next we will show that the solution is bounded.Boundedness:

with the assumption that x > a and l ¼ l 1 .

Here the integrating factor is e lt . Therefore after integration we get, N t ð Þ 6 x l þ ce Àlt . Now as t ! 1 we get,

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

System (1)-(3) admits two equilibria namely, the infection free equilibrium E 0 ¼ x l ; 0; 0 and the infected equilibrium

The basic reproduction number was calculated using the next generation matrix method [15] and the expression for R 0 is given by

8. Local dynamics and global dynamics of the model 1

The jacobian matrix of the system (1)-(3) at the infection free equilibrium E 0 is given by,

The characteristics equation is given by,

The first root of Eq.

Using the relation (9) the roots of the quadratic part of Eq. (8), are given by,

There are two cases that we need to consider here. Case I: When R 0 < 1 When R 0 < 1 there are further two subcases, (a):

Hence all the eigen values of the characteristics Eq. (8) are negative. Therefore the infection free equilibrium point E 0 is asymptotically stable. Sub-case (b): When A 2 þ 4 R 0 À 1 ð Þ D < 0 the eigenvalues of the quadratic part of Eq. (8) are complex conjugates with the negative real parts. Therefore we again have E 0 to be locally asymptotically stable.

Hence we conclude that E 0 is locally asymptotically stable provided R 0 < 1.

Case II: When R 0 > 1 For the case R 0 > 1, the characteristics Eq. (8) has two negative eigenvalues and one positive eigenvalue. Therefore whenever R 0 > 1 the infection free equilibrium E 0 is unstable.

With the definition of R 0 the infected equilibrium is given by,

Therefore the infected equilibrium exists only if R 0 > 1, otherwise E 1 will become negative which does not make sense.

The jacobian matrix of the system (1)-(3) is given by,

The characteristic equation of the jacobian J evaluated at E 1 is given by,

. Therefore if we substitute k ¼ Àkin Eq. (10) using Descartes rule of sign change we get all the roots of (10) to be negative. Hence we conclude that the infected equilibrium point E 1 exists and remains asymptotically stable provided R 0 > 1.

To establish the global stability of the infection free equilibrium E 0 we make use of the method discussed in Castillo-Chavez et al. [8] .

Theorem 8.1. Consider the following general system,

where X denotes the uninfected population compartments and Y denotes the infected population compartments including latent, infectious etc. Here the function G is such that it satisfies G X; 0 ð Þ¼0. Let U 0 ¼ X 0 ; 0 À Á denote the equilibrium point of the above general system. If the following two conditions are satisfied then the infection free equilibrium point U 0 is globally asymptotically stable for the above general system

Now we will prove the global stability of E 0 ¼ x l ; 0; 0 of system (1)-(3) by showing that system (1)-(3) can be written as the above general form and both the conditions A 1 and A 2 are satisfied.

Comparing the above general system (11) to the system (1)-(3) the functions F and G are given by

From the stability analysis of E 0 , we know that U 0 is locally asymptotically stable iff R 0 < 1. Clearly, we see that

The integrating factor is e lt and therefore after performing integration on the above Eq. (12) we get,

which is independent of c. This independency implies that X 0 ¼ x l 1 is globally asymptotically stable for the subsystem dS dt ¼ x À lS. So, the assumption A 1 is satisfied. Now, we will show that assumption A 2 holds. First, we will find the matrix A. As per the theorem,

Thus both the assumptions A 1 and A 2 are satisfied and therefore infection free equilibrium point E 0 is globally asymptotically stable provided R 0 < 1. Similarly The global stability of E 1 when R 0 > 1 can be proved by considering the Lyaponov function.

We now use the method given by Chavez and Song in [9] to do the bifurcation analysis.

Theorem 9.1. Consider a system,

Let 0 be the equilibrium point of the system such that f 0; / À Á ¼ 0; 8 / 2 R. Let the following conditions hold:

1. For the matrix A ¼ D X f 0; 0 À Á , zero is the simple eigenvalue and all other eigenvalues have negative real parts. 2. Corresponding to zero eigenvalue, matrix A has nonnegative right eigenvector, denoted as u and non-negative left eigenvectors, denoted as v.

Let f k be the k th component of f. Let a and b be defined as follows -

Then local dynamics of the system near the equilibrium point 0 is totally determined by the signs of a and b. Here are the following conclusions:

1. If a > 0 and b > 0, then whenever / < 0 with j/j ( 1, the equilibrium 0 is locally asymptotically stable, and moreover there exists a positive unstable equilibrium. However when 0 < / ( 1; 0 is an unstable equilibrium and there exists a negative and locally asymptotically stable equilibrium.

2. If a < 0; b < 0, then whenever / < 0 with j/j ( 1; 0 À is an unstable equilibrium whereas if 0 < / ( 1; 0 is locally asymptotically stable equilibrium and there exists a positive unstable equilibrium. 3. If a > 0; b < 0, then whenever / < 0 with j/j ( 1; 0 is an unstable equilibrium, and there exists a locally asymptotically stable negative equilibrium. However if 0 < / ( 1; 0 is stable, and a there appears a positive unstable equilibrium. 4. If a < 0; b > 0, then whenever / changes its value from negative to positive, the equilibrium 0 changes its stability from stable to unstable. Correspondingly a negative equilibrium, unstable in nature, becomes positive and locally asymptotically stable.

Applying the Theorem 3.2 to our Model 1: In our case, we have x ¼ S; I; B ð Þ2R 3 where x 1 ¼ S; x 2 ¼ I and x 3 ¼ V. Let us consider b (transmission rate of the infection) to be the bifurcation parameter.

We know that

can be written as follows:

The disease free equilibrium point E 0 is given by,

Þ denote the Jacobian matrix of the above system at the equilibrium point x à and R 0 ¼ 1. Now we see that, The characteristic polynomial of the above matrix is obtained as

Hence, we obtain the first eigenvalue of (14) as

The other eigenvalues k 2;3 of (14) are the solution the following equation,

substituting the expression for b à from (13) in (15) we get,

The eigen values of (16) are k 2 ¼ 0 and k 3 ¼ À p þ q þ l þ l 1 ð Þ Hence, the matrix D x f x à ; b à ð Þ has zero as its simple eigenvalue and all other eigenvalues with negative real parts. Thus, the condition 1 of the Theorem 3.2 is satisfied.

Next, for proving condition 2, we need to find the right and left eigenvectors of the zero eigenvalue (k 2 ). Let us denote the right and left eigenvectors by u and v respectively. To find u, we where u ¼ u 1 ; u 2 ; u 3 ð Þ T . As a result, we obtain the system of simultaneous equations as follows:

By choosing u 3 ¼ l in the above simultaneous Eq. (17)-(19) we obtain

Therefore, the right eigenvector of zero eigenvalue is given by

The simultaneous equations obtained thereby are as follows:

Àb

Therefore solving the above simultaneous Eqs. (20)- (22) we

Hence, the left eigenvector is given by

Now, we need to find a and b. As per the Theorem 3.2, a and b are given by

Expanding the summation in the expression for a, it reduces to

where partial derivatives are found at x à ; b à ð Þ. Since we know u and v we only need to find the partial derivatives in above expression. They are found to be

Substituting these partial derivatives along with u and v in the expression of a, we get,

Next, expanding the summation in the expression for b, we get,

We notice that condition (iv) of the theorem is satisfied. Hence, we conclude that the system undergoes bifurcation

Thus, we conclude that when R 0 < 1, there exists a unique disease free equilibrium which is globally asymptotically stable and negative infected equilibrium which is unstable. Since negative values of population is not practical, therefore we ignore it in this case. Further, as R 0 crosses unity from below, the disease free equilibrium point loses its stable nature and become unstable, the bifurcation point being at b ¼ b à implying R 0 ¼ 1 and there appears a positive locally asymptotically stable infected equilibrium point. There is an exchange of stability between disease free equilibrium and infected equilibrium at R 0 ¼ 1. Hence, a trans-critical bifurcation takes place at the break point b ¼ b à . Table 2 Parameter values for global stability of E 0 . Table 3 Parameter values for the stability of E 1 . Table 4 Values of parameters for transcritical bifurcation. 

We now numerically depict and verify the results obtained in Local and Global dynamics and Bifurcation Analysis sections. The theoretical results obtained are validated for a set of model parameters using MATLAB software. The values of x; l and l 1 ; a are approximated and chosen to be from [31, 16] respectively. The rest of the parameter values of the model are estimated minimizing the root mean square difference between the model predictive output and the experimental data from [17, 36] . All the parameter values of the model 1 are summarized in the following Table 1 .

The following table shows different values of parameters chosen so as to obtain the phase portraits of the system (1)-(3) for the case R 0 < 1. With the parameter values from Table 2 we have R 0 ¼ 0:77 < 1. Since R 0 < 1, the disease free equilibrium point, E 0 ¼ 20; 0; 0 ð Þis globally asymptotically stable for the system (1)-(3). The plots in Fig. 1 for different initial conditions depict the global stability of the infection free equilibrium E 0 of the system (1)-(3).

We know that the infected equilibrium E 1 exists only if R 0 > 1 and it is also locally asymptotically stable when R 0 > 1. For the parameter values in the following Table 3 the value of R 0 was calculated to be 1.929 and E 1 to be 25:9; 2:44; 1:85 ð Þ . 2 demonstrates that E 1 is locally asymptotically stable whenever R 0 > 1.

In this section, an exchange of stability between the disease free equilibrium and infected equilibrium as R 0 crosses unity from below is shown. Here the parameter x was varied in the interval 0; 5 ð Þ in steps of 0.1. The other parameters chosen were as given in the following Table 4 . For R 0 varying in the interval 0; 3 ð Þ, Fig. 3 depicts the occurence of the transcritical bifurcation at R 0 ¼ 1.

The characteristic of the disease Covid-19 are as follows:

(1) Target cells(monocytes) deplete to approximately 66 % to the original level (from 6x10 8 to 4x10 8 in severe cases) [36] .

(2) The median incubation period of SARS-Cov-2 is approximately 4-5 days and peak viremia occurs within 5-6 days of disease onset. [39] .

In this section we vary the parameters x ¼ d 1 þd 2 þd 3 þd 4 þd 5 þd 6 ð

Þand

Þ and do a two parameter heat plot to validate our model 1.

The parameters x and y were varied in the interval 0; 5 ð Þ and model 1 was able to reproduce both the above characteristics of the Covid-19. Other fixed parameter values were taken from the following Table 5 .

We choose the parameters x and y in the x-axis and y-axis respectively and vary them to check for reproduction of the characteristics (1) and (2).

In Fig. 4 (a) model 1 is able to recover the first characteristic (1) exactly for the parameter values in the region pointed by the arrow. In this region with yellow colour the final fraction of uninfected cells after infection lies between 60 À 70 ð Þ % and for the parameter values in this region the model 1 is able to recover the first characteristic (1).

The region with yellow colour in Fig. 4(b) pointed by the arrow is the region where model 1 is able to recover the second characteristics (2). From [39] the peak viremia occurs approximately during the second week of disease onset. In Fig. 4(b) time to peak viremia lies between 8-16 days and hence for the parameter values in this region model 1 is able to recover the second characteristics (2).

From the previous sections, it is clear that the infected cell population and the virus load die out when R 0 < 1. Therefore it is important to control the model parameters in a manner which will make R 0 less than one. Thus, determining the intervals in which the model parameters are sensitive becomes vital. As each parameter is varied in different intervals, the infected cell population, mean infected cell population and the mean square error are plotted with respect to time. These plots are used to determine the sensitivity of the parameter. The different intervals chosen are given in the following Table 6 . For all the plots in this section, the time scale is the following: x-axis: 10 units = 1 day, y-axis: 1 unit = 1 cell.

The results related to sensitivity of a, varied in three intervals as mentioned in Table 6 , are given in Fig. 5 . The plots of infected population for each varied value of the parameter a per interval, the mean infected population and the mean square error are used to determine the sensitivity. We conclude 

Interval

Step Size Other Parameters from these plots that the parameter a is sensitive in interval II and insensitive in I and III.

The results related to sensitivity of l 1 , varied in three intervals as mentioned in Table 6 , are given in Fig. 6 . The plots of infected population for each varied value of the parameter l 1 per interval, the mean infected population and the mean square error are used to determine the sensitivity. We conclude from these plots that the parameter l 1 is sensitive in interval I and II and insensitive in III. In similar lines, the sensitivity analysis is done for other parameters. The results are summarized in Table 7 . The corresponding plots are given in Appendix -A owing to the brevity of the manuscript.

In this section, we will formulate an optimal control problem for the model 2 with drug interventions as control. The controls to be considered are:

1. Drug intervention to modulate the immune response: The innate response mediated through NK cells, macrophages, neutrophils, mast cells, basophils, eosinophils etc. mainly acts to contain the pathogen. The acquired immune response mediated via activation of T and B lymphocytes targets the specific antigen. In severe Covid-19 the uncontrolled inflammatory innate response and impaired Table 6 . The plots depict the infected population for each varied value of the parameter a per interval along with the mean infected population and the mean square error in th.e same interval.

acquired immune response may lead to severe damage of the host cells locally as well as systematically. Therefore, u 11 t ð Þ incorporate the interventions which control the innate response and reduce number of infected cells, while u 12 t ð Þ incorporate the interventions towards improving T and B cell response in reducing viral load and interventions that decrease the viral load due to decrease in the infected cells.Therapeutic use of Natural Killer Cells for treatment of Covid-19 is already under clinical trial [34] . INF-1 and its subtypes such as INF-a was recommended by some workers also reflected in the guidelines by China Govt. for treatment of Covid-19 cases [3]. Also INF-a2b [13] was already registered for clinical trials for its therapeutic use in patients of Covid-19. The specific immune mediators such as IgG, IgM studied successfully and recommended in vitro experimentation. Also the T-cell response in the elimination of virus was extensively studied in SARS-CoV [4] . Chimeric Antigen Receptor T-Cell therapy was recom-mended whenever required during the treatment of Covid-19 patients [38] .Several other interventions such as Camostat mesylate, Nafamostat mesylate [43] , Chloroquine phosphate and hydroxychloroquine [42] , Cepharanthine (antiinflammatory compound)/selamectin (anti-helminthic)/ meoquine hydrochloride (used extensively for chemoprophylaxis in malaria) plays crucial role in reducing the viral load, were also evaluated and recommended by several authors. 2. Drug intervention to prevent viral replication: The pharmaceutical interventions targeting towards preventing viral replication which in turn reduces its multiplication rate. We use the control variable u 2 t ð Þ to denote this intervention.Several pharmaceutical agents interfering with the replication of virus such as Remdesivir [24] , Lopinavir/ritonavir [12] , Arbidol (Umifenovir) [14] , Favipiravir [6] are either already being used or recommended to use in treatment of Covid-19. Some naturally occurring immunomod- Table 6 . The plots depict the infected population for each varied value of the parameter l 1 per interval along with the mean infected population and the mean square error in th.e same interval.

ulators and phytochemicals such as Luteolin, Myricetin, Apigenin, Quercetin, Kaempferol, Baicalin, Wogonoside; which interfere with activation of NLRP3, also reported to reduce the viral replications and can be considered for treatment of Covid-19 especially in prevention of ARDS.

The set of all admissible controls is given by. U ¼ u 11 t ð Þ; u 12 t ð Þ; ð f u 2 t ð ÞÞ : u 11 t ð Þ 2 0; u 11 max ½ ; u 12 t ð Þ 2 0; u 12 max ½ ; u 2 t ð Þ 2 0; u 2 max ½ ; t 2 0; T ½ g Without medical interventions, u 11 ; u 12 ;; and u 2 are just constant parameters with u 2 ¼ 0. u 11 and u 12 have some value based on the inherent release of these cytokines and chemokines by the body. In the control problem, these variables are considered as functions of time. The upper bounds of control variables are based on the resource limitation based on avail-ability and the limit to which these drugs would be prescribed to the patients.

Since these antiviral drugs administered are not inherently present in the human body and being foreign particles to the body, there could be side effects once administered. Also cost of these drugs could be an issue that needs to be addressed.

Based on these we now propose and define the optimal control problem with the goal to reduce the cost functional

such that u ¼ u 11 t ð Þ; u 12 t ð Þ; u 2 t ð Þ ð Þ 2 Usubject to the system

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

is called the Lagrangian or the running cost.

Here, the cost function represents the number of infected cells and viral load throughout the observation period, and the side effects of the drug on the body. Effectively, we want to minimize the infected cells and the virus load in the body with the optimal medication that is also least harmful to the body. Since the drugs administered have multiple effects, the non-linearity for the control variables in the objective become justified [19] .

The admissible solution set for the Optimal Control Problem (23) and (24) is given by X ¼ S; I; V; u 11 ; u 12 ; u 2 ð Þ j S; I and V that satisfy 24 ð Þ8 u 2 U f g

We will show the existence of optimal control functions that minimize the cost functions within a finite time span 0; T ½ showing that we satisfy the conditions stated in Theorem 4.1 of [18] .

Theorem 16.1. There exists a 3-tuple of optimal controls u à 11 t ð Þ; u à 12 t ð Þ; u à 2 t ð Þ À Á in the set of admissible controls U such that the cost functional is minimized i.e.,

corresponding to the optimal control problem (23) and (24) .

Proof. In order to show the existence of optimal control functions, we will show that the following conditions are satisfied:

1. The solution set for the system (24) along with bounded controls must be non-empty, i:e., X -/. 2. U is closed and convex and system should be expressed linearly in terms of the control variables with coefficients that are functions of time and state variables.

3. The Lagrangian L should be convex on U and L S; I; V ; u 11 ; u 12 ; u 2 ð Þ P g u 11 ; u 12 ; u 2 ð Þ , where g u 11 ; u 12 ; u 2 ð Þ is a continuous function of control variables such that j u 11 ;u 12 ;u 2 ð Þ j À1 g u 11 ;u 12 ;u 2 ð Þ ! 1 whenever j u 11 ; u 12 ;u 2 ð Þ j!1 , where j:j is an l 2 0;T ð Þ norm. Now we will show that each of the conditions are satisfied:

1. From Positivity and boundedness of solutions of the system (24) all solutions are bounded for each bounded control variable in U.Also,the right hand side of the system (24) satisfies Lipschitz condition with respect to state variables.Hence, using the positivity and boundedness condition and the existence of solution from Picard-Lindelof Theorem [29] , we have satisfied condition 1. 2. U is closed and convex by definition. Also, the system (24) is clearly linear with respect to controls such that coefficients are only state variables or functions dependent on time. Hence condition 2 is satisfied. 3. Choosing g u 11 ; u 12 ; u 2 ð Þ¼c u 2 11 þ u 2 12 þ u 2 2 À Á such that c ¼ min A 1 ; A 2 ; A 3 f g , we can satisfy the condition 3.

Hence there exists a control 3-tuple u à 11 ; u à 12 ; u à 2 À Á 2 U that minimizes the cost function (23).

We will obtain the necessary conditions for optimal control functions using the Pontryagin's Maximum Principle [25] and also obtain the characteristics of the optimal controls.

The Hamiltonian for this problem is given by H S; I; V; u 11 ; u 12 ; u 2 ; k ð Þ :¼ L S; I; V; u 11 ; u 12 ; u 2 ð Þ þ k 1 dS dt þ k 2 dI dt þ k 3 dV dt

Here k = (k 1 ; k 2 ; k 3 ) is called co-state vector or adjoint vector. Now the Canonical equations that relate the state variables to the co-state variables are given by

Substituting the Hamiltonian value gives the canonical system

along with transversality conditions k 1 T ð Þ ¼ 0; k 2 T ð Þ ¼ 0; k 3 T ð Þ ¼ 0. Now, to obtain the optimal controls, we will use the Hamiltonian minimization condition @H @ui = 0, at u i ¼ u à i for i = 11, 12, 2.

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

In this section, we perform numerical simulations to understand the efficacy of multiple drug interventions. This is done by studying the effect of control on the dynamics of the system. These simulations also validate the theoretical results obtained in the previous section. The efficacy of various combinations of controls considered are:

1. Pharmaceutical substances which acts to improve the performance of host immune system towards reducing the number of infected cells and viral load. 2. Several drugs which prevents viral replication. 3. The therapeutic agents in combination of both the above types of mechanisms of intervention.

For our simulations, we have taken the total number of days as T ¼ 30. The parameter values considered are as follows: x = 10, l= 0.05, l 1 = 1.1, b = 0.005, a = 0.5.

We first solve the state system numerically using Fourth Order Runge-Kutta method in MATLAB without any interventions. We take the initial values of state variables to be S 0 ð Þ ¼ 3:2 Â 10 5 ; I 0 ð Þ ¼ 0; V 0 ð Þ ¼ 5:2 [5, 16] with the control parameters as constant values u 11 = 0.8866, u 12 = 0.56, u 2 = 0. Now, to simulate the system with controls, we use the Forward-Backward Sweep method stating with the initial values of controls and solve the state system forward in time. Follow- ing this we solve the adjoint state system backward in time due to the transversality conditions, using the optimal state variables and initial values of optimal control. Now, using the values of adjoint state variables, the values of optimal control are updated and with these updated control variables, we go through this process again. We continue this till the convergence criterion is met [25] . The positive weights chosen for objective coefficients are A 1 = 80, A 2 = 80, A 3 = 7480. A 3 is chosen high compared to A 1 and A 2 as the effort to stop the viral replication process is higher than the effort to enhance innate immune response as it is already existing in human body. Fig. 7 shows the change in the cell population S t ð Þ; I t ð Þ; V t ð Þ respectively with time. We observed that the susceptible cells reduce and the infected cells increase exponentially due to the increase in viral load over a period of time. Fig. 8a shows the behaviour of susceptible cells under all possible combinations of control. We observe that in the presence of immune boosting medication only, there is a slight reduction in the number of susceptible cells getting infected by the virus but when viral replication is prevented the impact is more. Even lesser damage is seen on susceptible cells when all control interventions are applied together.

In Fig. 8b the first frame shows how infected cells grow or decay with various combinations of controls. The increase in infected cells is slightly lesser in the presence of immunity boosters. There is a huge difference if viral replication preventing medicines are used. Here too, the best results are obtained only when all the controls are applied. The second frame in Fig. 8b gives detailed view of the cases when only antiviral replication medicines are used and all 3 controls are used. We see that when only control u 2 is used, the infected cells increase slowly, and reduce later but increase again after 20 days. This may be harmful to the patient. When all 3 controls are used, the infected cells reach peak around the 8th day and start decaying then on and even become nearly zero after 30 days. In Fig. 8c the first frame shows the viral load under all combinations of control interventions. When no medication is provided, there is exponential increase in the viral load but when only immunity boosting medication is provided (u 2 ¼ 0), there is a little reduction seen in the increase of viral load. The combinations of both these drug interventions show much better results. From the second frame in Fig. 8c it can be seen that when only anti-viral replication medication is provided (u 11 ¼ u 12 ¼ 0), the viral load becomes very less by the 5th day but due to the absence of immunity boosters, the viral load tends to increase around 25th day. The best results are shown when all the control interventions are administered together. The detailed view of the figure explains how the viral population tends to become nearly zero around 7th day and remains there throughout. These results are in line with the clinical findings discussed in [7] .

Based on the pathogenesis of Covid-19, two models are proposed. The first model deals with natural history and the course of infection while the second model incorporates the drug interventions. The results of these studies show that the disease system admits two steady states: one being the disease-free equilibrium and the other being the infected equilibrium. The dynamics of the system show that the disease takes it course to one of these states based on the reproduction number R 0 . Specifically, when R 0 < 1 the system tends to stabilize around the disease free equilibrium and when R 0 > 1 the system tends to stabilize around the infected equilibrium. The system also undergoes a trans-critical bifurcation at R 0 ¼ 1.

This result is in line with the conclusions made at the population level for Covid-19 [20] . From the sensitivity analysis was observed that the burst rate of virus particles and the natural death rate of the virus are sensitive parameters of the system. The proposed model 1 was validated using two-parameter heat plots that reproduce the characteristics of Covid-19.

Results from the optimal control studies suggest that the antiviral drugs that target on viral replication and the drugs that enhance the immune system response both reduce the infected cells and viral load when taken individually. In particular, the antiviral drugs that target viral replication seem to yield better results than the drugs that enhance the innate immune response. From Fig. 8c , it can be seen that on administering control intervention u 2 there is a substantial decrease in the viral load from day 2. This result validates the clinical findings in [7] which states that "Ivermectin, an FDAapproved anti-parasitic previously shown to have broadspectrum anti-viral activity in vitro, is an inhibitor of the causative virus (SARS-CoV-2), with a single addition to Vero-hSLAM cells 2 h post infection with SARS-CoV-2 able to effect 5000-fold reduction in viral RNA at 48 h."

When applied in combination, these drugs yield the best possible results. Hence, the optimal control strategy and the best drug regime would be the combination therapy involving one or more antiviral and one or more immunomodulating drugs administration which helps in patient's recovery faster. This is inline with the clinical observations in [1] which states that "In a randomized trial, people who received interferonbeta, ribavirin, and lopinavir/ritonavir had more-rapid clearance of virus and quicker clinical improvement than a control group who received just lopinavir/ritonavir.". For all the plots in this section, the time scale is the following:

x-axis: 10 units = 1 day, y-axis: 1 unit = 1 cell.

A.1. Parameter u 11 

Car t cell therapy during the covid-19 pandemic

Minimal within-host dengue models highlight the specific roles of the immune response in primary and secondary dengue infections

The fda-approved drug ivermectin inhibits the replication of sars-cov-2 in vitro

On the computation of reproduction number and its role in global stability

Dynamical models of tuberculosis and their applications

A mathematical model for simulating the phasebased transmissibility of a novel coronavirus

Detectable serum sars-cov-2 viral load (rnaaemia) is closely associated with drastically elevated interleukin 6 (il-6) level in critically ill covid-19 patients, medRxiv

Role of lopinavir/ritonavir in the treatment of sars: initial virological and clinical findings

Induction of pro-inflammatory cytokines (il-1 and il-6) and lung inflammation by coronavirus-19: antiinflammatory strategies

Arbidol combined with lpv/r versus lpv/r alone against corona virus disease 2019: A retrospective cohort study

The construction of next-generation matrices for compartmental epidemic models

-host modelling of covid-19 kinetics in humans, medRxiv

Virological assessment of hospitalized cases of coronavirus disease

Optimal control of an hiv immunology model

Modeling the dynamics of novel coronavirus (2019-ncov) with fractional derivative

A mathematical model of combined drug therapy of hiv infection

Early dynamics of transmission and control of covid-19: a mathematical modelling study

Molecular immune pathogenesis and diagnosis of covid-19

Rapid review for the anti-coronavirus effect of remdesivir

Calculus of variations and optimal control theory: a concise introduction

A conceptual model for the coronavirus disease 2019 (covid-19) outbreak in Wuhan, China with individual reaction and governmental action

Clinical and biochemical indexes from 2019-ncov infected patients linked to viral loads and lung injury

Genomic characterisation and epidemiology of 2019 novel coronavirus: implications for virus origins and receptor binding

The picard algorithm for ordinary differential equations in coq

Efavirenz or nevirapine in three-drug combination therapy with two nucleoside or nucleotide-reverse transcriptase inhibitors for initial treatment of hiv infection in antiretroviral-naı¨ve individuals

A micro-epidemic model for primary dengue infection

Modelling in vivo hiv dynamics under combined antiretroviral treatment

Multiple drug rescue therapy for hiv-infected individuals with prior virologic failure to multiple regimens

Immune responses in covid-19 and potential vaccines: Lessons learned from sars and mers epidemic, Asian Pac

Single cell rna sequencing of 13 human tissues identify cell types and receptors of human coronaviruses

Dysregulation of immune response in patients with covid-19 in

Koelle, minimal within-host dengue models highlight the specific roles of the immune response in primary and secondary dengue infections

Sars: systematic review of treatment effects

The trinity of covid-19: immunity, inflammation and intervention

Antiretroviral treatment of adult hiv infection: 2012 recommendations of the international antiviral society-usa panel

A review of sars-cov-2 and the ongoing clinical trials

Remdesivir and chloroquine effectively inhibit the recently emerged novel coronavirus (2019-ncov) in vitro

Identification of nafamostat as a potent inhibitor of middle east respiratory syndrome coronavirus s protein-mediated membrane fusion using the split-protein-based cell-cell fusion assay

A mathematical model for the novel coronavirus epidemic in wuhan, china

Sars coronavirus infection: pathology and pathogenesis of an emerging virus disease

The authors from SSSIHL acknowledge the support of SSSIHL administration for this work.

The outbreak of novel coronavirus in Wuhan, China marked the introduction of a virulent coronavirus into human society. On Feb. 11, 2020, the World Health Organization named novel corona viral pneumonia induced disease as Coronavirus disease , which is caused by Severe Acute Respiratory Syndrome coronavirus-2 (SARS-CoV-2). Although a lot of research is being done, effective approaches to treatment and epidemiological control are still lacking.In this context, the within-host mathematical modelling studies can be extremely helpful in understanding the natural history of this new disease, role of various interventions individually and in combination. The efficacy of the pharmaceutical interventions requiring animal and human trials are already underway. The results of drugs proven effective on the previous similar diseases such as SARS and MERS also available in hand. In this scenario, the clinical trials on various combination of drugs also opens an arena of research. Mathematical models are known to provide useful information in short period of time and more importantly in the non-invasive way. Therefore, at this juncture where several clinical trials on newer interventions are being done, the effect of these interventions in combination, felt the need of the hour to, provide optimal treatment options to the treating physicians. The current study revealed that, to control progression of the disease in an individual, interventions needs to be applied towards two important mechanisms; modulation of immune response and preventing the viral replication.The immune modulators mainly acts through either activation of innate or the acquired immune response. Among the available interventions, Natural Killer Cells, IFN-1 and its subtypes IFN-a, IFN-a2b, specific immune mediators such as IgG, IgM, Chimeric Antigen Receptor T-Cell therapy. The drugs preventing viral replication such as Remdesivir, Lopinavir/ritonavir, Arbidol (Umifenovir) Favipiravir are recommended for use in treatment of Covid-19.Conclusions from the optimal control studies suggest the following.

affect the way of spread of infection but they become effective only along with the drugs preventing viral replication. 2. Viral replication, when prevented, reduces the infected cells and the viral load drastically but the immune system needs to be sufficiently modulated if eradication of viral particle and complete cure is anticipated. 3. When all controls (interventions) are used effectively in combination, then the patient can be cured completely, with optimal adjustment of the dosage of drugs to maximize the efficacy with least adverse effects.The use of such combination interventions was already in use for other viral diseases such as for HIV [21, 30, 32, 33] . In the works [21, 32] the combination of reverse transcriptase inhibitor and protease inhibitor were tried with sustained levels of CD4+. Currently the combined chemotherapy named Highly Active Anti-Retroviral Therapy (HAART) is the most recommended chemotherapy for HIV infection [40] . In fact it was concluded in [33] that Multiple Drug Rescue Therapy (MDRT) induced a substantial antiviral response and in [30] it was stated that combination therapy had good benefits on initial treatment of HIV. Recently, this kind of combination intervention strategy is also being recommended for use in Covid-19 disease [12, 14] .This within-host modelling study involving the crucial biomarkers of Covid-19 is first of its kind for Covid-19 and the results will help treating physicians in decision making while treating the patients and also help researchers, epidemiologists in formulating appropriate policies and guidelines.

This work is an initial attempt to understand the basic dynamics of Covid-19 and its course. The consequences and the outcomes of the two main functions of antiviral drug interventions is modelled and discussed. Further research can be focused on incorporating the side effects of these drugs into the model. Future studies can also focus on the various other drug interventions along with their effects on Covid-19 dynamics.

The authors from SSSIHL and SSSHSS dedicate this paper to the founder chancellor of SSSIHL, Bhagawan Sri Sathya Sai Baba. The corresponding author also dedicates this paper to his loving elder brother D. A. C. Prakash who still lives in his heart and the first author also dedicates this paper to his loving father Purna Chhetri.

None.Appendix A. Sensitivity plots for other parameters A.5. Parameter u 12