Mathematical modelling for malaria under resistance and population movement

1. Introduction

Malaria is a hematoprotozoan parasitic infection transmitted by certain species of anopheline mosquitoes. Four species of plasmodium commonly infect to humans, but one, Plasmodium falciparum is the most lethal in humans, causing many deaths per year. Malaria also provides an unbalance that impairs the economic and social development of certain zones of the planet [17]. In reviewing history, control programs have been focused in two directions: control of the anopheles mosquito through removal of breeding sites, use of insecticides, prevention of contact with humans (by using of screens and bed nets), and use of antimalarial drug (or effective case management) [32]. Unfortunately, the implementation of these control mechanisms has not been entirely effective. Amongst the reasons we can mention: a) resistance of the malaria parasites to antimalarial drugs such as chloroquine and sulfadoxine–pyrimethamine. In this case, and from a mathematical point of view, Aneke in [5] describes the phenomenon of antimalarial drug resistance in a hyperendemic region by a model of ordinary differential equations (ODEs). Esteva et al. in [14] present a deterministic model for monitoring the impact of antimalarial drug resistance on the transmission dynamics of malaria in a human population. Tchuenche et al. in [29] formulate and analyze a mathematical model for malaria with treatment and three levels of resistance in humans incorporing both, sensitive and resistant strains of the parasites. Agusto in [1] formulates and analyzes a deterministic system of ODES for malaria transmission incorporating human movement, as well as the development of antimalarial drug resistance in a multipatch–type system. Other works to underline in this topic are [19], [6], [24]. b) The use of pyrethroid insecticides (a man-made pesticides similar to the natural pesticide pyrethrum) in malaria vector control. Here we can find the work of Luz et al. in [22] in which a model of the seasonal population dynamics of Aedes aegypti, both to assess the effectiveness of insecticide interventions on reducing adult mosquito abundance, and to predict evolutionary trajectories of insecticide resistance. In addition, Aldila et al. formulate and analyze a mathematical model for transmission of temephos resistance in Aedes aegypti population [2], meanwhile in the works [3], [16], the authors treat the insecticide resistance in general cases. c) The population migration problem. The movement of infected people or infected mosquitos from areas where malaria is still endemic to areas where the disease had been eradicated led to resurgence of the disease, and this situation also results in a increasing of resistance to insecticides and antimalarial drug [9]. With respect to migration problem, the works have been addressed through multipatch–type models see for instance [15], [26], [1]. Migration problems for dengue virus and other general epidemic models have been reviewed in [18], [8] and [20], [33], [7], [23], [10], respectively.

As far as we know, mathematical models considering resistance to antimalarial drug and insecticides and movement of populations simultaneously, as factors that hinder the malaria control, do not exist. This paper present a response to this situation, including numerical experiments that allow us to verify the feasibility of our theoretical results.

In this paper we propose two mathematical models for the malaria transmission dynamics and whose equations are based in [27]. More precisely, in the first model, we consider the interaction between humans and mosquitoes inside a patch when the human population is resistant to antimalarial drug and mosquitoes population is resistant to insecticides. Existence and stability conditions for the equilibrium solutions in terms of the local basic reproductive number are determined. For the second model, human–mosquitoes population movements in two patches is considered under the same conditions established in a patch and also following the ideas from [20]. Besides, by incorporating the use of antimalarial drugs and insecticides as control strategies, we formulate an optimal control problem for the disease.

2. One patch model

In this section we consider a single patch with a susceptible–infected–recovered (SIR) structure for humans and a susceptible–infected (SI) structure for mosquitoes. In order to present the complete model, we describe the dynamic equations that form our model as follows: let us denote as S_h (t), I_h (t) and R_h (t) the number of susceptible, infected, and recovered humans at time t, respectively. The total human population at time t is denoted by N_h (t) = Sh(t) + I_h (t) + R_h (t). Similarly, let us denote as S_v (t) and I_v (t) the number of susceptible, and infected mosquitoes at time t, respectively. The total mosquito population at time t is denoted by N_v (t) = S_v (t) + I_v (t).

Moreover, from [27], we define the force of infection for humans by β_hϵ , where β_h represents the probability of a human being infected by the bite of an infected mosquito, and ϵ represents the per capita biting rate mosquitoes. Similarly, we define the force of infection for mosquitoes as β_vϵ , where β_v represents the probability of infection of mosquito by contact with infected humans.

