Controlling malaria in a population accessing counterfeit antimalarial drugs

1.   Introduction

Malaria, although preventable and treatable, continues to be one of the most prevalent and lethal human infections worldwide ^[1]. Globally, 627,000 deaths were recorded in 2020, with Sub-Saharan Africa accounting for almost 96% of them ^[2]. Children under 5 years become the most vulnerable to infection and death after losing their maternal antibodies, which protect them during their first six months ^[3,4].

A massive deployment of effective prevention and treatment tools by the World Health Organization (WHO) in 2000 has led to a reduction of cases in the Americas and the Western Pacific Region, no indigenous cases in the European Region, a decline of 66% in all age mortality rates in Africa and 663 million cases being averted in Sub-Saharan Africa from 2001 to 2015 ^[5,6]. During this period, it is estimated that Sub-Saharan Africa saved USfanxiexian_myfh900 million in case management costs. In 2019, the governments of the endemic countries and their international partners invested funding estimated at USfanxiexian_myfh 3.0 billion for malaria control and elimination ^[2].

Early treatment of infected humans with effective antimalarial drugs gives complete recovery and prevents the transmission of infection to mosquitoes during a blood meal, recrudescence, severe malaria and even death ^[7]. Recrudescence is common in Plasmodium falciparum infection and may occur when parasites, which remain in the red blood cells after an episode of malaria, start multiplying and cause the recurrence of the clinical symptoms due to treatment failure in the patient ^[8]. P. falciparum is the most deadly species of malaria parasite globally and the most prevalent in Africa. In most Sub-Saharan countries, after P. falciparum resistance was identified, the treatment for uncomplicated malaria was changed to artemisinin-based combination therapy (ACT) ^[7,9]. WHO recommends that an entire course of highly effective ACT must be used by both semi-immune and non-immune malaria patients to have a complete cure from both sexual and asexual forms of the parasites and partial immunity ^[5,7].

The availability, distribution, trade and use of monotherapies and other substandard antimalarial drugs continue throughout most malaria-endemic countries ^{[10,11,12,13,14,15,16]}. It is estimated that about a third of antimalarial drugs that end up in Africa are counterfeit ^[13,15]. These antimalarial drugs may contain too few or too many required active ingredients and may fail to be adequately absorbed by the body ^[17]. Recently, evidence has shown that counterfeit antimalarial drugs pose a public health threat of prolonging the illness, incomplete recovery, treatment failure, recrudescence, severe disease, spreading drug resistance and asymptomatic infection carriers, resulting in more than 122,000 deaths of African children under five years annually ^{[7,12,15,17,18]}. According to ^[4,13], the use of these drugs may jeopardize the success made so far in controlling and eliminating malaria, particularly in Sub-Saharan Africa.

Several mathematical models have considered the effects of effective treatment ^{[19,20,21,22,23]}, recovery or immunity ^{[19,24,25,26,27,28]}, reinfection ^[27,29,30] and age-structure ^{[21,22,31,32]} on the dynamics of malaria. In particular, ^{[19,24,25,26,27,28,33]} incorporated into their models recovered or semi-immune humans, who are reservoirs of infection to mosquitoes. In addition, ^[27] assumed that the recovered humans could relapse. Evidence has shown that in most malaria-endemic areas in Sub-Saharan Africa, symptomatic humans often use counterfeit and effective drugs for treatment ^{[4,9,11,12,13,14,15,16]}. The aforementioned models did not consider the effect of both effective and counterfeit antimalarial drugs on the dynamics of malaria.

Many studies have also been carried out to examine the impact of control strategies on the transmission of malaria infection. The potential impact of personal protection, treatment and possible vaccination on the transmission dynamics of malaria was theoretically assessed in ^[20]. The optimal control theory has been successfully used in decisions concerning the cost minimization of several control intervention models after implementing Pontryagin's maximum principle ^[34]. The studies ^{[23,24,35,36,37,38]} applied optimal control theory to determine the impact of control strategies on the transmission dynamics of malaria. The studies investigated optimal strategies for controlling the spread of malaria disease using treatment, treated bednets and spraying of mosquito insecticide in models with mass action ^[24], standard ^[23,37] or non-linear ^[36] incidence rates and cost-effectiveness ^[23] of the controls. Other studies ^[38] used optimal control problem strategies to study how genetically modified mosquitoes should be introduced into the environment. The aforementioned studies did not consider optimal control in age-structured models, especially where the treatment can be either effective or ineffective due to using effective or counterfeit antimalarial drugs. Optimal control has, however, been applied to an age-structured model for HIV ^[39].

This paper seeks to answer the following question: What optimal control strategy best eliminates P. falciparum malaria infection in an age-structured population where counterfeit drug use persists? We do this by developing a deterministic model for malaria transmission incorporating the infant and adult populations, counterfeit drug use and three of the malaria control measures adopted by most endemic countries in Sub-Saharan Africa ^[4,40], namely, use of highly effective antimalarial drugs (HEAs), insecticide-treated bednets (ITNs) and indoor residual spraying (IRS). The paper is organized as follows. Section 2 briefly describes the formulation of the model and its basic properties. In Section 3, the dynamics of the model is presented. The analysis of the reproduction number is done in Section 4. In Section 5, the analysis of the optimal control problem is undertaken to find the conditions for optimal malaria control using Pontryagin's maximum principle, and the numerical simulations of the model are illustrated. The discussion of results and conclusion are presented in Sections 6 and 7, respectively.

2.   Model formulation

The transmission model for malaria in human and female anopheles mosquito populations is formulated with the total population sizes at time $t$ given by $N_{h}(t)$ and $N_{m}(t)$, respectively. The human population is divided into two mutually exclusive age subgroups: infants aged 0–5 years and adults above five. Each age subgroup is divided into susceptible, infectious, counterfeit antimalarial drug users and effective antimalarial drug users epidemiological classes. The mosquito population is divided into susceptible and infectious epidemiological classes. Figure 1 shows the flows between the classes. We let $S_{A} (S_{B})$, $I_{A} (I_{B})$, $U_{A} (U_{B})$ and $T_{A} (T_{B})$ represent infants (adults) who are susceptible, infectious, counterfeit drug users and effective drug users, respectively. For the mosquito population, $S_m$ and $I_m$ denote the susceptibles and infectives, respectively. The total human population $N_{h}(t)$ at a time $t$ is given by

The total mosquito population $N_{m}(t)$ at a given time $t$ is given by $N_{m}(t)$ = $S_{m}(t)+I_{m}(t)$.

Susceptible humans are those with no merozoites or gametocytes in their bodies. Infectious humans show clinical symptoms of malaria and can infect feeding mosquitoes. Users of effective antimalarial drugs recover from malaria either naturally or due to using an effective antimalarial drug. These humans have partial immunity and become susceptible when this immunity wanes. Counterfeit drug users are humans who are removed from the infectious class due to using counterfeit drugs to treat malaria. They are asymptomatic, can infect mosquitoes (at a lower rate compared to the infectious humans) and can recrudesce into the infectious class. Adult counterfeit drug users can become susceptible due to the high level of acquired immunity resulting from their repeated exposures to malaria ^[3]. The susceptible mosquitoes have no sporozoites in their body. Infective mosquitoes can infect humans during a blood meal. It is assumed that the infective stage of the mosquitoes ends with their death due to their short life cycle. Merozoites are the parasites released into the human bloodstream when a hepatic or erythrocytic schizont bursts, and gametocytes are the sexual stages of malaria parasites that infect anopheline mosquitoes when taken up during a blood meal. On the other hand, sporozoites are the motile malaria parasites inoculated by feeding female anopheline mosquitoes to invade the human hepatocytes ^[7].

The susceptible infants class, $S_A$, is generated by the recruitment of newborns at a birth rate $\alpha_h$ and the loss of post-treatment prophylaxis by effectively treated infants at a per capita rate $\phi_A$. This class is decreased when infants in it die naturally, mature into susceptible adults or get infected, at a per capita natural death rate $\mu_h$, per capita rate $\tau$ of maturing and per capita infection rate of infants $(1 - u_{1})\rho^{n}\lambda_A$, respectively. The control $u_{1}$ represents the effort of preventing malaria with insecticide-treated bednets (ITNs), so $(1-u_{1})$ describes the failure probability of this prevention effort. The daily survival probability of a human is assumed to be $\rho_{1} = e^{-\mu_h}$. The probability of survival of a human over the average latent period of length $n_{h}$ to be infectious is given by $\rho^{n} = e^{-\mu_h n_h}$. $\lambda_A$ is the force of infection in infants.

