International Journal of Mathematics and Mathematical Sciences
Volume 2018 (2018), Article ID 6180892, 11 pages
https://doi.org/10.1155/2018/6180892

Research Article

Analysis of Fractional Order Mathematical Model of Hematopoietic Stem Cell Gene-Based Therapy

Mohammad Imam Utoyo, Windarto, and Aminatus Sa’adah

Department of Mathematics, Faculty of Science and Technology, Universitas Airlangga, Indonesia

Correspondence should be addressed to Mohammad Imam Utoyo; m.i.utoyo@fst.unair.ac.id

Received 18 May 2018; Accepted 7 July 2018; Published 1 August 2018

Academic Editor: Hans Engler

Copyright © 2018 Mohammad Imam Utoyo et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Hematopoietic stem cell (HSC) has been discussed as a basis for gene-based therapy aiming to cure immune system infections, such as HIV. This therapy protects target cells from infections or specifying technic and immune responses to face virus by using genetically modified HSCs. A mathematical model approach could be used to predict the dynamics of HSC gene-based therapy of viral infections. In this paper, we present a fractional mathematical model of HSC gene-based therapy with the fractional order derivative . We determine the stability of fractional model equilibriums. Based on the model analysis, we obtained three equilibriums, namely, free virus equilibrium (FVE) , CTL-Exhaustion Equilibrium (CEE) , and control immune equilibrium (CIE) . Besides, we obtained Basic Reproduction Number that determines the existence and stability of the equilibriums. These three equilibriums will be conditionally locally asymptotically stable. We also analyze the sensitivity of parameters to determine the most influence parameter to the spread of therapy. Furthermore, we perform numerical simulations with variations of to illustrate the dynamical HSC gene-based therapy to virus-system immune interactions. Based on the numerical simulations, we obtained that HSC gene-based therapy can decrease the concentration of infected cells and increase the concentration of the immune cells.

1. Introduction

An immune system is body’s primary defense system that has a function to fight against microbes. Humoral immunity consists of innate immunity and adaptive immunity. The main principle of innate immunity is an initial defense which responds quickly against microbes. Meanwhile, adaptive immunity would be activated when innate immunity failed to eradicate microbes. Adaptive immunity consists of B lymphocytes cells (B cells) and T lymphocytes cells (T cells). T cells consist of CD4+ T cells and cytotoxic T lymphocytes (CTL and CD8+ T cells). Adaptive immunity has some capabilities such as specificity, diversity, and memory, so it can eradicate microbes effectively [1].

Viruses are one of the microorganisms that infect the immune system. A virus will be detected and then eradicated by innate immunity. However, there are several viruses that can escape from innate immunity and infect CD4+ T cells. Infected CD4+ T cells will produce cytokines to stimulate proliferation and differentiation of precursor CTL cells to become effector CTL cells that would eradicate infected CD4+ T cells [2].

Hematopoiesis is the process to derive blood cells that are located in bone marrow. The first stage of hematopoietic is a stem cell. Hematopoietic stem cells (HSCs) have an ability to multipotency and self-renewing. While in embryonic phase, HSCs that migrate to thymus would be differentiated to mature T cells [1]. Gene-based therapy has been discussed to cure immune impairing infections, such as HIV. The therapy’s concept is specifying technic and immune responses to face viruses by using genetically modified HSCs. In 2012, Kitchen et al. performed an in vivo experiment with HSCs gene-based therapy. Based on the experiment, engineered HSCs have the ability to establish functional antivirus responses that suppressed HIV replication. Hence, it will allow suppression of infected CD4+ T cells and prevent suppression of uninfected CD4+ T cells [3].

Korpusik [4] developed a mathematical model to study the influence of a constant influx of CTLs on the dynamic of the virus and immune system interactions, using a modification of the basic mathematical model of virus-induced impairment of help [2]. In addition, Korpusik and Kolev [5] developed a mathematical model to study the dynamical virus-immune system interactions after a single injection of CD8+ T cells derived from HSCs.

