Dynamical demeanour of SARS-CoV-2 virus undergoing immune response mechanism in COVID-19 pandemic.

COVID-19 is caused by the increase of SARS-CoV-2 viral load in the respiratory system. Epithelial cells in the human lower respiratory tract are the major target area of the SARS-CoV-2 viruses. To fight against the SARS-CoV-2 viral infection, innate and thereafter adaptive immune responses be activated which are stimulated by the infected epithelial cells. Strong immune response against the COVID-19 infection can lead to longer recovery time and less severe secondary complications. We proposed a target cell-limited mathematical model by considering a saturation term for SARS-CoV-2-infected epithelial cells loss reliant on infected cells level. The analytical findings reveal the conditions for which the system undergoes transcritical bifurcation and alternation of stability for the system around the steady states happens. Due to some external factors, while the viral reproduction rate exceeds its certain critical value, backward bifurcation and reinfection may take place and to inhibit these complicated epidemic states, host immune response, or immunopathology would play the essential role. Numerical simulation has been performed in support of the analytical findings.

Introduction

The interhuman transmission of Severe Acute Respiratory Syndrome Coronavirus 1 (SARS-CoV-1) caused drastic concussion on public health over the past decades. The recent coronavirus disease (COVID-19) is imposing extreme destruction in human civilization as about 210 countries around the world have been facing the COVID-19 pandemic cases [1]. The pathological agent of COVID-19 pandemic is Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2), a single-stranded RNA virus, belonging to the Coronaviridae family [2-4]. According to the worldometer data [5, 6], around 24 million people have already been infected with SARS-CoV-2 (and the number is still rising) among whom almost 7% of the COVID-19 patients have been died, and about 1% of the active cases are in critical condition [4]. The World Health Organization (WHO) posted several guidelines about non-pharmaceutical control techniques such as wearing masks, washing hands frequently, social distancing, and others [7, 8]. The countries and corresponding states have been enforced Standard Operating Procedure (SOP) in the public places and work spaces to maintain the COVID-19 guidelines recommended by WHO. Although at the recent time, the total death cases due to COVID-19 pandemic has been reduced to 3.4%, the global fatality rate is converting very rapidly. In this scenario, with the aim to control and overcome this COVID-19 pandemic, it is necessary to conduct more comprehensive research studies at the cellular level for better understandings.

To make accurate entry in host cells, SARS-CoV-2 utilizes angiotensin-converting enzyme 2 (ACE2) receptor which assists to turn the host cells more vulnerable [9]. Even though the expression of ACE2 is abundant in many epithelial cells like myocardial epithelial cells, kidney tubular epithelial cells, and gastrointestinal epithelial cells, type II alveolar epithelial cells of lungs contain the better copious expression of ACE2 and thus are considered as the major target cells of COVID-19 infection [10, 11]. The well-known traits of COVID-19 infection are common fever, cough, fatigue, breathlessness, loss of odor, and tastelessness. The individuals suffering with co-morbidity factors like cardiovascular issues, diabetes, liver diseases, and renal disorders are more endangered group of COVID-19 infection. While any pneumonia patient be affected by COVID-19 infection and, at the same time, that individual would face acute respiratory illness, multiple organ failure would take place at that case. From the previous waves of the COVID-19 pandemic, it has been detected that an outsized percentage of the COVID-19-tolerant experienced mild symptoms and became cured in consort with their own immunity; about 20% of them faced multiple organ failure and ultimately death [12, 13]. Lymphopenia, abnormal degradation in lymphocytes count (bellow the range ), has also been reported in severe COVID-19 patients [13-15].

After successful exposure to SARS-CoV-2, cell-mediated adaptive immunity (if innate immunity is not sufficient) be triggered through the activation and differentiation of T cells. The differentiated T cells benefit the cytokines, viz., INF-, IL-6, and IL-10, to be secreted [16] which are responsible to heavily stimulate the host immunity. Among different types of T cells, CD4T cells and CD8T cells have the most influence to fight against SARS-CoV-2 where CD4T cells assist to generate virus specific antibodies with the activation of T-dependent B cells and the role of CD8T cells is to neutralize the infected epithelial cells. It has been found that in a SARS-CoV-2 host, a large number of CD8T cells present where near about 80% of these cells are infiltrator and inflammatory in nature perceived in the interstitial pulmonary tract performing a crucial role in removing SARS-CoV-2 from the body. The depletion of CD4T cells is correlated with reduced conscription of lymphocytes and neutralizing the assembly of antibodies and cytokines, resulting in severe immune-mediated interstitial pneumonitis and delayed SARS-CoV-2 clearance [9, 17]. Researchers studied that there is a long-lasting and protracted response of T cells to the S and other structural proteins (including the proteins M and N), which provide sufficient knowledge to draft the SARS vaccine by combining viral structural proteins. These types of vaccine may provide a strong, efficient, and long-term response to the virus by memory T cells [18]. Moreover, the clinical trials examined that monoclonal antibody therapy is an effective intervention strategy with better antiviral response to SARS-CoV-2 [19].

