Modelling the dynamics of Trypanosoma rangeli and triatomine bug with logistic growth of vector and systemic transmission

1.   Introduction

Chagas disease, known as American trypanosomiasis, is a protozoan parasitic disease caused by Trypanosoma cruzi (T. cruzi). The disease was discovered firstly by Doctor Chagas in 1908 and it was named after that. Chagas disease is an illness that can cause serious consequences including heart disease and cardiomyopathy, and many people infected with Chagas disease may die due to these complications ^[1,2]. It is mainly prevalent in Central and South America, such as Argentina, Bolivia, Brazil, Chile, etc. About 13% of the Latin American population is at risk of T. cruzi infection ^[3]. Due to convenient transportation and the globalization, Chagas disease is spreading very widely in the world, see Figure 1(a), (b) for details. An estimated 8 million people are infected with Trypanosoma cruzi worldwide, mainly in Latin America where Chagas disease remains one of the biggest public health problems, causing incapacity in infected individuals and more than 10,000 deaths per year ^[4,5,6]. In particular, patients with Chagas disease may be coinfected with other epidemic diseases including HIV ^[7] and COVID-19. These patients are at risk of severe COVID-19 manifestations and should be a priority group to be vaccinated ^[8].

Trypanosoma cruzi, a protozoan parasite that parasitizes human and mammalian blood and tissue cells, can be transmitted by blood-sucking triatomine bugs to cause the symptoms of Chagas disease. It is spread mainly through the faeces of the infected blood-sucking triatomine bugs. These bugs usually live in the crevices of poorly built houses in rural or suburban areas. They hide during the day and come out at night to feed on human blood. They bite exposed areas of the skin, such as the face, and defecate near the bites. If one scratches on the bites, this leads to feces spreading to the sites of eye, mouth, or any skin break, and then the parasites enter the body, and eventually go into the heart, survive and proliferate inside ^[3,9]. Triatomine bugs, the vectors to transmit Chagas disease, have a relatively short life span, ranging from 4 to 14 months depending largely on species and environmental conditions. They suck the blood of vertebrates, especially mammals (such as dogs, bats, armadillos, squirrels, guinea pigs and humans), and then release feces on the skin of the bitten animal. The feeding time may take 10–30 minutes ^[10]. Once healthy triatomine bugs ingest with the parasites of Chagas disease, they will be infected quickly ^[11,12,13].

Trypanosoma rangeli (T. rangeli) is a kind of parasite which is pathogenic to some vector species including triatomine bug, it always influence the transmission dynamics of the infected bugs, so this kind of transmission behavior deserves further research. Rhodnius prolixus (R. prolixus) is one of the triatomine species and T. rangeli is one of the trypanosoma which is spread between hosts and triatomine vectors. Although T. rangeli can infect mammals through the same triatomines, it is not pathogenic to human. However, it is still important to study the transmission dynamics of T. rangeli because it shares soluble antigenic epitopes with T. cruzi and the crossed serological reactions affect the diagnosis of Chagas disease ^[14,15]. Both T. rangeli and T. cruzi have common hosts and triatomine vectors. The transmission dynamics of T. rangeli between hosts and triatomine bugs can affect the effectiveness of T. cruzi transmission ^[10]. Moreover, T. rangeli has pathogenic effect on triatomine bugs in the sense that the infection of T. rangeli can change the behavior of triatomine bugs, which alters the transmission of Chagas disease ^[14,16]. There are some studies showing the possible interaction between T. rangeli and T. cruzi ^[17]. Therefore, it is essential to study the transmission behavior of T. rangeli. T. rangeli's infection pattern is similar to that of T. cruzi. The healthy population will get infected if they bite the infected counterparts through systemic transmission, where T. rangeli parasite can enter and multiply at the common hosts' and vectors' bloodstream. In addition to this normal insect-host-insect transmission mode, there is another mode of T. rangeli transmission called insect-to-insect co-feeding transmission. Susceptible triatomine bugs can get infected if they are feeding with infected counterparts on the same hosts. This has been studied in some papers ^{[10,18,19,20,21,22]} and it differs from the transmission of T. cruzi.

There are a number of mathematical models that studied the transmission dynamics of Chagas disease, such as the different transmission routes of the interaction between hosts and vectors ^[2,23,24], the disease transmission in the host movement and host community composition ^{[9,16,25,26,27]}, the triatomine population with temporal or spatial variations ^{[28,29,30,31,32,33]}, and the optimal control of Chagas disease ^{[34,35,36,37,38]}. Recently, Wu et al. ^[10] formulated a new model by considering the Ricker's type growth of triatomine bugs and T. rangeli's pathogenic effect on triatomine bugs. However, the logistic growth of triatomine bugs is very common but was not investigated in the model. In this paper, we assume the generation rate of triatomine bugs follows the logistic growth instead of Ricker's type function and further study the dynamical behavior of triatomine-rangeli-host transmission. The new model is shown to have interesting dynamics, provide more insights into the interaction between triatomine bugs and T. rangeli, and may help to prevent and control the Chagas disease.

The paper is outlined as follows. In the next section, we propose the model with logistic growth and T. rangli's pathogenic effect on triamine bugs. In Section 3, the existence and stability of vector-free equilibrium, parasite-free equilibrium and parasite-positive equilibrium are considered. In Section 4, the bifurcation analysis including forward bifurcation and Hopf bifurcation is studied. Numerical simulations are also performed in Section 5. Conclusion and discussion are given in Section 6.

2.   Model

The model developed in the reference ^[10] is

