Saturated lysing efficiency of CD8⁺ cells induced monostable, bistable and oscillatory HIV kinetics

1.   Introduction

The human immunodeficiency viruses (HIV) are two species of Lenivirus that infect human immune cells (e.g., CD4⁺ T cells) and cause acquired immunodeficiency syndrome (AIDS) ^[1,2]. Combined antiretroviral therapy (cART) can transform HIV infection from a seeming death sentence to a manageable chronic disease ^[3]. However, cART can eventually progress to both HIV drug resistance and systemic side effects ^[4].

Human leukocyte antigens (HLA) class Ⅰ presents short HIV-derived peptides, which lead to different antigenic specificities of CD8⁺ T cell responses ^[5,6]. Variations in HLA genes lead to differences in how the infected cells alert the immune system, and ultimately differences in the disease outcomes ^[6,7]. Approximately 1% of infected patients (known as HIV elite controllers) efficiently control HIV infectivity to very low or undetectable levels and exhibit superior HIV-1-specific CD8⁺ cell responses ^[8,9]. An elite controller must carry a protective HLA allele; however, having protective HLA alleles may not always confer the individual's superior HIV infection outcome ^{[10,11,12,13]}.

Numerous mathematical models have been constructed to understand the role of the lysing efficiency of CD8⁺ cells on HIV ^[8,9,14], influenza ^[10,15,16], hepatitis B ^[6,17], and SARS-cov-2 ^[18] dynamics. However, among these models, the lysing efficiency of CD8⁺ cells is assumed to linearly increase with the CD8⁺ cells or the infected cell count, thereby ignoring saturation. Ganusov et al. ^[19] demonstrated that the lysing efficiency using law mass-action provides a better description of the killing than saturation; however, Myers et al. ^[20] suggested that the saturated model using the Hill function is a better candidate. The saturated lysing efficiency provides a smaller sum of squared error and better Akaike information criterion (ki) values compared to the unsaturated lysing efficiency modelled by the law of mass action ^[21]. Including saturation has been previously considered, albeit not in large numbers of studies. First, Althaus et al. ^[22] used the saturated lysing efficiency using the Hill function to indicate that the HIV kinetics converged to different viral loads according to varying efficacies of their CD8⁺ responses. Next, Burg et al. ^[23] illustrated that the HIV kinetics with the saturated lysing efficiency exhibited damped oscillations upon approaching the steady states. Crucially, the saturation effect always led to bistable behaviors, for example, neutrophil-induced bistable bacteria kinetics ^[24], antibody-induced bistable influenza kinetics ^[15], and antibiotic-induced bacteria kinetics ^[25,26]. Additionally, cytokines have been reported to induce bistable HIV kinetics ^[27,28,29]. Moreover, it is necessary to understand how the lysing rate induced by different HLA alleles, the saturation effect, and CD8⁺ cell count interplay to determine the HIV kinetics. Specifically, it is unclear whether CD8⁺ cells can lead to different HIV kinetics, including monostability, bistability, or periodic oscillations; if they can, then the parameter regions that correspond to the different HIV kinetics have not been explored.

In this work, we conduct a further investigation of our mathematical model that describes HIV dynamics with the saturated lysing efficiency of CD8⁺ cell and biologically relevant parameters ^[21]. After a discussion of the basis for our approach, we investigate the well-posedness of the proposed model, the existence of an equilibrium, analyze their local stability, and numerically determine which parameters can lead to monostability, bistability, and periodic oscillations. Interestingly, a low lysing rate and a low saturation effect can lead to periodic oscillations, and a low lysing rate and a high saturation effect can trigger monostability to a high viral titer; this explains how having non-protective HLA alleles always coincides with bad infection control. Moreover, a high lysing rate and a high saturation parameter results in bistability, in which a high viral inoculum size results in the maintenance of viral infectivity, but a small viral inoculum size count can lead to the eradication of viral infectivity; this explains why patients with protective HLA alleles are not always associated with HIV infection outcomes. Then, we present a bifurcation analysis of the model which leads to unexpected findings that CD8⁺ cell lead to oscillating HIV kinetics, monostable HIV kinetics, and bistable HIV kinetics. Finally, we discuss advantages and limitations of our approach and potential future work.

2.   Model development

The HIV dynamics with the semi-saturated lysing efficiency of CD8⁺ cell are captured using a mathematical framework (Reference ^[21], shown in Figure 1). Comparing the short time scale (days) used in this model with the 8–11 years life expectancy without cART and several decades with cART, the effector CD8⁺ cell concentration is assumed to be a constant. We describe the rate of change of susceptible CD4⁺ cell, infected CD4⁺ cells, and the HIV viral titer as follows:

where $t$ is time (days post infection), $T\left(t\right)$ refers to the susceptible CD4⁺ cell concentration, $I\left(t\right)$ refers to the productively infected CD4⁺ cell concentration, $V\left(t\right)$ refers to the free HIV virus concentration, and $E$ refers to the constant CD8⁺ cell concentration. Note that the interaction among the CD4⁺ cell, HIV, and CD8⁺ cells was considered spatially homogeneous, and they were well-mixed. The interpretation and value of all parameters are given in Table 1.

Figure 1.  Flow diagram for the interaction among CD4⁺ cells, HIV virions, and CD8⁺ cells.

DownLoad: Full-Size Img PowerPoint

Table 1.  Summary of HIV replication and CD8+ cell lysing parameter.

Parameter/Meaning	Value	Dimensions	Reference
$\lambda$, Birth rate of CD4+ susceptible cell	$295$	cells μL−1day−1	[30,31]
${\delta }_{T}$, Death rate of CD4+ susceptible cell	$0.18$	Day-1	[31]
$\beta$, HIV infectivity	$3.9\times {10}^{-3}$	μL virions−1day−1	[32]
${\delta }_{I}$, Death rate of active infected cell	$1.0$	Day−1	[33]
$\pi$, Burst size	$5.5\times {10}^{4}$	Virions cell−1day−1	[34]
${\delta }_{V}$, Virion clearance rate	$36$	Day-1	[35]
$\kappa$, CD8+ cell lysing rate	[0, 10]	μL−1cell-1day−1	[21]
$\gamma$, Saturation parameter controlling lysing efficiency with increase of infected cells	[0, 10]	μL−1cell-1	[21]
$\eta$, Saturation parameter controlling lysing efficiency with increase of CD8+ cell count	0	μL−1cell-1	[21]
$E$, CD8+ cell count	[0, 1600]	μL−1cells	[36,37]

---

 | Show Table

DownLoad: CSV