Respect to susceptible humans population, it is increasing due to recruitment at a constant rate Λ _h and by recovered humans from infection, which are represented by the term ωR_h . Simultaneously, this population decrease due to infection by contact with infected mosquitoes through the term β_hϵ S_h and by natural death through the term µ_hS_h . Thus, the ODE that represents the variation of the susceptible humans population is

where the symbol · corresponds to the derivative in time, i.e, = S_h (t). Now, respect to the infected humans population, it is treated with drug at a constant rate ξ ₁ θ ₁, where ξ ₁ is the drug efficacy and θ ₁ is the recovery rate due to the drug. Besides, the number of infected individuals resistant to the drug (by selective pressure) is ξ ₁ θ ₁ q ₁ I _h, where q ₁ ∈ [0, 1] represents the resistance acquisition ratio to the drug. Thus the term ξ ₁ θ ₁(1−q ₁)I_h represents the proportion of sensitive individuals to the drug. Additionally, a proportion of infected individuals recover spontaneously at a rate δ (by action of the immune system), others die from infection at a rate ρ and others from natural death at a rate µ_h . thus, the equation for the variation of the infected humans population is given by

Finally, in our model the recovered humans population increase by the action of the drug and by spontaneous recovery, and decrease as consequence of natural death and loss of immunity. Thus, the variation of the recovered humans population in time is described by

On the other hand, the description for the SI model is the following: the susceptible mosquitoes population is recruited at a constant rate Λ _v . It is diminished by infection due to contact with infected humans, which is described through the term β_vϵ S_v . Simultaneously, it is reduced due to natural death with a rate µ_v and by action of insecticides at a rate ξ ₂ θ ₂, where ξ ₂ represents the efficacy of insecticide and θ ₂ is the death of mosquitoes due to insecticides. The number of mosquitos resistant to the insecticides is ξ ₂ θ ₂ q ₂, where q ₂ ∈ [0, 1] represents the resistance acquisition ratio to the insecticides. Thus, the expression ξ ₂ θ ₂(1 − q ₂) represents the proportion of sensitive mosquitos to the insecticides. Then, the system describing the variation of the mosquitoes population in time is

In summary, from (1)-(2), our model for malaria under resistance in a patch is given by

where (N _h (0), N _v (0)) denotes a initial condition and N _h and N _v are vectors formed by S_h , I_h , R_h and S_v , I_v , respectively.

Remark 2.1. The novelty in this work involves the parameters ξ_i , θ_i and q_i with i = 1, 2. Their interpretation and values are given in Tables 1 and 2 from Section 2.2. A complete description and interpretation of the others parameters involved in the model (3) can be found in [27].

Now, a set of biological interest for the solutions of the system (3) is defined as follows:

The following lemma establishes the invariance property for Ω.

Lemma 2.2. For (N _h (0), N _v (0)) a non–negative initial condition, the system (3) has a unique solution and all state variables remain non–negative for all time t ≥ 0. Moreover, the set defined on (4) is positively invariant with respect to the system (3).

Proof. Since the vector field defined on the right side of (3) is continuously differentiable, the existence and uniqueness of the solutions is fulfilled. On the other hand,

Thus

Multiplying both sides of the above inequality by the integrating factor e^µhτ and integrating from 0 to t, we obtain that

from where

Similar calculation shows that N_v (t) → as t → ∞. Thus, the region Ω is positively invariant. This completes the proof.

2.1. Qualitative análisis

In this subsection, we first compute the local basic reproductive number associated to the system (3) . Afterward, conditions for existence and stability of the equilibrium solutions are developed.

Local basic reproductive number

It is well known that a disease–free equilibrium (DFE) is a steady state solution of a system where there is no disease, in our case, S_h = > 0, S_v = > 0, and all others variables I_h , I_v , R_h are zero. It will be denoted by E _0one = (, 0, 0, , 0), where

Since the basic reproductive number, commonly denoted by R ₀ (but in this case denoted by R _0one ) is the average number of secondary infective generated by a single infective during the curse of the infection in a whole susceptible population, it is a threshold for determining when an outbreak can occur, or when a disease remains endemic. Using the next generation operator method [30] on the system (3), the Jacobian matrices F _one and V _one evaluated in the DFE are given by

