Modelling intra-host competition between malaria parasites strains

Abstract

An intra-host epidemiological model is formulated for the co-infection of drug-sensitive and drug-resistant malaria parasites to examine the impact of aggressive treatment on the effect of competitive release in an infected host. The analysis of the existence of equilibrium and their stability of the model is conducted, and the results reveal that the intra-host competition and treatment play a key role in the prevalence of drug-resistant strains. The mathematical outcomes qualitatively match the experimental fact, that the rapid elimination of drug-resistant strains could promote the very evolution it is intended to retard.

Appendix

Proof

(Proof of Theorem 2) The Jacobi matrix J(E) of right hand side of (2) at E is

$$\begin{aligned} J(E)=\left( \begin{array}{ccc}-\beta _1 I_s-\beta _2 I_r-d_1&{}-\beta _1 S&{}-\beta _2 S\\ \alpha \beta _1 I_s+p \alpha \beta _2I_r&{}\alpha \beta _1 S -d_2-\mu &{}p \alpha \beta _2 S\\ (1-p)\alpha \beta _2 I_r&{}0&{}(1-p)\alpha \beta _2 S-d_3\end{array}\right) . \end{aligned}$$

At $P_0$, the eigenvalues of $J(P_0)$ are

$$\begin{aligned} -d_1,\quad \frac{\varLambda \alpha \beta _1 }{d_1}-d_2-\mu ,\quad \frac{\varLambda \alpha \beta _2(1-p) }{d_1}-d_3 \end{aligned}$$

These quantities are negative whenever $R_0<1$. At $\tilde{P}=(\tilde{S},\tilde{I}_s,0)$,

$$\begin{aligned} J(\tilde{P})=\left( \begin{array}{ccc}-\beta _1 \tilde{I}_s-d_1&{}-\beta _1 \tilde{S}&{}-\beta _2 \tilde{S}\\ \alpha \beta _1 \tilde{I}_s&{}0&{}p \alpha \beta _2 \tilde{S}\\ 0&{}0&{}(1-p)\alpha \beta _2 \tilde{S}-d_3\end{array}\right) . \end{aligned}$$

The stability of $E_1$ is determined by the eigenvalue $(1-p)\alpha \beta _2 \tilde{S}-d_3$, which is also equal to $\frac{(1-p)(d_2 +\mu )\beta _2}{\beta _1}-d_3$. If $R_2 < R_1$, then this quantity is always less than 0. For $\hat{P}$, one has

$$\begin{aligned} J(\hat{P})=\left( \begin{array}{ccc}-\frac{\varLambda }{\hat{S}}&{} -\beta _1\hat{S}&{}-\beta _2\hat{S}\\ \frac{d_2\hat{I}_s}{\hat{S}}&{}\alpha \beta _1 \hat{S} -d_2-\mu &{}p \alpha \beta _2 \hat{S}\\ (1-p)\alpha \beta _2 \hat{I}_r&{}0&{}0\end{array}\right) . \end{aligned}$$

The characteristic equation of $J(\hat{P})$ is given by

$$\begin{aligned} \lambda ^3+A\lambda ^2+B\lambda +C=0. \end{aligned}$$

where

$$\begin{aligned} \begin{aligned} A=&\,\,d_1R_2+\frac{T_1}{\beta _2},\quad C=\frac{d_1d_3(R_2-1)T_1}{\beta _2}\\ B=&\,\,\frac{d_1R_2T_1}{\beta _2}+\frac{p\beta _1d_1d_2d_3(R_2-1)}{(1-p)T_2}+\frac{d_1d_3(R_2-1)T_1}{T_2},\\ T_1=&\,\beta _2(d_2+\mu )-\beta _1d_3/(1-p),\quad T_2=\,\beta _2(d_2+\mu )-\beta _1d_3 \end{aligned} \end{aligned}$$

It follows from (3) that

$$\begin{aligned} \alpha \beta _1 \hat{S} -d_2-\mu =\frac{\beta _1d_3}{\beta _2(1-p)}-(d_2+\mu )<0. \end{aligned}$$

and therefore, $T_2>T_1>0$, $A>0$, $B>0$ and $C>0$. Since