In HIV kinetics with a constant CD8⁺ cell concentration, uninfected CD4⁺ cells are generated at homeostasis with a rate of $\lambda$ and die with a rate of ${\delta }_{T}$. Susceptible CD4⁺ cells become infected with a rate of $\beta$ upon exposure to free HIV virions $V$, and infected CD4⁺ cells release HIV virions at a burst size of $\pi$ (shown in Figure 1). HIV virions naturally degrade with a rate of ${\delta }_{V}$. CD8⁺ cells lyse actively infected cells with the saturated lysing efficiency $\frac{\kappa IE}{1+\gamma I+\eta E}$, in which the lysing efficiency of CD8⁺ cell increase but gradually converge to its maximum value as both the CD8⁺ cell and CD4⁺ cell counts increase. By fitting the virus inhibition assay data, the parameter controlling saturation effect in the lysing efficiency as CD8⁺ cell count $E$ increase $\eta$ is fitted to be zero; then, the lysing efficiency $\frac{\kappa IE}{1+\gamma I+\eta E}$ is degraded to $\frac{\kappa IE}{1+\gamma I}$ (shown in Figure 1, Reference ^[21]). $\kappa$ represents the CD8⁺ cell lysing rate and $\gamma$ controls the saturation effect in the lysing efficiency as the infected cells count $I$ increases (shown in Figure 1).

3.   The well-posedness of the system

Because system (1) proposes a biological phenomenon, solutions of system (1) are supposed to be positive and bounded with the initial conditions. First, it is easy to know that $\frac{dT}{dt}{\left.\mathrm{}\right|}_{T = 0} = \lambda > 0$, $\frac{dI}{dt}{\left.\mathrm{}\right|}_{I = 0} = \beta T\left(0\right)V\left(0\right) > 0$, for all $T\left(0\right)$, $V\left(0\right) > 0$, $\frac{dV}{dt}{\left.\mathrm{}\right|}_{V = 0} = \pi I\left(0\right) > 0$ for all $I\left(0\right) > 0$. Thus, the positivity is obvious. Next, we prove $\phi = T+I+\frac{{\delta }_{I}}{2\pi }V$. In this case, we have $\frac{d\phi }{dt} = \frac{dT}{dt}+\frac{dI}{dt}+\frac{{\delta }_{V}}{2\pi }\frac{dV}{dt} = \lambda -{\delta }_{T}T-\frac{\kappa IE}{1+\gamma I+\eta E}-\frac{{\delta }_{I}}{2}I-\frac{{{\delta }_{I}\delta }_{V}}{2\pi }V < \lambda -C\phi$, where $C = \mathrm{m}\mathrm{i}\mathrm{n}({\delta }_{T}, \frac{{\delta }_{I}}{2}, \frac{{{\delta }_{I}\delta }_{V}}{2\pi })$. Then, we have $0\le \phi \left(t\right)\le \theta$, where $\theta = \frac{\lambda }{C}$. Thus, boundedness is obvious.

4.   Existence of equilibrium

Setting the right-hand-side of system (1) to zero, three equilibria are found: (a) a solution with the maximum susceptible cells amount equaling its carrying capacity and viral titer equaling to zero; (b) a solution with a low susceptible cell count and a low viral titer; and (c) a solution with a low susceptible cell count and a high viral titer. The populations corresponding to such cases are as follows: (a) virus-free equilibrium ${E}_{0} = \left({T}_{0}, {I}_{0}, {V}_{0}\right) = (\frac{\lambda }{{\delta }_{T}}, \mathrm{0, 0})$; (b) the first non-trivial equilibrium ${E}_{1} = \left({T}_{1}, {I}_{1}, {V}_{1}\right) = \left({T}_{1}^{\mathrm{*}}{\delta }_{V}, \frac{1-{\delta }_{T}{\delta }_{V}{T}_{1}^{\mathrm{*}}}{\beta \pi {T}_{1}^{\mathrm{*}}}, \frac{1-{\delta }_{T}{\delta }_{V}{T}_{1}^{\mathrm{*}}}{\beta {\delta }_{V}{T}_{1}^{\mathrm{*}}}\right)$; and (c) the second non-trivial equilibrium ${E}_{2} = \left({T}_{2}, {I}_{2}, {V}_{2}\right) = \left({T}_{2}^{\mathrm{*}}{\delta }_{V}, \frac{1-{\delta }_{T}{\delta }_{V}{T}_{2}^{\mathrm{*}}}{\beta \pi {T}_{2}^{\mathrm{*}}}, \frac{1-{\delta }_{T}{\delta }_{V}{T}_{2}^{\mathrm{*}}}{\beta {\delta }_{V}{T}_{2}^{\mathrm{*}}}\right)$. ${T}_{\mathrm{1, 2}}^{\mathrm{*}} = \frac{-{k}_{1}\pm \sqrt{{\mathit{\Delta}} }}{2{k}_{2}}$ are the two solutions of bifurcation equation ${k}_{2}{T}^{2}+{k}_{1}T+{k}_{0} = 0$, and its discriminant is ${\mathit{\Delta}} = {k}_{1}^{2}-4{k}_{2}{k}_{0}$, where ${k}_{2} = \beta \pi (\beta \pi -{\delta }_{T}{\delta }_{V}\gamma)$, ${k}_{1} = \pi \left(\gamma \lambda -E\beta \kappa \right)-{\delta }_{I}(\beta \pi -{\delta }_{T}{\delta }_{V}\gamma)$, and ${k}_{0} = -{\delta }_{I}\gamma \lambda$. ${T}_{1}^{\mathrm{*}}$ and ${T}_{2}^{\mathrm{*}}$ cannot be simplified by the factor ($\beta \pi -{\delta }_{T}{\delta }_{V}\gamma$), because ${k}_{1} = \pi \left(\gamma \lambda -E\beta \kappa \right)-{\delta }_{I}(\beta \pi -{\delta }_{T}{\delta }_{V}\gamma)$ cannot be factored by ($\beta \pi -{\delta }_{T}{\delta }_{V}\gamma$). In this case, ${T}_{1}^{\mathrm{*}}$ and ${T}_{2}^{\mathrm{*}}$ cannot be simplified and cannot provide a simplified conditions/ expression for the existence of an equilibrium.

The Jacobian matrix of the system was given by the following:

where we discuss the character of the eigenvalues for the above equilibria below.

5.   Local stability analysis

5.1.   Stability of virus-free equilibrium ${E}_{0}=\left(\frac{\lambda }{{\delta }_{T}}, \mathrm{0, 0}\right)$