The infectious infants class, $I_A$, is generated by the per capita rate of acquiring infection $(1 - u_{1})\rho^{n}\lambda_A$ and per capita recrudescence rate, $\xi_A$, of infants who used counterfeit drugs. This class is reduced by the per capita natural death rate $\mu_h$, per capita disease-induced death rate $\delta_A$, per capita removal rate due to the use of counterfeit drugs $\eta_A$, per capita recovery rate due to the use of a highly effective antimalarial and natural recovery $u_{2}\gamma_A+\gamma_{1}$, and per capita rate of maturing into infectious adults $\tau$ of infants in it. The control $u_{2}$ represents the treatment efforts with highly effective antimalarials and involves diagnosing patients, administering drug intake and follow-up of drug management. The recovery rate in infants due to using a highly effective antimalarial drug is $\gamma_A$, and $\gamma_1$ is the natural recovery rate in infants.

The infant class of counterfeit drug users, $U_A$, increases at a per capita removal rate $\eta_A$ of infectious infants due to counterfeit drug use and decreases at a per capita recrudescence rate $\xi_A$ into infectious infants, per capita natural death rate $\mu_h$ and per capita rate of maturing into adults $\tau$ of infants who used counterfeit malaria drugs.

The infant class of effective drug users, $T_A$, increases with the per capita recovery rate $u_{2}\gamma_A+ \gamma_{1}$ of infectious infants. This class decreases due to the loss of post-treatment prophylaxis at a per capita rate $\phi_A$, natural death at a per capita natural death rate $\mu_h$ and maturing into adults at a per capita rate $\tau$ of infants in it.

The susceptible adults class, $S_B$, is generated when the susceptible infants mature above 5 years, recovered adults lose their partial immunity and post-treatment prophylaxis at a per capita rate $\phi_B$, and counterfeit drug users become susceptible at a per capita rate $\theta$. This class is decreased by natural death and infections, with, respectively, per capita natural death rate $\mu_{h}$ and per capita rate of acquiring infection from infectious mosquitoes $(1 - u_{1})\rho^{n}\lambda_B$, where $\lambda_B$ is the force of infection in adults.

The infectious adults class, $I_B$, increases when infectious infants mature above 5 years, and adults recrudesce at per capita rate $\xi_B$ and acquire infection at rate $(1 - u_{1})\rho^{n}\lambda_{B}$. The class is reduced at per capita natural death rate $\mu_h$, per capita disease-induced death rate $\delta_B$, per capita removal rate $\eta_B$ due to counterfeit drugs use and recovery at per capita rate $u_{2}\gamma_B+\gamma_{2}$ of adults in it. The recovery rate of adults due to using an effective malarial drug is $\gamma_B$, and $\gamma_2$ is the natural recovery rate of adults.

The adult class of counterfeit drug users, $U_B$, increases when infants who used counterfeit drugs mature above 5 years and when the infectious adults who used counterfeit drugs are removed at a per capita rate $\eta_B$. This class decreases by natural death, recrudescence and reverting to the susceptible class of adults in it, with, respectively, per capita natural death rate $\mu_h$, recrudescence rate $\xi_B$ and susceptibility of adults $\theta$.

The adult class of effective drug users, $T_B$, increases when infants in the infant class of effective drug users mature above 5 years and when the infectious adults recover at a per capita rate $u_{2}\gamma_B+\gamma_{2}$. This class decreases by natural death, loss of immunity and post-treatment prophylaxis of adults with, respectively, per capita natural death rate $\mu_h$ and per capita loss of partial immunity and post-treatment prophylaxis rate $\phi_B$.

Susceptible female Anopheles mosquitoes are recruited at a per capita rate $\alpha_m(1-\frac{N_{m}}{K})$, where $\alpha_{m}$ is the birth rate of mosquitoes, and $K$ is the environmental carrying capacity of pupa ^[41]. They die at a per capita natural death rate $\mu_{m}$ and a per capita indoor residual spraying (IRS) induced death rate $\varphi u_{3}$. The control $u_{3}$ represents the prevention effort with the IRS and includes spraying with insecticide and surveillance. $\varphi$ is the proportion of mortality induced by IRS. We assume that the mosquito population does not go extinct, and hence $r_{m}$ = $\alpha_{m}-\mu_{m}-\varphi u_{3} > 0$. Most of the analysis is redundant if $r_{m} < 0$. The mosquitoes become infected by humans in the infectious and counterfeit drug users class at a rate $(1 - u_{1})\rho^{m}\lambda_m$, where $\lambda_m$ is the force of infection rate in mosquitoes. The daily survival probability of a mosquito is assumed to be $\rho_{2} = e^{-\mu_{m}}$. The probability of survival of a mosquito over the average latent period of length $n_{m}$ and becoming infectious is given by $\rho^{m} = e^{-\mu_{m}n_m}$.

The infective mosquito population increases when mosquitoes acquire infection from humans at a per capita rate $(1 - u_{1})\rho^{m}\lambda_m$. It decreases with a per capita natural death rate $\mu_{m}$ and a per capita IRS-induced death rate $\varphi u_{3}$.

Differences have been observed in malaria transmission rates, degree of infectivity, recovery and immunity due to strains, exposure, locality, age and treatment ^{[3,4,7,42,43]}. In humans, the gametocyte rate and density ^[42], as well as the parasite rate ^[42,44], have been found to decrease with age. Thus, we assume that infants have higher infectiousness to mosquitoes than adults. The degree of infectivity in the human population is denoted by the parameter $\pi_{i}, i = 1, 2, 3$, where $\pi_{1}, \pi_{2}$ and $\pi_{3}$ are associated, respectively, with $I_{B}, U_{A}$ and $U_{B}$ such that $0 < \pi_{3} < \pi_{2} < \pi_{1} < 1$. Further, the transmission probabilities in infants, adults and mosquitoes are denoted by $\beta_A$, $\beta_B$ and $\beta_m$, respectively. The average numbers of mosquito bites per infant and adult per time unit are denoted by $a$ and $b$, respectively, and $c$ is the average number of bites given by a mosquito per time unit.

Usually, the recruitment of children is dependent on the adult population. At the same time, the children mature and move to the adult population. If the population of adults is much larger than that of children, we can assume that it is constant and not affected by the latter. In this case, the total birth rate $\alpha_{h}$ is just the per capita birth rate in the adult population $\alpha_{B}$ multiplied by $N_{B}$, $\alpha_{h} = \alpha_{B}N_{B}$.

The above formulations and assumptions give the following system of ordinary differential equations for the structured malaria model:

with initial conditions

where

The parameter values of the malaria control model are shown in Table 1 in Section 5. The age-structured models (2.1)–(2.10) is an extension of the malaria model developed in ^[45], by including differentiated susceptibility, infectivity and infectiousness of infant and adult populations and control measures to malaria.

2.1.   Basic properties of the structured model

The first step in showing that the malaria models (2.1)–(2.10) makes sense epidemiologically is to prove that the populations remain non-negative, that is, that all solutions of systems (2.1)–(2.10) with positive initial conditions will remain positive for all times $t > 0$. We define a feasible region $\Omega$, such that

Theorem 2.1. The feasible region $\Omega$ is positively invariant and attracting.

Proof. The right-hand side of models (2.1)–(2.10) is continuous with continuous partial derivatives in $\Omega$. Since $S_{A}(0) > 0$, $S_{B}(0) > 0$, the Picard theorem gives the existence of solutions at least on some (maximum) interval $[0, \omega)$. It can be seen that $S_{A}^{'}\geq0$ if $S_{A} = 0$, $S_{B}^{'}\geq0$ if $S_{B} = 0$, $I_{A}^{'}\geq0$ if $I_{A} = 0$, $I_{B}^{'}\geq0$ if $I_{B} = 0$, $U_{A}^{'}\geq0$ if $U_{A} = 0$, $U_{B}^{'}\geq0$ if $U_{B} = 0$, $T_{A}^{'}\geq0$ if $T_{A} = 0$, $T_{B}^{'}\geq0$ if $T_{B} = 0$, $S_{m}^{'}\geq0$ if $S_{m} = 0$, and $I_{m}^{'}\geq0$ if $I_{m} = 0$. Therefore, given the initial condition, the solutions $(S_{A}, S_{B}, I_{A}, I_{B}, U_{A}, U_{B}, T_{A}, T_{B}, S_{m}, I_{m})$ are nonnegative on $[0, \omega)$. We let $0 < \delta_h = min\{\delta_A, \delta_B\}$. Adding the first eight and last two equations of (2.1)–(2.10) gives, respectively,

that is,

and

Solving the inequality (2.13) and Eq (2.14) gives

and

Thus, we see that if $N_{h}(0) > 0$ and $N_{m}(0) > 0$, then neither $N_{h}$ nor $N_{m}$ can become $0$ at any finite time (in particular, on $[0, \omega)$). Then, for $0 < N_h(0)\leq \frac{\alpha_{h}}{\mu_h}$ and $0 < N_m(0)\leq\frac{r_{m}K}{\alpha_{m}}$, we see that on $[0, \omega),$

