Burden of cytokines storm on prognosis of SARS-CoV-2 infection through immune response: dynamic analysis and optimal control with immunomodulatory therapy.

Immune responses have a crucial role to play against SARS-CoV-2 virus as the adaptive and innate immune systems of the human body help restoring the body to a healthy stage by annihilating this deadly viral infection. Cytokines also play a significant role in modulating a balance between innate and adaptive immune responses but excess of it can have a detrimental affect on critically ill patients. Therefore, this paper is a novel attempt to formulate a within-host mathematical model showing the impact of cytokines storm on healthy cells. The dynamics of the system is analysed which involves basic reproduction number, steady state solutions and global dynamics for disease-free point and endemic equilibrium using geometric approach. Further, an optimal control problem is discussed considering immunomodulatory therapy (targeting cytokines signaling) as control using linear feedback control method to increase the level of healthy cells, which provides vitality for our system. Through numerical simulations, analytic solutions are validated followed by the curve-fit for the cytokines using real data and an optimization algorithm for optimal fit. Finally, sensitivity analysis for the basic reproduction number and the rate of change of healthy cells using Latin Hypercube Sampling method (LHS) is performed. Our finding suggests that immunomodulatory therapy (tocilizumab) can act as a key component to control cytokines storm for critically ill patients to restore the body to a healthy state.

Introduction

Since being declared as a global pandemic, the severe-acute-respiratory syndrome corona virus-2 has led to the death of millions. Even after numerous studies and researches, the pathophysiology of this deadly virus is far from being completely understood. Specially for people with other existing ailments, the SARS-CoV-2 virus sometimes makes recovery harder or even worse. Immune response has a major role to play in it and several biological and bio-physics literature are available for better understanding of the immune process in context to this disease. For instance, the evolution of nanotechnology and application of nanoparticle (NP) has also a significant contribution in influencing the immune cell behaviour. It is essential to have clarity on the interplay and accumulation of NPs by the ECMs (extracellular matrix). ECMs acts as a pool for signaling molecules including inflammatory mediators and has an important role in influencing inflammatory cytokines in inflamed tissues [1] and a mechanistic understanding can be seen in [2]. The cellular infection and virus replication activates the inflammatory response which leads to the release of pro inflammatory cytokines [3] and gives an amplified inflammatory response. These cytokines subsequently send signals and activate natural killers (NK) which are innate immune cells and the adaptive immune T-lymphocytes cells. T-cells further activate B- lymphocytes. These B-cells produce antibodies, namely IfM and IgG. Even though the protective functions of cytokines are ideal responses but there may be an extremely excessive response due to SARS-CoV-2 [4]. This uncontrolled response due to the infection is known as cytokines release syndrome or cytokines storm. Here, along with their task to kill the infected cells [5], the cytokines start attacking the healthy cells as well [6]. The reception of a particular signal from the cytokines can modify the usual behavior of immune responses. In Covid-19, exaggerated pro-inflammatory cytokines initiate different inflammatory signaling pathways via their receptors on immune cells which leads to complicated medical symptoms [7]. This can create overwhelmed ICU (intensive care unit) resources due to the critically-ill patients in need of constant care. Even vaccination may face a problem due to existing comorbidities or pharmacological treatments. The vaccine efficiency may also be effected due to cytokines storm syndrome. There might occur an overlap between the factors that causes further damage to the prognosis of the infection namely cytokines storm and those which may effect vaccine effectiveness [8]. Different treatment approaches and therapies targeting cytokines are currently being used to treat cytokines storm. An overview in [9] is provided for specific treatments for cytokines blockades. Certain cytokines inhibitors such as tocilizumab, baricitinib and anakinra which are sometimes used for the treatment provide with promising results. The use of cytokines inhibitor such as tocilizumab is relatively safe and effective [10]. The control effort limits the fall of healthy cells and tissues due to cytokines storm which may enable the system/body to fight the virus further. The mathematical knowledge used to describe, predict and control biological systems relies on the approaches of non-linear dynamics [11]. Similar to physics, for understanding the complex and nonlinear systems, within host cell mechanisms may too require some current interactive and interdisciplinary research through which mathematical models motivated by data can provide understandings of the cytokines immune response-virus mechanism. A few significant contributions of the researchers are as follows: In [12] a VARI (Vector Autoregressive Integrated) model of covid-19 cases in Indonesia and Singapore is build and development of covid-19 cases is analysed. In [13], STEM (Spatiotemporal Epidemiological Modeler) is used to analyse Covid-19 SEIR model.