The virus-free equilibrium ${E}_{0}$ corresponds to a failure of the HIV virion to infect a CD4⁺ cell due to the lysing effect of CD8⁺ cells. Evaluating the Jacobian at this equilibrium gives the following:

which gives rise to the characteristic equation, $\rho \left(\tau; {\delta }_{T}, {\delta }_{I}, {\delta }_{V}\right) = \frac{1}{{\delta }_{T}}\left(\tau +{\delta }_{T}\right)$$\left({{\delta }_{T}\tau }^{2}+\left(\kappa E{\delta }_{T}+{\delta }_{T}{\delta }_{I}+{{\delta }_{T}\delta }_{V}\right)\tau +\left({\kappa E\delta }_{T}{\delta }_{V}+{\delta }_{T}{\delta }_{I}{\delta }_{V}-\lambda \pi \beta \right)\right)$.

The overall stability of the virus-free equilibrium ${E}_{0}$ is dependent on the real part of the root of ${\tau }_{0B}$ and ${\tau }_{0C}$ of the quadratic factor, because the root ${\tau }_{0A} = {-\delta }_{T} < 0$ of the linear factor is negative. After calculating ${\tau }_{0B}$ and ${\tau }_{0C}$, we find that the equilibrium is either a stable node or a stable focus when ${\kappa E\delta }_{T}{\delta }_{V}+{\delta }_{T}{\delta }_{I}{\delta }_{V} > \lambda \pi \beta$. Because both the lysing rate $\kappa$ and the CD8⁺ cell count $E$ are bifurcation parameters, the real parts of the two eigenvalues are numerically analyzed in the numerical simulation section.

5.2.   Stability of two non-trivial equilibria ${E}_{1}=\left({T}_{1}, {I}_{1}, {V}_{1}\right)$ and ${E}_{2}=\left({T}_{2}, {I}_{2}, {V}_{2}\right)$

The model might admit two non-trivial equilibria, ${E}_{1} = \left({T}_{1}, {I}_{1}, {V}_{1}\right)$ and ${E}_{2} = \left({T}_{2}, {I}_{2}, {V}_{2}\right)$, where susceptible CD4⁺ cells coexist with the HIV virions. After substituting ${E}_{1} = \left({T}_{1}, {I}_{1}, {V}_{1}\right)$ and ${E}_{2} = \left({T}_{2}, {I}_{2}, {V}_{2}\right)$ into the Jacobian matrix, the characteristic equation for these two solutions, is given by $\rho \left(\tau \right) = {{\tau }^{3}+a}_{2}{\tau }^{2}+{a}_{1}\tau +{a}_{0}$. The interpretation of three coefficients—${a}_{2}$, ${a}_{1}$ and ${a}_{0}$—are provided in Appendix. For this cubic function, the Routh-Hurwitz criterion is used to deduce the parameter values that the three roots have negative real parts ^[38]. This criterion states that, given a general cubic of the form $\rho \left(\lambda \right) = {{\tau }^{3}+a}_{2}{\tau }^{2}+{a}_{1}\tau +{a}_{0}$, two conditions need to be simultaneously met for all roots to have negative real parts (i.e., (ⅰ) $\frac{{a}_{1}{a}_{2}-{a}_{0}}{{a}_{2}} > 0$ and (ⅱ) ${a}_{0} > 0$).

6.   Numerical simulation

In this section, we integrate the HIV replication parameters (shown in Table 1) into system (1) to investigate how the CD8⁺ cell lysing rate $\kappa$, the CD8⁺ cell count $E$, and the saturation parameter $\gamma$ determine the CD4⁺ cell and HIV kinetics. First, we understand how these three parameters determine the existence of the equilibria ${E}_{0}$, ${E}_{1}$, and ${E}_{2}$. Next, we examine the dynamical behaviors of the system on the parameter regions that correspond to a different number of equilibria, for example, the stability of an equilibrium and the existence of a periodic cycle. Conclusively, we demonstrate that the saturated lysing efficiency leads to unexpectedly complex dynamical behaviors, including monostability to either a high or low viral titer, bistability, and periodic oscillations.

6.1.   Existence of equilibria

The CD8⁺ cell counts range from $300\ cells/\mu l$ to $2100\ cells/\mu l$ before and during cART, respectively ^[36]. Before cART, the CD8⁺ cell counts approximately range from $600\ cells/\mu l$ to $1200\ cells/\mu l$; during cART, the CD8⁺ cell counts approximately range from $100\ cells/\mu l$ to $600\ cells/\mu l$. First, we fixed the CD8⁺ cell count as $E = 400$ (shown in Figure 2A), $E = 800$ (shown in Figure 2B), and $E = 1600$ (shown in Figure 2C) to investigate the role of the lysing rate $\kappa$ and the saturation parameter $\gamma$ on HIV kinetics. The parameter region of the lysing rate $\kappa$ and the saturation parameter $\gamma$ was divided into six regions according to the positivity of ${T}_{\mathrm{1, 2}}$, ${I}_{\mathrm{1, 2}}$, and ${V}_{\mathrm{1, 2}}$. First, one positive equilibrium ${E}_{2}$ exists when ${T}_{2} > 0$, ${I}_{2} > 0$, and ${V}_{2} > 0$ are in the green and yellow regions (shown in Figure 2A–C). ${T}_{1} > 0$, ${I}_{1} < 0$, and ${V}_{1} < 0$ are in the green region, but ${T}_{1} < 0$, ${I}_{1} < 0$, and ${V}_{1} < 0$ are in the yellow region (shown in Figure 2A–C). Therefore, regardless of the magnitude of the saturation parameter $\gamma$, the low lysing rate $\kappa$ leads to one positive equilibrium ${E}_{2}$.

Figure 2.  Regions representing the existence of the nonzero equilibrium A–F. In Panel A–C, the existence of the nonzero equilibrium changes with the saturation parameter $\eta$ and the CD8⁺ cell lysing rate $\kappa$ at a CD8⁺ cell count $E = 400, \ 800, \ 1600$. In Panel D–F, the existence of the nonzero equilibrium changes with the lysing rate $\kappa$ and the CD8⁺ cell count $E$ at the saturation parameter $\gamma = 3, \ 5, \ 8$. The green region represents ${T}_{2}$, ${I}_{2}$, ${V}_{2}$, and ${T}_{1}$ are positive, and ${I}_{1}$ and ${V}_{1}$ are negative; the yellow region represents ${T}_{2}$, ${I}_{2}$, and ${V}_{2}$ are positive, and ${T}_{1}$, ${I}_{1}$, and ${V}_{1}$ are negative; the purple region represents ${T}_{2}$, ${I}_{2}$, ${V}_{2}$, ${T}_{1}$, ${I}_{1}$, and ${V}_{1}$ are positive; the cyan region represents ${T}_{2}$, ${I}_{2}$, ${V}_{2}$, ${T}_{1}$, ${I}_{1}$, and ${V}_{1}$ do not exist; and the blue region represents ${T}_{2}$, ${I}_{2}$, ${V}_{2}$, ${T}_{1}$, ${I}_{1}$, and ${V}_{1}$ are negative, ${T}_{2}$ is positive, and ${I}_{2}$, ${V}_{2}$, ${T}_{1}$, ${I}_{1}$, and ${V}_{1}$ are negative.