Ordinary differential equation (ODE) system that consists of first-order differential equations can be generalized to fractional differential equation (FDE) system that consists of fractional order differential equations α, with a fractional order parameter [6]. In most biological systems, FDE are naturally connected to systems with memory [7]. In addition, memory effect has an important role in the disease spread. The presence of memory effects on past events will affect the disease spread in the future. The distance of memory effect indicates the history of disease spread. Thus, memory effects on the spread of an infectious disease can be investigated using fractional derivatives [8–13].

In this paper, we extend the ODE model from Korpusik and Kolev [14] into a fractional differential equation system model. We also determined equilibria and stability of the equilibria from the proposed model. Finally, we perform numerical simulation to support the mathematical model interpretation.

2. Mathematical Model Formulation

In this section, we proposed a mathematical model of HSCs gene-based therapy based on [5]. The model was constructed under the following assumptions:

(1)The model consists of five compartments, namely, concentration of uninfected CD4+ T cells , concentration of infected CD4+ T cells , concentration of precursor CTL , concentration of effector CTL , and concentration of CD8+ T lymphocytes derived from HSCs .(2)Uninfected CD4+ T cells are produced with constant rate.(3)Viral infections occur only in a human body.(4)Free virus particle only infected uninfected CD4+ T cells.(5)Effector CTL only killed infected CD4+ T cells.(6)Injected HSCs will be differentiated into a precursor CTL.(7)The proliferation rate of precursor CTL cells depends on susceptible CD4+ T cells, infected CD4+ T cells, and CTL precursor cells at the current time.

The transmission diagram of the model is shown in Figure 1.

Figure 1: HSCs gene-based therapy transmission diagram.

The basic model from Korpusik and Kolev [5] was given as follows:

Here, all parameters , and . The description of the parameter for the model is presented in Table 1.

Table 1: The description of parameters in the model.

Parameter
Description
Unit

Production rate of uninfected CD4+ T cells
Concentration cells
/time unit

Death rate of uninfected CD4+ T cells
1/ time unit

Infection rate
Concentration cells⁻¹/ time unit

Death rate of infected CD4+ T cells
1/ time unit

Effector CTL rate of killing infected CD4+ T cells
Concentration cells⁻¹/ time unit

Fractions of HSCs that successfully differentiated into functional CTL cells
-

Proliferation rate of precursor CTL
Concentration cells⁻²/ time unit

Fraction of precursor CTL cells that differentiated into effector CTL cells
-

Death rate of precursor CTL
1/ time unit

Death rate of effector CTL
1/ time unit

Differentiation rate of CD8+ T lymphocytes
1/ time unit

The differential equation (1a) describes the dynamic of uninfected CD4+ T cells concentration. Concentrations of CD4+ T cells increase by production of bone marrow and decreased caused by natural death rate and infected by free virus particle. The differential equation (1b) shows the dynamic of infected CD4+ T cells concentrations. Concentrations of infected CD4+ T cells increase, caused by infection of free virus particle. Infected CD4+ T cells will be decreased because of natural death rate and it was killed by effector CTL cells.

The differential equation (1c) describes the dynamic of precursor CTL cells concentrations. Concentrations of precursor CTL cells increased by CTL proliferation process and HSCs were differentiated into CD8+ T cells that are parts of precursor CTL cells. Concentrations of precursor CTL cells are decreased, caused by differentiating into effector CTL, and decreased by natural death rate. Equation (1d) shows the dynamic of concentrations of effector CTL cells at time. Concentrations of effector CTL cells are increased, caused by differentiating of precursor CTL into effector CTL phase, and decreased, caused by natural death rate. The last equation (1e) shows the dynamic of CD8+ T cells concentrations that derived from HSCs. Concentrations of CD8+ T cells are decreased, caused by successfully passing the thymic selection and differentiating into functional precursor CTL cells.

Next, we consider a fractional order model of (1a)-(1e) above. The fractional model corresponding to system (1a)-(1e) is as follows:

where the fractional order derivative . Fractional derivative of the model (2) is adopted from Caputo’s definition. The main advantages of Caputo approach are the initial values for fractional differential equations with the Caputo derivatives taking on the same form as for integer order differential equations [15]. The Caputo fractional derivative is defined as follows.

