A MODEL FOR TRANSMISSION OF PARTIAL RESISTANCE TO ANTI-MALARIAL DRUGS

Anti-malarial drug resistance has been identified in many regions for a long time. In this paper we formulate a mathematical model of the spread of anti-malarial drug resistance in the population. The model is suitable for malarial situations in developing countries. We consider the sensitive and resistant strains of malaria. There are two basic reproduction ratios corresponding to the strains. If the ratios corresponding to the infections of the sensitive and resistant strains are not equal and they are greater than one, then there exist two endemic non-coexistent equilibria. In the case where the two ratios are equal and they are greater than one, the coexistence of the sensitive and resistant strains exist in the population. It is shown here that the recovery rates of the infected host and the proportion of anti-malarial drug treatment play important roles in the spread of anti-malarial drug resistance. The interesting phenomena of “long-time” coexistence, which may explain the real situation in the field, could occur for long period of time when those parameters satisfy certain conditions. In regards to control strategy in the field, these results could give a good understanding of means of slowing down the spread of anti-malarial drug resistance.

1. Introduction.Malaria is a critical public health problem in many countries.This disease is caused by Plasmodium spp.parasites and transmitted to humans through the bites of infected female Anopheles mosquitoes.It is estimated that 350-500 million clinical malaria cases occur annually, with most of these cases caused by the infection of Plasmodium falciparum [21].
Unfortunately, Plasmodium parasites have been resistant to some anti-malarial drugs in many regions.What we mean by anti-malarial drug resistance is "the ability of a parasite strain to multiply or to survive in the presence of concentrations of a drug that normally destroys parasites of the same species or prevents their multiplication" [6].The spread of anti-malarial drug resistance in a population is a major problem for the control of malaria.Currently, it is generally believed that the widespread use of a new anti-malarial drug will inevitably be followed by the appearance of resistance in the parasite population.Therefore, it is very important to make a rational management of the use of anti-malarial drugs in order to slow down the spread of the resistance [20].
Ideally, there is a parasite screening in order to give proper treatments for malaria patients.Parasites that have been resistant to an anti-malarial drug can not be eliminated by the drug.So, the hosts infected by the resistant strain should not be given the drug.Unfortunately, it is too costly to do a parasite screening in many developing countries.In these countries, malaria patients are treated without a prior knowledge of which strain infects them, so hosts infected by the sensitive strain or the resistant one get the same treatments.
In some instances the parasites are not fully resistant to a drug, that is, they are in the intermediate stage of drug tolerance.If this case occurs, the parasites still can be eliminated by the drug [13].
Mathematical modeling of malaria transmission has a long history [18,1,3].Recently, some mathematical models on the spread of anti-malarial drug resistance have been proposed and studied.The models may help identify those parameters that are important for the spread and fixation of drug resistant alleles in the parasite population, and may provide insight into some approaches toward controlling the spread of drug resistance.The attempts to understand the spread of anti-malarial drug resistance in a population using mathematical models can be classified into a population genetic approach, such as in [9,12,17] and an epidemiological approach, such as in [2,5,14,16].
In this paper we propose a new model to explore the effect of the proportion of treatment to the spread of the resistance.We think that the model is more appropriate for malarial situations in developing countries.The proposed model in this paper is similar, but different from the model of [14] in handling the treatment scenario and incorporation of the full vector-host dynamics.In [14], the resistant strain is assumed to be fully resistant to the anti-malarial drug, so treatment is given only to the hosts infected by the sensitive one.Here, we consider the intermediate stage of drug tolerance and the hosts infected by the resistant one also get treatment.
2. Model formulation.It is known that vertical transmission of the malaria parasite cannot occur in the vectors [8].In the host population, the infected mothers pass passive immunity to their babies.This immunity is transient, that is, it is effective only for the first three to six months of life.Although the infants could have immunity, they are not protected against infection, but the parasite density is reduced and the length of parasitaemia is shortened [19].
We assume the following things.The host and vector populations are constant for all time.The latent period of infected hosts is ignored.There is no superinfection, in the sense that the infected hosts can not have a new infection during the period of illness.The anti-malarial drug resistance occurs to a single drug and is believed to arise as a consequence of mutation(s) in a single gene of the parasite.The parasite population in an infected host contains the sensitive allele a and the resistant allele A. The resistant allele is partially resistant to the drug.Some fraction f a (f A ) of the sensitive-(resistant-) infected host subpopulations is under treatment.It is also assumed that the virulence of the parasite is somewhat negligible, so that the death rate is not exclusively associated with the disease.
Let the variables S, Ẽ, Ĩ, and R denote the susceptible, exposed, infected and recovered (temporary immune) subpopulations respectively and the indices H, V , a, A, U and T denote host, vector, sensitive-infected (in the sense that the frequency of sensitive allele a is greater than of the resistant allele A), resistant-infected (in the sense that the frequency of resistant allele A is greater than of the sensitive allele a), untreated and treated respectively.
We use the transmission scheme of anti-malarial drug resistance as in Figure 1 below.The description of the parameters can be seen in Table 1.   1 and their description.
For the biological interest, the parameters and the variables are non-negative.We assume that γ HaT > γ HaU , since recovery rates should be faster for the sensitiveinfected hosts under treatment than for those without treatment.We also assume that γ HAT ≤ γ HaT , since the drug should kill more sensitive strains than the resistant ones.For the resistant-infected recovery rates, we assume that γ HAU < γ HAT , where the resistant strain is partially resistant to the drug, so the drug still has the ability to kill the resistant strain.
Using the above assumptions and the transmission diagram, we formulate the following eleven-dimensional system of differential equations: where j ∈ {a, A}, ÑH = SH + RH + j∈{a,A} k∈{T ,U } ĨHjk and ÑV = SV + ẼV a + ẼV A + ĨV a + ĨV A .
Model ( 1) is well-posed in the non-negative region R 11 ≥0 because the vector field on the boundary does not point to the exterior.So, if it is given an initial condition in the region, then the solution is defined for all time t ≥ 0 and remains in the region.Moreover, the solution is bounded since the host and vector populations are constant for all time.
Since the host and vector populations are constant, we can eliminate the equations for RH and SV in model (1).Next, we scale model (1) using transformations