DownLoad: Full-Size Img PowerPoint

Next, two positive equilibria, namely ${E}_{1}$ and ${E}_{2}$, exist when ${T}_{2} > 0$, ${I}_{2} > 0$, ${V}_{2} > 0$, ${T}_{1} > 0$, ${I}_{1} > 0$, ${V}_{1} > 0$ are in purple region (shown in Figure 2A–C); more specifically, the high lysing rate $\kappa$ and the high saturation parameter $\gamma$ always lead to the two equilibria ${E}_{1}$ and ${E}_{2}$. Finally, no positive equilibrium and only the virus-free equilibrium ${E}_{0}$ exists in the cyan, blue, and red regions (shown in Figure 2A–C). Due to a negative discriminant of the quadratic bifurcation equation (Eq (2)), ${T}_{\mathrm{1, 2}}$, ${I}_{\mathrm{1, 2}}$, and ${V}_{\mathrm{1, 2}}$ did not have real solutions in the cyan region. Because ${T}_{2} < 0$, ${I}_{2} < 0$, ${V}_{2} < 0$, ${T}_{1} < 0$, ${I}_{1} < 0$, and ${V}_{1} < 0$ were in the blue region, no positive equilibrium existed. Though ${T}_{2} > 0$ but ${I}_{1} < 0$, ${V}_{1} < 0$, ${T}_{2} < 0$, ${I}_{2} < 0$, ${V}_{2} < 0$, and ${T}_{1} < 0$ were in the red region, no positive equilibrium existed either. Thus, the high lysing rate $\kappa$ and the low saturation parameter $\gamma$ did not lead to the existence of a positive equilibrium. Moreover, the parameter region which corresponded to one positive equilibrium (${E}_{2}$, green and yellow regions) and two positive equilibria (${E}_{1}$ and ${E}_{2}$, purple region) decreased with the CD8⁺ cell count; the parameter region which corresponded to no positive equilibrium (${E}_{0}$, cyan, blue, and red regions) increased with the CD8⁺ cell count (shown in Figure 2A–C). Biologically, the higher the CD8⁺ cell counts were, the higher the possibility that good infection outcomes were achieved.

Next, we fixed the saturation parameter $\mathrm{\gamma } = 3$, $\mathrm{\gamma } = 5$, and $\mathrm{\gamma } = 8$, to investigate the role of the CD8⁺ cell lysing rate $\kappa$ and the CD8⁺ cell count $E$ on HIV kinetics. For the saturation parameter $\eta = 3$, one positive equilibrium ${E}_{2}$ existed in the yellow region, and only the virus-free equilibrium ${E}_{0}$ existed in the red region (shown in Figure 2D); specifically, one positive equilibrium ${E}_{2}$ existed if either 1) the lysing rate $\kappa$ was low and the CD8⁺ cell count $E$ was high, or 2) the lysing rate $\kappa$ was high and the CD8⁺ cell count $E$ was low. Otherwise, only the virus-free equilibrium ${E}_{0}$ existed. For the saturation parameter $\eta = 5$ and $\eta = 8$, one positive equilibrium ${E}_{2}$ existed in the green region, two positive equilibria (${E}_{2}$ and ${E}_{1}$) existed in the purple region, and only the virus-free equilibrium ${E}_{0}$ existed in the cyan and blue regions (shown in Figure 2E, F); specifically, one positive equilibrium ${E}_{2}$ existed if either 1) the lysing rate $\kappa$ was low and the CD8⁺ cell count $E$ was high, or 2) the lysing rate $\kappa$ was low and the CD8⁺ cell count $E$ was high. Only the virus-free equilibrium ${E}_{0}$ existed if the lysing rate $\kappa$ and the CD8⁺ cell count $E$ were high. Otherwise, the two positive equilibria ${E}_{2}$ and ${E}_{1}$ existed. Additionally, the parameter region which corresponded to two positive equilibria (${E}_{2}$ and ${E}_{1}$) and one positive equilibrium (${E}_{2}$) increased with the saturation parameter $\gamma$. The parameter region which corresponded to the virus-free equilibrium (${E}_{0}$) decreased with the saturation parameter $\gamma$.

In the following context, we used a CD8⁺ cell count $E = 800$ as an example to investigate the dynamical behaviors of HIV dynamics on parameter region which correspond to a different number of equilibria, for example, monostable, bistable, and oscillating HIV kinetics (shown in Figure 3A). The CD8⁺ cell induced monostable HIV kinetics implies that the viral titer converges to an equilibrium independent of the initial conditions (shown in Figure 3B). The CD8⁺ cell induced oscillating HIV kinetics implies that the viral titer exhibits a periodical behavior independent of the initial conditions (shown in Figure 3C). The CD8⁺ cell induced bistable HIV kinetics implies that high viral inoculums remained infectious, and low viral inoculums led to an inhibition of the viral infectivity (shown in Figure 3D).

Figure 3.  Regions representing the existence of equilibria and three possible dynamical behaviors. A: Region representing the existence of equilibria ${E}_{0}$, ${E}_{1}$, and ${E}_{2}$ with the saturation parameter $\eta$ and the lysing rate $\kappa$ at a CD8⁺ cell count $E = 800$. B: The effector cell induced monostable HIV kinetics implies that the viral titer converges to an equilibrium independent of the initial conditions. C: The effector cell induced oscillating HIV kinetics implies that the viral titer exhibits a periodical behavior independent of the initial conditions. D: The effector cell induced bistable HIV kinetics implies that the high viral inoculums remain infectious, and a low viral inoculum leads to an inhibition of the viral infectivity.

DownLoad: Full-Size Img PowerPoint

6.2.   Local analysis of equilibrium