A kinetic model for Covid-19 describing the dynamics of the variation SIR cases in [14] based on a lag logistic equation. SIR and SEIR epidemiological models are generalised to situations with anomalous kinetics [15] with the aim to describe a single epidemiological peak. For within host mechanics a full density functional quantum mechanical (DFT) simulations is described of Mpro( SARS-CoV-2 main protease) in complex in [16] to obtain absolute ligand binding energies with various ligands. And in [17] within host dynamics of a SARS-CoV-2 system is discussed but it has not considered the impact of cytokines storm among critically ill patient. A few more literature on the work done for within host-SARS-CoV-2 related systems can be seen in [18-20]. Still within-host modeling and analysis for SARS-CoV-2 virus have not been explored enough compared to the need to understand the inner working of our bodies. Therefore, we will aim to analyse and understand a within-host model of the SARS-CoV-2 virus.

Goal and structure of the study

With the aforesaid information, it is clear that cytokines storm can impact the system leading to systemic hyper-inflammation. Thus, in this paper, our aim is to study a system showing the impact of cytokines storm. A healthy cell-immune responses model is proposed to simulate the effects of the cytokines storm. We will propose an SARS-CoV-2 model consisting of healthy cells, cytokines, NK-cells, T-cells, B-cells and antibodies, with the impact of cytokines storm in Sect. 2. The basic reproduction number, existence of equilibrium points, feasible region and the global stability of disease-free and endemic equilibrium points are calculated in Sect. 3 with their detailed proof in the Appendix. In Sect. 4, we will formulate an optimal control problem using immunomodulatory agent and discuss it using LQR method. In Sect. 5, we will perform numerical simulation for validation of our analytic results, curve fit for real data of cytokines, curve-fit using an algorithm based on initial guess approach for the system and local sensitivity analysis for reproduction number and endemic equilibrium point. We shall also study global sensitivity analysis using PRCC (Global sensitivity analysis method) to determine the sensitive parameters which are important for transmission of infection and healthy cells. Thus, we shall proceed with our study in the following sections to address the following research questions as how the cytokine storm impacts the disease.The novelty of the paper lies in the concept of developing a within-host model in the presence of cytokines storm and its effect on healthy cells. The paper encapsulates dynamical analysis of the system, model fitting to real cytokines data and extending the model to an optimal control problem considering immunomodulatory therapy as control.

What is the impact of cytokines storm on disease dynamics for a SARS-CoV-2 patient?

Is the activation of the ideal immune responses of the cytokines required for suppression of transmission of infection (basic reproduction number)?

Cytokines inhibitors are considered relatively safe and effective to support and build a stronger system (healthy cells). What is the impact of using immunomodulatory therapy on the system during cytokines storm ?

Model formulation

We have formulated a model referring [17] but in the formulation of our model, we have explicitly focused on critically ill patient. H(t) is the concentration of healthy cells which are uninfected, C(t) is the concentration of cytokines, I(t) denote the concentration of T cells which are infected, N(t) are the concentration of natural killer cells, T(t) are the concentration of T-lymphocytes cells, B(t) are the concentration of B-lymphocytes cells, A(t) are the concentration of antibodies and V(t) are the concentration of the virus at time t. The formulation of the equations are done as follows:The parameters are explained in the Table 1 and a schematic diagram of the model has been shown in Fig. 1 for better understanding.

Table 1

Parameter description