To design effective treatment strategies that can target both SARS-CoV-2 and the infected epithelial cells, we need to understand the intermediate relationships among uninfected epithelial cells, SARS-CoV-2 virus particles, and host immune response. The infection processes of SARS-CoV-1, SARS-CoV-2, and MARS-CoV are almost identical. The functional interaction analysis between the host immune response and other HCoVs, combined with evolutionary sequence analysis of SARS-CoV-2, may guide to come up with new treatment and preventive interventions. Mathematical modeling may be an efficient tool that will aid us to recognize the within-host interactions in COVID-19 infection. Nadim et al. [20] discussed the short-term dynamics and preventive approaches for COVID-19 infection. In this study, they have shown the declining inclination of new COVID-19 cases through establishing effective management of quarantined persons. Volpert et al. [21] have attempted to measure the efficiency of quarantine strategies by means of a new mathematical modeling of COVID-19 infection progression. In this study, it is observed that the peak of infection (maximum of daily cases) is accomplished about 10 days after the commencement of confined process and the authors suggested to implement stricter measures as radically obligatory methods. Mondal et al. [14] presented a seven compartmental SEIQR type model to explore the COVID-19 disease progression and the effect of pharmaceutical and non-pharmaceutical interventions as control input and there effects in reducing the number of the infected population.

Even though there are a few mathematical modeling highlighting the transmission dynamics of COVID-19 infection at the population level [22-26], yet the intrahost viral dynamics of COVID-19 infection has not been investigated on a large scale. Mathematical modeling with real data assists to extensively explore dynamical aspects of any infection under cellular level [27-29]. Tang et al. [30] proposed a four-dimensional model depicting the basic virus and host immune response dynamics incorporating the concentration of DPP4 receptors for MARS-CoV infection. Chatterjee and Basir [31] explored a mathematical model indicating the dynamical behaviors of epithelial cells during SARS-CoV-2 infection in the presence of CTL response and they studied the function of the ACE2 receptor concluding the fact that the proper dose of immunostimulant drug will aid to reduce COVID-19 infection. Hernandez et al. [32] proposed a model to examine the cellular level dynamics and T-cell responses against the viral replication of SARS-CoV-2 during the COVID-19 infection. Wang et al. [33] evaluated the effect of the several pathogenic characteristics of SARS-CoV-2 in viral dynamics and host immune response, and from this study, it is observed that anti-inflammatory treatment strategies or combination of antiviral drugs with interferon is effective in reducing the plateau phase in viral load and shortening the recovery time. Chatterjee and Basir [13] formulated a mathematical model to examine the consequences of adaptive immune response to the viral mutation in controlling disease transmission and the effect of the combined antiviral drug therapy on the model dynamics. Chatterjee et al. [34, 35] proposed a set of fractional differential equations model in cellular level accounting the lytic and non-lytic effects of immune response in the kinetics of the model exploring the effect of a commonly used antiviral drug in COVID-19 treatment along applying an optimal control-theoretic approach.

In this mathematical study, through our proposed model, the intermediate ambiences of the intrahost immune response on SARS-CoV-2 viral dynamics during COVID-19 pandemic has been calibrated. The overall study is designed in the following manner: in Sect. 2, the structure of our proposed compartmental model is presented. Section 3 is build up with the qualitative characteristics of the model such as non-negativity and boundedness of all the solutions of the system. In Sect. 4, an epidemic subsystem of our proposed model is considered. In the consecutive subsections, all the possible equilibrium points executed by the subsystem and the basic reproduction number of this subsystem are computed. The stability of the system is analyzed, and existence of backward bifurcation is studied. Section 5 is accomplished with numerical simulation to examine whether the analytical results are aligned with the numerical findings or not. Finally, we discuss about our overall analytical and numerical findings in Sect. 6.