where the population is divided into four compartments: susceptible and infected competent hosts, susceptible and infected triatomine bugs, denoted by $S_{h}, I_{h}, S_{v}, I_{v}$ in order. $\Lambda_{h}$ is the constant recruitment rate of susceptible competent host per unit time. The transmission rate from infected bugs to susceptible competent hosts is denoted by $\tilde{\beta}_{h} = \frac{ba}{N_{c}+\alpha N_{q}}$, where $b$ is the transmission probability from infected bugs to susceptible competent hosts per bite, $a$ is the number of bites per triatomine bug per unit time, $\alpha$ is the biting preference of quasi-competent hosts to competent hosts, $N_{c}$ is the total number of competent hosts, and $N_{q}$ is the total number of quasi-competent hosts. The transmission rate from infected competent hosts to susceptible bugs is denoted by $\tilde{\beta}_{v} = \frac{ca}{N_{c}+\alpha N_{q}}$, where $c$ is the transmission probability from infected hosts to susceptible triatomine bug per bite. The total infection rate through co-feeding transmission between susceptible and infected bugs is $\beta_{c}$, which is transmitted by both the competent and quasi-competent hosts. Here $\beta_{c} = \frac{1}{\delta}\beta_{ch}N_{c}(\frac{a}{N_{c}+\alpha N_{q}}\frac{\tau_{1}}{\omega})^{2}+\frac{1}{\delta}\beta_{cq}N_{q}(\frac{a\alpha}{N_{c}+\alpha N_{q}}\frac{\tau_{2}}{\omega})^{2}$, where $\omega$ is the unit time, $\delta$ is the ratio of night to unit time, $\beta_{ch}$ and $\beta_{cq}$ are the transmission rates from infected bugs to susceptible bugs on an average competent host and quasi-competent host during night, respectively. $\tau_{1}$ and $\tau_{2}$ are the feeding times per bite on a competent host and a quasi-competent host, respectively. The Ricker's type function $b(x) = rxe^{-\sigma x}$ was chosen to model the reproduction rate of R. prolixus. Integrating the pathogenic effect, the growth rate of triatomine bugs is modeled as $r(S_{v}+\theta I_{v})e^{-\sigma(S_{v}+I_{v})}$, where $r$ is the maximal number of offsprings that a triatomine bug can produce per unit time, $\theta \in [0, 1]$ is the reproduction reduction of bugs due to the pathogenic effect of T. rangeli on bugs, $\sigma$ is the density-dependency strength measuring the reproduction of bugs. $\mu_{1}$ and $\mu_{2}$ are the natural death rates of competent hosts and triatomine bugs, respectively. $d$ is the death rate of infected vectors induced by pathogenic effect.

Model (2.1) includes the systemic and co-feeding transmission routes among vectors and hosts. Two thresholds were derived to study the dynamical behavior of this model ^[10]. Sustained oscillations were found numerically by changing the parameters $d$ and $\theta$. Furthermore, the oscillation amplitude is larger if $d$ is larger or $\theta$ is smaller.

In this paper, we assume the generation rate of triatomine bugs follows the logistic growth instead of Ricker's type function. Here we only consider the systemic transmission, then model (2.1) can be changed to

where $\beta_{h} = \frac{ba}{N_{h}}, \beta_{v} = \frac{ca}{N_{h}}, \; N_{h}$ is the total number of hosts. The infection rate of susceptible competent hosts after bugs' biting at time $t$ is $\beta_{h}I_{v}(t)S_{h}(t)$, and the infection rate of susceptible vectors after bugs' biting at time $t$ is $\beta_{v}S_{v}(t)I_{h} (t)$. In the logistic growth $r(S_{v}+\theta I_{v})(1-\frac{S_{v}+I_{v}}{K})$ of triatomine bugs, $K$ is the carrying capacity. The other parameters are the same as those in model (2.1). All the parameters are non-negative and their biological meanings and ranges are given in Table 1.

Denote the total population of competent hosts by $N_{h} = S_{h}+I_{h}$. Adding the first and second equations leads to

It follows that

Thus, in the limiting system we have $S_{h}(t) = N_{h}-I_{h}(t)$.

Accordingly, system (2.2) can be reduced to the following three-dimensional limiting system:

It is easy to know that the feasible region of system (2.3) is

3.   Existence and stability of equilibria

3.1.   Existence of equilibria

Let the right-hand side of the equations of system (2.3) be zero. There are one vector-free equilibrium $E_{0}(0, 0, 0)$ and one parasite-free equilibrium $E_{S}(0, \frac{K(r-\mu_{2})}{r}, 0)$. We will derive two thresholds $R_{v}$ and $R_{0}$ to study the dynamic behavior of system (2.3) where $R_{v}$ is the triatomine bug basic reproduction number and $R_{0}$ is the T. rangeli basic reproduction number.

The Jacobian matrix of system (2.3) at the vector-free equilibrium $E_{0}(0, 0, 0)$ is

The eigenvalues of the Jacobian matrix $J(E_{0})$ at $E_{0}(0, 0, 0)$ are $-\mu_{1}, -(d+\mu_{2})$ and $r-\mu_{2}$, respectively. Let

We will show that it provides a threshold to determine the persistence or extinction of the vector population.

For the parasite-free equilibrium $E_{S}(0, S_{v}^{0}, 0)$ with $S_{v}^{0} = \frac{K(r-\mu_{2})}{r}$ to be biologically feasible, we need $r-\mu_{2} > 0$, namely, $R_{v} > 1$. Next, we calculate the T. rangeli basic reproduction number $R_{0}$ of system (2.3). Using the method in ref. ^[40], we have

Thus, the T. rangeli basic reproduction number of system (2.3), given by the spectral radius of the next generation matrix, is

where $S_{v}^{0} = \frac{K(r-\mu_{2})}{r}$.