Parameter	Dimension	Description
Δ	Cells ml^-1 day^-1	Rate of production for healthy cells
μ_H	day^-1	Per capita death rate of healthy cells
Υ	ml RNA copies^-1 day^-1	Rate at which the healthy are infected
p_1	ml cells^-1 day^-1	Rate at which the healthy cells are reduced due to N
p_2	ml cells^-1 day^-1	Rate at which the healthy cells are reduced due to T
s_1	day^-1	Rate at which virus is produced due to infected cells
s_2	day^-1	Rate at which cytokines are produced
μ_I	day^-1	Per capita death rate of infected cells
β	RNA copies^-1	Strength of immunosupresion
p_3	ml cells^-1 day^-1	Rate at which the infected cells are reduced due to T
μ_V	day^-1	Per capita death rate of the virus
p_4	ml cells^-1 day^-1	Rate at which the infected cells are reduced due to N
p_5	ml cells^-1 day^-1	Rate at which the virus are reduced due to C
p_6	ml cells^-1 day^-1	Rate at which the virus are reduced due to A
r	day^-1	Rate of activation of natural killer cells
k	day^-1	Production rate of antibodies
μ_N	day^-1	Per capita death rate of Natural killer cells
μ_C	day^-1	Per capita death rate of cytokines
μ_T	day^-1	Per capita death rate of T-cells
μ_B	day^-1	Per capita death rate of B-cells
μ_A	day^-1	Per capita death rate of antibodies
η_1	cells^-1 day^-1	Rate of activation of T cells
η_2	cells^-1 day^-1	Rate of activation of B cells
p_7	ml cells^-1 day^-1	Rate at which the antibodies are reduced due to V

Fig. 1

Schematic diagram depicting the interaction among all the compartments

Healthy cell H(t) gets infected due to virus at the rate . Along with it, due to cytokines storm, we assume that the natural killers and T-cells also negatively effect H-cells at rate and respectively [4, 6] which was not considered by [17], is the natural death rate. Hence, the evolution equation is:

The infected cell increases due to interaction of health cells with the virus. As an immune response N and T proceed to kill these infected cells at the rate and and is the natural death rate of I, giving rise to the equation:

New virus particles are produced by I at a rate . Also, the antibodies and cytokines work to reduce the viral load at a rate and which gives the evolution equation:

Cytokines are activated due to cellular infection and viral replication at the rate along with immuno-suppression at the rate and natural death rate making the equation of the following form:

These cytokines subsequently activate natural killers at the rate r and these natural killers help in reducing the infected cells. These natural killers also naturally decay at the rate . Hence, the equation is as below:

Cytokines activate the adaptive immune cells T-lymphocytes at the rate with natural decay at a rate making the equation look like:

T-cells further activate B-lymphocytes at the rate with as the natural death rate of B-lymphocytes. Therefore, the evolution equation becomes:

Finally, these B-cells produce antibodies at the rate k. Critically/severely ill patients of COVID-19 develop SARS-Cov-2 specific antibodies [21, 22]. Thus, at the time of cytokines storm for a critically ill patient, antibodies exist without delay and these antibodies are neutralized by the virus at the rate .

Dynamical analysis

Basic reproduction number and existence of equilibrium points

The basic reproduction number is given bywith the proof provided in Appendix A.

The system has following four equilibrium points:

Disease free equilibrium point :

where

Equilibrium point without immune response : where and which exists for .

Equilibrium point without adaptive immune response:

and assuming , we have provided , , , and for we get a cubic equation : where , , and .

Now and for : . if where i.e. provided . Then by Descartes’ rule of sign there exist one positive root of .

Equilibrium point with adaptive immune response:

and assuming , we have ,

. Now , if where i.e. provided . And if and if .

Theorem 1

The feasible region is and is positively invariant.

The detailed proof for Theorem 1 is provided in Appendix A.

Global stability for disease-free equilibrium and endemic equilibrium point

Theorem 2

If i.e if then the DFE will be globally asymptotically stable(GAS)in .

is globally asymptotically stable for .

We shall use geometric approach method [23] to establish global stability for endemic equilibrium and the proof is mentioned in the Appendix A.

The next section would emphasis especially on control of cytokines through immunomodulatory therapy in specific tocilizumab treatment which would protect healthy cells. The positive impact of tocilizumab treatment which supports our study is also shown in [7], where COVID-19 patients in the ICU across hospitals in the United States showed that the risk of mortality was lower in patients who received tocilizumab treatment (by inhibiting cytokines signaling) immediately after ICU admission than in those who did not receive early tocilizumab intervention. As healthy cells are the part and parcel of our body it will give vitality for the system.

Optimisation model

Immunomodulatory therapy as control by linear feedback

