Fractional-order dynamics of Chagas-HIV epidemic model with different fractional operators

1.   Introduction

American trypanosomiasis, commonly known as Chagas disease, is a fatal disease transmitted by biting insects called "kissing bugs". The Trypanosoma cruzi parasite is the cause of Chagas disease. People who have Chagas are unaware of the disease for many years, as there are no symptoms. After several years, some of the people who have Chagas face heart damage, which consequently leads to sudden death. In Latin America, no other parasitic disease kills as many people as Chagas does. Almost 6 million people have been infected from Chagas, and per year, 173,000 new cases appear. In Latin America, it is endemic in 21 countries, and 9,490 people have been killed in 2019. Unfortunately, countries in which Chagas is not endemic have recently been affected due to the migration of people from the endemic regions^[1,2], which results in the contact of Chagas-infected migrants with HIV-infected people, and hence a co-infection occurs ^[3,4]. It has been observed that vertical transplacental transmission, contaminated blood products, and blood transfusions are the transmission means of both diseases. Intravenous drug users are a high prevalence group for Chagas, when they use contaminated needles. The prevalence of co-infection is significantly higher in this group ^[5].

A serious situation develops in patients when a co-infection of Chagas disease with HIV infection occurs, which results in high fatality rates^[6,7]. It has been reported that after diagnosis, it takes 10–20 days to fatal evolution. Myocarditis, meningoencephalitis and cutaneous lesions are the clinical manifestations which may present in coinfected patients ^[8,9,10]. Early diagnosis with highly active antiretroviral treatment can only improve the survival. In non-endemic regions, the professionals rarely suspect the HIV-positive patients of having Chagas.

While a number of clinical studies have been reported on the Chagas epidemic in the setting of advanced AIDS ^{[11,12,13,14]}, not much has been analyzed in the setting of the Chagas epidemic and HIV co-infection despite its contingency and clinical consequence. A number of mathematical models have also been studied to simulate the spread and control of either the Chagas epidemic or advanced AIDS ^{[15,16,17,18,19]}.

In this paper, we use fractional derivatives to simulate the proposed disease structure. Fractional calculus is a modification of classical calculus. To simulate a mathematical model of any proposed phenomena related to a real-world problem, the fractional-order derivatives are most suitable because they provide a higher degree of accuracy and are the best fit to capture the memory effects and spanning nature. Epidemic models defined in the sense of fractional derivatives provide more information about the disease dynamics as compared to the integer-order models ^[20,21]. Various types of fractional derivatives in the forms of different kernel properties have been proposed by mathematicians ^[23,24] and have been used to simulate different problems ^[25,26,27]. For example, in ^[28], some researchers proposed a malaria model with Caputo-Fabrizio and Atangana-Baleanu derivatives. Some recent studies on mathematical modeling of infectious diseases using fractional derivatives are ^{[29,30,31,32,33]}.

The organization of our paper is as follows: Section 2 deals with the feasible region, basic reproduction number and the equilibrium point of the proposed Chagas and HIV co-infection model. Section 3 is concerned with the local stability of the proposed model. In Section 4, the deterministic model is then fractionalized by using the Atangana-Baleanu fractional derivative in the Caputo sense. Also, some preliminaries are available regarding different fractional derivatives in the said section. Existence and uniqueness are determined in Section 5. Section 6 deals with the Ulam-Hyers stability of our model. Section 7 depicts some simulations carried out by a new numerical scheme proposed by Atangana and Seda. Section 8 is devoted to justifying the concluding remarks along with the further scope of the study.

2.   Model formulation

In the model, we split the total human population size $\mathbb{N}$ into five different classes: Chagas susceptible class $\mathbb{S}_C$, class of HIV susceptible $\mathbb{S}_H$, Chagas infectious class $\mathbb{I}_C$, class of HIV infectious $\mathbb{I}_H$ and the class of infected population with both Chagas and HIV $\mathbb{I}_{CH}$. So, the complete human population size is specified by