And

Thus, the next generator operator of model (3) is given by

It follows that the local basic reproduction number of the system (3), denoted by R _0one ,

Existence of endemic equilibria

In this subsection, conditions for existence of endemic equilibria of the model (3) are studied. First of all, the existence of the DFE, denoted by E _0one , is guaranteed as consequence of the previous subsection. Now, in order to analyze the endemic equilibria of the model (3) we consider the solutions to the algebraic equation system

Let us define

Thus, after some algebraic manipulations of the system (7), we obtain the following expressions for S_h , R_h , S_v and I_v in terms of I_h :

and the following cuadratic equation for I_h :

From (11), we have that the coefficients a and b are non–negatives, while c ≥ 0 if ≤ 1, otherwise c < 0. Thus, the polynomial P(I_h ) = + bI_h + c has only one sign change, and by the Descartes’ rule of sign [4] it has one or zero positive roots. This result is summarized in the following theorem.

Theorem 2.3. For the model (3) the DFE contained in Ω always exists. Additionally,

1. If R _{0 one} ≤ 1, there are not endemic equilibria.

2. If R _0one > 1, there exist one endemic equilibrium.

Stability análisis

In this subsection, we proof the stability of the equilibrium solutions of the system (3) given on Theorem 2.3. First, using the linearization of the system (3) at the DFE, we proof it local stability, which is determined by the sign of the real part of the eigenvalues of the Jacobian matrix denoted by J(E _{0 one} ), which is given by

It is easy to determine the eigenvalues of J(E _{0 one} ), namely: η ₁ = −µ_h , η ₂ = −(ω + µ_h ) and η ₃ = − [ξ ₂ θ ₂(1 − q ₂) + µ_v ], while the others eigenvalues are given by the roots of the quadratic equation

where

From above, the coefficients a0 and a1 are positive, while the sign of the coefficient a ₂ depends of R _{0 one} . From the Routh–Hurwitz criterion [13], we can guarantee that the quadratic equation (13) has roots with negative real part if, and only if, its coefficients are positive and the following determinants called minors of Hurwitz are positive:

We verify that ∆₁ > 0 and ∆₂ > 0 if, and only if, R _{0 one} ≤ 1. In consequence, when R _{0 one} ≤ 1 the DFE is a locally asymptotically stable (LAS) equilibrium point of the system (3).

Now, we are going to proof the stability of the endemic equilibrium of the system (3). For this end, we use results based on the center manifold theory described in [11] to show that the system (3) exhibits a forward bifurcation when R _0one = 1, or equivalently, when

The eigenvalues of the Jacobian matrix given in (12) evaluated in (E _{0 one} , β^∗ ) are: 0 and −µ_h , −(ω + µ_h ), − [ξ ₂ θ ₂(1 − q ₂) + µ_v ] and

where the last four have negative real part. In consequence, in β^∗ , the DFE is a non– hyperbolic equilibrium. Let W = (w ₁, w ₂, w ₃, w ₄, w ₅) ^T a right eigenvector associated to the zero eigenvalue, which satisfies J(E _0one , β ^∗ )W = 0W = 0, or equivalently,

The vectorial form for the solutions of above linear system is given by

where the parameter α is defined in (8). Similarly, a left eigenvector V = (v ₁, v ₂, v ₃, v ₄, v ₅) of the matrix J(E _0one , β^∗ ) associated to the zero eigenvalue satisfies V J(E _0one , β^∗ ) = 0V = 0, or equivalently, v ₁ = v ₃ = v ₄ = 0 and

from here,

The values for w ₅ and v ₅ such that W · V = 1, are

Thus, the coefficients and given on Theorem 4.1 from [11],

can be explicitly computed as follows. Let us denote as f_i , i = 1, ..., 5 to the scalar functions of the right hand of the system (3), and x ₁ = S_h , x ₂ = I_h , x ₃ = R_h , x ₄ = S_v , x ₅ = I_v . The coefficients w_p and v_p , with p = 1, 2, ...5 of (18), represent the components of the eigenvectors W and V defined on (15) and (16), respectively. After some calculations we have that the second order partial derivatives evaluated in (E _0one , β^∗ ) are given by