First, the parameter region which corresponds to only the virus-free equilibrium ${E}_{0}$ existing consists of a red square ($\kappa > 1.2$ and $0 < \eta < 3.3$) and a cyan and blue trapezoid region ($\kappa > 1.2$ and $3.3 < \eta < 10$, shown in Figure 4A). In this parameter region, the virus-free equilibrium ${E}_{0}$ is a locally asymptotically stable node due to three negative eigenvalues: ${\lambda }_{0A}$, ${\lambda }_{0B}$ and ${\lambda }_{0C}$. Obviously, ${\lambda }_{0A} = {-\delta }_{T}$ is negative due to a positive death rate of the susceptible cells, and ${\lambda }_{0B}$ and ${\lambda }_{0C}$ are also negative in the corresponding parameter region (shown in Figure 4B, C).

Figure 4.  Bifurcation diagram, two eigenvalues on the virus-free equilibrium ${E}_{0}$, and typical numerical simulation of susceptible cells, infected cells, and viral titer. A: Bifurcation diagram with the saturation parameter $\eta$ and the lysing rate $\kappa$ corresponding to the case that no positive equilibrium exists. B and C: Two eigenvalues on the virus-free equilibrium ${E}_{0}$ are negative on the parameter region where no positive equilibrium exists. D–F: Susceptible cell kinetics, infected cells kinetics, and viral titer kinetics converge to a high susceptible cell count, a low infected cell count, and a low viral titer independent of the initial infected cell count for the fixed parameter $\kappa = 4$ and $\gamma = 4$. Blue, cyan, red, black, purple, and green curves represent the susceptible cell kinetics, infected cells kinetics, and viral titer kinetics with the initial infected cell counts 101,101.5,102,102.5,103, and 103.5 in Panel D–F. Initial susceptible cell $S\left(0\right) = {10}^{2}$, and initial virial titer $I\left(0\right) = {10}^{0}$.

DownLoad: Full-Size Img PowerPoint

For the fixed lysing parameters $\kappa = 4$ and $\eta = 4$ in the blue parameter region, the model converges to a high susceptible cell count, a low infected cell count, and a low virus titer independent of the initial conditions (shown in Figure 4D–F), because the three eigenvalues at the virus-free equilibrium ${E}_{0} = (1638, \ \ \ \ 0, \ 0)$ are $-0.18$, $-3212$, and $-24.9$. Because the infected cell count $I\left(t\right)$ is positively correlated with the viral titer $V\left(t\right)$, we provide the susceptible cell count $S\left(t\right)$ and the viral titre $V\left(t\right)$ in following text.

Next, the parameter region which corresponds to one positive equilibrium ${E}_{2}$ existing (${T}_{2} > 0$, ${I}_{2} > 0$, ${V}_{2} > 0$, ${T}_{1} > 0$, ${I}_{1} < 0$ and ${V}_{1} < 0$, green region in Figure 3A) ranges from $\kappa = 0$ to $\kappa = 1.2$ and from $\eta = 3.3$ to $\eta = 10$ (shown in Figure 5A). Unexpectedly, the system converges to a low susceptible cell count and a high viral titer in a majority of the parameter region (green region in Figure 5A), but also exhibits oscillating viral kinetics in a minority of the parameter region (black region, low concern in Figure 5A). In this parameter region, the virus-free equilibrium ${E}_{0}$ is a saddle point due to two negative eigenvalues ${\tau }_{0A} = -{\delta }_{T} < 0$ and ${\tau }_{0B} < 0$ (shown in Figure 5B) and one positive eigenvalue ${\tau }_{0C} > 0$ (shown in Figure 5C). Due to $\frac{{a}_{1}{a}_{2}-{a}_{0}}{{a}_{2}} > 0$ and ${a}_{0} > 0$ in a majority of the parameter region (shown in Figure 5D, E), one positive equilibrium ${E}_{2}$ is locally asymptotically stable by the Routh–Hurwitz criterion (green region in Figure 5A); similarly, equilibrium ${E}_{2}$ is unstable in a minority of the parameter region (black region in Figure 5A).

Figure 5.  Bifurcation diagram with one positive equilibrium ${E}_{2}$, two eigenvalues on virus-free equilibrium ${E}_{0}$, and a discriminant for the Routh–Hurwitz criterion and constant term in the characteristic equation on the positive equilibrium ${E}_{2}$. A: Bifurcation diagram with the saturation parameter $\eta$ and the lysing rate $\kappa$ corresponding to the case that one positive equilibrium ${E}_{2}$ exists. Green region represents that the system converges to a high viral titer and low susceptible cell counts; black region represents the system exhibiting oscillating viral kinetics. B and C: One eigenvalue on the virus-free equilibrium ${E}_{0}$ is negative on the parameter region that one positive equilibrium ${E}_{2}$ exists, and another is positive. D and E: Discriminant for Routh–Hurwitz criterion and constant term in the characteristic equation on the positive equilibrium ${E}_{2}$.

DownLoad: Full-Size Img PowerPoint

Numerically, for the fixed lysing parameters $\kappa = 1$ and $\eta = 5$, the system provides two equilibria: ${E}_{0} = (1638, \ 0, \ 0)$ and ${E}_{2} = (3.67, \ 134.5, \ \mathrm{205,825})$. Three eigenvalues at the equilibrium ${E}_{0}$ are ${\tau }_{0A} = -0.18$, ${\tau }_{0B} = -844$, and ${\tau }_{0C} = 7$; three eigenvalues at the equilibrium ${E}_{2}$ are ${\tau }_{2A} = -82$, ${\tau }_{2B} = -34$, and ${\tau }_{2C} = -1$. Thus, the system monotonically converges to a low susceptible cell count and a high virus titer independent of the initial conditions (shown in Figure 6A, C, E).

Figure 6.  Typic numeric simulations of system (1) demonstrating oscillating and monostable viral kinetics. In Panel A, C, and E, system (1) converges to a lower susceptible cell count and a high viral titer for the fixed parameters $\kappa = 1$ and $\eta = 5$. In Panel B, D, and F, system (1) exhibits an oscillating susceptible cell and viral kinetics for the fixed parameters $\kappa = 1.18$ and $\eta = 3.4$. Blue, cyan, red, black, purple, and green curves represent the susceptible cell kinetics, infected cells kinetics, and viral titer kinetics with the initial infected cell counts 101,101.5,102,102.5,103, and 103.5. Initial susceptible cell $S\left(0\right) = {10}^{2}$, and initial virial titer $I\left(0\right) = {10}^{0}$.

DownLoad: Full-Size Img PowerPoint