and thus, by the nonnegativity, the solution ($S_{A}$, $S_{B}$, $I_{A}$, $I_{B}$, $U_{A}$, $U_{B}$, $T_{A}$, $T_{B}$, $S_{m}$, $I_{m}$) is a priori bounded, and hence it is defined for all $t\geq0$ and stays in $\Omega$. Thus, the region is positively invariant. Further, if $N_{h}(0) > \frac{\alpha_{h}}{\mu_h}$ or $N_{m}(0) > \frac{r_{m}K}{\alpha_{m}}$, then they are also isolated from zero, and $\limsup_{t\rightarrow \infty}N_{h}(t)\leq\frac{\alpha_{h}}{\mu_h}$, $\lim_{t\rightarrow \infty} N_{m}(t) = \frac{r_{m}K}{\alpha_{m}}$. Hence, $N_{h}(t)$ either becomes smaller than $\frac{\alpha_{h}}{\mu_h}$ or approaches $\frac{\alpha_{h}}{\mu_h}$, while $N_{m}(t)$ converges to $\frac{r_{m}K}{\alpha_{m}}$. Therefore, the region attracts the solutions in $\Re_{+}^{10}$ with $N_{h}(0) > \frac{\alpha_{h}}{\mu_h}$ and $N_{m}(0) > \frac{r_{m}K}{\alpha_{m}}$. □

Since the region $\Omega$ is positively invariant and attracting, it is sufficient to consider the dynamics of the flow generated by the model in $\Omega$.

3.   Model dynamics

3.1.   Disease-free equilibrium

The disease-free equilibrium (DFE) of the models (2.1)–(2.10) is given by

The local stability of $E_0$ is established by applying the next generation matrix method ^[46] to (2.1)–(2.10). It follows that the basic reproduction number, $R_0$, of the age-structured malaria systems (2.1)–(2.10) is given by $R_0$ = $\rho(FV^{-1})$ where $\rho$ represents the spectral radius. The associated matrices $F$, $V$ and $FV^{-1}$ are given in Appendix A.1. This way, we obtain

where

and

We note that $d_{1}d_{3}-\eta_{A}\xi_{A} > 0$ and $d_{2}d_{4}-\eta_{B}\xi_{B} > 0$.

The result below follows from Theorem 2 of ^[46].

Lemma 1. The DFE, $E_{0}$, of (2.1)–(2.10) is locally asymptotically stable if $R_{0} < 1$ and unstable if $R_{0} > 1$.

3.2.   Existence of endemic equilibrium

The conditions for the existence of an endemic equilibrium $E_{e}$ for the model (2.1)-(2.10) are established in this section. Here, $E_{e} = (S^{**}_A, S^{**}_B, I^{**}_A, I^{**}_B, U^{**}_A, U^{**}_B, T^{**}_A, T^{**}_B, S^{**}_m, I^{**}_m),$ and the forces of infection $\lambda^{*}_{A}, \lambda^{*}_{B}$ and $\lambda^{*}_{A}$ are non-zero. Expressing $S^{**}_{A}$, $S^{**}_{B}$, $I^{**}_{A}$, $I^{**}_{B}$, $U^{**}_{A}$, $U^{**}_{B}$, $T^{**}_{A}$, $T^{**}_{B}$, $S^{**}_{m}$, $I^{**}_{m}$, $\lambda^{*}_{m}$ and $\lambda^{*}_{B}$ in terms of $\lambda^{*}_{A}$, we get

where $\sigma; B_{1}; C_{i}\quad \text{and}\quad F_{i}, i = 1, 2, 3, 4, 5, 6, 7;$ $D_{i}, i = 1, 2, 3, 4, 5;$ $E_{i}, i = 1, 2, 3, 4, 5, 6, 7, 8, 9, 10;$ $G_{i}, i = 1, 2, 3;$ $J_{i}, i = 1, 2, 3;$ $K_{i}, i = 0, 1, 2, 3, 4, 5, 6, 7, 8;$ $L_{i}, i = 1, 2, 3, 4, 5, 6$ are shown in Appendices A.2 and A.3.

Substituting $I^{**}_{m}, I^{**}_{A}$ and $I^{**}_{B}$ into the expression for $\lambda^{*}_{A}$, we obtain the following equation for $\lambda^{*}_{A}$:

where $P_{i}, i = 0, 1, 2, 3, 4, 5, 6,$ are given in Appendix A.3.

Lemma 2. The malaria models (2.1)–(2.10) has at least one endemic equilibrium, $E_{e},$ when $R_{0} > 1$.

Proof. From Eq (3.5) we see that $f(\lambda^{*}_{A})$ is continuous on $[0, \infty)$. We have that $\lim_{\lambda^{*}_{A}\rightarrow \infty}$ $f(\lambda^{*}_{A}) = \infty$ (since $P_{6} > 0$) and $f(0) = P_{0} < 0$, when $R_{0} > 1$. This implies that (3.5) admits at least one positive solution $\lambda^{*}_{A} > 0$. □

Lemma 2 ensures that there exists at least one endemic equilibrium provided $R_{0} > 1$. There is also the possibility of multiple endemic equilibria when $R_{0} > 1$ or $R_{0} < 1,$ as $f(\lambda^{*}_{A})$ is a sextic polynomial. The maximum number of positive roots of Eq (3.5) can be identified using Descartes's rules of signs. It can be determined that when $R_{0} > 1$, $f(\lambda^{*}_{A})$ will have one, three or five positive roots. On the other hand, when $R_{0} < 1$, there will be zero, two, four or six positive roots. The possibility of the existence of multiple endemic equilibria when $R_{0} < 1$ suggests that a backward bifurcation may occur ^{[29,30,36,45]}.

3.3.   Backward bifurcation

We explore the possibility and establish the conditions that ensure the existence of backward bifurcation in the systems (2.1)–(2.10) using the center manifold theory (Theorem 4.1, ^[47]). To do so, we introduce new variables and consider $\beta_{m}$ as a bifurcation parameter with $\beta_{m} = \beta_{m}^*$ corresponding to $R_0 = 1$.

Let $x_{1} = S_{A}, x_{2} = S_{A}, x_{3} = I_{A}, x_{4} = I_{B}, x_{5} = U_{A}, x_{6} = U_{B}, x_{7} = T_{A}, x_{8} = T_{B}, x_{9} = S_{m}, x_{10} = I_{m}$ so that $N_{h} = x_{1}+x_{2}+x_{3}+x_{4}+x_{5}+x_{6}+x_{7}+x_{8}$ and $N_{m} = x_{9}+x_{10}$. Let $X = (x_{1}, x_{2}, x_{3}, x_{4}, x_{5}, x_{6}, x_{7}, x_{8}, x_{9}, x_{10})^{T}$ and $\frac{dX}{dt} = (f_{1}, f_{2}, f_{3}, f_{4}, f_{5}, f_{6}, f_{7}, f_{8}, f_{9}, f_{10})^{T}$ so that (2.1)–(2.10) can be written in the form

with $H_{1} = \frac{c\beta_m (x_{3}+\pi_{1}x_{4}+\pi_{2}x_{5}+\pi_{3}x_{6})}{\sum_{i = 1}^{8} x_{i}}$ and $H_{2} = (1-u_{1})e^{-\mu_{m} n_{m}}H_{1}+\mu_{m}+\varphi u_{3}.$

The bifurcation parameter $\beta^*$ for which $R_0 = 1$ is given by

where $H_{3} = a\beta_{A}\mu_{h}[(d_{2}d_{4}-\eta_{B}\xi_{B})(d_{3}+\pi_{2}\eta_{A})+\tau(\pi_{1}(d_{3}d_{4}+\eta_{A}\xi_{B})+\pi_{3}(\eta_{B}d_{3}+\eta_{A}d_{2}))],$ and $H_{4} = b\beta_{B}\tau(d_{1}d_{3}-\eta_{A}\xi_{A})(\pi_{1}d_{4}+\pi_{3}\eta_{B})$.

The Jacobian matrix, $J_{E_{0}}$, of the transformed systems (2.1)–(2.10) evaluated at the DFE, $E_{0}$, with $\beta_{m} = \beta^{*}_{m}$, is given in Appendix A.2, where, in addition, explicit expressions are provided for the right and left eigenvectors, $w = (w_{1}, w_{2}, w_{3}, w_{4}, w_{5}, w_{6}, w_{7}, w_{8}, w_{9}, w_{10})^{T}$ and $v = (v_{1}, v_{2}, v_{3}, v_{4}, v_{5}, v_{6}, v_{7}, v_{8}, v_{9}, v_{10}),$ respectively, corresponding to the zero eigenvalue of $J_{E_{0}}$. The associated non-zero partial derivatives of $f$ evaluated at the DFE are also listed in Appendix A.2.