For any parasite-positive equilibrium $E^{*} = (I_{h}^{*}, S_{v}^{*}, I_{v}^{*})$ of system (2.3), its elements satisfy

and $I_{h}^{*}$ is the positive root of the following equation:

where

Let $\Delta = B^{2}-4AC$. We have

If $\Delta\geq0$, then the equation $f(I_{h}^{*}) = 0$ may have two roots, which are denoted by

From (3.2), we can see that $S_{v}^{*} > 0, I_{v}^{*} > 0$ as long as $I_{h}^{*} > 0$. Therefore, we study the existence of the parasite-positive equilibria in the following cases:

1) $R_{0} = 1$, namely, $a^{2}bcK(r-\mu_{2}) = rN_{h}\mu_{1}(d+\mu_{2})$.

In this case, we have $B > 0, C = 0, \Delta > 0, A = \frac{ac\mu_{1}^{2}}{r-\mu_{2}}((d+\mu_{2})(r\theta-d-\mu_{2})+acr \theta (r-\mu_{2}))$. By Vieta theorem, if $A > 0$, then $I_{h1}^{*} < 0, I_{h2}^{*} = 0$. Thus, system (2.3) has no parasite-positive equilibrium. If $A < 0$, then $I_{h1}^{*} = 0, I_{h2}^{*} > 0$, and system (2.3) has a unique parasite-positive equilibrium $E_{2}(I_{h2}^{*}, S_{v}^{*}(I_{h2}^{*}), I_{v}^{*}(I_{h2}^{*}))$. If $A = 0$, then Eq (3.3) has a zero root, i.e., there is no parasite-positive equilibrium of system (2.3).

2) $R_{0} > 1$, namely, $a^{2}bcK(r-\mu_{2}) > rN_{h}\mu_{1}(d+\mu_{2})$.

In this case, we have $B > 0, C < 0$. If $A > 0$, then $\Delta > 0$, $I_{h1}^{*} < 0, I_{h2}^{*} > 0$, and system (2.3) has a unique parasite-positive equilibrium $E_{2}(I_{h2}^{*}, S_{v}^{*}(I_{h2}^{*}), I_{v}^{*}(I_{h2}^{*}))$. If $A < 0, \Delta > 0$, then $I_{h1}^{*} > 0, I_{h2}^{*} > 0$, and system (2.3) has two parasite-positive equilibria $E_{1}(I_{h1}^{*}, S_{v}^{*}(I_{h1}^{*}), I_{v}^{*}(I_{h1}^{*}))$ and $E_{2}(I_{h2}^{*}, S_{v}^{*}(I_{h2}^{*}), I_{v}^{*}(I_{h2}^{*}))$. If $A < 0, \Delta = 0$, then $E_{1} = E_{2},$ which means that there is a parasite-positive equilibrium of multiplicity 2. If $A < 0, \Delta < 0$, then there is no parasite-positive equilibrium. If $A = 0$, there exists only one root of Eq (3.3) and the root is positive, i.e., one parasite-positive equilibrium $E^{*}(I_{h}^{*}, S_{v}^{*}(I_{h}^{*}), I_{v}^{*}(I_{h}^{*}))$ of system (2.3).

3) $R_{0} < 1$, namely, $a^{2}bcK(r-\mu_{2}) < rN_{h}\mu_{1}(d+\mu_{2})$.

In this case, we have $B > 0, C > 0$. If $A < 0$, then $\Delta > 0$, $I_{h1}^{*} < 0, I_{h2}^{*} > 0$, and system (2.3) has a unique parasite-positive equilibrium $E_{2}(I_{h2}^{*}, S_{v}^{*}(I_{h2}^{*}), I_{v}^{*}(I_{h2}^{*}))$. If $A\geq0$, there is no positive root of (3.3), i.e., there is no parasite-positive equilibrium of system (2.3).

We summarize the results as follows:

Theorem 3.1 For system (2.3), we have the following results on the existence of equilibria.

1) The vector-free equilibrium $E_0(0, 0, 0)$ always exists. The parasite-free equilibrium $E_S(0, \frac{K(r-\mu_2)}{r}, 0)$ exists if and only if $R_v > 1$.

2) When $R_{0}\leq1$, there is a unique parasite-positive equilibrium $E_{2}$ if $A < 0$; Otherwise, there is no parasite-parasite-positive equilibrium.

3) When $R_{0} > 1$, there is a unique parasite-positive equilibrium if $A\geq0$, and there are two parasite-positive equilibria $E_{1}$ and $E_{2}$ if $A < 0, \Delta > 0$, and the two equilibria coalesce to $E$ if and only if $A < 0, \Delta = 0$.

3.2.   Stability of equilibria

3.2.1.   Stability of vector-free and parasite-free equilibria

Theorem 3.2 The vector-free equilibrium $E_{0}(0, 0, 0)$ of system (2.3) is globally asymptotically stable if $R_{v} < 1$ and unstable if $R_{v} > 1$.

Proof. At the vector-free equilibrium $E_{0}(0, 0, 0)$, the eigenvalues of the Jacobian matrix of system (2.3) are $-\mu_{2}, -(d+\mu_{2})$ and $r-\mu_{2}$. If $R_{v} = \frac{r}{\mu_{2}} > 1$, namely, $r > \mu_{2}$, then there are two negative eigenvalues and one positive eigenvalue, i.e., $E_{0}$ is unstable. If $R_{v} < 1$, then all the eigenvalues are real and negative, which indicates that $E_{0}$ is locally asymptotically stable for $R_{v} < 1$.

Let $N_{v} = S_{v}+I_{v}$. We have

Thus, we have