Similarly, for the fixed lysing parameters $\kappa = 1.18$ and $\eta = 3.4$, the system provides two equilibria: ${E}_{0} = (1638, \ 0, \ 0)$ and ${E}_{2} = (28.74, \ 16.92, \ \mathrm{25,853})$. Three eigenvalues at the equilibrium ${E}_{0}$ are ${\tau }_{0A} = -0.18$, ${\tau }_{0B} = -982$, and ${\tau }_{0C} = 1$; three eigenvalues at the equilibrium ${E}_{2}$ are ${\tau }_{2A} = -51$, ${\tau }_{2B} = 1.9+1.84i$, and ${\tau }_{2C} = 1.9-1.84i$. Due to well-posedness, the model exhibits an oscillating susceptible cell and viral kinetics independent of the initial conditions (shown in Figure 6B, D, F).

The parameter region which corresponds to one positive equilibrium ${E}_{2}$ existing (${T}_{2} > 0$, ${I}_{2} > 0$, ${V}_{2} > 0$, ${T}_{1} < 0$, ${I}_{1} < 0$ and ${V}_{1} < 0$) ranges from $\kappa = 0$ to $\kappa = 1.2$ and from $\eta = 0$ to $\eta = 3.3$ (shown in Figure 7A). System (1) converges to a low susceptible cell count and a high viral titer (yellow region, shown in Figure 7A), and exhibits oscillating viral kinetics (black region, shown in Figure 7A). The parameter region which corresponds to monostable virus kinetics approximately equals to that which corresponds to oscillating viral kinetics (shown in Figure 7A). In this parameter region, the virus-free equilibrium ${E}_{0}$ is a saddle point due to two negative eigenvalues ${\tau }_{0A} = -{\delta }_{T} < 0$ and ${\tau }_{0B} < 0$ (shown in Figure 7B) and one positive eigenvalue ${\tau }_{0C} > 0$ (shown in Figure 7C). According to $\frac{{a}_{1}{a}_{2}-{a}_{0}}{{a}_{2}} > 0$ and ${a}_{0} > 0$ (shown in Figure 7D, E), one positive equilibrium ${E}_{2}$ is locally asymptotically stable by the Routh–Hurwitz criterion in the yellow region (shown in Figure 7A); similarly, equilibrium ${E}_{2}$ is unstable in the black region (shown in Figure 7A).

Figure 7.  Bifurcation diagram with one positive equilibrium ${E}_{2}$, two eigenvalues on virus-free equilibrium ${E}_{0}$, and discriminant for the Routh–Hurwitz criterion and constant term in the characteristic equation on the positive equilibrium ${E}_{2}$. A: Bifurcation diagram with the saturation parameter $\eta$ and the lysing rate $\kappa$ corresponds to the case that one positive equilibrium ${E}_{2}$ exists. Yellow region represents that the system converges to a high viral titer and low susceptible cell counts; black region represents the system exhibiting oscillating viral kinetics. B and C: One eigenvalue on the virus-free equilibrium ${E}_{0}$ is negative on the parameter region that one positive equilibrium exists, and another is positive. D and E: Discriminant for the Routh–Hurwitz criterion and the constant term in the characteristic equation on the positive equilibrium ${E}_{2}$.

DownLoad: Full-Size Img PowerPoint

Numerically, for the fixed lysing parameters $\kappa = 0.5$ and $\gamma = 2$, the model provides two equilibria: ${E}_{0} = (1638, \ 0, \ 0)$ and ${E}_{2} = (5.189, \ 95.1119, \ \mathrm{145,310})$. Three eigenvalues at the equilibrium ${E}_{0}$ are ${\tau }_{0A} = -0.18$, ${\tau }_{0B} = -480$, and ${\tau }_{0C} = 43$; three eigenvalues at the equilibrium ${E}_{2}$ are ${\tau }_{2A} = -61$, ${\tau }_{2B} = -31$, and ${\tau }_{2C} = -1$. Thus, the system monotonically converges to a low susceptible cell count and a high virus titer independent of the initial conditions (shown in Figure 8A, C, E).

Figure 8.  Typic numeric simulations of system (1) demonstrating oscillating and monostable viral kinetics. In Panel A, C, and E, system (1) converges to a lower susceptible cell count and a high viral titer for the fixed parameters $\kappa = 0.5$ and $\gamma = 2$. In Panel B, D, and F, system (1) exhibits an oscillating susceptible cell and viral kinetics for the fixed parameters $\kappa = 0.5$ and $\gamma = 0.7$. Blue, cyan, red, black, purple, and green curves represent the susceptible cell kinetics, infected cells kinetics, and viral titer kinetics with the initial infected cell counts 101,101.5,102,102.5,103, and 103.5. Initial susceptible cell $S\left(0\right) = {10}^{2}$, and initial virial titer $I\left(0\right) = {10}^{0}$.

DownLoad: Full-Size Img PowerPoint

Similarly, for the fixed lysing parameters $\kappa = 0.5$ and $\gamma = 0.7$, the model provides two equilibria: ${E}_{0} = (1638, \ 0, \ 0)$ and ${E}_{2} = (415.2985, \ 0.8901, \ 1359.8)$. Three eigenvalues at the equilibrium ${E}_{0}$ are ${\tau }_{0A} = -0.18$, ${\tau }_{0B} = -480$, and ${\tau }_{0C} = 43$; three eigenvalues at the equilibrium ${E}_{2}$ are ${\tau }_{2A} = -205$, ${\tau }_{2B} = 15$, and ${\tau }_{2C} = 0.7$. Due to well-posedness, the model exhibits an oscillating susceptible cell count and viral kinetics independent of the initial conditions (shown in Figure 8B, D, F). Compared to the lysing parameters $\kappa = 1.18$ and $\eta = 3.4$, the lysing parameters $\kappa = 0.5$ and $\gamma = 0.7$ provide a higher frequency of the viral titer and the susceptible cell count.