To design a feedback law that may stabilize the equilibrium point or to increase the healthy cell count to a certain desired level, we shall use optimal control theory [24] and we have the dynamic system form with control as follows:We can use an immunomodulatory agent [25] as a control strategy such as tocilizumab with the aim to suppress the cytokines storm i.e to control the response of the storm. We assume that the immunomodulatory agents can directly target cytokines. With the goal of cytokines storm control strategy, we aim to increase the level of healthy cell population and maintain at level ( i.e increased desired level of H to be less than the level at DFE ) by . is the designed healthy cell population. The equations satisfied by the desired positive steady state with control are given below as:From the above equations, we obtain the control variable ,When controls the desired steady state of the system i.e (), it may be unstable. Therefore, u (the feedback control) can be made such that the desired state will be asymptotically stable. Defining new variables as follows:where . We shall now substitute Eq. (13) into Eqs. (10) and admitting (11) we have:Hence we are able to obtain the error system as below:matrices A and B are given asTo find the feedback control u we shall use the Theorem 1 of [26] for the following:

Theorem 3

[24, 26]. For any matrix P andfunction l(y) is positive definite at the neighbourhood of the origin (0, 0, 0, 0, 0, 0, 0, 0).

The detailed proof for above Theorem 3 is provided in Appendix A.

Numerical simulations

In this section, we shall now analyse the system based on certain parametric values mentioned in Table 2 and run simulations using MATLAB. A few parametric values have been taken from [17, 27–31] and the fitted values are obtained from cytokines load real data of RT-PCR positive patient COVID-19 patients [32]. The observations in [32] were indicative of a strong relation between the cytokines level and severity of symptoms/organ damage. We used the real cytokines levels data from [32] which were observed in patients(admitted to the Mount Sinai Health System, New York) in relation to SARS-CoV-2 PCR status(positive). The remaining parameters are further adjusted so as the model describes an uncomplicated SARS-CoV-2 infection system

Table 2

Parameter values

Parameter	Value	Source	Parameter	Value	Source
Δ	8000	[31]	μ_H	0.14	[17]
Υ	0.9	Assumed	p_1	0.0000000001	Assumed
p_2	0.0000000001	Assumed	s_1	0.002	Assumed
s_2	0.0005	Fitted	μ_I	0.003	[17]
β	0.5	Fitted	p_3	0.001	[27]
μ_V	0.002	[17]	p_4	0.000574	[17]
p_5	(0–1)	[17]	p_6	(0–1)	[17]
r	0.52	[17]	k	0.2	[28]
μ_N	0.07	[30]	μ_C	0.7	Fitted
μ_T	1	[17]	μ_B	0.2	[29]
μ_A	0.07	[17]	η_1	0.9	Assumed
η_2	0.01	[17]	p_7	0.0000003	[17]

Equilibrium points

Disease free equilibrium point: For and rest parameters as in Table 2 we get as seen in Fig. 2. The global stability condition for i.e if : is met as per the analytical solution).

Fig. 2

Dynamics of equilibrium points

Equilibrium point without immune response: For and rest parameters as in Table 2 we get . For the existence as seen before in Sect. 3.1, i.e is validated numerically.

Equilibrium point without adaptive immune response: For and rest parameters as in Table 2 we get . if i.e which is in line with the existence of equilibrium.

Equilibrium point with adaptive immune response: For and rest parameters as in Table 2 we get if i.e is validated numerically for the existence of equilibrium.

Model validation: curve-fit of cytokines level C(t) to real data

Cytokines load data was obtained from [32]. For our analysis we shall be using a data set of patients (1422 independent samples of patients admitted to the Mount Sinai Health System in New York) with respiratory symptoms for their IL-6 (cytokines) level whose RT-PCR came positive. The data was collected with the help of online software [33] and fitted using Matlab. The fitting is shown in Fig. 3 (The value for both the plots are in real scenario). It can be seen that the cytokines level for our model solution is close to the actual numbers from real data, Thus, proving validity of our model.

Fig. 3

Fitting model solution for cytokines to RT-PCR positive patient with respiratory symptoms real data

Curve-fit for system using optimization algorithm