Thus, the mathematical model for defining the proposed structure of Chagas and/or HIV infections is derived by the following system ^[34]:

The parameters and the classes used in the above proposed model are specified in Tables 1 and 2.

2.1.   Basic reproduction number $\mathbb{R}_0$

The region

is positively invariant for (2.2), and all solutions of $(\mathbb{S}_C, \mathbb{S}_H, \mathbb{I}_C, \mathbb{I}_H, \mathbb{I}_{CH}) \in \mathbb{R^5_+}$ remain in $\Delta$ for all $t > 0$.

The disease-free equilibrium (DFE) denoted by $E^0$ is

Let $(\mathbb{I}_C, \mathbb{I}_H, \mathbb{I}_{CH})$ be our infected compartment, and then it follows from system (2.2) that

the Jacobian matrix of the model is

Arranging $J$ such that $J = F-V$, we get

where the elements in the matrix F constitute the new infection terms

The matrix $V$ represents the exchange of infection from one compartment to another. Therefore, the next generation matrix defined by $FV^{-1}$ is

Thus, $\mathbb{R}_0$, which is the dominant eigenvalue of matrix $FV^{-1}$, is obtained as

where $\mathbb{R}_{0C}$ and $\mathbb{R}_{0H}$ are the reproduction numbers for Chagas and HIV, respectively.

3.   Local stability

Theorem 3.1. The disease-free equilibrium $E^0$ of the system (2.2) is locally asymptotically stable if $\mathbb{R}_0 < 1$.

Proof. The Jacobian matrix of the system (2.2) at $E^0$ is given by

where $a_{11} = (\mu+\delta_1+\rho_C-\gamma_1 \frac{\alpha \sigma}{\mu})$ and $a_{22} = (\mu+\delta_2+\gamma_4+\rho_H-\gamma_2 \frac{(1-\alpha)\sigma}{\mu})$.

The characteristic equation of $J^{\mid0\mid}$ is given by

The three eigenvalues of the characteristic equation of $J^{\mid0\mid}$ are negative, i.e., $\lambda_1 = -\mu$, $\lambda_2 = -\mu,$ and $\lambda_5 = -(\mu+\delta_3+\dot{\theta})$.

The eigen value $\lambda_3 = -(\mu+\delta_1+\rho_C-\gamma_1 \frac{\alpha \sigma}{\mu})$ is negative if $\mu+\delta_1+\rho_C-\gamma_1 \frac{\alpha \sigma}{\mu} > 0$, that is

and by definition

Similarly, $\lambda_4 = -(\mu+\delta_2+\gamma_4+\rho_H-\gamma_2 \frac{(1-\alpha)\sigma}{\mu})$ is negative if $(\mu+\delta_2+\gamma_4+\rho_H-\gamma_2 \frac{(1-\alpha)\sigma}{\mu}) > 0$, that is

and by definition

So, $\mathbb{R}_0 = max(\mathbb{R}_{0C}, \mathbb{R}_{0H})$ implies that $\mathbb{R}_0 < 1$.

This shows that the DFE point is asymptotically stable if $\mathbb{R}_0 < 1$.

4.   Preliminaries

Here, we present some definitions of integral and differential operators, starting with the Caputo fractional derivative

the Caputo-Fabrizio fractional derivative

and the Atangana-Baleanu-Caputo (ABC) fractional derivative

The fractal-fractional derivatives with power-law kernel, exponential decay kernel, and Mittag-Leffler kernel are given by

where

The fractal-fractional integrals with power-law, exponential decay and Mittag-Leffler kernels, respectively, are given below:

4.1.   Fractional order model

To the best of our knowledge, no one has yet considered the fractional-order Chagas-HIV epidemic model in the sense of the ABC fractional derivative. Therefore, getting motivation from the above-given model (2.2), we consider a fractional-order Chagas-HIV model using the ABC fractional derivative ^[26], which is as follows:

with the initial sizes of given classes

5.   Existence and uniqueness of solution for the ABC model

First, we derive the existence and uniqueness of the solution with respect to the Atangana-Baleanu-Caputo derivative for the system (4.7). Consider a continuous real-valued function denoted by $B(J)$ associated to the supremum-norm characteristic, is a Banach space on $J = [0, b]$ and $P = B(J) \times B(J) \times B(J) \times B(J) \times B(J)$ with norm $\left\|\left(\mathbb{S}_{C}, \mathbb{S}_{H}, \mathbb{I}_{C}, \mathbb{I}_{H}, \mathbb{I}_{CH} \right)\right\| = \|\mathbb{S}_{C}\|+\|\mathbb{S}_{H}\|+\left\|\mathbb{I}_{C}\right\|+\left\|\mathbb{I}_{H}\right\|+\|\mathbb{I}_{CH}\|$, where $\|\mathbb{S}_{C}\| = \sup _{t \in J}|\mathbb{S}_{C}(t)|$, $\|\mathbb{S}_{H}\| = \sup _{t \in j}|\mathbb{S}_{H}(t)|$, $\left\|\mathbb{I}_{C}\right\| = \sup _{t \in j}|\mathbb{I}_{C}(t)|$, $\left\|\mathbb{I}_{H}\right\| = \sup _{t \in j}|\mathbb{I}_{H}(t)|$, $\left\|\mathbb{I}_{CH}\right\| = \sup _{t \in j}|\mathbb{I}_{CH}(t)|$. Using the Atangana-Baleanu-Caputo fractional integral operator on both sides of Eq (4.7), we get

Now, the definition (4.3) leads us to

where

The symbols $\mathfrak{K}_{1}$, $\mathfrak{K}_{2}$, $\mathfrak{K}_{3}$, $\mathfrak{K}_{4}$ and $\mathfrak{K}_{5}$ have to hold for the Lipschitz condition only if $\mathbb{S}_{C}(t)$, $\mathbb{S}_{H}(t)$, $\mathbb{I}_{C}(t)$, $\mathbb{I}_{H}(t)$ and $\mathbb{I}_{CH}(t)$ possess an upper bound. Taking that $\mathbb{S}_{C}(t)$ and $\mathbb{S}_{C}^*(t)$ are couple functions, we get

Taking into account $\eta_{1}: = \left\|-\left(\left(\gamma_1(\mathbb{I}_C + \pi_1 \mathbb{I}_{CH})\right)+\mu \right)\right\|$, one reaches

Also, we can get

where

which shows that the Lipschitz condition holds. Continuing in a recursive manner, the expressions in (5.2) yield

Together with $\mathbb{S}_{C0}(t) = \mathbb{S}_{C}(0)$, $\mathbb{S}_{H0}(t) = \mathbb{S}_{H}(0)$, $\mathbb{I}_{C0}(t) = \mathbb{I}_{C}(0)$, $\mathbb{I}_{H0}(t) = \mathbb{I}_{H}(0)$ and $\mathbb{I}_{CH0}(t) = \mathbb{I}_{CH}(0)$, differences of consecutive terms yield