with any feasible initial solution $N_{v}(0) = S_{v}(0)+I_{v}(0)$ when $R_{v} < 1$. That is to say, the solutions of $S_v$ and $I_v$ with any feasible initial conditions will tend to zeroes if $R_v < 1.$ For subsystem of $I_{h}$, it is cooperative with the positive invariance set $[0, N_{h}]$. The vector-free equilibrium $E_0$ is unique for system (2.3) if $R_{v} < 1$. From this, we know that $E_{0}$ is globally asymptotically stable if $R_{v} < 1$.

Theorem 3.3 The parasite-free equilibrium $E_{S}(0, \frac{K(r-\mu_{2})}{r}, 0)$ of system (2.3) is

1) a saddle-node point when $R_{0} = 1$;

2) unstable when $R_{0} > 1$;

3) locally asymptotically stable when $R_{0} < 1$.

Proof. At the parasite-free equilibrium $E_{S}(0, S_{v}^{0}, 0)$, where $S_{v}^{0} = \frac{K(r-\mu_{2})}{r}$, the Jacobian matrix of system (2.3) is

The corresponding characteristic polynomial of (3.4) is

where

The eigenvalues of system (2.3) at $E_{S}$ are the roots of $P(\lambda) = 0$ and denoted by $\lambda_{1}, \lambda_{2},$ and $\lambda_{3}$.

Let $\Delta_{0} = b_{0}^{2}-4c_{0}$, i.e.,

Then the eigenvalues of (3.5) are

Since $r > \mu_{2}$, we have $\Delta_{0} > 0$, which means that $\lambda_{2}, \lambda_{3}$ are real.

When $R_{0} = 1$, we have $c_{0} = 0, b_{0} = \sqrt{\Delta_{0}}$. Therefore, $\lambda_{1} = \mu_{2}-r < 0, \lambda_{2} = -b_{0} = -d-\mu_{1}-\mu_{2} < 0$ and $\lambda_{3} = 0$. Because of the zero eigenvalue, we need to further study the type of $E_{S}$.

We let $x = I_{h}, y = S_{v}-S_{v}^{0}, z = I_{v}$ and shift the equilibrium to the origin. System (2.3) becomes

Then we make the following transformations

where

This leads to the following system

We know that $E_{S}$ is a saddle-node point when $R_{0} = 1$. Moreover, if $R_{0} > 1$, then $\lambda_{1} < 0, \lambda_{2} < 0, \lambda_{3} > 0$, i.e., $E_{S}$ is unstable; if $R_{0} < 1$, then $\lambda_{1} < 0, \lambda_{2} < 0, \lambda_{3} < 0$, i.e., $E_{S}$ is locally asymptotically stable.

Based on the above analysis, we conclude that the parasite-free equilibrium $E_S$ of system (2.3) is a saddle-node point when $R_{0} = 1,$ unstable when $R_{0} > 1,$ and locally asymptotically stable when $R_{0} < 1$ ^[41].

3.2.2.   Stability of the unique parasite-positive equilibrium

In this section, we will study the global stability of the unique parasite-positive equilibrium $E_{1}$ by Li-Muldowney global-stability criterion ^[42] when $R_{0} > 1$, $R_{v} > 1$ and $A\geq0$.

Let $|\cdot|$ denote a vector norm in $R^{n}$ and the induced matrix norm in $R^{n\times n}$, the space of all $n\times n$ matrices. For matrix $A$ in $R^{n\times n}$, the Lozinsk$\tilde{\rm i}$ measure or the logarithmic norm of $A$ with respect to $|\cdot|$ ^[43] is

Let $y(t)$ be a solution of linear differential equation

where $A(t)$ is $m\times m$ matrix-valued continuous function. For $t\geq t_{0}$, we have

Let $B$ be an $n\times n$ matrix. The second additive compound matrix of $B$, denoted by $B^{[2]}$, is an $(^{n}_{2})\times (^{n}_{2})$ matrix. For instance, if $B = (b_{ij})$ is a $3\times 3$ matrix, then

Consider the following autonomous system

where $f:\Omega\rightarrow R^{n}, \Omega\subset R^{n}$ is an open set, simply connected, and $f\in C^{1}(\Omega)$. Let $x(t, x_{0})$ be the solution of system (3.8) such that $x(0, x_{0}) = x_{0}$. Suppose $x^{*}$ is an equilibrium of system (3.8), i.e., $f(x^{*}) = 0$. A set $K$ is said to be absorbing in $\Omega$ for system (3.8) if $x(t, K_{1})\subset K$ for each compact set $K_{1}\subset \Omega$ and sufficiently large $t$. Assume the following assumptions hold:

$(H_{1})$ System (3.8) has a unique equilibrium point $x^{*}$ in $\Omega$.

$(H_{2})$ System (3.8) has a compact absorbing set $K\subset \Omega$.

Let $Q:\Omega\mapsto Q(x)$ be $(^{n}_{2})\times (^{n}_{2})$ matrix-valued function with its inverse $Q^{-1}(x)$. Let $\mu$ be a Lozinsk$\tilde{\rm i}$ measure on $R^{N \times N}$, where $N = (^{n}_{2})$. Define

where

the matrix $Q_{f}$ is obtained by replacing each entry $q_{ij}$ of $Q$ by its derivative in the direction of $f$, $(q_{ij})f$, and $J^{[2]}$ is the second additive compound matrix of the Jacobian matrix $J$ of system (3.8).

Lemma 3.1 ^[42] Assume that $\Omega$ is simply connected and assumptions $(H_{1})$ and $(H_{2})$ hold. Then, the unique equilibrium $x^{*}$ of system (3.8) is globally asymptotically stable in $\Omega$ if there exist a function $Q$ and a Lozinsk$\tilde{\rm i}$ measure $\mu$ such that $\bar{q_{2}} < 0$.