In the above expressions we did not consider the zero and cross partial derivatives. Ad[1]ditionally, the second order partial derivatives with respect to the bifurcation parameter β^∗ evaluated in E _0one are all zero, except

Thus, the coefficients and given on (18) can be expressed by

From (19) we have that > 0, while the sign of a˜ depends of the sign of w ₂, v ₂(w ₂ + w ₃) and w ₂(w ₁ + w ₃). From (15) and (16) we verify that w ₂ ≥ 0, v ₂(w ₂ + w ₃) ≥ 0, and

Thus, from [11, Theorem 4.1], the endemic equilibrium is LAS when R _0one > 1, which suggest the global stability of the DFE. The previous results are summarized in the following theorem.

Theorem 2.4. If R _0one ≤ 1 the DFE is LAS in Ω, and the endemic equilibrium is unestable. If R _0one > 1, the DFE becomes an unstable hyperbolic equilibrium point, and the endemic equilibrium is LAS in Ω.

Figure 1 shows the bifurcation diagram.

Theorem 2.5. If ≤ 1, then the DFE is globally asymptotically stable (GAS) in Ω.

Proof. From Theorem 2.4, when ≤ 1, the DF E is LAS in Ω. Let (N _h (t), N _v (t)) a positive solution of the system (3); then, by Lemma 2.2 it satisfies

We prove the existence of a Lyapunov function for the traslated system = f(y + E _0one ) − f(E _0one ) = F(y), where f is the vectorial field defined from right hand of the system (3) and y = 0 is a trivial solution of the system = F(y). Let us consider the function

and let

The function V defined in (21) satisfies the following properties:

(P1) V (, 0, 0, , 0) = V (E0 _one ) = V ∗ (0) = 0.

(P2) V > 0 ∀(, ) E _0one in Ω (V is positive definite).

(P3) The orbital derivative of V along the trajectories of (3) is negative definite. In fact,

Thus, the DFE is globally stable in Ω. To verify its global asymptotic stability, let us consider △ = {(N _h , Nv) : V˙ = 0}. Then △ ⊂ {(N _h , N _v ) : I_h = 0}. Let △′ ⊂ △ be the biggest invariant set with respect to (3) and (N _h , N _v ) a solution of (3) in △′ ; then, (N _h , N _v ) is defined and is bounded ∀t ∈ and I_h (t) = 0 in △′ for all t. Replacing this value in the system (3) we obtain that R_h (t) = I_v (t) = 0 for all t, while, from the first and fourth equation of (3) we obtain that S_h ( _t ) = Λ _h /µ_h = and S_v (t) = Λ _v / (ξ ₂ θ ₂(1 − q ₂) + µ_v ) = . Thus, △′ = {E _0one }, and from the Lasalle invariance principle [31] E _{0 one} is GAS in Ω.

2.2. Numerical experiments

In this subsection, we validate our theoretical results with numerical experiments. To this end, we take data from rural areas of Tumaco (Colombia) reported in [27] and make some numerical simulations. For the values of the parameters corresponding to insecticides, we assume that the fumigation is done with two pyrethroids insecticides (deltamethrin and cyfluthrin) according to the recommendations of Palomino et al. in [25]. Pyrethroids insecticides are a special chemicals class of active ingredients found in many of the modern insecticides used by pest management professionals. Due to the low concentrations in which these products are applied, a constant safety of use and a decrease in the toxic impact on vector control have been achieved. For the values of the parameters corresponding to the drug, we assume that the infected patients are treated with artemisinin–based combination therapy (ACT) according to the recomendations of Smith in [28]. Artemisinin (also called qinghaosu), is an antimalarial drug derived from the sweet wormwood plant: Artemisia annua. Fast acting artemisinin–based compounds are combined with other drugs, for example, lumefantrine, mefloquine, amodiaquine, sulfadoxine/pyrimethamine, piperaquine and chlorproguanil/dapsone. The artemisinin derivatives include dihydroartemisinin, artesunate and artemether [28]. Tables 1 and 2 show the values of the parameters corresponding to the drugs and insecticides supply, respectively.