Finally, the parameter region which corresponds to two positive equilibria ${E}_{1}$ and ${E}_{2}$ existing (${T}_{2} > 0$, ${I}_{2} > 0$, ${V}_{2} > 0$, ${T}_{1} > 0$, ${I}_{1} > 0$ and ${V}_{1} > 0$) ranges from $\kappa = 1.2$ to $\kappa = 3.5$ and from $\eta = 3.3$ to $\eta = 10$ in a triangle region (shown in Figure 9A). System (1) exhibits a bistable behavior in the purple triangle region, and a monostable lower viral titer in the green region (shown in Figure 9A). The parameter region which corresponds to the bistable viral kinetics is much larger than that which corresponds to monostable viral kinetics (shown in Figure 9A). In this parameter region, the virus-free equilibrium ${E}_{0}$ is a stable node point due to three negative eigenvalues: ${\lambda }_{0A} = -{\delta }_{T} < 0$, ${\lambda }_{0B} < 0$ (shown in Figure 9B) and ${\lambda }_{0C} < 0$ (shown in Figure 7C). Due to $\frac{{a}_{1}{a}_{2}-{a}_{0}}{{a}_{2}} < 0$ and ${a}_{0} < 0$ (shown in Figure 9D, E), the positive equilibrium ${E}_{1}$ is unstable by the Routh–Hurwitz criterion. Similarly, the positive equilibrium ${E}_{2}$ is stable when $\frac{{a}_{1}{a}_{2}-{a}_{0}}{{a}_{2}} > 0$ and ${a}_{0} > 0$ (purple region in Figure 9A, F, G), and is unstable when $\frac{{a}_{1}{a}_{2}-{a}_{0}}{{a}_{2}} < 0$ and ${a}_{0} > 0$ (green region in Figure 9A, F, G).

Figure 9.  Bifurcation diagram with two positive equilibria ${E}_{1}$ and ${E}_{2}$, two eigenvalues on virus-free equilibrium ${E}_{0}$, and discriminant for the Routh–Hurwitz criterion and constant term in the characteristic equation on the positive equilibria ${E}_{1}$ and ${E}_{2}$. A: Bifurcation diagram with the saturation parameter $\eta$ and the lysing rate $\kappa$ corresponds to the case that two positive equilibria, namely ${E}_{1}$ and ${E}_{2}$, exist. Purple region represents that the system exhibits a bistable susceptible cell and viral kinetics; green region represents the system converging to a high susceptible cell count and a lower viral titer. B and C: Two eigenvalues on the virus-free equilibrium ${E}_{0}$ are negative. D and E: Discriminant for the Routh–Hurwitz criterion and the constant term in the characteristic equation on the positive equilibrium ${E}_{1}$ are negative. F and G: Discriminant for the Routh–Hurwitz criterion and the constant term in the characteristic equation on the positive equilibrium ${E}_{2}$.

DownLoad: Full-Size Img PowerPoint

Numerically, for the fixed lysing parameters $\kappa = 2$ and $\eta = 7$, the model provides three equilibria: ${E}_{0} = (1638, \ 0, \ 0)$, ${E}_{1} = (694.293, \ 0.411, \ 627.93)$, and ${E}_{2} = (7.5159, \ 65.5726, \ \mathrm{100,180})$. Three eigenvalues at the equilibrium ${E}_{0}$ are ${\tau }_{0A} = -0.18$, ${\tau }_{0B} = -1623$, and ${\tau }_{0C} = -13$; three eigenvalues at the equilibrium ${E}_{1}$ are ${\tau }_{1A} = -198$, ${\tau }_{1B} = 55$, and ${\tau }_{1C} = -0.1$; three eigenvalues at the equilibrium ${E}_{2}$ are ${\tau }_{2A} = -50$, ${\tau }_{2B} = -24$, and ${\tau }_{2C} = -1$. Thus, the model exhibits a bistable behavior, namely large viral inoculum sizes lead to the maintenance of viral infectivity and small viral inoculum sizes lead to an inhibition of the viral infectivity (shown in Figure 10A, C).

Figure 10.  Typic numeric simulations of system (1) demonstrating bistable and monostable viral kinetics. In Panel A and C, system (1) converges to bistable susceptible cells and viral kinetics for the fixed parameters $\kappa = 2$ and $\eta = 7$. In Panel B and D, system (1) exhibits monostable kinetics to a high susceptible cell count and a lower viral titer for the fixed parameters $\kappa = 2.4$ and $\eta = 7$. Blue, cyan, red, purple, and green curves represent the susceptible cell kinetics, infected cells kinetics, and viral titer kinetics with the viral inoculum sizes 102.9,103.2,103.5,103.8, and 104.1. Initial susceptible cell $S\left(0\right) = {10}^{2}$, and initial virial titer $I\left(0\right) = {10}^{0}$.

DownLoad: Full-Size Img PowerPoint

Moreover, for the fixed lysing parameters $\kappa = 2.4$ and $\eta = 7$, the model provides three equilibria: ${E}_{0} = (1638, \ 0, \ 0)$, ${E}_{1} = (192.9549, \ 2.2638, \ 3458.6)$, and ${E}_{2} = (27.0437, \ 18.0055, \ \mathrm{27,508})$. Three eigenvalues at the equilibrium ${E}_{0}$ are ${\tau }_{0A} = -0.18$, ${\tau }_{0B} = -1939$, and ${\tau }_{0C} = -17$; three eigenvalues at the equilibrium ${E}_{1}$ are ${\tau }_{1A} = -88$, ${\tau }_{1B} = 431$, and ${\tau }_{1C} = -0.1$; three eigenvalues at the equilibrium ${E}_{2}$ are ${\tau }_{2A} = -50$, ${\tau }_{2B} = 1.39+2.15\mathrm{i}$, and ${\tau }_{2C} = 1.39-2.15\mathrm{i}$. The model monotonically converges to a high susceptible cell count and a low virus titer independent of the initial condition (shown in Figure 10B, D).

7.   Discussion

The model analyzed in this work displayed several interesting and unexpected features, both from mathematical and biological perspectives. First, a range of possible dynamical outcomes based on the value of the model parameters were found. Numerous nontrivial dynamical behaviors were discovered, with the presence of an important system equilibrium (i.e., bistability and periodic oscillation) that was described by a Jacobian matrix and the Routh–Hurwitz criterion. The lysing rate, CD8⁺ cell count, and saturation effect that controlled the lysing efficiency as the infected cell count increased were combined to determine the HIV kinetics. The HIV kinetics with the CD8⁺ cell exhibited unexpectedly complex behaviors that were both mathematically and biologically interesting, for example, monostability, bistability, and periodic oscillations depending on the initial conditions.

The model provides a few insights into the interactions among CD4⁺ cells, HIV, and CD8⁺ cells with the saturated lysing efficiency. One of the main limitations of the present approach is that CD8⁺ cell concentration is assumed to be a constant. Upon analysis (Figures 4, 6, 8, and 10), we found that the viral titer and the susceptible cell kinetics rapidly converged to stable steady states in several days or a period of periodic oscillations that lasted between ten to fifty days. Because the life expectancy of patients without cART is around ten-fifteen years and is more than several decades with cART, it is reasonable to accept the assumption that the CD8⁺ cell count is constant during a couple of days.