Definition 1 (see [15]). Let , , and . Caputo fractional derivative , with fractional order , of function is defined by

where is the gamma function.

3. Stability Analysis

In this section, we study stability of the equilibriums of the fractional order model (2) above. We begin by computed the basic reproduction number of model (2). Basic reproduction number is defined as the average number of new cases of an infection caused by one typical infected individual, in a population consisting of susceptible only [16]. If , then the infections will die out, while if , then there is an epidemic case [17].

The stability theorem on fractional order system was given in the following theorem.

Theorem 2 (see [18]). Consider a nonlinear fractional order system

where , , and . The equilibrium points of system (4) are calculated by solving equation . The equilibrium points are locally asymptotically stable if all the eigenvalues of the Jacobian matrix evaluated at the equilibrium points satisfy the following condition:

The equilibriums of the fractional mathematical model in (2) satisfy the following equations:

The fractional model in (2) has three equilibriums, namely, virus-free equilibrium , CTL-Exhaustion Equilibrium , and control immune equilibrium . The virus-free equilibrium is condition when there is CD4+ T cells in human body ), no infected CD4+ T cells , and CTL precursor is not activated ). The virus-free equilibrium of model (2) above is given by .
Basic reproduction number is an important parameter in epidemiological cases. Basic reproduction number is defined by average secondary infections caused by one primary infection in susceptible population. In this paper, we use Next Generation Matrix (NGM) that developed by [19] to determine . By using the Next Generation method, we obtained the basic reproduction number . Stability of virus-free equilibrium is presented in Theorem 3.

Theorem 3. The virus-free equilibrium of model (2) is locally asymptotically stable if and only if .

Proof. The Jacobi matrix that evaluated at the virus-free equilibrium of the model (2) above is as follows:

The eigenvalues of the Jacobi matrix are , , , , and . Hence, we have . Meanwhile, if and only if If the condition is satisfied, then it is clear that all the eigenvalues satisfy the condition , for . Hence, the virus-free equilibrium of model (2) is locally asymptotically stable if .

We continue with the stability analysis of the CTL-Exhaustion Equilibrium of the model (2). The CTL-Exhaustion Equilibrium is condition when there is CD4+ T cells in human body ), there are infected CD4+ T cells , but CTL precursor is not activated yet . The CTL-Exhaustion Equilibrium of the model in (2) is given by . This equilibrium will exist whenever .

Theorem 4. Let and . The CTL-Exhaustion Equilibrium of the model in (2) is locally asymptotically stable if and only if .

Proof. The Jacobi matrix evaluated at the free virus equilibrium of the model in (2) above is as follows:

The eigenvalues of the Jacobi matrix above are , , , and is root of polynomial , where , and . We have . Then we determined the argument of and . Based on the condition of existence of that , we obtained . By using Routh-Hurwitz criterion, we obtained .
Since all parameters are assumed to have positive value, the third eigenvalue is a real number. We obtained that if and only if the condition is fulfilled. It is clear that all the eigenvalues satisfy the condition , for . Hence, the CTL-Exhaustion Equilibrium of model (2) is locally asymptotically stable if and only if or . This completes the proof

Last, we analyze the stability of the control immune equilibrium . The control immune equilibrium is the condition when there are CD4+ T cells in human body , there is infected CD4+ T cells , and CTL precursor is activated so that can give the immune response. The control immune equilibrium is given by , where

The equilibriums will exist if the conditions , and are satisfied.

The Jacobi matrix of the model in (2) that evaluated in the equilibrium point is

By using the control immune equilibrium condition, the Jacobi matrix could be simplified into

From the Jacobi matrix above, we get the polynomial characteristic equations as follows:

where , , and .

Based on the characteristic equations (16) above, we obtain eigen values , , and is the roots of the following equation:

The argument of the first and second eigen values is , . Next, we will determine the argument of the third, fourth, and fifth eigen values. The third, fourth, and fifth eigen values are roots of polynomial with degree of three (17), satisfying the following conditions:

(i)(ii)(iii)

Then, we will analyze two conditions when and . If , then we have Theorem 5. On the other hand, if then we obtain Theorem 6 below.