We have the following theorem for our model.

Theorem 3.4 Assume that $R_{0} > 1$, $R_{v} > 1$ and $A\geq0$. The unique parasite-positive equilibrium $E_{1}$ of system (2.3) is globally asymptotically stable if $\sigma = K^{2}\beta_{v}+\tilde{v} < 0$, where

Proof. When $R_{0} > 1$, $R_{v} > 1$ and $A\geq0$, it is easy to show the uniqueness of the parasite-positive equilibrium $E_{1}$. By Theorems 3.2 and 3.3, we know that the vector-free equilibrium $E_{0}$ and the parasite-free equilibrium $E_{S}$ are unstable when $R_{v} > 1$ and $R_{0} > 1$. It can also be checked that the conditions $(H_{1})$ and $(H_{2})$ are satisfied. Next we show the global stability of the unique parasite-positive equilibrium $E_{1}(I_{h}, S_{v}, I_{v})$.

The Jacobian matrix of system (2.3) at $E_{1}$ is

where $J_{22} = r(1-\frac{I_{v}+S_{v}}{K})-I_{h}\beta_{v}-\mu_{2}-\frac{r(S_{v}+\theta I_{v})}{K}, J_{23} = r\theta(1-\frac{I_{v}+S_{v}}{K})-\frac{r(S_{v}+\theta I_{v})}{K}.$ The second additive compound matrix $J^{[2]}$ is

Let $P(I_{h}, S_{v}, I_{v}) =$ diag$(1, \frac{S_{v}}{I_{v}}, \frac{S_{v}}{I_{v}})$. Then $P_f =$ diag$(0, \frac{\dot{S_{v}}}{I_{v}}-\frac{S_{v}\dot{I_{v}}}{I_{v}^2}, \frac{\dot{S_{v}}}{I_{v}}-\frac{S_{v}\dot{I_{v}}}{I_{v}^2})$. We have $P_{f}P^{-1} =$ diag$(0, \frac{S_{v}^{'}}{S_{v}}-\frac{I_{v}^{'}}{I_{v}}, \frac{S_{v}^{'}}{S_{v}}-\frac{I_{v}^{'}}{I_{v}}),$ and let $B = P_{f}P^{-1}+PJ^{[2]}P^{-1}$.

From system (2.3), we obtain

Straightforward calculations yield

where

Take the norm $|(I_{h}, S_{v}, I_{v})| = \max\{|I_{h}|, |S_{v}|+|I_{v}|\}$ in $R^{3}$. $\mu(\cdot)$ is the Lozinsk$\tilde{\rm i}$ measure with the vector norm ^[44]. We have

where $|B_{12}|, |B_{21}|$ are the matrix norms with respect to $l_{1}$ vector norm.

Calculations show that

Thus, we have

where $\tilde{v}, \; v_1, \; v_2$ are defined by (3.9).