$$\begin{aligned} \begin{aligned} D&:=AB-C\\&=\frac{d_1R_2T_1A}{\beta _2}+\frac{d_1d_3(R_2-1)}{T_2\beta _2(1-p)}\left[ p\beta _1\beta _2d_1d_2R_2 +\beta _2d_1R_2T_1(1-p)+p\beta _1T_1(d_2-d_3)\right] \end{aligned} \end{aligned}$$

the Routh-Hurwitz criterion tells that $\hat{P}$ is locally asymptotically stable if and only if $D>0$. $\square$

Proof

(Proof of Theorem 4) The Jacobi matrix J(E) of right hand side of (2) is

$$\begin{aligned} J(E)=\left( \begin{array}{ccc}-\beta _1 I_s-\beta _2 I_r-d_1&{}-\beta _1 S&{}-\beta _2 S\\ \alpha \beta _1 I_s+p \alpha \beta _2I_r&{}\alpha \beta _1 S -\gamma _1 I_r-d_2-\mu &{}p \alpha \beta _2 S-\gamma _1 I_s\\ (1-p)\alpha \beta _2 I_r&{}-\gamma _2 I_r&{}(1-p)\alpha \beta _2 S-\gamma _2 I_s-d_3\end{array}\right) . \end{aligned}$$

At $P_0$, the eigenvalues of $J(P_0)$ are

$$\begin{aligned} -d_1,\quad \frac{\alpha \beta _1 \varLambda }{d_1}-d_2-\mu ,\quad \frac{(1-p)\alpha \beta _2 \varLambda }{d_1}-d_3 \end{aligned}$$

which are negative whenever $R_0<1$. At $\tilde{P}=(\tilde{S},\tilde{I}_s,0)$,

$$\begin{aligned} J(\tilde{P})=\left( \begin{array}{ccc}-\beta _1 \tilde{I}_s-d_1&{}-\beta _1 \tilde{S}&{}-\beta _2 \tilde{S}\\ \alpha \beta _1 \tilde{I}_s&{}0&{}p \alpha \beta _2 \tilde{S}-\gamma _1 \tilde{I}_s\\ 0&{}0&{}(1-p)\alpha \beta _2 \tilde{S}-\gamma _2 \tilde{I}_s-d_3\end{array}\right) . \end{aligned}$$

The stability of $\tilde{P}$ is determined by the eigenvalue $(1-p)\alpha \beta _2 \tilde{S}-\gamma _2 \tilde{I}_s-d_3$, which is equal to $\frac{(1-p)d_2\beta _2}{\beta _1}-\frac{\gamma _2(\varLambda \alpha \beta _1-d_1(d_2+\mu ))}{\beta _1(d_2+\mu )}-d_3$. Recall that $d_2<d_3$ and $\beta _1>\beta _2$ due to fitness cost of resistant strain. Then, this quantity is always less than 0.

At $P^*=(S^*,I_s^*,I_r^*)$, we have

$$\begin{aligned} J(E_2)=\left( \begin{array}{ccc}-\beta _1 I_s^*-\beta _2 I_r^*-d_1&{}-\beta _1 S^*&{}-\beta _2 S^*\\ \alpha \beta _1 I_s^*+p \alpha \beta _2I_r^*&{}\alpha \beta _1 S^* -\gamma _1 I_r^*-d_2-\mu &{}p \alpha \beta _2 S^*-\gamma _1 I_s^*\\ (1-p)\alpha \beta _2 I_r^*&{}-\gamma _2 I_r^*&{}(1-p)\alpha \beta _2 S^*-\gamma _2 I_s^*-d_3\end{array}\right) . \end{aligned}$$

The eigenvalues of $J(P^*)$ are given by zeros of a third order polynomial. However, we will not provide the detailed expression of characteristic equation for its complexity. The expression of $J(P^*)$ is enough for carrying out numerical simulations.

In order to show the occurrence of fixed point bifurcation, we change the order of equations of system (1), and let $I_r(t)=x, I_s(t)=y, S(t)=z$. Then,