Using the approach in ^[47], we have the expressions for $a_1$ and $b_1$ as

and

where

and

with $g_{2} < 0$ if $w_{1}+w_{2} > w_{3}+w_{4}+w_{5}+w_{6}+w_{7}+w_{8}$ or $g_{2}\geq 0$ otherwise.

Clearly, $b_1 > 0$, and hence it follows from (Theorem 4.1, ^[47]) that the direction of the bifurcation of the transformed systems (3.6)–(3.15) at $R_{0} = 1$ is determined by the sign of $a_{1}:$ if $a_{1} > 0$, then it will be backward, and if $a_{1} < 0$, then it will be forward.

To determine the sign of $a_{1}$, we expand the terms in parentheses in (3.16) and obtain

We note that $a_{1} > 0$ if

The above result is summarized below.

Theorem 3.1. The models (2.1)–(2.10) exhibits backward bifurcation at $R_0 = 1$ when (3.19) holds.

Remark 1. We note that in this particular case, the existence of backward bifurcation can be proved in a more elementary way. Consider again

We considered a bifurcation parameter $\beta_m$, such that $R_0$ and all coefficients $P_i$ can be expressed as functions of this parameter. In particular, we defined $\beta_m^*$ by $R_0(\beta_m^*) = 1$. Here, $P_0(\beta_m) = K_hE_7F_7 (1-R_0^2(\beta_m))$ ($K_hE_7F_7$ are independent of $\beta_m$) and $P_1 = F_7(K_hE_6-G_3) + K_hE_7F_6-ua\beta_A (E_7F_2 + E_6F_3).$ We see that the equation

has an isolated single root $\lambda_A^* = \lambda_A(\beta_m^*) = 0$ as soon $P_1(\beta_m^*)\neq 0$. Using the implicit function theorem, for $\beta_m$ close to $\beta_m^*$ (that is, $R_0$ close to $1$), there are differentiable solutions $\lambda_A(\beta_m)$, satisfying

Indeed,

Differentiating (3.21) with respect to $\beta_m$, we obtain

Since, by (3.2)–(3.4), $R_0$ is an increasing function of $\beta_m$, if we assume that $P_1(\beta_m^*) < 0$, then $\frac{d \lambda_A}{d \beta_m}(\beta_m^*) < 0$, that is, $\beta_m\mapsto \lambda_A(\beta_m)$ is decreasing (in some neighborhood of $\beta_m = \beta_m^*$). Since $\lambda_A(\beta_m^*) = 0$, we see that $\lambda_A(\beta_m) > 0$ whenever $\beta_m < \beta_m^*$ is sufficiently close to $\beta_m^*$, that is, whenever $R_0 < 1$ is sufficiently close to $1$. For such $\beta_m < \beta_m^*$ we have $f(0, \beta_m) > 0$, $f(\lambda_A(\beta_m), \beta_m) = 0,$ and, from the above, we see that $\left.\frac{\partial f(\lambda_A^*, \beta_m)}{\partial \lambda_A^*}\right|_{\lambda_A^* = \lambda_A(\beta_m)} < 0.$ Therefore, $f(\lambda_A^*, \beta_m) < 0$ for some $\lambda_A^* > \lambda_A(\beta_m)$. Since $\lim\limits_{\lambda_A^*\to\infty}f(\lambda_A^*, \beta_m) = \infty$, it follows that there is another positive solution to $f(\lambda_A^*, \beta_m) = 0$, and consequently, we have backward bifurcation.

Previous malaria studies ^[29,30,45] have shown the existence of backward bifurcation, with some attributing this to the use of the standard incidence function instead of mass action, and also to partial immunity, repeated infection and a high disease death rate. Our result in Theorem 3.1 indicates that increasing the birth rate of the population above the threshold ($\alpha^{*}_{h}$) can remove the backward bifurcation.

4.   Analysis of the reproduction number

The formula for $R_{0}$ in (3.2)–(3.4) can also be expressed in the form (4.1) given below to provide some insight into the adult, infant and mosquito malaria transmission and control. When a newly infected mosquito is introduced into a completely susceptible human population, during its average infectious period, $\frac{1}{(\mu_m+\varphi u_{3})}$, it will infect humans (infants and adults) who are not using ITNs at the rate $\frac{(1-u_{1})a\beta_{A}S^{*}_{A}}{N^{*}_{h}}+\frac{(1-u_{1})b\beta_{B}S^{*}_{B}}{N^{*}_{h}}$. Humans survive the latent period with probability $e^{-\mu_{h}n_{h}}$ and become infectious. Thus, the total number of infants and adults who become infectious due to this mosquito during its entire infectious period is approximately equal to $R_{I_{mh}} = \frac{e^{-\mu_{h}n_{h}}(1-u_{1})(a\beta_{A}\mu_{h}+b\beta_{B}\tau)}{(\mu_m+\varphi u_{3})(\tau+\mu_{h})}$.

Similarly, for humans, during the average infant's (adult's) infectious period, $\frac{d_{3}}{d_{1}d_{3}-\eta_{A}\xi_{A}}$ $\left(\frac{d_{4}}{d_{2}d_{4}-\eta_{B}\xi_{B}}\right)$, the infant (adult) will infect mosquitoes at the rate $\frac{(1-u_{1})c\beta_{m}S^{*}_{m}}{N^{*}_{h}}$ $\left(\frac{(1-u_{1})\pi_{1}c\beta_{m}S^{*}_{m}}{N^{*}_{h}}\right)$. The mosquitoes survive the latent period with probability $e^{-\mu_{m}n_{m}}$ and become infective, and hence the total number of mosquitoes who become infective due to this infant (adult) during the infectious period is approximately equal to $R_{I_{Am}} = \frac{e^{-\mu_{m}n_{m}}(1-u_{1})c\beta_{m}d_{3}\mu_{h}r_{m}K}{\alpha_{h}\alpha_{m}(d_{1}d_{3}-\eta_{A}\xi_{A})}$ $\left(R_{I_{Bm}} = \frac{e^{-\mu_{m}n_{m}}(1-u_{1})c\beta_{m}\pi_{1}d_{4}\mu_{h}r_{m}K}{\alpha_{h}\alpha_{m}(d_{2}d_{4}-\eta_{B}\xi_{B})}\right)$.

Additionally, when an infant (adult) who used a counterfeit drug is infectious with a degree $\pi_{2} (\pi_{3})$, then he/she has a probability $\frac{\eta_{A}d_{3}}{d_{1}d_{3}-\eta_{A}\xi_{A}}$ $\left(\frac{\eta_{B}d_{4}}{d_{2}d_{4}-\eta_{B}\xi_{B}}\right)$ of surviving the infectious period using the counterfeit drug, and $\frac{1}{d_{3}}$ $\left(\frac{1}{d_{4}}\right)$ is the average time of using the drug. This infant (adult) will infect mosquitoes at the rate $\frac{(1-u_{1})\pi_{2}c\beta_{m}S^{*}_{m}}{N^{*}_{h}}$ $\left(\frac{(1-u_{1})\pi_{3}c\beta_{m}S^{*}_{m}}{N^{*}_{h}}\right)$. The infected mosquitoes survive the latent period with probability $e^{-\mu_{m}n_{m}}$ and then become infective. Hence, the total number of mosquitoes that become infective due to the infant (adult) counterfeit drug user is approximately equal to $R_{U_{Am}} = \frac{e^{-\mu_{m}n_{m}}(1-u_{1})c\beta_{m}\pi_{2}\eta_{A}\mu_{h}r_{m}K}{\alpha_{h}\alpha_{m}(d_{1}d_{3}-\eta_{A}\xi_{A})}$ $\left(R_{U_{Bm}} = \frac{e^{-\mu_{m}n_{m}}(1-u_{1})c\beta_{m}\pi_{3}\eta_{B}\mu_{h}r_{m}K}{\alpha_{h}\alpha_{m}(d_{2}d_{4}-\eta_{B}\xi_{B})}\right)$.