The compartmental model

The mathematical modeling in epidemiology advantages us to make out the basic dynamics of any communicable disease and, in general cases, target cell-limited mathematical models are formed to study the viral dynamics. In this research work, to observe the dynamics of SARS-CoV-2 and the role of host immune response in COVID-19 infection, we are aimed to propose a target cell-limited model consisting of five populations, viz.

: the uninfected and susceptible target cells of COVID-19 infection; these are the surface epithelial cells with ACE-2 receptor located at the respiratory tracts including lungs, nasal, and trachea/bronchial tissues,

: the SARS-CoV-2-infected virus-producing epithelial cells,

V: the SARS-CoV-2 virus particles,

A: the antibody response, and

C: the CTL response.

Our proposed five-dimensional and ODE compartmental mathematical model is as follows:with non-negative initial conditionThe first equation of the epidemic system (1) expresses the dynamics of uninfected epithelial cells () and the second equation shows the dynamics of the infected epithelial cells . The replication of SARS-CoV-2 virus (V(t)) is considered in the third equation of the system (1) where the COVID-19 infection promotes endothelins on several organs as a direct consequence of viral involvement. The fourth and fifth equations of the model (1) represent the kinetics of the antibody response and the CTL response, respectively, to combat against the COVID-19 infection.

The constant recruitment rate of the uninfected epithelial cells is denoted as  (). During the host cell entry of SARS-CoV-2, the coronaviruses first bind to a cell surface receptor for viral attachment [36, 37] and let (), be the rate at which SARS-CoV-2 infects the susceptible epithelial cells. After the commencement of infection virions are secreted from the infected epithelial cell [36, 37] and let p () be the rate at which new virus are produced. Virus particles are removed at a rate (). The uninfected epithelial cells expire at a rate () due to natural apoptosis and the infected cells are removed from the system at a rate () as a result of cytopathic viral effects and immune response.

Cytokines (IFNs) are too essential to inhibit the viral replication and in modulating the downstream effects of the host immune response. Specific cytokines activate natural killer cells (NK) to fight against the infected cells. It has been observed that SARS-CoV-2 often targets the JaK/STAT pathway (a pathway which is activated by the IFN signaling cascades) to reduce the production of IFNs. We consider this immune suppressing mechanism of SARS-CoV-2 as , a functional form of the declinature in the production of cytokines. The adaptive immunity of the host (specifically the T cells and B cells) be activated through the cytokines to produce antigen-specific antibody response. B cells mainly secrete IgM and IgG antibodies which are released from blood and lymph fluid benefit to neutralize the viral particles. We extend the target cell-limited version of our proposed COVID-19 infection system by incorporating the antibody response A(t) with depletion rate (). We consider the expansion of CTL response mediated by the infected epithelial cells in a functional form as and at the rate (), the infected epithelial cells would be neutralized by CTL response. Here, the constants and stand for the simulation rate of antibody response and CTL response, respectively. The constants and describe the half-maximal simulation thresholds. The antibody response loses its productivity at the rate () and the CTL response loses its productivity at the rate (). The system will undergo the direct viral cytopathicity above normal target death rate, while , for , and for , cytopathic effect of SARS-CoV-2 would be absent in the system. All the system parameters are non-negative. The host–pathogen interaction between epithelial cells and SARS-CoV-2 and the effect of host immune response (both antibody and CTL responses) in this interaction process are explained through Fig. 1.

Fig. 1

Graphical explanation regarding the effects of antibody response and CTL response in COVID-19 infection for the intrahost mathematical model (1)

Qualitative characteristics of the model

In this section, we study two qualitative characteristics—positivity and boundedness of the state variables of our proposed model system (1). These properties are essential for the epidemic system (1) to be biologically convincing.

Positivity

Theorem 1

For all the solution trajectories of the system (1) along with the initial condition (2), the state variables remain positive, i.e., the conditions , , , , and hold for all time window .

Proof

The fourth and fifth equations of the epidemic system (1) may be reformed asOn integration of the above inequality (3), we achieveThe third equation of the system Eq. (1) may be written asand integrating the Eq. (4), we see thatNext, the second equation of the epidemic system (1) could be rewritten asThus, integrating Eq. (5), we haveIn the similar fashion, from the first equation of the epidemic system (1), we may obtainOn integration of Eq. (6), we getTherefore, all the five state variables of the epidemic system (1) are positive, i.e., all the solution trajectories of the system (1) together with the non-negative initial condition (2) remain positive for all time window .