$$\begin{aligned} \left\{ \begin{aligned} x'=&\,(1-p)\alpha \beta _2 x z-\gamma _2 x y-d_3 x\\ y'=&\,\alpha \beta _1 y z+p\alpha \beta _2 x z-\gamma _1 x y-(d_2+\mu ) y\\ z'=&\,\varLambda -\beta _1 y z-\beta _2 x z-d_1 z \end{aligned} \right. \end{aligned}$$

(9)

Shifting the fixed point $P_0=(\frac{\varLambda }{d_1},0,0)$ through ${\overline{x}}=x,{\overline{y}}=y,\overline{z}=z-\frac{\varLambda }{d_1}$, we have

$$\begin{aligned} \left\{ \begin{aligned} {\overline{x}}'=&\,((1-p)\alpha \beta _2 \frac{\varLambda }{d_1}-d_3) {\overline{x}}+(1-p)\alpha \beta _2 {\overline{x}} {\overline{z}}-\gamma _2 {\overline{x}} {\overline{y}}\\ {\overline{y}}'=&\,p\alpha \beta _2 \frac{\varLambda }{d_1} {\overline{x}}+ (\alpha \beta _1 \frac{\varLambda }{d_1}-(d_2+\mu )){\overline{y}}+ \alpha \beta _1 {\overline{y}} {\overline{z}}+p\alpha \beta _2 {\overline{x}} {\overline{z}}- \gamma _1 {\overline{x}} {\overline{y}}\\ {\overline{z}}'=&\,-\beta _2 \frac{\varLambda }{d_1} {\overline{x}}-\beta _1 \frac{\varLambda }{d_1} {\overline{y}} -d_1 {\overline{z}}-\beta _1 {\overline{y}} {\overline{z}}-\beta _2 {\overline{x}} {\overline{z}} \end{aligned} \right. \end{aligned}$$

(10)

Let $d_3^*=(1-p)\alpha \beta _2 \frac{\varLambda }{d_1}$ and $\xi =d_3^*-d_3$ be the bifurcation parameter. The Jacobian matrix of system (10) at (0, 0, 0) when $\xi =0$ is

$$\begin{aligned} J\big |_{(0,0,0)}=\left( \begin{array}{ccc}0&{}0&{}0\\ \frac{p \alpha \beta _2 \varLambda }{d_1}&{}\frac{\alpha \beta _1 \varLambda }{d_1}-(d_2+\mu )&{}0\\ -\frac{\beta _2 \varLambda }{d_1}&{}-\frac{\beta _1 \varLambda }{d_1}&{}-d_1\end{array}\right) \end{aligned}$$

Accordingly, its eigenvalues are $\lambda _1=0, \lambda _2=\frac{\alpha \beta _1 \varLambda }{d_1}-(d_2+\mu ), \lambda _3=-d_1$, and the corresponding eigenvectors are

$$\begin{aligned} \left( \begin{array}{c}1\\ a\\ b\end{array}\right) , \left( \begin{array}{c}0\\ 1\\ c\end{array}\right) , \left( \begin{array}{c}0\\ 0\\ 1\end{array}\right) \end{aligned}$$

with $a=\frac{p \alpha \beta _2 \varLambda }{d_1 (d_2+\mu )-\alpha \beta _1 \varLambda }, b=-\frac{1}{d_1}\left( \frac{\beta _2 \varLambda (d_1 (d_2+\mu )-\alpha \beta _1 \varLambda )+ \varLambda ^2 \alpha \beta _1 \beta _2 p}{d_1(d_1 (d_2+\mu ) - \alpha \beta _1 \varLambda )}\right) , c=\frac{\beta _1 \varLambda }{d_1 (d_2+\mu )-\alpha \beta _1 \varLambda -d_1^2}$.

Make the following change of variables

$$\begin{aligned} \left( \begin{array}{c}{\overline{x}}\\ {\overline{y}}\\ \overline{z}\end{array}\right) =T\left( \begin{array}{c}u\\ v\\ w\end{array}\right) \end{aligned}$$

(11)

where T is the invertible matrix whose columns are the three eigenvectors obtained above. Then,

$$\begin{aligned} \left( \begin{array}{c}u\\ v\\ w\end{array}\right) '=\left( \begin{array}{ccc}0&{}0&{}0\\ 0&{}\frac{\alpha \beta _1 \varLambda }{d_1}-d_2&{}0\\ 0&{}0&{}-d_1\end{array}\right) \left( \begin{array}{c}u\\ v\\ w\end{array}\right) + \left( \begin{array}{c}g_1(u,v,w,\xi )\\ g_2(u,v,w,\xi )\\ g_3(u,v,w,\xi )\end{array}\right) \end{aligned}$$

where

$$\begin{aligned} \begin{aligned} g_1(u,v,w,\xi )&=\xi u+((1-p)\alpha \beta _2 b -a \gamma _2)u^2+((1-p)\alpha \beta _2 b -\gamma _2)uv+(1-p)\alpha \beta _2 u w\\ g_2(u,v,w,\xi )&=-a \xi u+(-(1-p)\alpha \beta _2 a b + a^2 \gamma _2+\alpha \beta _1 a b + p \alpha \beta _2 b - a \gamma _1)u^2\\&\quad +(-(1-p)\alpha \beta _2 a b + a \gamma _2 + \alpha \beta _1 a c + \alpha \beta _1 b + p \alpha \beta _2 b c - \gamma _1)uv\\&\quad +(-a(1-p)\alpha \beta _2 + \alpha \beta _1 a + p \alpha \beta _2)uw+\alpha \beta _1 c v^2 + \alpha \beta _1 v w \\ g_3(u,v,w,\xi )&=(2ac-b)\xi u+ ((ac-b) ((1-p)\alpha \beta _2 b - a \gamma _2){-}c(\alpha \beta _1 a b {+}p \alpha \beta _2 b {-} a \gamma _1)\\&\quad -(\beta _1 a b + \beta _2 b))u^2 +((ac-b)((1-p) \alpha \beta _2 b - \gamma _2) \\&\quad -c(\alpha \beta _1 a c+ \alpha \beta _1 b + p \alpha \beta _2 b c -\gamma _1)-(\beta _1 a c + \beta _1 b + \beta _2 c))u v\\&\quad +((ac-b)((1-p) \alpha \beta _2)-c(\alpha \beta _1 a+p \alpha \beta _2)-(\beta _1 a+\beta _2))u w \\&\quad -(c^2 \alpha \beta _1 + c \beta _1)v^2-(c \alpha \beta _1+\beta _1)v w \end{aligned} \end{aligned}$$

Denote

$$\begin{aligned}&A=0,~B=\left( \begin{array}{cc}\lambda _2&{}0\\ 0&{}\lambda _3\end{array}\right) , ~f(u,v,w,\xi )\\&\quad =g_1(u,v,w,\xi ),~~g(u,v,w,\xi )= \left( \begin{array}{c}g_2(u,v,w,\xi )\\ g_3(u,v,w,\xi )\end{array}\right) \end{aligned}$$

Assume the center manifold takes the form of

$$\begin{aligned} v= & {} h_1(u,\xi )=a_1 u^2+a_2 u \xi + a_3 \xi ^2 + o(3)\\ w= & {} h_2(u,\xi )=b_1 u^2+b_2 u \xi + b_3 \xi ^2 + o(3) \end{aligned}$$

It then follows that

$$\begin{aligned} N(h)\triangleq D_u h(u,\xi )(A u +f(u, h_1, h_2 \xi ))-B h -g(u, h_1, h_2, \xi )=0. \end{aligned}$$

After some symbolic manipulations, we have

$$\begin{aligned} \begin{aligned} a_1&=- \frac{1}{\lambda _2}(-(1-p)\alpha \beta _2 a b + a^2 \gamma _2 + \alpha \beta _1 a b + p \alpha \beta _2 b - a \gamma _1)\\ a_2&=\frac{a}{\lambda _2},\quad a_3=0\\ b_1&= - \frac{1}{\lambda _3}((ac-b)((1-p)\alpha \beta _2 b - a \gamma _2)-c(\alpha \beta _1 a b +p \alpha \beta _2 b+ a \gamma _1)-(\beta _1 a b + \beta _2 b))\\ b_2&= - \frac{2ac-b}{\lambda _3}, \quad b_3=0 \end{aligned} \end{aligned}$$

Consequently, the restricted equation of (1) on the center manifold is given by

$$\begin{aligned} u'=\xi u+r u^2 +(a_1 s+t b_1)u^3+(a_2 s + t b_2)u^2 \xi \end{aligned}$$

where

$$\begin{aligned} r=(1-p)\alpha \beta _2 b - a \gamma _2, \quad s=(1-p)\alpha \beta _2 b - \gamma _2,\quad t=(1-p)\alpha \beta _2 b. \end{aligned}$$

$\square$