Finally, when an infant who is infectious (or used a counterfeit drug) matures to an adult during the infectious period, the total number of mosquitoes that become infective due to it is approximately $R_{I_{ABm}} = \frac{e^{-\mu_{m}n_{m}}(1-u_{1})c\beta_{m}\pi_{1}\tau\mu_{h} r_{m}K(d_{3}d_{4}+\eta_{A}\xi_{B})}{\alpha_{h}\alpha_{m}(d_{1}d_{3}-\eta_{A}\xi_{A})(d_{2}d_{4}-\eta_{B}\xi_{B})}$ $\left(R_{U_{ABm}} = \frac{e^{-\mu_{m}n_{m}}(1-u_{1})c\beta_{m}\pi_{3}\tau\mu_{h} r_{m}K(\eta_{B}d_{3}+\eta_{A}d_{2})}{\alpha_{h}\alpha_{m}(d_{1}d_{3}-\eta_{A}\xi_{A})(d_{2}d_{4}-\eta_{B}\xi_{B})}\right)$. Therefore,

represents the average number of secondary infections in a completely susceptible human population resulting from introducing one infective human who generates infections in a fully susceptible mosquito population ^[25].

4.1.   Effects of control strategies on $R_{0}$

The effective use of the control measures by infants leads to a reduction in $R_{0}$ since all local reproduction numbers $R_{I_{Am}}$, $R_{U_{Am}}$, $R_{I_{ABm}}$ and $R_{U_{ABm}}$ will decrease. Similarly, adults using control measures will reduce $R_{0}$ because $R_{I_{Bm}}$ and $R_{U_{Bm}}$ will also decrease. We analyze the impact of insecticide-treated bednets (ITNs) on $R_{0}$ when highly effective antimalarial drugs (HEAs) and indoor residual spraying (IRS) are absent. Let us find the difference between the reproduction number with no control measure and with only ITNs as control, and the partial derivative with respect to $u_{1}$ of the reproduction number when ITNs are used. By (3.2)–(3.4), we obtain

and

We find that $R_{0}|_{u_{1} = u_{2} = u_{3} = 0} - R_{0}|_{u_{1}\neq0, u_{2} = u_{3} = 0} > 0$ and $\frac{\partial R_{0}|_{u_{1}\neq0, u_{2} = u_{3} = 0}}{\partial u_{1}} < 0$ for all $u_{1}\in [0, 1)$. Hence, if humans sleep under the ITNs, contact between mosquitoes and humans will become difficult, reducing transmission to and from mosquitoes and hence $R_0$.

On the other hand, treating infectious humans (infants and adults) will decrease the average length of the infection period and transmission period, causing decreases in $R_{0}$ and the disease incidence. The duration of the infection can be reduced if a highly effective antimalarial drug is used ^[7,48] and also if the recovery rate due to using an effective antimalarial drug is increased ^[19,45]. Even though using counterfeit drugs for treatment reduces disease symptoms and may prevent death, it creates and increases a reservoir of infection. This is because humans who use such drugs can still transmit parasites (gametocytes) to feeding mosquitoes and also can recrudesce. Hence, if more humans opt for counterfeit drugs, then it will be difficult to control malaria transmission. On the other hand, if both infants and adults use only highly effective antimalarial drugs for the treatment of malaria, then $R_{U_{Am}}$, $R_{U_{Bm}}$ and $R_{U_{ABm}}$ will equal zero, and $R_{0}$ will be vastly reduced.

Lastly, the effects of the IRS on $R_{0}$ are considered. We let $R_{I_{mh}}$ denote the local reproduction number from mosquitoes to humans. Taking the difference between this reproduction number without control and with the IRS control only, and the partial derivative with respect to $u_{3}$ of this reproduction number with the IRS control, yields

and

Since $R_{I_{mh}}|_{u_{1} = u_{3} = 0} - R_{I_{mh}}|_{u_{1} = 0, u_{3}\neq0} > 0$ and $\frac{\partial R_{I_{mh}}}{\partial u_{3}} < 0$ for all $\varphi > 0$ and $u_{3} > 0$, using IRS also reduces $R_{0}$.

5.   Application of optimal control

This section uses control theory to derive optimal prevention and treatment methods to curtail malaria infection in an endemic area while minimizing the implementation cost. We apply Pontryagin's maximum principle to determine the necessary conditions for the optimal control of the age-structured model. The controls used are based on treatment and preventive tools adopted by most endemic countries in Sub-Saharan Africa ^[4,40].

5.1.   Analysis of optimal control

Together with the malaria models (2.1)–(2.10), we consider an optimal control problem with the objective (cost) functional given by

where $t_{f}$ is the final time, and the coefficients $z_{1}$, $z_{2}$, $z_{3},$ $z_{4}$ and $z_{5}$ represent, respectively, the positive weight constants of the infectious infants and adults, infants and adults who used counterfeit drugs, and the total mosquito population. On the other hand, the coefficients $Y_{1}$, $Y_{2}$ and $Y_{3}$ are positive constant weights for prevention efforts with insecticide-treated bednets (ITNs), treatment efforts with highly effective antimalarial drugs (HEAs) and prevention effort with indoor residual spraying (IRS), respectively, which regularize the optimal control. The terms $z_{1}I_{A}, z_{2}I_{B}, z_{3}U_{A}, z_{4}U_{B}$ and $z_{5}N_{m}$ are the costs of infection in infants, adults and in the total mosquito population. It is assumed that the cost of prevention and treatment, given in the quadratic form in (5.1), that is, $Y_{1}{u_{1}}^{2}, Y_{3}{u_{3}}^{2}$ and $Y_{2}{u_{2}}^{2},$ are the costs of prevention with ITNs, IRS and treatment with HEAs, respectively. The cost of prevention is associated with the purchase and use of ITNs and IRS. Similarly, the cost of treatment is the cost associated with the diagnosis or medical examination, HEAs and hospitalization.

Note that for bounded Lebesgue measurable controls and non-negative initial conditions, non-negative bounded solutions to the state system exist ^[49]. Our goal is to minimize the number of infectious humans, counterfeit drug users, the total number of mosquitoes and the cost of implementing the controls $u_{1}(t), u_{2}(t)$ and $u_{3}(t)$. Hence, we seek to find optimal controls ${u_{1}}^{*}, {u_{2}}^{*}$ and ${u_{3}}^{*}$ such that

where the control set is $\Gamma = \{(u_{1}, u_{2}, u_{3})$ such that $u_{1}, u_{2}, u_{3}$ are Lebesgue measurable with $0\leq u_{1}\leq1$, $0\leq u_{2}\leq 1$, $0\leq u_{3}\leq1$ for $t\in [0, t_{f}]\}$.

The necessary conditions that an optimal solution must satisfy come from the Pontryagin maximum principle ^[34]. This principle converts (2.1)–(2.10) and (5.1) into a problem of minimizing pointwise the Hamiltonian $H$, given below, with respect to $u_{1}, u_{2}$ and $u_{3}$.

where $S_{1} = e^{-\mu_h n_h}(1-u_{1})\frac{a\beta_{A}}{N_{h}}I_{m}, \quad S_{3} = e^{-\mu_h n_h}(1-u_{1})\frac{b\beta_{B}}{N_{h}}I_{m}, \quad S_{5} = e^{-\mu_{m}n_{m}}(1-u_{1})\frac{c\beta_{m}}{N_{h}}(I_{A}+\pi_{1} I_{B}+\pi_{2} U_{A}+\pi_{3} U_{B}),$ and $\lambda_{S_{A}}$, $\lambda_{S_{B}}$, $\lambda_{I_{A}}$, $\lambda_{I_{B}}$, $\lambda_{U_{A}}$, $\lambda_{U_{B}}$, $\lambda_{T_{A}}$, $\lambda_{T_{B}}$, $\lambda_{S_{m}}$ and $\lambda_{I_{m}}$ are the adjoint variables.

Theorem 5.1. Given the solutions $S_{A}$, $S_{B}$, $I_{A}$, $I_{B}$, $U_{A}$, $U_{B}$, $T_{A}$, $T_{B}$, $S_{m}$ and $I_{m}$ of the state system (2.1)–(2.10) and optimal controls $u^{*}_{1}$, $u^{*}_{2}$, $u^{*}_{3}$ that minimize $J(u_{1}, u_{2}, u_{3})$ over $\Gamma$, there exist adjoint variables $\lambda_{S_{A}}$, $\lambda_{S_{B}}$, $\lambda_{I_{A}}$, $\lambda_{I_{B}}$, $\lambda_{U_{A}}$, $\lambda_{U_{B}}$, $\lambda_{T_{A}}$, $\lambda_{T_{B}}$, $\lambda_{S_{m}}$ and $\lambda_{I_{m}}$ satisfying

where

with transversality conditions

for $i$ = $S_{A}, S_{B}, I_{A}, I_{B}, U_{A}, U_{B}, T_{A}, T_{B}, S_{m}$ and $I_{m}$. Furthermore, the controls $u^{*}_{1}, u^{*}_{2}, u^{*}_{3}$ satisfy the optimality condition