According to $S_{v}+I_{v}\leq K$, we have $\mu_{B}\leq\frac{S_{v}^{'}}{S_{v}}+K^{2}\beta_{v}+\tilde{v},$ that is

where $\sigma = K^{2}\beta_{v}+\tilde{v}.$

Along each solution $(I_{h}, S_{v}, I_{v})\subset D$ of system (2.3), then for $t > \bar{t},$ we have

which means that $\bar{q_{2}}\leq \sigma < 0.$ This completes the proof.

T. rangeli may not be pathogenic to every triatomine species ^[17]. Thus, we also study the dynamics of system (2.3) in the absence of pathogenic effect on triatomine bugs, i.e., $d = 0$ and $\theta = 1.$ This allows us to compare the obtained results with and without the pathogenic effect.

Theorem 3.5 Assume $R_{v} > 1$ and $R_{0} > 1$. In the absence of pathogenic effect on triatomine bugs, namely, $\theta = 1$ and $d = 0$, system (2.3) admits a unique parasite-positive equilibrium $E^{*} = (I_{h}^{*}, S_{v}^{*}, I_{v}^{*})$, which is locally asymptotically stable when $d^{*} > 0$, unstable when $d^{*} < 0$, where $d^*$ is defined by (3.13).

Proof. In the case of $\theta = 1, d = 0$, system (2.3) becomes

According to Theorem 3.1, the system (3.10) has a unique parasite-positive equilibrium $E^{*} = (I_{h}^{*}, S_{v}^{*}, I_{v}^{*})$ when $R_{0} > 1$. Adding the second and third equations of system (3.10), we have

Letting the right-hand side of equation (3.11) be equal to zero, we obtain a unique parasite-positive equilibrium $N_{v}^{*} = S_{v}^{0}$ if $R_{v} > 1$, and a unique zero equilibrium which is globally asymptotically stable if $R_{v}\leq1$.

When $R_{v} > 1$, the limiting system of system (3.10) can be reduced to

System (3.12) has an unstable parasite-free equilibrium $(0, 0)$. When $R_{0} > 1$, there is a unique parasite-positive equilibrium $E_{1} = (I_{h}^{*}, I_{v}^{*})$, where $S_v^*$ is replaced by $S_v^0-I_v^*$ for susceptible vectors at equilibrium. The Jacobian matrix of system (3.12) at $E_{1}$ is

It is easy to know that the trace of $J(E_{1})$ is negative. The determinant of $J(E_{1})$ is

The eigenvalues of Jacobian matrix $J(E_{1})$ have negative real parts when $d^{*} > 0$. When $d^{*} < 0$, one of the eigenvalues of Jacobian matrix $J(E_{1})$ has a negative real part and the other has a positive real part. Therefore, the equilibrium $E_{1}(I_{h}^{*}, I_{v}^{*})$ of system (2.3) is locally asymptotically stable if $d^{*} > 0$, unstable if $d^{*} < 0$. That is to say, the parasite-positive equilibrium $E^{*} = (I_{h}^{*}, S_{v}^{*}, I_{v}^{*})$ of system (2.3) is locally asymptotically stable if $d^{*} > 0$, unstable if $d^{*} < 0$.

Denote

Obviously, both $f$ and $g: \mathbb{R}^2_+\rightarrow \mathbb{R}$ are continuously differentiable maps. We have $\frac{\partial f}{\partial I_v} = \beta_{h}(N_{h}-I_{h})-\mu_{1}I_{h}\geq0,$ and $\frac{\partial g}{\partial I_h} = \beta_{v}(S_{v}^{0}-I_{v})-\mu_{2}I_{v}\geq0.$ The system is cooperative in a domain $\mathbb{D} = \{(I_n, I_v)\in \mathbb{R}^2: I_h\in [0, N_h], I_v\in [0, S_v^0]\}.$ System (3.12) has a parasite-free equilibrium $(0, 0)$ and a unique parasite-positive equilibrium $E_1.$ According to Theorem 3.2.2 in ^[45], $E_1$ is globally attractive. Thus, the parasite-positive equilibrium $E^*$ of sub-system (3.12) is globally asymptotically stable. This shows that the pathogenic effect may cause the system to be unstable and be responsible for the occurrence of sustained oscillations.

4.   Bifurcation analysis

In this section, we will analyze the existence of forward bifurcation and Hopf bifurcation of system (2.3).

Theorem 4.1 System (2.3) exhibits a forward bifurcation from $E_{S}(0, \frac{K(r-\mu_{2})}{r}, 0)$ when $R_{0} = 1$. Furthermore, no backward bifurcation occurs.

Proof. The Jacobian matrix of system (2.3) at $E_{S}(0, \frac{K(r-\mu_{2})}{r}, 0)$ is

Its eigenvalues are

where

Because all parameter values are non-negative, we know that $\lambda_{1}$ is always negative.

If $R_{0} = 1$, then

Substituting $c = c^{*}$ into $\lambda_{2}$ and $\lambda_{3}$, we have $\lambda_{2} < 0, \lambda_{3} = 0$. Also, the parasite-free equilibrium $E_{S}$ is locally stable when $c < c^{*}$, and unstable when $c > c^{*}$. Therefore, $c = c^{*}$ is a bifurcation value.

We obtain a right eigenvector $u$ and a left eigenvector $\bar{v}$ associated with the zero eigenvalue, where

By the orthogonal condition $< u, \bar{v} > = 1$, we get

By the transformation

and noticing that system (2.3) has the form $\frac{dx}{dt} = f$, where $x = (x_{1}, x_{2}, x_{3})^{T}$ and $f = (f_{1}, f_{2}, f_{3})^{T}$, we have

The formula of the bifurcation coefficient in system (4.1) at $E_{S}$ is:

Since $\bar{v}_{2} = 0$, what we need to consider are the cross derivatives of $f_{1}$ and $f_{3}$ in system (4.1) at the equilibrium $E_{S}$. We obtain some non-zero terms

Now we calculate the values of $\bar{a}$ and $\bar{b}$.

From ^[41], we know that the local dynamical behavior of system (4.1) at equilibrium $E_{S}$ is determined by the signs of $\bar{a}$ and $\bar{b}$. From the above calculation, we have $\bar{a} < 0$ and $\bar{b} > 0$. Thus, there exists a locally asymptotically stable endemic equilibrium of system (4.1) showing a forward bifurcation near the equilibrium $E_{S}.$

Remark 4.2 We conclude that no backward bifurcation occurs for system (4.1). The forward bifurcation of system (4.1) is shown by Figure 7 in Section 5.

For any parasite-positive equilibria $E^{*} = (I_{h}^{*}, S_{v}^{*}, I_{v}^{*})$, the Jacobian matrix of system (2.3) at $E^{*}$ is

where

The corresponding characteristic polynomial is

where

$J_{22}^*, J_{23}^*$ are defined by (4.2), $S_{v}^{*}$ and $I_{v}^{*}$ are defined in (3.2) and the coordinate of $I_{h}^{*}$ is a positive root of (3.3).

From the above calculation and reference ^[46], we have the following theorem.

Theorem 4.3 The parasite-positive equilibrium $E^{*}$ of system (2.3) undergoes

1) a static bifurcation if $c_{1}|_{E^{*}} = 0$ and $\Delta_{1, 2}|_{E^{*}} > 0$;

2) a Hopf bifurcation if $\Delta_{2} = 0, \frac{d\Delta_{2}}{d(Bif.)}\neq0, \Delta_{1}|_{E^{*}} > 0$ and $c_{1}|_{E^{*}} > 0$.

Here $a_{1}, b_{1}, c_{1}$ are the coefficients of the characteristic polynomial (4.3), $d(Bif.)$ is the differentiation of the bifurcation parameter, and $\Delta_{1} = a_{1}, \Delta_{2} = a_{1}b_{1}-c_{1}$ are the first and the second Hurwitz arguments, respectively.

5.   Numerical results

In this section, we conduct numerical simulations for system (2.3) by using Matlab and Auto07P ^[47]. We choose the following parameter values,

which were also used in the reference ^[10]. It is easy to calculate that there are one vector-free equilibrium, one parasite-free equilibrium, and one parasite-positive equilibrium point with the above parameter values.