Boundedness

Theorem 2

All the positive solution trajectories in of the system (1) together with the non-negative initial condition (2) are uniformly bounded in the positively invariant and attracting region (defined in the proof).

Adding the first two equations of the epidemic system (1), we may writewhere we assume that . From the above expression, it is implied thatwhere . Now, taking on both sides of the above equation, we getFrom the third equation of the system Eq. (1), it is followed that:and thus, we getThus, taking on both sides, we may obtainIn a similar way, from the fourth equation of the system (1), we may writeand accordingly, we haveAs a result, taking on the both sides of the above equation, we obtainFinally, from the last equation, it may be seen thatand thus, the above expression is leading toTherefore, assembling the above computations, it is notable that all the positive solutions initiating in of the system (1) together with the non-negative initial condition (2) are uniformly bounded in the following positively invariant and attracting region:Hence, the system (1) is well posed and epidemically realistic.

Consideration of subsystem

The data on the CTL response and antibody response are not available properly due to a large number of parameters that cannot be currently estimated. Therefore, we simplify the model (1) using a quasi-steady state approximation, assuming that the dynamics of antibody and CTL stimulation is faster than the time course of acute SARS-CoV-2 infection [32–34, 38]. This derives a model similar to the target cell-limited model, given belowwith the above-mentioned initial condition (2) of the model system (1). Here, we considered the terms and as () and () defining the destruction rate of infected epithelial cells through productive CTL response and antibody response, respectively.

Equilibria of the subsystem

The epidemic subsystem (11) comprises two equilibrium points: (i) infection-free equilibrium point: , where there is no COVID-19 infection in the system (11) (exists whatever the circumstances) and (ii) endemic equilibrium point: for enduring COVID-19 infection in the system.

Basic reproduction number

Basic reproduction number of an epidemic system is the expected number of secondary infections directly propagated from the case where it is assumed that primary classes are all susceptible to infection. We focus on to determine the basic reproduction number of the epidemic system (11), using the next-generation matrix method [39] at the infection-free equilibrium point . The epidemic system (11) is containing two infected classes and V. Following the approach of [14], the two matrices required in computation of basic reproduction number are and at , and they are defined as follows:Thus, we get the next-generation matrix of the system (11) asLet us consider be the eigenvalues of the matrix and the spectral radius of the matrix is defined as . The basic reproduction number () of the system (11) is the spectral radius of the next-generation matrix and is obtained asBasic reproduction number () is the crucial fact in analyzing the kinetics of the epidemic system (11).

Existence conditions of endemic equilibrium

In this subsection, we compute the components of the endemic equilibrium point and at , solving the system Eq. (11), we getprovided the fact that . The value of could be obtained from the following quartic equation:whereFrom the above quartic Eq. (12), it is noticeable that the coefficient is positive irrespective of any condition and the coefficient would be positive whether and . Again, we observe that the coefficient would be positive if and . Now, it is noted that the coefficient be positive with the condition that along with and the coefficient would be positive (or negative) for (or ). Therefore, summarizing these conditions, we can conclude that the existence and positivity of the endemic equilibrium point of the epidemic system (11) be possible if (i) , (ii) , and (iii) depend on if less than unity or greater than unity.

Next, the possible number of positive real roots of the quartic Eq. (12) depending upon is arranged in the tabular form (in Table 1).

Table 1

Possible number of positive real roots of Eq. (12)

Case	σ_1	σ_2	σ_3	σ_4	σ_5	R_0	Changes	No. of possible
							of sign	positive real roots
1	+	+	+	+	+	R_0<1	0	0
	+	+	+	+	-	R_0>1	1	1
2	+	+	+	-	+	R_0<1	2	0,2
	+	+	+	-	-	R_0>1	1	1
3	+	+	-	-	+	R_0<1	2	0,2
	+	+	-	-	-	R_0>1	1	1
4	+	+	-	+	+	R_0<1	2	0,2
	+	+	-	+	-	R_0>1	3	1,3
5	+	-	-	-	+	R_0<1	2	0,2
	+	-	-	-	-	R_0>1	1	1
6	+	-	-	+	+	R_0<1	2	0,2
	+	-	-	+	-	R_0>1	3	1,3
7	+	-	+	-	+	R_0<1	4	0,2,4
	+	-	+	-	-	R_0>1	3	1,3