Proof. Corollary 4.1 in ^[50] gives the existence of an optimal control on account of the convexity of the integrand of $J$ with respect to $u_{1}, u_{2}, u_{3}$, a priori boundedness of the state solutions and the Lipschitz property of the state system with respect to the state variables. The differential equations governing the adjoint variables are obtained by differentiating the Hamiltonian, evaluated at the optimal control. Then, the adjoint system can be written as

for $i$ = $S_{A}, S_{B}, I_{A}, I_{B}, U_{A}, U_{B}, T_{A}, T_{B}, S_{m}$ and $I_{m}$. To obtain the optimality conditions (5.5), we also differentiate the Hamiltonian with respect to $(u_{1}, u_{2}, u_{3})$ $\in \Gamma$ and set it equal to zero. Thus,

Solving for the optimal control, we obtain

By standard control arguments involving the controls' bounds, we conclude that

where

□

The optimality system consists of the state systems (2.1)–(2.10) with initial conditions, the adjoint system (5.4) with transversality conditions (5.4) and the optimality condition (5.5).

5.2.   Numerical simulation

As mentioned in the introduction, our next objective is to investigate numerically the optimal control strategies that can eliminate malaria infection in the age-structured population using two main scenarios, namely, the population without counterfeit drug use, i.e., $\eta_{A} = \eta_{B} = \xi_{A} = \xi_{B} = \theta = 0$, and the population with counterfeit drugs. For the first scenario, we use the following control strategies: (ⅰ.) insecticide-treated bednets only, $(u_{1}, 0, 0)$, (ⅱ.) highly effective antimalarial drugs only, $(0, u_{2}, 0)$, (ⅲ.) indoor residual spraying only, $(0, 0, u_{3})$, (ⅳ.) insecticide-treated bednets and highly effective antimalarial drugs only, $(u_{1}, u_{2}, 0)$, (ⅴ.) insecticide-treated bednets and indoor residual spraying only, $(u_{1}, 0, u_{3})$, (ⅵ.) highly effective antimalarial drugs and indoor residual spraying, $(0, u_{2}, u_{3})$ and (ⅶ.) insecticide-treated bednets, highly effective antimalarial drugs and indoor residual spraying, $(u_{1}, u_{2}, u_{3})$. For the second scenario, we shall use the best control strategy results from the first scenario and incorporate the effects of counterfeit drug use considering the following sub-categories: (ⅰ) high removal rate and high recrudescence rate using $\eta_{A} = \frac{1}{7}$, $\eta_{B} = \frac{1}{7}$, $\xi_{A} = \frac{1}{7}$, $\xi_{B} = \frac{1}{7}$, (ⅱ) high removal rate and low recrudescence rate using $\eta_{A} = \frac{1}{7}$, $\eta_{B} = \frac{1}{7}$, $\xi_{A} = \frac{1}{30}$, $\xi_{B} = \frac{1}{30}$, (ⅲ) low removal rate and high recrudescence rate using $\eta_{A} = \frac{1}{30}$, $\eta_{B} = \frac{1}{30}$, $\xi_{A} = \frac{1}{7}$, $\xi_{B} = \frac{1}{7}$ and (ⅳ) low removal rate and low recrudescence rate using $\eta_{A} = \frac{1}{30}$, $\eta_{B} = \frac{1}{30}$, $\xi_{A} = \frac{1}{30}$, $\xi_{B} = \frac{1}{30}$. A fourth-order Runge-Kutta iterative scheme is used to solve the optimality system. The state system is solved forward in time with initial conditions $(S_{A}, S_{B}, I_{A}, I_{B}, U_{A}, U_{B}, T_{A}, T_{B}, S_{m}, I_{m})$ = $(508, 6505,100,300, 0, 0, 0, 0, 14583,300),$ and the adjoint system is solved backwards in time with transversality conditions (5.4). The initial total populations are estimated using (2.15) and the parameter values in Table 1. At first, the state system and the adjoint system are solved with an initial guess for the control $(u_{1}(t), u_{2}(t), u_{3}(t)) = (0, 0, 0)$. The control functions are updated using the optimality conditions given by (5.5), and the process is repeated. This iterative process terminates when the values converge sufficiently: the differences between the current state, adjoint and control values and the previous state, adjoint and control values are negligibly small ^[51]. The current values are then taken as the solution. The numerical values of the parameters used for solving the optimality system to obtain the optimal solution are given in Table 1. Most of the parameter values are taken from the literature on Ghana. Other parameter values not directly found in the literature were estimated using the assumptions made during the model formulation and following literature indications. The following weight constants were used: $Y_{1} = 20, Y_{2} = 40, Y_{3} = 25$ and $z_1 = 30, z_2 = 25, z_3 = 20, z_4 = 15, z_5 = 20$. The control is applied over $5$ years ($1825$ days).

5.3.   Simulation results

5.3.1.   The control strategies without counterfeit antimalarial drug use

We investigate the effect of each of the seven control strategies mentioned above on the dynamics of the population when no counterfeit antimalarial drugs are used (i.e., $\eta_{A} = \eta_{B} = \xi_{A} = \xi_{B} = \theta = 0$).

1) Insecticide-treated bednets only:

The strategy considers insecticide-treated bednets ($u_{1}$) as the only control. The highly effective antimalarial drugs ($u_{2}$) and indoor residual spraying ($u_{3}$) are set to zero to optimize the objective function (5.1). The benefits of using ITNs are seen in the reduction of the number of infectious individuals in both mosquito and human populations, as compared to when there is no control, as shown in Figure 3(a)–(c). Also, the number of susceptible mosquitoes is higher than when no control is used (Figure 3(d)). The control profile in Figure 3(e) indicates that, in this case, the coverage of insecticide-treated bednets ($u_{1}$) should be maintained at $75\%$ for the entire duration of the intervention.

2) Highly effective antimalarial drugs only:

Here, only highly effective antimalarial drugs $(u_{2})$ are used to optimize the objective function (5.1), while the insecticide-treated bednets $(u_{1})$ and indoor residual spraying $(u_{3})$ are absent. From the results in Figure 4(a)–(c), we observe a significant decrease in the numbers of the infectious human and mosquito populations as compared to when no control is used. The results also show an increase in the susceptible mosquito population. (Figure 4(d)). From the control profile shown in Figure 4(e), we observe that strategy $(u_{2})$ should be maintained at coverage of $75\%$ for 148 days, then gradually reduced to 25% by day 350 and maintained at this level for the remaining intervention period.

3) Indoor residual spraying only:

We use indoor residual spraying $(u_{3})$ to optimize the objective function (5.1), while insecticide-treated bednets $(u_{1})$ and highly effective antimalarial drugs $(u_{2})$ are absent. The results reveal a reduction in the number of infectious humans $I_{A}$, $I_{B}$ and mosquitoes $I_{A}$ when the control is applied as compared to the case without control, as observed in Figure 5(a)–(c). The graph in Figure 5(d) also shows an increase in the susceptible mosquito population. We observe from the control profile shown in Figure 5(e) that the IRS $(u_{3})$ only strategy should be maintained at coverage of $75\%$ for the total duration of the intervention.

4) Insecticide-treated bednets and highly effective antimalarial drugs:

We use the strategy which has both insecticide-treated nets $(u_{1})$ and highly effective antimalarial drugs $(u_{2})$ as controls while setting the indoor residual spraying to zero $(u_3 = 0)$ and optimize the objective function (5.1). Figure 6(a)–(d) shows a significant decrease in the number of infectious humans and mosquitoes and a considerable increase in the number of susceptible mosquitoes due to the application of this control. Its profile in Figure 6(e) also reveals that to be optimal, the highly effective antimalarial drugs coverage should be at 75% for 74 days, and then gradually be reduced to 25% on day 230 for the rest of the duration, while the insecticide-treated net's coverage should begin at 44%, then increase to 45% from day 91 to day 1735 and finally be reduced gradually to 25% by the end of the intervention.

5) Insecticide-treated bednets and indoor residual spraying:

When the insecticide-treated bednets $(u_1)$ and indoor residual spraying $(u_3)$ are used to optimize the objective function (5.1) without highly effective antimalarial drugs $(u_2 = 0)$, we observe a reduced number of infectious humans $I_{A}$, $I_{B}$ and mosquitoes $I_{m}$ and an increased number of susceptible mosquitoes (Figure 7(a)–(d)) as compared to the populations with no control. The control profile in Figure 7(e) shows that for this strategy to be optimal, the insecticide-treated bednets and the indoor residual spraying must both be maintained at 75% coverage for the total duration of the intervention.

6) Highly effective antimalarial drugs and indoor residual spraying:

Here, we use highly effective antimalarial drugs $(u_2)$ and indoor residual spraying $(u_3)$ while setting the control on insecticide-treated bednets to zero $(u_1 = 0)$ to optimize the objective function (5.1). As shown in Figure 8(a)–(d), with this control, the infectious populations $I_{A}$, $I_{B}$ and $I_{m}$ are significantly smaller, and the susceptible mosquitoes population $S_{m}$ is larger, as compared to the populations without control. The control profile in Figure 8(e) shows that the coverage of indoor residual spraying should be maintained at 75% for the whole duration of the intervention. At the same time, administration of highly effective antimalarial drugs should begin with 75% coverage for 61 days before being reduced to 25% by day 132 and maintained at this level for the remaining treatment time.

7) Insecticide-treated bednets, highly effective antimalarial drugs and indoor residual spraying:

We explore the combination of all three controls by using insecticide-treated bednets $(u_1)$, highly effective antimalarial drugs $(u_2)$ and indoor residual spraying $(u_3)$ to optimize the objective function (5.1). The benefits of this strategy are a significant decrease in the populations of infectious humans and mosquitoes and a significant increase in the number of susceptible mosquitoes in comparison to when no controls are used (Figure 9(a)–(d)). The control profile indicates that for the strategy to be optimal, the indoor residual spraying coverage at 75% should be maintained for the entire intervention period, the insecticide-treated bednets should be applied with coverage of 32% between day 1 to day 1769 and then reduced gradually to 25%, and, lastly, highly effective antimalarial drugs should be administered with coverage of 75% for the first 47 days and then reduced to 25% by day 105 and maintained at this level to the end of the period (Figure 9(e)).

5.3.2.   The effect of counterfeit drugs on the best control strategy

In this section, we examine the impact of four possible cases resulting from counterfeit drug use on the performance of the best control strategies discussed in Section 5.3.1, i.e., the insecticide-treated bednets, highly effective antimalarial drugs and indoor residual spraying.

1) Counterfeit drug with high removal rate and high recrudescence rate:

As the control we use insecticide-treated bednets $(u_1)$, highly effective antimalarial drugs $(u_2)$ and indoor residual spraying $(u_3)$ and set the per capita removal rates $\eta_{A} = \frac{1}{7}$ and $\eta_{B} = \frac{1}{7}$ and the per capita recrudescence rate $\xi_{A} = \frac{1}{7}$ and $\xi_{B} = \frac{1}{7}$ for both infants and adults. Figures 10(a)–(g) show significant increases in the populations of susceptible humans $S_{A}$ and $S_{B}$ and significant decreases in the numbers of infectious individuals $I_{A}$, $I_{B}$, $I_{m}$ and counterfeit antimalarial drug users $U_{A}$, $U_{B}$ when the control is used. The control profile in Figure 10(h) shows that maintaining indoor residual spraying coverage at 75% for the whole intervention period, highly effective antimalarial drugs coverage at 75% for the first 94 days and gradually reduced to 25% by day 191, and, lastly, the insecticide-treated bednets coverage at 36%, reduced to 35% on day 1744 and then gradually to 25%, is optimal.

2) Counterfeit drug with high removal and low recrudescence rates:

Similarly, we use the insecticide-treated bednets $(u_1)$, highly effective antimalarial drugs $(u_2)$ and indoor residual spraying $(u_3)$ control to optimize the objective function (5.1) and take the per capita removal rates $\eta_{A} = \frac{1}{7}$ and $\eta_{B} = \frac{1}{7}$ and the per capita recrudescence rates $\xi_{A} = \frac{1}{30}$ and $\xi_{B} = \frac{1}{30}$ for both infants and adults. We observe in Figure 11(a)–(g) that the control strategy causes a significant difference in the populations. Namely, the number of susceptibles increases while that of the infectious and counterfeit drug users decreases, as compared to when no control is used. The control profile in Figure 11(h) shows that it is optimal to use the indoor residual spraying coverage at 75% for the entire intervention period, the insecticide-treated bednets coverage at 60%, reduced to 55% by day 27, increased to 57% until day 1532 and again gradually reduced to 25%, and, lastly, starting with highly effective antimalarial drugs coverage at 75% for 258 days and then gradually reducing it to 25% on day 575 and maintaining this level to the end of the period.

3) Counterfeit drug with low removal and high recrudescence rates:

We use insecticide-treated bednets $(u_1)$, highly effective antimalarial drugs $(u_2)$ and indoor residual spraying $(u_3)$ to optimize the objective function (5.1) and take the per capita removal rates $\eta_{A} = \frac{1}{30}$, $\eta_{B} = \frac{1}{30}$ and the per capita recrudescence rates $\xi_{A} = \frac{1}{7}$, $\xi_{B} = \frac{1}{7}$. The control produces an increase in the susceptible population compared to when no control is applied (Figure 12(a), (b)). Also, the infectious and counterfeit antimalarial drug users populations are observed to decrease (Figure 12(c)–(g)). Figure 12(h) illustrates that indoor residual spraying with coverage of 75% for the entire intervention period, insecticide-treated bednets coverage at 33% at the beginning, then reduced and maintained at 32% for 1794 days and reduced again to 25%, and, lastly, highly effective antimalarial drugs coverage at 75% for the first 61 days and reduced to 25% by day 127 is optimal for disease control.

4) Counterfeit drug with low removal and low recrudescence rates:

Finally, we use insecticide-treated bednets $(u_1)$, highly effective antimalarial drugs $(u_2)$ and indoor residual spraying $(u_3)$ to optimize the objective function (5.1) and take the per capita removal rates $\eta_{A} = \frac{1}{30}$ and $\eta_{B} = \frac{1}{30}$ and the per capita recrudescence rates $\xi_{A} = \frac{1}{30}$ and $\xi_{B} = \frac{1}{30}$. The control produces significant decreases in the infectious classes $I_{A}$, $I_{B}$, $I_{m}$ and also in the number of counterfeit antimalarial users in the classes $U_{A}$, $U_{B}$ (Figure 13(c)–(g)) in comparison with the populations without control. The control profile indicates that maintaining the indoor residual spraying coverage at 75% for the entire intervention period, beginning with the insecticide-treated bednets coverage of 37%, then reducing it and maintaining it at 36% from day 4 to day 1735 and reducing it to 25% by the end of the intervention period, and lastly, administering the highly effective antimalarial drugs with coverage of 75% for 127 days before reducing it to 25% by day 308 is optimal for disease control (Figure 13(h)).

5.3.3.   Comparison of control strategies

The simulation results from all seven control strategies without counterfeit drugs (in Figures 3–9) reveal that all the strategies reduce the populations of infectious humans and infectious mosquitoes and increase the susceptible populations.

When counterfeit drug use is absent in the population, we observe that any strategy with highly effective antimalarial drugs achieves a 100% reduction in the infectious population. Still, any strategy without highly effective antimalarial drugs reduces the numbers of infectious adults and infectious mosquitoes significantly more than the number of infectious infants. Also, a strategy with only one prevention control measure, either IRS or ITNs, may be insufficient to reduce the disease burden in humans and to control malaria unless they are combined with the other controls.

Comparing the four strategies with highly effective antimalarial drugs, we observe similar benefits of the 100% reduction in the infectious populations (Figures 4, 6, 8 and 9). However, the four strategies differ in the increase in the susceptible mosquito population and the control profiles. The two strategies, (Ⅰ) highly effective antimalarial drugs and indoor residual spraying and (Ⅱ) insecticide-treated bednets, highly effective antimalarial drugs and indoor residual spraying, give the smallest increases in the susceptible mosquitoes numbers, which is a part of the objective of our optimal control problem. From the control profiles, the duration of the coverage of highly effective antimalarial drugs in (Ⅰ) must be longer than in (Ⅱ) (Figures 8(e) and 9(e)).

We observe that using the different sub-categories of counterfeit drugs without control significantly reduces the population of susceptible infants (in Figures 10–13). However, when the combination of insecticide-treated bednets, highly effective antimalarial drugs and indoor residual spraying is used as control, we observe a considerable increase in that population. We also observe that a combination of insecticide-treated bednets, highly effective antimalarial drugs and indoor residual spraying achieves a 100% reduction in the infectious population even when counterfeit drug use persists. From the control profiles in Figures 9(e), 10–13(h), we observe that the required intervention periods with high coverage are more extended when there is counterfeit drug use. This is particularly visible when applying highly effective antimalarial drugs or insecticide-treated bednets as controls. On the other hand, the rates of removal or recrudescence of counterfeit drugs seem to have no impact on the increases in the susceptible population.

Ranking the effect of the four subcategories of counterfeit drugs on the performances of the best control strategy, we see that when the counterfeit drugs used have a high removal rate and low recrudescence rate, then the control strategy has to be applied for the most extended period to achieve the elimination; this is followed by counterfeit drugs with low removal and recrudescence rates, then counterfeit drugs with high removal and recrudescence rates and, finally, counterfeit drugs with low removal rate and high recrudescence rate. Thus, counterfeit drugs with low removal and high recrudescence rates are the easiest to deal with using the considered control strategies (Figure 12(h)).