For the fixed CD8⁺ cell count, the low lysing rate led to monostability to a high viral titer with a high magnitude of the saturation effect (shown in Figures 5A and 6A, C) and periodic oscillations with a low magnitude of the saturation effect (shown in Figures 7A and 8A, C); this explains why patients with non-protective HLA alleles always lead to bad infection outcomes. A high CD8⁺ cell lysing rate induced to bistability (shown in Figures 9A and 10A, C); this explains why patients with protective HLA alleles sometimes exhibited good infection outcomes but with some bad infection outcomes. This theoretical result may challenge the idea that the lysing efficiency of the CD8⁺ cells is used as a predictor of good or bad infection outcomes, because the HIV replication kinetics and the CD8⁺ cell lysing kinetics combined to determine the infection outcomes. Moreover, the lysing efficiency can be used as a predictor of infection outcomes if the HIV kinetics only exhibits monostable kinetics independent of the initial conditions.

The HIV infection outcomes with the effector CD8⁺ cell is strongly associated with the CD8⁺ lysing rate and the effector CD8⁺ cell concentration (shown in Figures 5A and 9A). For a given saturation parameter, we have two possible scenarios: 1) Figure 5A consists of monostability and periodic oscillations; and 2) Figure 9A consists of bistable viral kinetics. For scenarios 1), slightly increasing the lysing rate above the critical value can change the monostability to a high viral titer and to periodic oscillations. For scenarios 2), slightly increasing the lysing rate above the critical values changes the bistable viral kinetics to monostability to a low viral titer; then, using cART pushes the infected CD4⁺ cell count and the viral titer under the critical values, thereby inducing good infection outcomes.

Reference ^[28] proposed that immunosuppressed viral infection exhibits bistable HIV kinetics through a saddle–node bifurcation. Reference ^[27] used monotonic and nonmonotonic immune responses of CD8⁺ cell to display bistable HIV kinetics through a saddle–node bifurcation. Here, we found that the semi-saturated lysing model led to bistable HIV kinetics through a saddle–node bifurcation and oscillatory HIV kinetics though the Hopf bifurcation.

In the future modelling works, it is worthwhile to use the Monod-Haldane function to the delineate immune response stimulated by the HIV titer ^[27]. For the non-protective HLA allele, it is meaningful to determine which CD8⁺ cell count interval corresponds to monostable viral kinetics and bistable viral kinetics. It is meaningful to investigate whether using cART pushes the infected CD4⁺ cell count and viral titer under critical values for the protective HLA allele carrier, thereby inducing good infection outcomes. In a similar framework, it is interesting to understand how CAR (Chimeric Antigen Receptor) -T cell and CAR-NK therapy affect HIV infections.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The author declares there is no conflict of interest.

Appendix

A1. Coefficients in characteristic polynomial

In Section Stability of two positive equilibria ${E}_{1} = \left({T}_{1}, {I}_{1}, {V}_{1}\right)$ and ${E}_{2} = \left({T}_{2}, {I}_{2}, {V}_{2}\right)$, we provided characteristic polynomial on equilibria ${E}_{1} = \left({T}_{1}, {I}_{1}, {V}_{1}\right)$ and ${E}_{2} = \left({T}_{2}, {I}_{2}, {V}_{2}\right)$, $\rho \left(\tau \right) = {{\tau }^{3}+a}_{2}{\tau }^{2}+{a}_{1}\tau +{a}_{0}$. Here, we provide closed-form expression of three coefficients in characteristic polynomial, ${a}_{2}$, ${a}_{1}$ and ${a}_{0}$.

First, the closed-form expression of coefficient ${a}_{2}$ is ${a}_{2} = \frac{1}{{\left(1+\gamma I\right)}^{2}}({I}^{2}{V}^{2}\beta {\gamma }^{2}+{I}^{2}{\gamma }^{2}{\delta }_{T}$+${I}^{2}{\gamma }^{2}{\delta }_{I}+{I}^{2}{\gamma }^{2}{\delta }_{V}+2IV\beta \gamma$+$2I{\delta }_{T}\gamma +2I{\delta }_{I}\gamma +2I{\delta }_{V}\gamma +\kappa E+V\beta +{\delta }_{T}+{\delta }_{I}+{\delta }_{V})$. Next, the closed-form expression of coefficient ${a}_{1}$ is ${a}_{1} = \frac{1}{{\left(1+\gamma I\right)}^{2}}(-{I}^{2}T\beta {\gamma }^{2}{\pi }_{1}+{I}^{2}V\beta {\delta }_{I}{\gamma }^{2}+{I}^{2}V\beta {\delta }_{V}{\gamma }^{2}$+${I}^{2}{\delta }_{T}{\delta }_{I}{\gamma }^{2}+{I}^{2}{\delta }_{T}{\delta }_{V}{\gamma }^{2}+{I}^{2}{\delta }_{I}{\delta }_{V}{\gamma }^{2}-2IV\beta \gamma \pi$+$2IV\beta {\delta }_{I}\gamma +2IV\beta {\delta }_{V}+EV\beta \kappa +2I{\delta }_{T}{\delta }_{I}\gamma$+$2I{\delta }_{T}{\delta }_{V}\gamma +2I{\delta }_{I}{\delta }_{V}\gamma +E{\delta }_{T}\kappa +E{\delta }_{V}\kappa -\pi T\beta +V\beta {\delta }_{I}+V\beta {\delta }_{V}+{\delta }_{T}{\delta }_{I}+{\delta }_{T}{\delta }_{V}+{\delta }_{I}{\delta }_{V})$. Next, the closed-form expression of coefficient ${a}_{0}$ is ${a}_{0} = \frac{1}{{\left(1+\gamma I\right)}^{2}}(-{I}^{2}T\beta {\delta }_{T}\gamma {\pi }_{1}+{I}^{2}V\beta {\delta }_{I}{\delta }_{V}{\gamma }^{2}$+$I{\delta }_{T}{\delta }_{I}{\delta }_{V}{\gamma }^{2}-2IT\beta {\delta }_{T}\gamma \pi +2IV\beta {\delta }_{I}{\delta }_{V}\gamma$+$EV\beta {\delta }_{V}\kappa +2I{\delta }_{T}{\delta }_{I}{\delta }_{V}\gamma +E{\delta }_{T}{\delta }_{V}\kappa -T\beta {\delta }_{T}\pi +V\beta {\delta }_{I}{\delta }_{V}+{\delta }_{T}{\delta }_{I}{\delta }_{V})$, $T$, $I$ and $V$ are value in ${E}_{1} = \left({T}_{1}, {I}_{1}, {V}_{1}\right)$ and ${E}_{2} = \left({T}_{2}, {I}_{2}, {V}_{2}\right)$.