We study the effect of the number of competent hosts $N_h$ and the carrying capacity of triatomine bugs $K$ on the basic reproduction number of T. rangeli $R_0$. $R_0$ is proportional to $K.$ Thus, $R_0$ increases as $K$ increases. In particular, $R_0$ increases with an increasing number of hosts and the carrying capacity of bugs $K$. This is shown in Figure 2(a). $R_0$ depends on $\beta_h,$ which is an important parameter and its expression is a combination of $N_h$ and $K.$ We find that the slopes of the curves increase when the number of competent hosts decreases. Because $R_0$ is related to $\beta_h$ and $K,$ the transmission rate from infected bugs to susceptible hosts is inversely proportional to the number of hosts, as shown in Figure 2(b). Figure 3(a) shows the occurrence of sustained oscillation as the parameter values are defined in (5.1) while Figure 3(b) illustrates the stability of the parasite-positive equilibrium in Theorem 3.4.

Next, we consider the influence of the maximal number of offsprings that a triatomine bug produces per unit time $(r)$, the carrying capacity of bugs ($K$), the pathogenic effect ($d$), and the transmission probability from infected hosts to susceptible triatomine bug per bite $(c)$ in model (2.3) by using one-parameter and two-parameter bifurcation analysis.

We start with the maximal number of offsprings that a triatomine bug can produce per unit time. We choose $r$ as the primary bifurcation parameter and keep the other parameters fixed as in (5.1). The one-parameter bifurcation diagram is shown in Figure 4. There exist two transcritical bifurcation points $TC_1(0, 0, 0)$ and $TC_2(0, 3.10847, 0)$ when $r = 8.3\times10^{-3}$ and $r = 8.326\times10^{-3},$ one supercritical Hopf bifurcation point $HB_{1}(1.94275, 3.12364, 0.135572)$ when $r = 9.39283\times 10^{-3}$, and one supercritical Hopf bifurcation point $HB_{2}(8.18863, 3.17343, 0.580539)$ when $r = 1.23408\times 10^{-2}$. The numbers of infected competent hosts, susceptible triatomine bugs and infected triatomine bugs will increase gradually when $0 < r < 9.39283\times 10^{-3}$ and $r > 1.23408\times 10^{-2}.$ There is an unstable interval in which supercritical Hopf bifurcation occurs when $9.39283\times 10^{-3} \leq r\leq1.23408\times 10^{-2}.$ The red solid curve represents the stable limit cycle branch bifurcating from the supercritical Hopf bifurcation point, which indicates the appearance and the disappearance of stable limit cycle with the increase of the parameter $r$. Thus, when the maximal number of offsprings of susceptible triatomine bugs increases, the number of infected competent hosts and the number of infected triatomine bugs will increase out of the Hopf bifurcation interval $r\in[9.39283\times 10^{-3}, 1.23408\times 10^{-2}]$. Further, if we use $\mu_1$ and $r$ as the primary bifurcation parameters, then we obtain a two-parameter Hopf bifurcation curve, where the line $\mu_1 = 0.0025$ corresponds to the parameter values for Hopf bifurcation shown in Figure 4(a)–(c). The limit cycle branch connects the two Hopf bifurcation points and the period of limit cycle is finite. From Figure 4(d), we find that there are one or two Hopf bifurcation points when $0\leq \mu_{1}\leq0.0055$, and no Hopf bifurcation occurs when $\mu_{1} > 0.0055$.

The carrying capacity of triatomine bugs is the number that the ecosystem can sustainably support. It may vary due to many factors. We also use the carrying capacity $K$ as the primary bifurcation parameter. With $\theta = 0.3$, we obtain the one-parameter bifurcation diagram for system (2.3) shown in Figure 5. There are one transcritical bifurcation point $TC(0, 3.10847, 0)$ and one supercritical Hopf bifurcation point $HB(33.9734, 3.39698, 2.57824)$ when $K = 4.45926$ and $K = 5.45626\times 10^2,$ respectively. The amplitudes and periods of limit cycles bifurcating from $HB$ become larger as $K$ increases. When the carrying capacity of susceptible triatomine bugs increases, there will be always a supercritical Hopf bifurcation point, i.e., all the state variables will vary periodically, and all the competent hosts, susceptible triatomine bugs and infected triatomine bugs will coexist. In addition, if we use $r$ and $d$ as the primary two bifurcation parameters, then we obtain the two-parameter Hopf bifurcation curves which are all closed curves when $K = 800, 1000, 1200,$ respectively. This also indicates that no Bogdanov-Takens bifurcation of the parasite-positive equilibrium occurs for system (2.3).

Next, we use $d$ as the primary bifurcation parameter to study the influence of the pathogenic effect in system (2.3). The parameter $\theta$ is set to 0.2 and other parameters are fixed as in (5.1). We obtain the one-parameter bifurcation diagram, shown in Figure 6(a)–(c). There are two supercritical Hopf bifurcation points $HB_{1}(27.3612, 8.51017, 2.03959)$ and $HB_{2}(36.7308, 1.90856, 2.80866)$ when $d = 0.0756105$ and $d = 0.0100453$, respectively. Thus, when the death rate of infected vectors increases due to the strong pathogenic effect, the number of the competent hosts and infected triatomine bugs will decrease, and the number of susceptible triatomine bugs will increase except an unstable interval for the occurrence of Hopf bifurcation. Two-parameter Hopf bifurcation curve is also given to illustrate the occurrence of stable limit cycles (Figure 6(d)), where the red dotted line is $K = 1000.$ Hopf bifurcation will always occur when $K\geq252.498$.

