Stability Analysis and Cauchy Matrix of a Mathematical Model of Hepatitis B Virus with Control on Immune System near Neighborhood of Equilibrium Free Point

: Mathematical models are useful tools to describe the dynamics of infection and predict the role of possible drug combinations. In this paper, we present an analysis of a hepatitis B virus (HBV) model including cytotoxic T lymphocytes (CTL) and antibody responses, under distributed feedback control, expressed as an integral form to predict the effect of a combination treatment with interleukin-2 (IL-2). The method presented in this paper is based on the symmetry properties of Cauchy matrices C ( t , s ) , which allow us to construct and analyze the stability of corresponding integro-differential systems.


Introduction
The hepatitis B virus (HBV) disease is defined as the detection of the HBV surface antigen (HBsAg) on two occasions measured at least six months apart. HBV represents a huge problem for public health, increasing the risk of cirrhosis and hepatocellular carcinoma in the population [1].
HBV is a virus that belongs to the family of the Herpesviridae. It needs to integrate its DNA with the host cell DNA to survive, forming a covalently closed circular DNA, enhancing DNA replication and cellular division. For this reason, chronic HBV infection accounts for at least 50% of hepatocellular carcinoma cases [2], showing a high mortality rate.
Vaccination strategies have reduced the appearance of all of these consequences, although it is possible to lose the antibody protection against the virus and a vaccine booster would be required to avoid disease onset [3,4]. Unfortunately, it is impossible to check for this antibody protection without any serological control and often the early stages of the infection are completely asymptomatic, leading to a need for severe treatments [5].
Current therapy for chronic HBV infection is based on two strategies: interferons and nucleoside analogue [6]. Interferons (α, β, and γ) are cytokines that are endogenously produced by immune system cells in response to viral infections. They have antiviral effects, creating a barrier around the cell that reduces the rates of infected cell formation and virus replication [7]. For this reason, injectable formulations of pegylated interferon alfa (pegylated interferon alpha-2a and alpha-2b) are available for HBV therapy.
Five nucleoside analogues are approved in the United States: lamivudine, adefovir, entecavir, tenofovir disoproxil, and tenofovir alafenamide [8]. Their role is to block virus replication by blocking the reverse transcriptase.
Even with all of these therapeutical options, many resistance processes occur and therapeutic goals are rarely achieved [9][10][11][12]. Possible drug improvement could be obtained with the integration of interleukin-2 (IL-2) and new clinical trials are evaluating its impact on the combination therapy (NCT02360592, NCT00451984). There is much evidence supporting this hypothesis: in fact, increased IL-2 levels seem to stimulate the activities of NK cells and CD8 lymphocytes, which are responsible for viral clearing [13]. Mathematical models are a useful tool for elucidating and representing molecular mechanisms and try to abstract and predict new biological insights. Moreover, they are a tool for studying many diseases (diabetes, various types of cancer) and responses to therapy [14]. In recent years, many mathematical approaches have been used to effectively solve this problem, such as delay differential equations (DDEs) [15] or fractional order models [16]. Previous models have tried to predict the effectiveness of antiviral therapy, but without considering the important role of the immune system [17][18][19]. This approximation does not allow for realistic predictions, ignoring many crucial factors such as the lifespan of HBV-infected cells that varies greatly due to the strength of the anti-HBV CTL response [20].
Long et al. have tried to overcome these limits in their study [5], obtaining more reliable results, but without considering any treatment options. A complete model was presented in [21], where immune system and drug response were both taken into account. In this study, we propose a mathematical model for HBV infection, using a system of functional differential equations. This approach has already proven its efficacy in biology and medicine [22,23]. We have taken into account the role of the immune system and the two different treatments (interferon and nucleoside analogues), and we have included a control function to optimize the role of a possible IL-2 co-treatment. Yousdi et al. presented a global analysis of a hepatitis C virus (HCV) model with CTL and antibody responses [1]. The model consisted of five non-linear equations as follows: (1) Each of the five variables describes the state of the virus-host system over time. These variables are: uninfected cells, infected cells, free virus numbers, antibody response and CTL response. Those are denoted respectively by: X (cells/mL), Y (cells/mL), V (IU/mL), W (IU/mL), and Z (cells/mL). X cells are produced at a rate r, die at a rate dX and are infected by virus at a rate βVX. The infected cells, Y, grow at the rate of βVX, die (naturally) at a rate aY, and are killed by the CTLs at a rate pYZ. Free viruses V are produced by infected cells Y at a rate kY, decay at a rate uV and are neutralized by antibodies at a rate qVW. Antibodies W develop in response to free viruses at a rate gvW, and decay at a rate hW. The number of CTLs, Z, expands in response to infected cells at a rate cYZ and decays in the absence of infection at a rate bZ.
The dynamics of both HBV and HCV are similar in their essence. Both HCV and HBV are viruses belonging to the Herpesviridae family [24]. Hence, we can use a similar model to describe each one of them, using the five variables described above. We can adjust the mathematical model of Yousdi et al. to our study by adapting the correct course of treatment to each disease.