ÑV
, where j ∈ {a, A}, k ∈ {T, U }. Thus, we obtain the following nine-dimensional system of differential equations: is the scaled infection rate and σ Hjk = δ H + γ Hjk is the deathadjusted recovery rate.Model (2) has the region of biological interest In the following we will study model (2), because the reduced model inherits all properties of model (1).
3. Model analysis.Basic reproduction ratio is an important threshold in mathematical epidemiology.This threshold is given by the spectral radius of a next generation matrix K of a model.Using the method in [7], the next generation matrix K of model ( 2) is given by The element k ij of K denotes the expected number of new infective with index i, caused by an infective with index j.The indices {1, 2, 3, 4, 5, 6} of i and j correspond to {I HaT , I HaU , I V a , I HAT , I HAU , I V A } respectively.For our model, there are two basic reproduction ratios, one is for the infection caused by the sensitive strain, and the other one is caused by the resistant strain.
Let us define parameters The square root of R 0a and R 0A are the basic reproduction ratios of the sensitive and resistant infection respectively.These ratios represent the expected numbers of secondary cases of the sensitive (resistant) infection per primary case of the sensitive (resistant) infection in a 'virgin' population during the infectious period of primary case [10].The ratios also determine whether the sensitive strain or the resistant one will eventually fix in population after an initial infection occurs.The parameter R 0a is similar to the basic reproduction ratio R S in [14].Here, the basic reproduction ratio of the sensitive infection uses a square root sign, because the model incorporates the vector explicitly.So, it gives a finer result.
The parameter R 0j is a strictly decreasing function with respect to the treatment fraction f j , where j ∈ {a, A}.The increase (decrease) of f j will decrease (increase) R 0j linearly.Hence, increasing treatment fraction f j can be used as a strategy to decrease the parameter R 0j .
Let us define r 0a = βHa λV τa δV σHaU (δV +τa) .The square root of r 0a is the basic reproduction ratio for malaria infection in the absence of a resistant strain and treatment in the population.Moreover, we have the relation From Equation ( 4) we obtain that the treatment reduces the basic reproduction ratio of malaria infection in the absence of a resistant strain and treatment.In the absence of a resistant strain, there is a critical sensitive-treatment proportion The critical parameter ( 5) is useful in disease eradication, that is, when f a is greater than f crit a , then the disease will eventually disappear from the population, because R 0a < 1.In a case where the population is only infected initially by a resistant strain, a similar argument holds.
Model (2) has three boundary equilibria, these are, the disease-free equilibrium E d = (1, 0, 0, 0, 0, 0, 0, 0, 0), the sensitive infection equilibrium E a = (S a H , I a HaT , I a HaU , 0, 0, E a V a , I a V a , 0, 0), and the resistant infection equilibrium E A = (S A H , 0, 0, , where The equilibrium E d always exists.Meanwhile, the equilibria E a and E A exist if and only if R 0a and R 0A is greater than one respectively.Moreover, if E j exists, we have the relation .
The last equation gives the relation between the square of the basic reproduction ratio and the equilibrium states of the susceptible host and vector subpopulations.
We obtain that the greater the square of basic reproduction ratio, the less the susceptible host or vector subpopulations in the equilibrium state.
For R 0a = R 0A > 1 there is a one-dimensional manifold of interior equilibria.The manifold contains endemic coexistent equilibria.If we parameterize this manifold with m, which corresponds to I V A , then one of the equilibria is , where When the coexistent equilibrium E c exists, it depends on the parameters and initial condition.
The parameter m satisfies .
Moreover, the lower and upper bounds of m is given by Consider the limiting values of parameter m.If m = 0, then the equilibrium E c coincides with the equilibrium E a .And, if then the equilibrium E c is the equilibrium E A .So, the manifold connects the equilibria E a and E A .This will be confirmed with numerical simulation in the next section.
The following proposition gives the stability criteria of the disease-free equilibrium Proof.The local stability of the disease-free equilibrium E d is a corollary of Theorem 2 in [11].
In proving the global asymptotical stability of It can be seen that V (x) ≥ 0 for x ∈ Ω.Moreover, V (x) = 0 if and only if x ∈ {(S H , 0, 0, 0, 0, 0, 0, 0, 0) : 0 ≤ S H ≤ 1}.The last set is a positive invariant set under the flow generated by the vector field of model (2).In this positive invariant set, we have dSH dt = (δ H + η) (1 − S H ) and S H → 1 as t → ∞.Moreover, the derivative of V along solution curves of model ( 2) is given by which is less or equal to zero in Ω.
Let M be a subset of Ω, where M satisfies V = 0. Thus, the set M is defined by where j ∈ {a, A}.
From the inspection of model ( 2), it can be seen that {E d } is the largest invariant set under the flow generated by the vector field of model ( 2) that is contained in M .Therefore, by Theorem VIII of [15] the equilibrium E d is globally asymptotically stable in Ω for R 0a , R 0A ≤ 1.
From the biological point of view, Proposition 1 implies that malaria disease caused by the sensitive and resistant strains will eventually disappear from any initial size of population when R 0a , R 0A ≤ 1.
Let the equilibrium E a = (S a H , I a HaT , I a HaU , 0, 0, E a V a , I a V a , 0, 0) exists, that is if R 0a > 1.By using the coordinate system (S H , I HaT , I HaU , E V a , I V a , I HAT , I HAU , E V A , I V A ), the linearization of model ( 2) at point E a gives the Jacobian matrix where . The stability of equilibrium E a is determined by the eigenvalues of matrices A 1 and A 2 .
The eigenvalues of matrix A 1 are somewhat difficult to obtain analytically, because they relate to a quintic characteristic polynomial.But, numerical simulation with the parameter values which we consider indicates that all eigenvalues of the matrix have negative real parts if R 0a > 1.
Next, the matrix −A 2 is an M-matrix.The real part of all eigenvalues of matrix −A 2 is positive if and only if det(−A 2 ) > 0 (see [4]).Furthermore, all eigenvalues of A 2 have negative real parts if and only if det(A 2 ) > 0. The determinant of matrix Thus, if R 0a > R 0A , then the equilibrium E a is locally asymptotically stable.And, it is unstable if R 0a < R 0A .Hence, we obtain the following property.
Suppose that the equilibrium E a exists.If R 0a > R 0A , then the equilibrium E a is locally asymptotically stable.Furthermore, if R 0a < R 0A , then equilibrium E a is unstable.
By observing the symmetry property that exists between the coordinate of equilibria E a and E A , we can get a similar property as the following.Suppose that the equilibrium E A exists.If R 0a < R 0A , then the equilibrium E A is locally asymptotically stable.Moreover, if R 0a > R 0A , then the equilibrium E A is unstable.
For R 0a = R 0A > 1, the stability of equilibrium E c is difficult to obtain analytically.The linearization of model ( 2) at equilibrium E c produces a nine-degree characteristic polynomial.But, numerical exploration with parameter values which we consider indicates that the equilibrium E c is locally asymptotically stable.This can be seen in the next section.
The behavior of the equilibria can be summarized on a bifurcation diagram as depicted in Figure 2.
Figure 3 illustrates the typical dynamics of the infected host subpopulations.In case R 0A > R 0a (the left figure), we see that the resistant strain fixes (the resistantinfected host subpopulations approach some positive numbers) and the sensitive one vanishes in the population after a long period of time.Furthermore, for R 0a > R 0A (the right figure) the sensitive strain fixes and the resistant one eventually declines to zero in the population.In Figure 4 we give the typical dynamics of the infected host subpopulations for R 0A = R 0a > 1.The figure shows the projection of dynamics of model (2) on planes I HaT − I HaU and I HAT − I HAU .For this simulation, we use some fixed parameters and eleven initial conditions.We observe in this case that both of the resistant and sensitive strains approach some positive numbers, that is, both strains coexist in the population.This result confirms the existence of one-dimensional manifold of equilibria which connects the sensitive infection equilibrium E a (represented by the • sign) and the resistant-infection equilibrium E A (represented by the sign).It also confirms the stability of the manifold.
Let γ HAT = γHAU m , where m is the intermediate stage of drug tolerance and 0 < m ≤ 1.The greater of m, the more tolerant the parasites to the drug, so the drug kills less parasites.If m = 1, the parasites has been fully tolerant to the drug, so the drug kills no parasites.Figure 6 also illustrates the effect of varying the treatment fraction f A and the stage of drug tolerance m on the basic reproduction ratio √ R 0A .The left graph illustrates the effect in the three-dimensional graph.Meanwhile, the right graph is the contour plot of the left graph.From the right graph, we observe that for any fixed ratio √ R 0A , increasing f A will increase parameter m.It means that the parasites are more tolerant to the drug, so the drug kills less parasites.
Figure 7 illustrates the effect of varying f a or f A on the infected host population around the coexistent case.The initial condition for these simulations are the coexistent equilibrium, except for the evolution of coexistent case.From this result, varying f a or f A can destroy the existence of coexistent equilibrium.It should be careful in changing f a or f A , because it will determine which strain will eventually fix in the population.

Conclusion.
In this paper we formulate a deterministic model for the spread of anti-malarial drug resistance in a population.Here, we consider that the Plasmodium spp.parasite is partially resistant to an anti-malarial drug.The model incorporates the vector population explicitly in the model.In the model, both of the host and vector populations are constant for all time.The model is suitable with malarial situations in developing countries, where there is no parasite screening prior to a treatment.It also allows some treatment fractions for the infective subpopulations.Although the model is constructed for malaria disease, it can be used for some other diseases with some adaptations.
For the model, we obtain parameters R 0a and R 0A as in equation (3).These parameters correspond to the basic reproduction ratios, one for the sensitive strain infection ( √ R 0a ), the other one for the resistant strain infection ( √ R 0A ).These ratios determine the existence and the stability of the equilibria of the model.The stability of the equilibrium corresponds to which malarial strain will eventually fix in the population.
When both R 0a and R 0A are less than or equal to one, both of the strains can be eradicated from any initial size of population.If one of them is greater than one, then the strain which has the greater ratio will fix in the population.The coexistence of the sensitive and resistant strains in the long term can happen if and only if the two ratios are equal and they are greater than one.If this coexistent case occurs, the initial condition determines the equilibrium state in the long term.
Figure 6 illustrates the effect of varying the resistant-treatment fraction and the intermediate stage of drug tolerance on the basic reproduction ratio √ R 0A .We observe that for any fixed ratio √ R 0A , increasing the fraction will increase the intermediate stage, so the drug kills less parasites.support from DIKTI through BPPS scholarship and after that from Universitas Indonesia during the preparation of the manuscript.He wants to thank the editorial team of MBE for their proof reading on the manuscript.He also thanks Desiree and Huberta for their continuous support.

Figure 1 .
Figure 1.A transmission diagram of an anti-malarial drug resistance.

Figure 2 . 4 .
Figure 2. Bifurcation diagram.The scripts g, s, and u stand for globally asymptotically stable, locally asymptotically stable, and unstable respectively.

Figure 5 Figure 5 .Figure 6 .
Figure 5.The basic reproduction ratio √ R 0A as the intermediate stage of drug tolerance m and treatment fraction f A vary.

Figure 7 .
Figure 7. Dynamics of the infected host population.The solid curve corresponds to the coexistent case (f a = 0.9, f A ≈ 0.98613).

Table 1 .
Parameters used in Figure