It is vital to observe that $\mathbb{S}_{Cn}(t) = \sum_{i = 0}^{n} \mathbb{I}_{\mathbb{S}_{C}, i}(t)$, $\mathbb{S}_{Hn}(t) = \sum_{i = 0}^{n} \mathbb{I}_{\mathbb{S}_{H}, i}(t)$, $\mathbb{I}_{Cn}(t) = \sum_{i = 0}^{n} \mathbb{I}_{\mathbb{I}_{C}, i}(t)$, $\mathbb{I}_{{H}_{n}}(t) = \sum_{i = 0}^{n} \mathbb{I}_{\mathbb{I}_{H}, i}(t)$, $\mathbb{I}_{CHn}(t) = \sum_{i = 0}^{n} \mathbb{I}_{\mathbb{I}_{CH}, i}(t)$. Additionally, by using Eqs (5.5), (5.6) and considering that $\mathbb{I}_{\mathbb{S}_{C}, n-1}(t) = \mathbb{S}_{Cn-1}(t)-\mathbb{S}_{Cn-2}(t)$, $\mathbb{I}_{\mathbb{S}_{H}, n-1}(t) = \mathbb{S}_{Hn-1}(t)-\mathbb{S}_{Hn-2}(t)$, $\mathbb{I}_{\mathbb{I}_{C}, n-1}(t) = \mathbb{I}_{Cn-1}(t)-\mathbb{I}_{Cn-2}(t)$, $\mathbb{I}_{\mathbb{I}_{H}, n-1}(t) = \mathbb{I}_{Hn-1}(t)-\mathbb{I}_{Hn-2}(t)$, $\mathbb{I}_{\mathbb{I}_{CH}, n-1}(t) = \mathbb{I}_{CHn-1}(t)-\mathbb{I}_{CHn-2}(t)$, we reach

Theorem 5.1. Assume that the following condition holds:

Then, (4.7) has a unique solution for $t \in [0, b]$.

Proof. It is given that $\mathbb{S}_{C}(t)$, $\mathbb{S}_{H}(t)$, $\mathbb{I}_{C}(t)$, $\mathbb{I}_{H}(t)$ and $\mathbb{I}_{CH}(t)$ are bounded functions. In advance, as can be observed from Eqs (5.5) and (5.6), the symbols $\mathfrak{K}_{1}$, $\mathfrak{K}_{2}$, $\mathfrak{K}_{3}$, $\mathfrak{K}_{4}$ and $\mathfrak{K}_{5}$ hold for the Lipschitz condition. Therefore, using Eq (5.9), together with a recursive hypothesis, we derive

Thus, we can observe that sequences exist that satisfy $\left\|\mathbb{I}_{\mathbb{S}_{C}, n}(t)\right\| \rightarrow 0$, $\left\|\mathbb{I}_{\mathbb{S}_{H}, n}(t)\right\| \rightarrow 0$, $\left\|\mathbb{I}_{\mathbb{I}_{C}, n}(t)\right\| \rightarrow 0$, $\left\|\mathbb{I}_{\mathbb{I}_{H}, n}(t)\right\| \rightarrow 0$, $\left\|\mathbb{I}_{\mathbb{I}_{CH}, n}(t)\right\| \rightarrow 0,$ as $n \rightarrow \infty$. Moreover, from Eq (5.11) and applying the triangle inequality, for any $k$, we have

with $Z_{i} = \frac{1-\Xi}{B(\Xi)} \eta_{i} + \frac{\Xi}{B(\Xi) \Gamma(\Xi)} b^{\Xi} \eta_{i} < 1$ by hypothesis. Similarly, we can prove the existence of a unique solution for the proposed model in terms of other fractional derivatives.

6.   Hyers-Ulam stability

Definition 6.1. ^[23] The ABC fractional integral model proposed by Eq (5.2) is called Hyers-Ulam stable if there exist constants $\zeta_i > 0, \; i \in \textbf{N}^5$ satisfying the following: For every $\gamma_i > 0, \; i \in \textbf{N}^5$, when

there exist $(\dot{\mathbb{S}_{C}}(t), \dot{\mathbb{S}_{H}}(t), \dot{\mathbb{I}_{C}}(t), \dot{\mathbb{I}_{H}}(t), \dot{\mathbb{I}_{CH}}(t))$ satisfying

such that

Theorem 6.1. With the definition of $J$, the given model of fractional order (5.1) is Hyers-Ulam stable.

Proof. Using Theorem 5.1, the given ABC fractional system (5.1) contains a unique solution $({\mathbb{S}_{C}}(t), {\mathbb{S}_{H}}(t)$, ${\mathbb{I}_{C}}(t), {\mathbb{I}_{H}}(t), {\mathbb{I}_{CH}}(t))$ satisfying the equations of system (5.2). Then, we have