Theorem 5. If , then the immune control equilibrium of model (2) is unstable.

Proof. Let ; we obtained . Since , then based on the third condition (iii) above we have . There are two possibilities, either all the eigen values of the polynomial (17) are real number or one of the eigen values is real number and two other eigen values are complex numbers.

(1)Let all the eigen values be real number. Because of , then there are two possibilities. First, all the eigen values are real positive number. If all the eigen values are real positive numbers, then it is clear that . Second, one of the eigen values is real positive number and two others are real negative number. Let , then .(2)Let one of the eigen values be real number, , and two others eigen values are complex numbers, , , where . Because and , then we obtained . As a result, .

From the two conditions above, we find that the immune control equilibrium of the model (2) is unstable whenever .

Theorem 6. The immune control equilibrium of the model (2) is asymptotically stable if for some fractional order .

Proof. Let we consider the polynomial

where , , and .
Let . Hence we obtained . If , then based on the third condition (iii) above we have .

(i)Let the polynomial characteristic (18) above have complex number roots. Let one of the eigen values be real number, , and two other eigen values are complex number, , , where . Since and , then we obtained . As a result, . In addition, we can found a fractional order value such that the arguments of the third and fourth eigen values satisfy and . As a result, the immune control equilibrium of model (2) is asymptotically stable for some fractional order (ii)Let all the eigen values of the polynomial characteristic (18) above , , be real numbers. Then we will determine the conditions such that all roots of the polynomial (18) above are either negative or complex roots with negative real parts by using Routh-Hurwitz criterion. Based on Routh-Hurwitz criteria, the polynomial (18) will have real negative or real negative parts roots if it satisfies , , and . Hence, we obtained the following conditions:

(a)Since all parameters are assumed to be positive and , then it is clear that .(b)The coefficient will have positive value if is positive.(c)The coefficient will have positive value if it satisfies .(d)It is clear that the coefficients satisfy .

Based on [20], for these conditions are sufficient but not necessary. Therefore, the argument of will satisfy if the condition is satisfied. As a result, the immune control equilibrium of the model (2) is asymptotically stable. This completes the proof.

4. Sensitivity Analysis

In this section, we will analyze the sensitivity of parameters of the model (2) above. Sensitivity analysis is used to determine the relative importance of model parameters to disease transmission and prevalence based on the sensitivity index of basic reproduction number . The calculation of sensitivity index of follows the approach in Chitnis [21].

Definition 7 (see [21]). The normalized forward sensitivity index of a variable, , that depends differentially on a parameter, , is defined as

Based on Definition 7, the sensitivity indices of with respect to each parameter such as , and can be computed in the same way as (19). The sensitivity indices of are given in Table 2. Based on Table 2, it can be seen that parameters and positively affect the rate of model changed. Respectively, the parameters and negatively affect the rate of model changed. But, we did not know what is the relative importance of model parameters of . So, we also analyze the sensitivity parameters of .

Table 2: The sensitivity indices of .

Parameter
Sensitivity indices

1

1

-1

-1

Based on Definition 7, the sensitivity indices of with respect to each parameter such as , and can be computed in the same way as (19). The sensitivity indices of parameter are ; we next substituted the value of parameters of Table 4 so we obtain . The sensitivity indices of are given in Table 3.

Table 3: The sensitivity indices of .

Parameter
Sensitivity indices

0.945

1

-0.994

-1

-0.049

-0.006

0.055

Table 4: Parameters value of the model.

Parameter
Description
Value
Source

Production rate of uninfected CD4+ T cells
0.2
[14]

Death rate of uninfected CD4+ T cells
0.05
Assumption

Infection rate
0.15
[14].

Death rate of infected CD4+ T cells
0.8
[14]

Effector CTL rate of killing infected CD4+ T cells
0.016
[14]

Fractions of HSCs that successfully differentiated into functional CTL cells
0.01
Assumption

Proliferation rate of precursor CTL
0.1
Assumption

Fraction of precursor CTL cells that differentiated into effector CTL cells
0.1
Assumption

Death rate of precursor CTL
0.05
[14]