The transmission rate from infected competent hosts to susceptible bugs ($\beta_v$) depends on the number of bites per triatomine bug per unit time $(a)$ and the transmission probability from infected hosts to susceptible triatomine bugs per bite $(c).$ Since the two parameters $a$ and $c$ have a similar role in $\beta_v$, for simplicity, we only use $c$ as the primary bifurcation parameter and fix $\theta = 0.3$. We obtain the one-parameter bifurcation diagram (Figure 7(a)–(c)). There is a transcritical bifurcation point $TC(0, 6.9708\times 10^{2}, 0)$ when $c = 2.18504\times 10^{-2},$ which confirms the forward bifurcation as shown in Theorem 3.1. The numbers of infected hosts and infected vectors will increase firstly, then decrease dramatically to a low level. The appearance of the limit cycle indicates that T. rangeli parasites will persist with the increase of the transmission probability from infected hosts to susceptible triatomine bugs per bite. Also, there is a supercritical Hopf bifurcation point $HB(44.6167, 4.57113, 348737)$ when $c = 0.375043.$ The amplitudes and periods of limit cycles bifurcating from $HB$ become larger as $c$ increases. From Figure 7, we can see that all state variables coexist once the number of bites of per triatomine bug per unit time is greater than $0.865432.$ Therefore, T. rangeli parasites will always persist if the infection rates of susceptible triatomine bug and infected triatomine bug increase. The two-parameter $(c \ vs.\ K)$ Hopf bifurcation curve of system (2.3) is given in Figure 7(d), which tells the relationship of the carrying capacity and the transmission probability from infected hosts to susceptible triatomine bug per bite. There is only one supercritical Hopf bifurcation occurring for system (2.3).

From the above analysis, we know that the dynamics of triatomine bugs in system (2.3) are similar to HIV dynamics in ^[41,48]. System (2.3) undergoes a forward bifurcation instead of a backward bifurcation. However, from the limit cycle branches in Figures 5(a)–(c) and 7(a)–(c), we conclude that the T. rangeli parasites relevant to Chagas disease will persist due to the sustained oscillations from Hopf bifurcation when the carrying capacity $(K)$ or the transmission probability $(c)$ increases. Thus, it is challenging to eliminate T. rangeli parasites, a sister trypanosoma to T. cruzi and commonly causing the misdiagnosis of the Chagas disease.

Comparing with the dynamics of model (2.3) and revisiting model (2.1) numerically, we find that model (2.1) with Ricker's type function in ^[10] also doesn't undergo Bagdanov-Takens bifurcation at the parasite-positive equilibrium. However, model (2.1) has periodic-doubling bifurcation of limit cycles for $r,$ the maximal number of offsprings that a triatomine bug can produce per unit time. This indicates the occurrence of chaos for model (2.1), which differs from the dynamics of model (2.3).

6.   Discussion and conclusions

In this paper, we have formulated a model with a logistic growth of of triatmine bugs to study the dynamics of infected competent hosts, susceptible and infected triatomine bugs. The existence and stability of the vector-free equilibrium, parasite-free equilibrium and the parasite-positive equilibrium are studied. The direction of transcritical bifurcation and Hopf bifurcation is also investigated. Numerical simulations are conducted to illustrate and expand the theoretical results.

For many infectious disease models, the disease-free equilibrium usually loses its stability when $R_0$ increases to cross one, which results in a bifurcation where a curve of endemic equilibria emerges. The direction of this bifurcation is forward if the endemic curve occurs when $R_0$ is slightly above 1 and there is no parasite-positive equilibrium near the disease-free equilibrium for $R_0 < 1.$ In contrast, the bifurcation is backward if the bifurcating equilibrium occurs when $R_0 < 1.$ Basically, a backward bifurcation implies the occurrence of multiple endemic equilibria and the coexistence of a stable endemic equilibrium with a stable disease-free equilibrium. Thus, a forward bifurcation indicates that the infectious disease may be cured while it is not easy to eliminate the disease when a backward bifurcation occurs. Model (2.3) only goes through the forward bifurcation. However, this doesn't mean that the disease would be eliminated. The persistence of sustained oscillations of T. rangeli makes the disease eradication challenging.

One-parameter bifurcation diagrams for the parameters $r, K, d, c$ are showed, respectively. Oscillations always persists as $K$ or $c$ increases. This shows that the Trypanosoma rangeli will always exist in all the population and it is difficult to be eliminated. Bagdanov-Takens bifurcation of the parasite-positive equilibrium doesn't occur in models (2.1) and (2.3). Using numerical simulations, we also find that model (2.1) with Ricker's type function growth rather than the logistic function growth could go though periodic-doubling bifurcation of limit cycles when the maximal number of offsprings that a triatomine bug can produce per unit time $(r)$ increases. This suggests the emergence of chaos for model (2.1).

The pathogenic effect to the triatomine bug population is also studied. This may provide some critical insights for the prevention and control of Chagas disease. Although T. rangeli is not pathogenic to human, it is pathogenic to triatomine bugs, which have a great effect on the dynamics of T. cruzi population and triatomine bug population. It can cause the model to lose stability and undergo Hopf bifurcation. Moreover, the amplitudes of oscillations become larger as the parameters $K$ or $c$ increases. It is worth noting that we illustrate the non-existence of Bagdanov-Takens bifurcation of the parasite-positive equilibrium numerically without providing analytical results due to too many parameters in model (2.3).

This work provides more information that improves our understanding of the complexity of host-parasite-vector interactions.

Acknowledgments

The authors are very grateful to Professor Jianhong Wu for his helpful suggestions. This work was partially supported by the National Natural Science Foundation of China (No. 11671114) and NSF of Zhejiang (LY20A010002). X. Wu is supported by NSF of China (No.12071300). L. Rong is supported by the NSF grant DMS-1950254.

Conflict of interest

The authors declare there is no conflict of interest.