8	+	-	+	+	+	R_0<1	2	0,2
	+	-	+	+	-	R_0>1	3	1,3

Assembling the possible cases enlisted in Table 1, we may construct the following lemma:

Lemma 1

In consort with the conditions that , and , the epidemic system (11)

Executes unique positive endemic equilibrium point if and the case 1, case 2, case 3, or case 5 hold; occurrence of forward bifurcation is possible in these cases;

May possess one or more than one endemic equilibria if and the case 4, case 6, case 7, or case 8 hold; possible chance of hysteresis is associated with these cases;

May possess two endemic equilibria if the case 2, case 3, case 4, case 5, case 6, case 7, or case 8 hold; possible occurrence of backward bifurcation is possible in these cases;

May possess four endemic equilibria if and the case 7 holds; possible emergence of backward bifurcation happens in this case, and finally, (v). may possess no endemic equilibrium if and the case 1 holds; this case is linked to forward bifurcation occurrence.

Stability of the epidemic system

In this subsection, we are aimed to study the stability of the epidemic system (11) around both the infection-free and endemic equilibrium points, respectively.

Theorem 3

The epidemic system (11) is locally asymptotically stable around the infection-free equilibrium point on condition that .

To study the stability, first, we have to compute the Jacobian matrix of the system (11) at the infection-free equilibrium which is given as follows:From the above Jacobian matrix , it is observed that out of three eigenvalues of , one eigenvalue is and the rest two eigenvalues are calculated from the following quadratic characteristic equation (with respect to the eigenvalue ):To be stable, the rest two eigenvalues must be negative or be contained of negative real parts and that would be possible if the constant coefficient of the above quadratic Eq. (13) be non-negative. It is notable that the characteristic Eq. (13) will possess negative roots or roots having negative real parts only if .

Therefore, the epidemic system is locally asymptotically stable around the infection-free equilibrium point only if .

Theorem 4

The epidemic system (11) is globally asymptotically stable around the infection-free equilibrium point on condition that .

To study the global stability of the system (11) around the infection-free equilibrium , we apply the comparison theorem of [39]. In this respect, now, we rewrite the infected classes and V) of the system (11) with the help of the next-generation matrix method (discussed in the previous subsection) as follows:From the previous subsection, it could be observed that in , and consequentlyThus, the above linearized inequality system is expressing the fact that this system would be stable if the spectral radius of the next-generation matrix is less than unity, i.e., [40]. Next, using the standard comparison theorem [39], we obtainThus, returning back to the system Eq. (11), substituting the values and , it is obtained that as . Therefore, it is followed that:Hence, the system (11) is globally asymptotically stable around the infection-free equilibrium point in the region on condition that .

Theorem 5

The epidemic system (11) is locally asymptotically stable around the endemic equilibrium point on the conditions that

, ,

, and

, provided that .

To analyze the local stability of the epidemic system (11) around the endemic equilibrium , first, we compute the Jacobian matrix of the system (11) at which is calculated as follows:Now, the characteristic equation of the Jacobian matrix (with respect to the eigenvalue ) is given aswhere , , and . Using the well-known Routh–Hurwitz criterion for stability, we find that (i) , , , and (ii) , and also (iii) , while it is provided that , i.e., the Routh–Hurwitz criteria for stability be satisfied for the characteristic Eq. (14). Thus, the system (11) is locally asymptotically stable around the endemic equilibrium point .

Existence of transcritical bifurcation

In this subsection, we study the existence of transcritical bifurcation about the infection-free equilibrium point for the system (11) using Sotomayor’s theorem [41, 42].

Theorem 6

Assuming as bifurcation threshold, the epidemic system (11) experiences transcritical bifurcation around the infection-free equilibrium point .

To investigate the existence of transcritical bifurcation, first, we detect the eigenvalues exhibited by the Jacobian matrix of the system (11) for and evaluated at the infection-free equilibrium point . In this aspect, we observe that at , the Jacobian matrix possesses one simple zero eigenvalue. We arbitrarily choose as bifurcation parameter. Now associated with the zero eigenvalue of , we consider right eigenvector and left eigenvector .

Now, we rearrange the system (11) in the following form:where we consider . Differentiating the above system with respect to , we getNow, at the infection-free equilibrium point It is observed that . Therefore, the first condition of Sotomayor’s theorem is satisfied.