Figure 2 shows the behavior of human and mosquito populations when the patients are treated with ACT and the mosquitoes are fumigated with cyfluthrin and deltamethrin, respectively. In Figure 2, the solutions tend to an endemic equilibrium and = 2.15 (left), while (Figure 2 (right)) the solutions tend to the DFE and = 0.0012. In fact, given that cyfluthrin is an insecticide with less efficacy than deltamethrin, its application generates greater resistance hindering the disease control.

In Figure 3 we consider the effects of resistance in the population dynamics. In Figures 3 (a) and (b) we assume that there is no resistance (q ₁ = q ₂ = 0). Then, when the fumigation is done with cyflutrin, = 1.41 and the solutions tend to the endemic equilibrium (8136, 227, 1541, 250, 32), which evidences a considerable reduction in the persistence of the infection, while if the fumigation is done with deltamethrin, = 0.00095 and the solutions tend to the DFE. In Figures 3 (c) and (d) we assume total resistance (q ₁ = q ₂ = 1). Then, when the fumigation is done with cyflutrin, = 2724.4 and the solutions tend to the endemic equilibrium (0, 789, 0, 0, 4699), which evidences that after the first 30 days, all individuals (humans and mosquitoes) will be infected, while if the fumigation is done with deltamethrin, = 264.8 and the solutions tend to the endemic equilibrium (6, 824, 0, 195, 4617), which evidences a persistence of the infection.

3. Two patch model

In this section we model the malaria transmission dynamics between humans and mosquitoes within a patch and their spatial dispersal between two patches. Within a single patch, our model is defined by the equations (3), where the subscripts 1 and 2 refers to patch 1 and patch 2, respectively. The patches are coupled via the resident budgeting time matrix R = [λ _ij ]_2×2 for i, j = 1, 2 as in [20]. Here λ _ij α_ij + β_ji , being α_ij the probability of a human from patch i is visiting the patch j and β_ji the probability of a mosquito from patch j, is visiting the patch i. Some authors prefer not to consider the mobility of mosquitoes due to yours short life cycle (less than two weeks without captivity), in which case we assume β_ji = 0. Each λ _ij is a constant in [0, 1] and λ _ij = 1 for i = 1, 2. In this model we include bi–directional motion as in [20], that is, a susceptible human (mosquito) in patch i can be infected by an infected mosquito (human) from patch i as well as by an infected mosquito (human) from patch j who is visiting the patch i. Thus, the dynamic in two patches are represented through the following system of nonlinear ODEs:

where (N _h1 (0), N _v1 (0), N _h2 (0), N _v2 (0)) denotes an initial condition. Let us define N _H (t) = N _h1 (t) + N _h2 (t), N _V (t) = N _v1 (t) + N _v2 (t) and

A set of biological interest for the solutions of the system (22) is

Note that the proof of invariance of can be be made using the results of Lemma 2.2.

3.3. Global basic reproductive number and numerical experiments

In this subsection, we first compute the global basic reproductive number associated to the system (22). Then, we obtain numerical experiments to generate an application of the mathematical model (22) using data from [27]. Let us denote as E₀ = (, 0, 0, , 0, , 0, 0, , 0) with

to the DFE associated to the system (22).

Using a similar procedure to the presented in subsection 2.1 for

and V = (V1|V2), where

we get the following expression for the global basic reproductive number:

where

Considering the uncopling system (that is, λ₁₁ = 1 and λ₂₂ = 1) in (26), we obtain the local basic reproductive number for each patch given in (6).

In what follows, we make some numerical experiments. For this purpose, we are going to consider the following hypothesis: (a) patch 1 and patch 2 represent rural areas (RA) and urban areas (UA) from the municipality of Tumaco (Colombia), as in [27], respectively. (b) The epidemiological outbreak begins in RA and the individuals in UA acquire the infection due to the coupling between the two patches. Therefore (unless otherwise stated), the initial condition will be S_{h 1} (0) = 100000, I_{h 1} (0) = 30000, R _h1 (0) = 20000, S _v1 (0) = 50000, I _v1 (0) = 10000, S _h2 (0) = 100000, S _v2 (0) = 50000 and all others zero. (c) Mosquitos are fumigated only with cyflutrin (data from Table 2). (d) The infected patient are treated with ACT (data from Table 1). (d) The resistance acquisition ratio in RA is higher than UA due to in RA individuals are continuously exposed to the parasite, that is, q ₁₁ = 0.1, q ₁₂ = 0.09, q ₂₁ = 0.05 and q ₂₂ = 0.04. Besides, we will consider the following coupling scenarios poposed by Lee et al. in [20]:

(S1) Uncoupled: when there are no visits between patches, that is, λ₁₁ = λ₂₂ = 1 and others are equal to zero.

(S2) Weakly–coupling: small values for λ₁₂ and λ₂₁.

(S3) Strongly–coupling: when visitors from patch 2 spend quite an amount of time in patch 1, that is, λ₂₂ < λ₁₁.

Table 3 shows the values of the parameters in the residence–time matrix considering different scenarios of coupling.

Figure 4 shows the behavior of the solutions when the system (22) is uncoupled. If the disease begins in patch 1, the disease does not spread to patch 2.

Figure 5 shows the behavior of the humans and mosquitoes populations in patches 1 and 2, respectively, considering weakly–coupling. Here, the disease is spread from patch 1 to patch 2 during the first 50 days, then the disease is eliminated in patch 2, and remains at low load in patch 1. On the other ha hand, the strongly–coupling scenario is illustrated in Figure 6. Here, the disease is spread from patch 1 to patch 2 during the first 100 days and the infection persists in both patches. After 100 days, the disease is eliminated in both patches.

3.4. Local sensitivity analysis of parameters

In this subsection we determine the sensitivity indices of the parameters to the R ₀, considering strongly– coupling and data from [27] . The sensitivity indices are computed through the normalized forward sensitivity index [12], which allow us to measure the relative change of the variable R ₀ when a parameter changes. When the variable is a differentiable function of the parameter, the sensitivity index may be alternatively defined using partial derivatives [12]. If we denote the variable as u which depends on a parameter p, the sensitivity index is defined by . Given the explicit formula for R ₀ in (26), we determine an analytical expression for the sensitivity indices of R ₀ with respect to each parameter that comprise it. In Table 4 we show the values of the sensitivity indices, where P1 and P2 mean patch 1 and patch 2, respectively.

From Table 4, in both rural (patch 1) and urban (patch 2) areas, R ₀ is more sensitive to the parameters corresponding to recovery rate due to the drug θ _1i with i = 1, 2 and death rate due to the insecticides θ _2i with i = 1, 2. An interpretation of these indices is given as follows: in RA, given that Γ _θ11 = −0.54, increasing (or decreasing) θ ₁₁ in 10% implies that R ₀ decreases (or increases) in 5.4%. An analogous reasoning can be made for the others sensitivity indices. The information provided by the sensitivity indices to the R ₀, will be used in the next section, in which we will propose some control strategies for the malaria disease.

4. Optimal control problema

In this section an optimal control problem applied to the model (22) is formulated. Here, we are going to consider that the parameters corresponding to recovery rate due to the drug and death rate due to the insecticides θ_ij with i, j = 1, 2 will be the controls, therefore they will be functions depending on time. The first objective will be to minimize a performance index or cost function by the use of drugs and insecticides. For this purpose, we assume that θ _i1 with i = 1, 2 and θ _j2 with j = 1, 2 are the controls by drugs and insecticides, respectively, which assume values between 0 and 1, where θ_ij = 0 is assumed if the use of drugs (or insecticides) is ineffective and θ_ij = 1 if the use of drugs (or insecticides) is completely effective, that is, all individuals recover with medication and all mosquitoes die with insecticides. In this sense, for i and j fixed, the control variable θ_ij (t) provides information about amount of drug or insecticides that must be supplied at time t.

The second objective will be to minimize the number of infected humans and infected mosquitoes in each patch. For this purpose, the following performance index or cost function is considered:

where θ = (θ ₁₁, θ ₂₁, θ ₁₂, θ ₂₂) is the vector of controls, c ₁ and c ₂ represent social costs, which depend on the number of individuals with malaria and the number of mosquitoes with the parasite, and defines the absolute costs associated with the control strategies, such as, implementation, ordering, distribution, marketing, among others. For calculation purposes, we will denote the integrand of the performance index given on (27) as

where X represents the vector of states.

With the above considerations, the following control problem is formulated:

In above the formulation, we assume an initial time t ₀ = 0, a final time T fixed which represents the implementation time of the control strategies, free dynamic variables X ₁ in the final time, and the initial condition X ₀ being a non–trivial equilibrium of the system (22). Additionally, we assume that the controls are in a set of admissible controls U which contains to all Lebesgue measurables functions with values in the interval [0, 1] and t ∈ [0, T].

4.5. Existence of an optimal control

In this section, we use the classic existence theorem proposed by Lenhart and Workman [21] to prove the existence of an optimal control θ ^∗ for the formulation (29). Let U = [0, 1]⁴ be the set where θ assumes its values (set of controls), and f(t, X, θ) the state equations of the right side of (29). To guarantee the existence of optimal controls, hypotheses (H1) to (H5) from [21] must be verified, that is,

(H1) (a) |f(t, 0, 0)| ≤ C,

(b) |f _X(t, X, θ)| ≤ C(1 + |θ|),

(c) |f _θ(t, X, θ)| ≤ C.

(H2) The set of controls U is convex.

(H3) f(t, X, θ) = α(t, X) + β(t, X)θ.

(H4) The integrand of the performance index f ₀(t, X, θ) defined in (28) is convex for θ ∈ U.

(H5) f ₀(t, X, θ) ≥ c ₁|θ|^b − c ₂ with c ₁ > 0 and b > 1.

We will prove the hypothesis (H1)(a) and (H5), since the others are obvious. To this purpose, we need the following results.

Lemma 4.1. Let

Then,

where Λ _H , Λ _V , µ_H and µ_V are defined on (23), and f_θ (t, X, θ) is the matrix obtained by differentiating the state equations of the right side of the system (29) with respect to θ, whic is given by

Proof. Computing the Euclidean norm of matrix (32), we obtain

Lemma 4.2. The integrand of the performance index satisfies

Proof.

Remark 4.3. Note that the hypothesis (H5) is satisfied by taking b = 2, c ₂ = 0 and c ₁ = 1/2 min {d ₁, d ₂, d ₃, d ₄} in the last expression of (33).

4.6. Deduction of an optimal solution

In this subsection, the Pontryagin Principle for bounded controls [21] is used to compute the optimal controls of the problem (29). First, let us observe that the Hamiltonian associated to (29) is given by

where Z = (z ₁, z ₂, ..., z ₁₀) is the vector of adjoint variables which determine the adjoint system, and

The adjoint system and the state equations of (29) define the optimal system. The main result of this section is summarized in the following theorem.

Theorem 4.4. There are an optimal solution X ^∗ (t) that minimize J in [0, T], and an adjoint vector of adjoint functions Z such that

with transversality condition Z(t) = 0 and the following characterization of the controls:

Proof. The Pontryagin Principle guarantees the existence of the vector of adjoint variables Z whose components satisfy

Thus, the derivatives of the adjoint variables are

Replacing the derivatives of H with respect to the state equations in above equalities, we obtain the system (35). Additionally, the optimality conditions for the Hamiltonian are given by

so that

In consequence, satisfies

or equivalently,

Using a similar reasoning for , and , we obtain the characterization (36). Therefore the proof is complete

5. Numerical experiments

In this section we present some numerical simulations associated with the implementation of drugs and insecticides as control strategies, as well as their effects on the infected individuals under uncoupled and strongly–coupling scenarios. For the simulations, we use the forward-backward sweep method proposed by Lenhart and Workman [21]. The implementation time of the control strategies will be approximately 10 days, which is the duration of a malaria treatment. The values of the relative weights associated with the control, will be those of Table 7 from [27].

Figure 7 shows the behavior of the infected individuals in patches 1 and 2 under uncoupled scenario. Due to in this scenario, the disease only remains in patch 1 and does not spread to patch 2, the density of infected individuals decreases with control in patch 1 and the effects of the controls in patch 2 are not necessary.

In Figure 8 we can see the behavior of infected individuals in patches 1 and 2 under strongly–coupling scenario. Here, the infection decreases with control in both patches, but the efforts are greater in patch 1 than in patch 2.

In both cases, uncoupled and strongly–coupling scenario, the effects of the controls are highly effective and fast to eliminate the disease in patch 1, while in patch 2 the elimi[1]nation depends of the coupling scenario.

6. Discusión