Modified Model of HBV
Let us now discuss of the novelty of our approach. In this paper we consider the efficiency of the interferon therapy and the efficiency of therapy with nucleoside analogues. These are respectively denoted by η and ε, and their values are between 0 and 1. Due to their influence on the therapy, they directly effect the growth or the decay of the infected cells. Taking this into consideration, we obtain: ( Additionally, we wish to introduce the distributed control function U(t), which represents the impact of the IL-2 therapy, as the following integral form: Integro-differential systems were studied in the early stages of the development of theoretical population dynamics. Due to the growth of the population, and the various changes that occurred, Volterra included a functionals of an integral type, which later became classic differential models of population dynamics and mathematical ecology, such as the famous predator-prey system of Volterra. We chose to employ the Volterra integral, a type of control function and stability analysis that can be found in [25][26][27][28][29][30][31].
Since the role of IL-2 is still not very well-characterized and the literature is controversial, we have decided to consider that it is produced mainly by the CTLs, and we applied it to the Z equations due to its role in increasing CTL response and T lymphocyte duplication and activation [32]. Therefore we obtain: We wish to prove that our model is exponential stable. This will help us predict the patient's disease and advancement, and will allow us to properly adjust the therapeutic process to each subject. Yousfi at al. [1] proved that their model is stable using Lyapunov functions near the neighborhood of the equilibria disease-free point.
From a medical point of view, our desire is to support the immune system when transitioning from a state of infectious disease to a stable state that indicates an almost healthy subject. After obtaining this stable state, we take into consideration that every patient has her/his own individual conditions. As suggested in [27], we use a simple method of analysis and estimation based on a reduction of integro-differential systems to ones of ordinary differential equations.
The goal of this paper is to demonstrate a new approach of co-treatment for HBV using a distributed control in the model of the HBV disease. A model of infectious diseases is proposed based on the analysis of the integro-differential systems studied in [27].
Substituting Equations (2)-(4) in (1), the resulting modified model is as follows: Using the reduction method we pass from a integro-differential system (5) to an ordinary differential system (6):

Exponential Stability of System (6) in the Neighborhood of the Equilibrium Free Point
It is clear that: is an equilibrium free point of the system (6). We wish to study the disease free equilibria, since it represents a patient who is either a healthy subject or a recovered one. From this stationary state, it is possible to deduce the others.
The linearized system (6) near the neighborhood of point (7) is as follows: and the corresponding homogeneous system is: We construct the Jacobian matrix of system (8): We substitute our equilibrium point (7) into the Jacobian matrix and obtain the matrix of the coefficients of linearized system (8): The eigenvalues of (9) are: Let us denote: Theorem 1. If all the coefficients of system (6) are positive, and η and ε are parameters defined between 0 to 1, and if inequalities (11) and (12) are fulfilled, then system (6) is exponentially stable.
Proof. It follows that λ 1 and λ 2 are real and negative. Additionally, it is clear that d(a + u) > 0, so that λ 4 is real and negative. λ 3 will be real and negative if inequality (11) is fulfilled. Additionally, b is a real and positive coefficient, and α and D denote our cotreatment, and hence they are positive as well, so that (α − b) 2 + 4D > 0. As so, λ 6 is also real and negative. λ 5 will be real and negative if inequality (12) is fulfilled. (8) Let us construct the Cauchy matrix of system (8) (see, for example, [33]). Let us denote:

Constructing the Cauchy Matrix of System
, , According to (10) and (13) we can write the eigenvalues of (9) as follows: Using Maple, we obtain the eigenvectors of the matrix (9): Let us define a new matrix whose columns are the eigenvectors of (15).
Using Maple, we obtain the inverse matrix of (16): Let us now write the Cauchy matrix C(t, s) = C i (t, s) i=1,...,6 of the system (8). The Cauchy matrix can be written as: where: V i (t) are the eigenvectors we calculated in (15), and λ i are the eigenvalues we calculated in (14) while i = 1, ..., 6. Let us define: Substituting Equations (14), (15), (19), and (20) in (18), we can write the columns of the Cauchy matrix as follows: According to the properties of the Cauchy matrix, C(s, s) = I, where I is the identity (6 × 6) matrix, and it follows that: Let us define vector B i : (16), and b i are as follows: Substituting (22) and (24) in (23) gives us: Taking into account Equations (21) and (25), columns C i (t, s) 1≤i≤6 of the Cauchy matrix C(t, s) of system (8) are as follows:

Exponential Stability of a System with an Uncertain Coefficient
Let us write a system with an uncertain coefficient based on the work we have done in Section 4. Consider the following system with an uncertain coefficient: The uncertain coefficient ∆D(t) can be the result of individual conditions of the patient, because the assimilation of a drug in the body of different patients can have different rates of influence. We assumed that ∆D(t) is essentially a bounded function. Therefore, it is not expected to effect the stability of the system. As before, this system of integro-differential equations can be reduced to the following system of ordinary differential equations: According to the theory presented in [30], we can write our system as follows: where: The general solution of the auxiliary system X − A X = Z can be represented in the following form (see [34]): Finally, we have Z = (∆A(t)) X + F(t), and we obtain: t 0 C(t, s) Z(s)ds = F(t) + (∆A(t))C(t, 0) X(0), which can be written in the operator form: where Ω : ∞ is the space of six vector-functions with essentially bounded components. That leads to: Below we estimate the norm of the operator Ω. It is clear that: Let us denote: and ∆D * = ess sup t≥0 |∆D(t)|.
Therefore we obtain: According to (27): It is clear from the estimates of the elements of the Cauchy matrix that Q j ≤ Q * j for 1 ≤ j ≤ 6. Theorem 2. Let the assumption of Theorem 1 be fulfilled, c 13 = c 14 and the inequality max 1≤i≤6 {Q * i } < 1 be true. Then system (26) is exponential stable.
Proof. The proof follows from the estimates of (27).

Conclusions
In this paper, we presented a mathematical model for HBV infection, considering a possible co-treatment with the standard of care and IL-2.
Our main goal was to gain a better estimation of a patient's disease using our cotreatment approach. In the introduction we have reviewed the known model throughout its variables and understood the medical role of each parameter and variable. Then, in Section 2 we have expanded our observation regarding the uninfected and infected cells, respectively X and Y variables, and, furthermore, accounted for the CTL response to our supporting treatment in the Z equation.
In Section 3 we obtained results about exponential stability of an ODE reduced model near the neighborhood of the disease-free equilibria.
In Section 4 we built a Cauchy matrix of linearized system (8).
In Section 5, we have expanded our observation by adding more depth and complexity to it. It allowed us to gain a more realistic illustration of a person's condition. Each patient has his/her own individual conditions, and we accounted for them in ∆D(t), and wished to study the stability of the new modified model (26) we arrived at.
We used a known method to research this model, based on our previous work in Sections 3 and 4. The model illustrates a biological disease, and as so, all of the parameters and constants of the system must be real, positive, and bounded integers. As mentioned in [1], the parameters play a crucial role in characterizing the stability equilibrium point. A validation of the mathematical model presented in the paper could be achieved in the future through a series of clinical experiments that will test the level of antigens against HBV disease among patients who received the co-treatment we suggested. Moreover, this will provide us with information regarding the individual reactions of the human body to the co-treatment. As denoted in Section 5, the parameter ∆D(t) represents these variations between subjects and allows us to evaluate the effect of the treatment. Additionally, the sensitivity of different patients' reactions can be different, and it can be variable in time.
Funding: This research received no external funding.