Taking $\gamma_i = \mathbb{I}_i$, $\Delta_i = \frac{1-\Xi}{B(\Xi)}+\frac{\Xi}{B(\Xi) \Gamma(\Xi)}$, this implies

Following the same procedure, we have

From the results of Eqs (6.10) and (6.11), the AB fractional integral model (5.2) is Hyers-Ulam stable, and consequently the AB-fractional order model (5.1) is Hyers-Ulam stable. This ends the proof.

7.   Numerical scheme

Now, we derive a numerical scheme for our model. We shall start with the Caputo-Fabrizio fractional derivative, and this will be followed by the Caputo and Atangana-Baleanu fractional derivatives. Finally, we will solve the models with fractal-fractional derivatives. So, the Caputo-Fabrizio model is given by

For simplicity, we write the above equation as follows:

After applying the fractional integral with exponential kernel and putting Newton polynomials into these equations, we can solve our model as follows:

We have the following numerical scheme for the Mittag-Leffler case:

where

Finally, we have the following numerical approximation with the Caputo derivative:

Now, we consider our model with fractal-fractional operators. We start with the Caputo-Fabrizio fractal-fractional derivative

After applying the fractal-fractional integral with exponential kernel, we have the following scheme for this model:

For the Mittag-Leffler kernel, we have the following numerical scheme:

For the power-law kernel, we can get the following numerical scheme:

7.1.   Numerical results and discussion

This portion is devoted to the numerical simulation results based on the above numerical schemes for the proposed HIV and Chagas disease model in the given fractional derivative senses. The numerical simulations are performed using the following parameter values: $\alpha = 0.425$; $\sigma = 0.4898$; $\mu = 0.0012$; $\gamma_1 = 0.03$; $\gamma_2 = 0.09$; $\gamma_3 = 0.05$; $\gamma_4 = 0.004$; $\rho_c = 0.2$; $\rho_h = 0.05$; $\pi_1 = 0.002$; $\pi_2 = 0.005$; $\delta_1 = 0.01$; $\delta_2 = 0.1$; $\delta_3 = 0.05$; $\theta_h = 0.003$; $\theta_d = 0.076$. Variations in the given model classes at various fractional-order values can be observed from the family of Figures 1–6.

8.   Conclusions

In this research, we have reformulated Chagas and HIV co-infection models by using various fractional derivatives. The basic reproductive number $R_0$ of the aforementioned model is established along with the feasible region and disease-free equilibrium point $E^0$. It has been proven that the $E^0$ is locally asymptotically stable when $R_0$ is less than one. The model is then fractionalized by using the Atangana-Baleanu fractional derivative in the Caputo sense. The existence and uniqueness of the solution along with Ulam-Hyers stability have been established. Finally, the model has been solved by a new numerical scheme, which is generated using Newton polynomials and is considered more accurate than one using Lagrange polynomials which are used to generate in the Adam Bashforth method. Furthermore, the same model is numerically solved by taking some other fractional derivatives, like Caputo-Fabrizio, Caputo, and fractal-fractional with power-law, exponential-decay and Mittag-Leffler kernels.

In conclusion, the epidemiological profiles of AIDS and Chagas disease over the past years might facilitate the approximation of both diseases, thereby increasing the possibility of reactivation of Chagas disease in HIV patients. As patients are diagnosed after several years and some time in the last stages, its lethality rate is very high. Such cases have a big risk of not being correctly diagnosed, therefore adding to the lethality and seriousness of the epidemic. Serological tests for Chagas, even when they are indeterminate, should be taken fully into account, and the relevant parasitological tests must be done.

Acknowledgements

This research project was supported by Thailand Science Research and Innovation (TSRI), Basic Research Fund, fiscal year 2022, under project number FRB650048/0164.

Conflict of interest

The authors declare no conflicts of interest.