We have used python script to fit our model based on a system of non-linear differential equations for curve-fitting. Following our aim to curve fit our model based on the parameters(endemic equilibrium) we have referred or assumed, we have used standard algorithm curve_fit. We shall do the fitting by numerical approximation [34]. As parameter estimation is a special optimization problem, it would use an algorithm to minimise the squared error as objective and it’s performance is based on the initial guess approach. For an initial value problem, an algorithm can be used to approximate the solution for a given set of initial conditions and parameter values. This set of values are taken as initial guess by the algorithm. The algorithm used is the Levenberg–Marquardt algorithm which chooses the parameters and initial values to minimise error or noise [35]. However in Fig. 4a we do not get an adequate desired fit. Therefore, we use an optimizer(based of Nelder-Mean algorithm [36]) and random search algorithm which generates random initial guesses and chooses the best approximation. We see that this algorithm is quite efficient for finding an optimal fit. In Fig. 4b we can see that this approach helped generate a good fit for our model. The regression line(for least square estimation [37]) shows the predicted response based of the random search and algorithm. Thus, the estimated parameters by the algorithm are like our input parameters.

Fig. 4

Curve fit using a Levenberg–Marquardt algorithm and b Nelder-Mean algorithm

Sensitivity analysis for

In this section, we shall discuss the sensitivity analysis of with respect to parameters. We shall follow the method as mentioned in [38] to find the sensitivity index. We shall start with the sensitivity analysis for with respect to its parameters k:Using above, the normalised forward sensitivity indices with respect to the parameters we get:We see that are all positive and are negative. This implies that for say example , we decrease/increase by certain then will decrease/increase by same percentage. And for we decrease/increase by certain then will increase/decrease by same percentage and here is a highly sensitive parameter. And we can see for our system that if decreases then infection increases, and if increases then the healthy cells get infected more. If we increase the parameters by then increases/decreases by . In Fig. 5 For disease free point we can see that increase in and by will increase by . And increase in and by will decrease by . ( and will behave in similar manners and thus have not been added to the figure.) Similarly for endemic equilibrium point in Fig. 6 we see the same trend for the change in behavior of . Thus, we can take note of the important parameters one needs to keep in check in order to control the transmission of infection. These will be of help to maintain either for disease free system (restrict from increasing beyond 1) or endemic equilibrium system (suppress to bring down below 1).

Fig. 5

Behaviour of due to change in parameters for disease-free equilibrium point: the graphs report comparison between original and increase/decrease in parametric value

Fig. 6

Behaviour of due to change in parameters for endemic equilibrium point: the graphs report comparison between original and increase/decrease in parametric value

Uncertainty analysis of

PRCC (partial rank correlation coefficient) [39] is one technique which help’s us quantify the uncertainty for any model. For our basic reproduction number, we shall use PRCC to identify and quantify how the parameters (input)’s uncertainty may impact the transmission of infection. We have considered the output as for sample size . We consider the six parameters from (9) and have chosen normal distribution for the parameters. We find the PRCC values using Matlab using , , , , , .

For some selected combinations of sensitive parameters we have also found surface diagrams such as with and in Fig. 7.

Fig. 7

Phase plane of with the parameters

We get the PRCC values for our input parameters which can be seen in Fig. 8. We obtain following indexes for the parameters: and . The graphs shows that is positively correlated to and with maximum impact of which means that as the rate at which virus is produced due to infected cells are increased, the transmission of infection will increase as well. The effect of the parameter will bring about an opposite change in the transmission of infection as it is negatively correlated. As the death rate of virus increases, the infection transmission will decrease. Further, and satisfies the same negative correlation. Since the value parameter is close to 1, it indicates a strong correlation to change in . Thus, this hints at the need in suppressing the virus to decrease further transmission of infection through . Therefore, this emphasises the activation of the immune response of the body to fight the virus.

Fig. 8

PRCC: sensitivity indexes of

Uncertainty analysis of rate of change of H

In the similar manner, we shall use PRCC to identify and quantify how the parameters uncertainty may impact the rate of change of healthy cells in a SARS-COV-2 infected patient. For sample size we consider the five parameters and have chosen normal distribution for them. We find the PRCC values using Matlab using , , , , .