Death rate of effector CTL
0.05
[14]

Differentiation rate of CD8+ T cells
0.005
Assumption

Based on Table 3, it can be seen that if production rate of uninfected CD4+ T cells is increased (decreased) about 10%, then will increase (decrease) about 10%. Respectively, if death rate of infected CD4+ T cells is increased (decreased) about 10%, then will decrease (increase) about 9.94%, and the same for other parameters.

5. Numerical Simulation

In this section, we present numerical simulations to show the dynamics of the model at free virus condition, CTL-exhaustion condition, and control immune condition using MATLAB R2009a. The initial conditions of simulations of free virus condition are , and . The parameter values are presented in Table 4. Based on Table 4, we have that basic reproduction number R₀ is that shows free virus condition or there is no spread of viral infection. This simulation used variations of fractional order and time interval for 2500 time units.

The result of simulation of free virus condition can be seen in Figure 2. Based on Figure 2, greater the fractional order leads to the increase of concentrations of uninfected CD4+ T cells but concentrations of infected CD4+ T cells are decreasing that show free virus conditions. Besides, concentrations of precursor CTL cells are decreasing caused by the decrease of CD4+ T cells. That also implies the decrease of effector CTL. Meanwhile, concentrations of CD8+ T cells are decreasing because they were differentiated into precursor CTL cells.

Figure 2: The dynamics of mathematical model of HSCs gene-based therapy for free virus conditions.

Next, we will interpret the simulation of CTL-exhaustion condition. The initial conditions of simulations of CTL-exhaustion condition are , and . The parameter values are presented in Table 4 except for , and , which are , , and . Therefore, we have basic reproduction number and which is condition when there is spread of viral infections in the human body. This simulation used variations of fractional order and time interval for 2500 time units.

Based on Figure 3, it can be seen that greater the fractional order leads to the decrease of concentrations of uninfected CD4+ T cells but concentrations of infected CD4+ T cells are increasing. That implies the conditions of viral infections. But, concentrations of effector CTL cells are decreasing. That implies the CTL-exhaustion conditions where CTL are not activated yet so they cannot against the viral infections. Meanwhile, concentrations of CD8+ T cells were decreased because they differentiated into precursor CTL cells.

Figure 3: The dynamics of mathematical model of HSCs gene-based therapy for CTL-exhaustion conditions.

Last, we will show the simulation of control immune conditions. The initial conditions of simulations of CTL-exhaustion condition are , and . The parameter values are showed in Table 4 except for , that is, . Therefore, we have basic reproduction number and which is condition when there is spread of viral infections in the human body. This simulation used variations of fractional order and time interval for 2500 time units. Based on Figure 4, it can be seen that greater the fractional order leads to the increase of concentrations of uninfected CD4+ T cells but concentrations of infected CD4+ T cells are decreasing. Meanwhile, concentrations of precursor and effector CTL cells are increasing, which means the precursor and effector CTL cells are activated so they can fight against the viral infections. That implies control immune conditions. Meanwhile, concentrations of CD8+ T cells were decreased because they differentiated into precursor CTL cells.

Figure 4: The dynamics of fractional order mathematical model of HSCs gene-based therapy for control immune conditions.

6. Conclusions

We were discussed about the stability analysis of fractional order mathematical model of HSCs gene-based therapy. The model has three equilibriums points. The existence and stability of the equilibriums were achieved. Based on the numerical simulations, we conclude that HSCs gene-based therapy has potential to reduce the viral infections because it can decrease the concentration of infection cells and increase the concentration of the immune cells.

Data Availability

No data were used to support this study.

Conflicts of Interest

The authors state that there are no conflicts of interest concerning the preparation and publication of the manuscript.

References

R. Nairn and M. Helbert, Immunology for Medical Students, Mosby Elsevier, Philadelphia, 2nd edition, 2007.

D. Wodarz, Killer cell dynamics, vol. 32 of Interdisciplinary Applied Mathematics, Springer-Verlag, New York, 2007.