6.   Discussion

A deterministic model for P. falciparum malaria infection in a structured human population that uses counterfeit drugs was developed. The model incorporates the infant and adult populations, users of counterfeit drugs and three malaria control measures adopted by most endemic countries in Sub-Saharan Africa, namely, highly effective antimalarial drugs (HEAs), insecticide-treated bednets (ITNs) and indoor residual spraying (IRS). The model was developed to comprehensively understand disease transmission dynamics in children under five years, adults and mosquitoes and identify the best control strategy for disease elimination. We performed a standard analysis of the reproduction number $R_{0}$ and the equilibria. Our model revealed the possible existence of a backward bifurcation.

The $R_{0}$ and bifurcation analysis indicate that to control malaria transmission from the infectious population to others, it is necessary to reduce the value of $R_{0}$ below 1. However, this approach may not be sufficient for elimination when backward bifurcation is present. We found that a backward bifurcation occurs in our model when the birth rate of humans is lower than a determined threshold. Further, the analysis of $R_{0}$ showed that the control measures negatively impact malaria transmission. Thus, increasing coverage of the three control measures in the endemic setting through easy access, affordable pricing and individual adherence will reduce $R_{0}$ and hence the malaria burden.

Also, the optimal control analysis revealed that using any strategy without counterfeit drugs reduces the size of the infectious populations while increasing the susceptible populations, especially infants. For infectious diseases, when the infectious population is reduced, there is a reduction in the disease and the value of $R_{0}$. We also found that the ITNs-only or IRS-only strategy was less effective than the highly effective antimalarial drugs-only strategy for reducing disease in humans and mosquitoes. Thus, the ITNs-only and IRS-only strategies are only effective when used together or with highly effective antimalarial drugs. This result is supported by the evidence from previous studies ^[6,7,26], which showed that in an endemic region, where P. falciparum malaria infection is common, intervention strategies that involve highly effective antimalarial drugs are very effective in reducing the disease burden. In particular, our results suggest that the combination of highly effective antimalarial drugs and indoor residual spraying, and the combination of insecticide-treated bednets, highly effective antimalarial drugs and indoor residual spraying, have the best impact on malaria control and elimination. These strategies eliminate the infection from the population of infants, adults and mosquitoes, giving the highest increase in the susceptible human population and the smallest increase in the susceptible mosquito population. The study in ^[24] also confirmed that the combination of the three controls was the most effective, and the study in ^[23] even proved that the combination of highly effective antimalarial drugs and indoor residual spraying was the most cost-effective. Considering which of these two strategies should be adopted must be guided by the availability of resources, possible levels of coverage and their cost.

Finally, we observed that using the combination of insecticide-treated bednets, highly effective antimalarial drugs and indoor residual spraying as a control strategy, where counterfeit drugs persist, requires an increase in the duration of the coverage, especially for highly effective antimalarial drugs and insecticide-treated bednets, to achieve elimination. This increase is the highest when counterfeit drugs have high removal and low recrudescence rates. Since the counterfeit drug users remain asymptomatically infectious and may not seek further treatment, if many people in a community are using counterfeit drugs as an alternative to highly effective antimalarial drugs, then more resources, effort, time and funds will be required to control and possibly eliminate the disease ^[2]. Efforts to determine the level of use of counterfeit drugs are thus necessary to decide how an intervention can be properly implemented, especially when resources are scarce ^[58].

Our results exposed the negative impact of counterfeit drugs, even when the best control strategies were used. Hence, we conjecture that the effects will worsen when other strategies are used. The counterfeit drug effects could also be further compounded when factors such as the effects of environment, weather and movement of both humans and mosquitoes ^{[25,33,54,59,60]} on the disease dynamics and control are considered.

7.   Conclusions

The control and possible elimination of malaria can be achieved if the existing control strategies of insecticide-treated bednets, vector control methods and prompt treatment of infected humans with highly effective antimalarial drugs are accessible, affordable and duly implemented. This will also relieve infants under five years of the burden of the disease. However, the widespread use of counterfeit drugs may jeopardize the timely achievement of the malaria elimination milestones, especially in endemic areas.

Acknowledgments

The authors sincerely thank the three anonymous reviewers who provided insightful suggestions and to Prof. Lamb for helping in the final proofreading. We gratefully acknowledge all the help provided during the editorial process, especially the financial support.

Conflict of interest

The authors declare there is no conflict of interest.

Appendix

A.1

The associated non-negative matrix $F$ and the non-singular M-matrix $V$ of (2.1)–(2.10) are, respectively, given by

where

Further calculation gives $V^{-1}$, $FV^{-1}$ and $\rho(FV^{-1})$ as follows:

and

with

$Z_{1} = \frac{q_{1}}{\mu_m+\varphi u_{3}},$ $Z_2 = \frac{q_{2}}{\mu_m+\varphi u_{3}},$ $Z_3 = \frac{q_3d_3+q_{3}\pi_{2}\eta_A}{d_1d_3-\xi_A\eta_A}+\frac{q_3\pi_{1}\tau(d_3d_4+\xi_B\eta_A)+q_3\pi_{3}\tau(d_2\eta_A+d_3\eta_B)}{(d_1d_3-\xi_A\eta_A)(d_2d_4-\xi_B\eta_B)},$

$Z_4 = \frac{q_3\pi_1d_{4}+q_{3}\pi_{3}\eta_B}{d_2d_4-\xi_B\eta_B},$ $Z_5 = \frac{q_3\xi_{A}+q_{3}\pi_{2}}{d_1d_3-\xi_A\eta_A}+\frac{q_3\pi_{1}\tau(d_1\xi_B+d_4\xi_A)+q_3\pi_{3}\tau(d_1d_2+\xi_A\eta_B))}{(d_1d_3-\xi_A\eta_A)(d_2d_4-\xi_B\eta_B)}$

and $Z_6 = \frac{q_3\pi_1d_{2}+q_{3}\pi_{3}\eta_B}{d_2d_4-\xi_B\eta_B}.$

The spectral radius of $FV^{-1} = \sqrt{Z_1Z_3+Z_2Z_4}.$

A.2

The Jacobian matrix of the transformed systems (3.6)–(3.15) evaluated at the DFE, $E_{0},$ with $\beta_{m} = \beta^{*}_{m}$, is given by

where $r_{m} = \alpha_{m}-\mu_{m}-\varphi u_{3}$, $J_1 = \tau+\mu_h$, $J_2 = \phi_{A}+\tau+\mu_h$, $J_3 = \phi_{B}+\mu_h$.

The right eigenvector corresponding to the zero eigenvalue of $J_{E_0}$ is given by

where

and

with $B_{1} = d_{1}d_{3}-\xi_{A}\eta_{A}$, $B_{2} = d_{2}d_{4}-\xi_{B}\eta_{B}$, $B_{4} = d_{3}d_{4}+\xi_{B}\eta_{A}$, $B_{5} = d_{3}\eta_{B}+d_{2}\eta_{A}$, $B_{6} = (u_{2}\gamma_{A}+\gamma_{1})d_{3}B_{2}(J_{3}\phi_{A}+J_{1}\phi_{B})+J_{1}J_{2}(J_{3}\theta B_{5}+\phi_{B}(u_{2}\gamma_{B}+\gamma_{2})B_{4})$, $B_{7} = J_{3}\eta_{B}\theta+\phi_{B}(u_{2}\gamma_{B}+\gamma_{2})d_{4}$ and $B_{8} = q_{1}q_{3}[B_{2}(d_{3}+\pi_{2}\eta_{A})+\tau(\pi_{1}B_{4}+\pi_{3}B_{5})]+q_{2}q_{3}B_{1}(\pi_{1}d_{4}+\pi_{3}\eta_{B})$.

We note that $B_{6}-J_{2}J_{3}B_{1}B_{2} < 0$, $B_{7}-J_{3}B_{2} < 0,$ and $B_{8} > (\alpha_{m}-2r_{m})B_{1}B_{2}$.

Similarly, the left eigenvector corresponding to the zero eigenvalue of $J_{E_{0}}$ is given by

where

For the transformed systems (3.6)–(3.15), the associated non-zero partial derivatives of $f$ (evaluated at the DFE), which are needed in the computation of $a_{1}$ and $b_{1}$, are given by

where $m = c(1-u_{1})e^{-\mu_{m}n_{m}}$. For $i = 1, 2, 7, 8,$

A.3

The parameters in the endemic equilibrium, $E_{e},$ and (3.5) are defined as follows:

The coefficients of (3.5) are defined as