We get the PRCC values for our input parameters which can be seen in Fig. 9. We obtain following indexes for the parameters: and . The PRCC indices satisfies the same set of correlation as per the assumptions of formulation of rate of change of H i.e is positively correlated and the rest parameters are all negatively correlated. The value and parameters are close to , it indicates a strong correlation to the rate of change in H. We see that and are the parameters related to rate at which the healthy cells are reduced due to N and T respectively. This shows that this particular immune response which is due to the cytokines storm has an strong impact on the healthy cells. Thus, the exaggerated response of cytokines has an strongly negative impact on the healthy cells.

Fig. 9

PRCC: sensitivity indexes of rate of change of H

System with control

Now we shall numerically go for stabilising the system at desired steady state (10) with to avoid elimination of healthy cells and tissues due to cytokines storm. We aim to increase the healthy cell to get a system where the body is still able to fight the foreign entities at time of hyper-inflammatory immune response in critically ill patients. We calculated and from (11). From (12) we get . Matrix obtained for the system is as follows:AssumingUsing LQR commands in Matlab, from the Riccati equation we obtainWe shall investigate the local definiteness of l(y) (18) using graphic representation. Figure 10 shows the local definiteness of the function l(y). We have shown the local definiteness interaction between select variables as the rest variables have similar graphic representation as well.

Fig. 10

3D sections of l(y) at the origin

Finally, we get the form for optimal strategy as follows:For the system with control we get as in Fig. 11a. We can see that due to the control the population of healthy cells increases(Fig. 11b) as compared to the H population in Fig. 2d. Thus, the use of cytokines storm control strategy can help to steer clear of exhaustion of healthy cells and tissues with the aim to resume the host’s homeostasis. Thus, this can be an approach for severe COVID-19 patients to counteract cytokines storm.

Fig. 11

.

.

Conclusion

Though the vaccine drive has started in many countries, the transmission of SARS-CoV-2 infection is still continuing. Cytokines storm plays an important role as a cause for ARDS and MOF, thus leading to disease aggravation. Cytokines storms have been detected in Covid-19 critically ill patients, as seen in some clinical studies. Due to an unbalanced immune response, a strong cytokine storm may be extremely damaging to the patients, who may require intensive-care support. Thus, to save the lives of patients, one needs efficient strategies to suppress the cytokines storm.

The formulation of our complex model is with the aim to understand the immunopathological mechanisms of cytokines storm due to SARS-CoV-2 virus. Model was fitted to cytokines real data of SARS-CoV-2 PCR status (positive) patients with respiratory symptoms. To address our research questions, we studied the dynamic analysis of our system.

The observations from our analysis and numerical simulations brought the following suggestive measures which could be taken:

Since, the cytokines storm based immune response incorporation in the system has brought upon undesired effect as it reduces the healthy cells and tissues the body leading to ARDS or MOF which makes the patient’s recovery becomes tougher and tougher, our analysis give support to the study of development and use of treatments and vaccines for large populations in view of cytokines storm. By focusing on the impact of the uncontrolled immune response, decisions can be taken to tackle the critically ill patients specially in overwhelmed ICU’s more effectively.

The uncertainty analysis of suggests the need for the activation of immune response to suppress the virus which would further result in decrease in transmission of infection rate .

The global sensitivity analysis projects that and parameters with their PRCC indices as and respectively are strongly negatively correlated to the rate of change in H. These parameters which are related to rate at which the healthy cells are reduced due to N and T respectively showed that due to the cytokines storm this particular immune response has a strong impact on the healthy cells and may reduce the healthy cells.

The population of healthy cells increases to 953.4 as compared to the 534.1 in with the incorporation of cytokines storm control strategy. Immunomodulatory agent(such as tocilizumab) was considered as control strategy to control the response of the storm. To design a feedback law that may increase the healthy cell count to a certain desired level, we used optimal control theory. Therefore, working towards strategies that aim at the importance of functional features of cytokines storm-specific response for both vaccinated and non-vaccinated population may be the need of the hour. Further, to tackle the problematic exaggerated immune response due to cytokines storm, we may need to focus on the treatments leading to enhancement of the healthy cells by cytokines inhibition which may further require the understanding of the single interaction parameter with the healthy cells H. Cytokines storm control strategy was discussed with the aim to avoid healthy cells from depleting. This can help to ward off from exhaustion of healthy cells and tissues for severe COVID-19 patients due to cytokines storm. Healthy cells and tissues will be able to provide a system of great vitality for the immune system to work on to suppress the infection.