Now, differentiating with respect to X, we computeand at Thus, we obtainand consequently, . Therefore, the second condition of the Sotomayor’s theorem is satisfied.

Againand accordingly, . Therefore, the third condition of Sotomayor’s theorem holds.

Hence, all the three conditions of Sotomayor’s theorem hold, and thus, transcritical bifurcation occurs at the infection-free equilibrium point taking as bifurcation threshold implying the fact that stability of alters from stability to instability if crosses unity.

Analysis of backward bifurcation

In the previous subsection, in case of studying existence of endemic equilibrium of the epidemic system (11), from the Lemma 1, it is observed that if the condition (iii) or the condition (iv) holds, possible occurrence of backward bifurcation emerges. In the phenomena, backward bifurcation for any epidemic system suggests the co-existence of more than one endemic equilibrium even if which implies the persistence of infection in spite of the value of the basic reproduction number less than unity. To analyze the occurrence of the backward bifurcation in the system (11), we focus on the results of Castillo-Chavez and Song [43].

To apply the results of [43], first, we rewrite the system (11) in vector form asconsidering the state variables as , , and . Thus, the system (11) may be written asNow, arbitrarily choosing as bifurcation parameter, at , we observe that . Thus, for , the variational matrix of the system (11) around the infection-free equilibrium point, executes one simple zero eigenvalue and the rest two eigenvalues are negative and . Therefore, we can apply the center manifold theorem [44], and in this regard, let us consider the right eigenvector and left eigenvector corresponding to the zero eigenvalue of are and , respectively.

Now, according to [43], we compute the value of two constants a and b linked with the determination of bifurcation asAccording to [43], the sign of the constants a and b will determine the possibility if backward bifurcation occurs or not for the system (11). It is observed that b is always positive in any circumstances. For the existence of backward bifurcation with bifurcation threshold , the sign of the constant a must be positive and which would be possible if . Therefore, the epidemic system (11) experiences backward bifurcation whenever and . The appearance of backward bifurcation creates difficulties in control of any epidemic [40]. Thus, only the reduction of the basic reproduction number less than unity cannot eliminate an infection from an epidemic system. The productivity of immune response, especially the efficiency of the antibody response (), has to be increased to its maximal level to control the reproduction of SARS-CoV-2 and the transmission of the COVID-19 infection.

Numerical simulation

In this section, we are aimed to numerically study our proposed intrahost COVID-19 mathematical model to gain vivid insights about the host–pathogen (SARS-CoV-2) correspondence and how host immune response influences this interrelationship in COVID-19 infection. Figure 2 represents the best-fitted model for the viral load data of SARS-CoV-2-infected patients collected from [32–34, 45] and fitted each with our proposed model. All parameters of the best fits are listed in Table 2.

Fig. 2

Best fits of model (1) for SARS-CoV-2 to the viral load data continuous line are the simulation based on (1) and the red circles represent the data from [32–34, 45]

Table 2

Definition and default value of the parameters of the epidemic system (11)

Parameters	Biological meaning	Estimated	Sources
		Mean value
Π	Recruitment rate of epithelial cells	5	[33]
μ_1	Expiry rate uninfected epithelial cells	0.2	Estimated
μ_2	Death rate of infected epithelial cells	0.189	[33]
β	Rate of infection	0.0001	[32, 33]
p	Growth rate of virus in cells	70	[32, 33]
μ_3	Virus clearance rate	0.1	[33]
α_1	Rate of antibody response from immune cells	0.4	Estimated
α_2	Rate of CTL response from immune cells	0.4	Estimated
θ_1	Half maximal simulation for antibody response	0.02	Estimated
θ_2	Half maximal simulation for CTL response	0.1	Estimated

The left panel of Fig. 3 describes the numerical simulation of the epidemic system (11) whenever and the right panel of Fig. 3 shows the existence of two endemic equilibria if in spite of the condition . The left panel of Fig. 4 indicates the occurrence of forward bifurcation and the right panel of Fig. 4 describes the occurrence of backward bifurcation of system (11) using the Theorem 6 and the results from backward bifurcation analysis where the parameter values are as listed in Table 2. Here, in the left panel, the red curve indicates the stable infection-free equilibrium and the blue curve indicates the stable endemic equilibrium point; in the right panel, the blue line indicates the existence of two endemic equilibrium and the red curve indicates the infection-free equilibrium with . It is observed that when , no endemic equilibrium point exists and the infection-free equilibrium is stable in this case. However, a hysteresis loop observes when , and unexpectedly, new endemic equilibrium point appears in spite of the existence of infection-free equilibrium. Also, in the region where and , the infection-free state remains stable. When , the large endemic equilibrium exists, but the smaller one vanishes. Simultaneously, the infection-free state loses its stability.