In this work, we model the malaria transmission dynamics, considering three factors that hinder its control: resistance to drugs, resistance to insecticides and population movement. To illustrate the above factors, we divide our work into two mathematical models. (a) A mathematical model in a patch under the hypothesis that the parasites are resistant to the drugs, and the mosquitoes are resistant to the insecticides. In this first model, we make a qualitative analysis of the solutions of the system, which reveal the existence of a forward bifurcation and the global stability of the DFE. From the biological point of view, the existence of a forward bifurcation indicates that the disease can be controlled by keeping the local R _{0 one} below of one. Since the expression for R _{0 one} given on (6) depends directly on the resistance acquisition ratios q ₁ and q ₂, then at lower levels of resistance acquisition, the value of R _0one decreases, which implies that the infection levels decrease. On the other hand, since R _0one depends inversely on the effects of the drugs and insecticides, then an increase in the recovery rate humans due to drugs and the death of mosquitoes due insecticides, implies a decrease of R _0one and therefore the burden of infection. The numerical experiments for this first model corroborate the theoretical results. Here, we assume that the infected patients are treated with ACT (artemisinin–based combination therapy) and to contrast the fumigation of mosquitoes with deltamethrin and cifluthrin, where the first insecticide is more effective than the second one. With total resistance to the drugs and insecticides (q ₁ = q ₂ = 1), we verify that the burden of infection persists regardless of the type of drug and insecticide used, while without resistance (q ₁ = q ₂ = 0), the burden of infection decreases with the use of deltamethrin and is maintained at low levels with the use of cifluthrin. These results are alarms in public health, because despite the pharmaceutical industry is taking care day after day to create new drugs and new insecticides, if the phenomenon of resistance acquisition is not counteracted, the problem of malaria control will be increasingly difficult, and in some cases impossible.

(b) For the model in two patches, we consider the same hypotheses of the model in a single patch, and additionally, movement of populations between two patches. For this case, we determine the global basic reproductive number R ₀, and through numerical experiments, we illustrate the behavior of the solutions when the infection starts in the patch 1 (rural areas of Tumaco [27]), and under three coupling scenarios: (1) uncoupled scenario. When there is no movement between patches, the infection remains endemic in the patch 1 and does not spread to the patch 2. (2) Weakly–coupling. If the probabilities of visiting between both patches are low, the disease is endemic in the patch 1 and remains at a very low load in the patch 2. (3) Strongly–coupling. If the probabilities of visiting between both patches is high, the disease remains endemic in both patches. These results corroborate the phenomenon of reinfection in areas where malaria has been eradicated and is not endemic, as is the case of urban malaria. Here, a new alarm in public health is created, because if malaria has been completely eradicated in a sector and is not endemic there, the movement of humans (or mosquitoes) from endemic areas can activate the infection alarm again.

Finally, using results of a local sensitivity analysis of parameters to the global R ₀, we formulated an optimal control problem by using of drugs and insecticides as control strategies. The results of the theoretical and numerical analysis of the optimal control problem reveal that under uncoupled scenario, the control is effective and necessary in patch 1 but not in patch 2, while under strongly–coupling, greater efforts are required to control the disease in patch 1 than in patch 2.

An open problem through this research is to incorporate prophylaxis as a control strategy for the disease, that is, patient education campaigns both in the use of drugs and in the use of insecticides. In this way, the resistance phenomenon will be mitigated and the control campaigns for the disease will be more effective and less expensive.

Acknowledgments

The first author is supported by the Fondecyt Postdoctoral Grant N 3180100 and ANID-Basal Project FB0008 (C. Montoya).The second one is supported by a scholarship of the Fundación Ceiba, Colombia (J-P. Romero-Leiton).

References

Author notes

cdmontoy@mat.uc.cl

Additional information

To cite this article: C. Montoya, J.P. Romero-Leiton, Mathematical modelling for malaria under resistance and population movement, Rev. Integr. temas mat. 38 (2020), No. 2, 133-163. doi: 10.18273/revint.v38n2-2020006

Revista Integración
ISSN: 0120-419X
Vol. 38
Num. 2
Año. 2020

Mathematical modelling for malaria under resistance and population movement

CristhianJhoana P. MontoyaRomero–Leiton
Pontificia Universidad Católica de ChileUniversidad CesmagDatambiente,ChileColombiaColombia