S. G. Kitchen, B. R. Levin, G. Bristol et al., “In vivo suppression of HIV by antigen specific T cells derived from engineered hematopoietic stem cells,” PLoS Pathogens, vol. 8, no. 4, Article ID e1002649, 2012.

A. Korpusik, “Hematopoietic stem cell based therapy of immunosuppressive viral infection - Numerical simulations,” Biocybernetics and Biomedical Engineering, vol. 34, no. 2, pp. 125–131, 2014.

A. Korpusik and M. Kolev, “Single injection of CD8+ T lymphocytes derived from hematopoietic stem cells - Mathematical and numerical insights,” BioSystems, vol. 144, pp. 46–54, 2016.

S. Das and P. K. Gupta, “A mathematical model on fractional Lotka-Volterra equations,” Journal of Theoretical Biology, vol. 277, pp. 1–6, 2011.

E. Ahmed and A. S. Elgazzar, “On fractional order differential equations model for nonlocal epidemics,” Physica A: Statistical Mechanics and its Applications, vol. 379, no. 2, pp. 607–614, 2007.

F. Fatmawati, E. M. Shaiful, and M. I. Utoyo, “A Fractional-Order Model for HIV Dynamics in a Two-Sex Population,” International Journal of Mathematics and Mathematical Sciences, vol. 2018, Article ID 6801475, 11 pages, 2018.

J. Huo and H. Zhao, “Dynamical analysis of a fractional {SIR} model with birth and death on heterogeneous complex networks,” Physica A: Statistical Mechanics and its Applications, vol. 448, pp. 41–56, 2016.

J. Huo, H. Zhao, and L. Zhu, “The effect of vaccines on backward bifurcation in a fractional order HIV model,” Nonlinear Analysis: Real World Applications, vol. 26, pp. 289–305, 2015.

C. M. A. Pinto and A. R. M. Carvalho, “The HIV/TB coinfection severity in the presence of TB multi-drug resistant strains,” Ecological Complexity, vol. 32, pp. 1–20, 2017.

M. Saeedian, M. Khalighi, N. Azimi-Tafreshi, G. R. Jafari, and M. Ausloos, “Memory effects on epidemic evolution: the susceptible-infected-recovered epidemic model,” Physical Review E: Statistical, Nonlinear, and Soft Matter Physics, vol. 95, no. 2, 022409, 9 pages, 2017.

T. Sardar, S. Rana, and J. Chattopadhyay, “A mathematical model of dengue transmission with memory,” Communications in Nonlinear Science and Numerical Simulation, vol. 22, no. 1–3, pp. 511–525, 2015.

G. Huang, Y. Takeuchi, and A. Korobeinikov, “HIV evolution and progression of the infection to AIDS,” Journal of Theoretical Biology, vol. 307, pp. 149–159, 2012.

S. Z. Rida and A. A. Arafa, “New method for solving linear fractional differential equations,” International Journal of Differential Equations, Art. ID 814132, 8 pages, 2011.

O. Diekmann, J. A. P. Heesterbeek, and M. G. Roberts, “The construction of next-generation matrices for compartmental epidemic models,” Journal of the Royal Society Interface, vol. 7, no. 47, pp. 873–885, 2010.

F. Brauer and C. Castillo-Chávez, Mathematical Models in Population Biology and Epidemiology, vol. 40 of Texts in applied mathematics, Springer, New York, NY, USA, 2001.

M. S. Tavazoei and M. Haeri, “Chaotic attractors in incommensurate fractional order systems,” Physica D: Nonlinear Phenomena, vol. 237, no. 20, pp. 2628–2637, 2008.

P. van den Driessche and J. Watmough, “Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission,” Mathematical Biosciences, vol. 180, pp. 29–48, 2002.

E. Ahmed, A. M. A. El-Sayed, and H. A. A. El-Saka, “On some Routh-Hurwitz conditions for fractional order differential equations and their applications in Lorenz, Rössler, Chua and Chen systems,” Physics Letters A, vol. 358, no. 1, pp. 1–4, 2006.

N. Chitnis, J. M. Hyman, and J. M. Cushing, “Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model,” Bulletin of Mathematical Biology, vol. 70, no. 5, pp. 1272–1296, 2008.