Fig. 3

Left panel: the numerical simulation of the system (1). Right panel: the figure represents the existence of endemic equilibrium point

Fig. 4

Forward (left panel) and backward (right panel) bifurcation diagrams of system (1). Red curves represent infection-free equilibrium point and blue curves represent endemic equilibria

The solutions of the epidemic system (11) for the different values of non-lytic saturation factor are represented by Fig. 5. Here, we have observed that in the case of an increasing value of , the viral load as well as the infection level increase. Thus, the non-lytic effects of both the antibody and CTL responses play a pivotal role to control the infection and viral load. The solutions of the system (11) for the different values of infection rate () are depicted in Fig. 6. These numerical findings show that decreasing the value of viral transmission rate due to non-lytic immune responses protects the epithelial cells in case of SARS-CoV-2 infection.

Fig. 5

The solution trajectories of the epidemic system (11) for different values of with the set of parameters as in Table 2

Fig. 6

The solution trajectories of the epidemic system (11) for different values of with the set of parameters as enlisted in Table 2

To explore the sensitivity of with respect to the associated parameter variables, we used the Latin Hypercube Sampling and partial rank correlation coefficients (PRCCs) methods. Latin Hypercube sampling is a statistical sampling method and PRCCs rank each parameter by the effect it has on the outcome when all other parameters are kept at median values [46]. Figure 7(left panel) illustrates the degree of sensitivity of each parameters on . suggest that increases with the increasing values of the corresponding parameters, whereas mean the decrease of with decreasing values of the parameters. In this figure, it is clearly observed that the parameters , p and have positive effects, whereas , , and have negative effects on the outcome. Variations of against each model parameter are illustrated in Fig. 7(right panel), where all other parameters are kept at their sample values. Without any intervention, remains above 1. If the growth rate, disease transmission rate, and the virions production rate would be increased by a factor of 100, then the average is lowered but still above 1. If the death rates would be decreased 1% of its sample value, then the average is below 1 along with the upper quartile value below 1. This suggests that eradication is likely to take place. Figure 8 shows the local sensitivities of for all model parameters. The comparison among the non-normalization, half-normalization, and full-normalization techniques provides the sensitivity of some model critical parameters. From the comparison, it can be concluded that the half-normalization technique is more appropriate to spot the model critical parameters as compared with other techniques.

Fig. 7

Left panel: degree of sensitivity of each parameter on . Right panel: box-plots illustrating the variation in

Fig. 8

Local sensitivity analysis with non-normalization, half-normalization, and full-normalization techniques of with respect to all parameters in computational simulations using MATLAB

Discussion and conclusion

Addressing the importance of analyzing viral dynamics of SARS-CoV-2 with the aim to mitigate the COVID-19 infection, we formulated a five-dimensional deterministic intrahost mathematical model. In this study, some existing research works have been reviewed which are based on target cell-limited mathematical modeling of invaded SARS-CoV-2 in human respiratory tract, though in most of those studies, the role of cell-mediated immune response has not been considered. In our proposed model, the potency of host immune response against the interaction between SARS-CoV-2 and epithelial cells of our respiratory system is highlighted which is noteworthy to understand the characteristics of the COVID-19 pandemic. We investigated our model both analytically and numerically. The model parameters corresponding with host–pathogen interaction has been estimated from real sources of data using MATLAB software which is another novelty of the study and it may assist the researchers dealing with mathematical modeling to develop extension on the existing models. The sensitivity analysis of the model parameters associated with the basic reproduction number has been performed using three techniques: non-normalization, half-normalization, and full-normalization, and this sensitivity analysis points out the most critical model parameters indulged in COVID-19 infection. The existence criteria for backward bifurcation in the proposed epidemic system have been analyzed and it has been observed that only the reduction of the basic reproduction number below unity is not sufficient to draw declination in COVID-19 infection in the presence of the backward bifurcation phenomena. In this regard, we should pay more attention to the role of host immune response or immunopathology (specifically to the antibody response) in influencing the clinical infection outcome. This research work will advantage to explore the possible immunotherapeutic strategies to fight against the COVID-19 infection.