A mathematical model for malaria disease dynamics with vaccination and infected immigrants

1.   Introduction

Malaria is a deadly mosquito-borne disease caused by Plasmodium parasites transmitted by Anopheles female mosquitoes between humans. Once an infected mosquito bites a human, the parasites multiply in the host's liver, destroying red blood cells during infection ^[1]. Parasites that infect humans are Plasmodium falciparum (P. falciparum), Plasmodium vivax (P. vivax), Plasmodium ovale (P. ovale), Plasmodium knowlesi (P. knowlesi) and Plasmodium malariae (P. malariae) ^[2]. P. falciparum and P. vivax are unicellular protozoan parasites of humans, and are the most important Plasmodium species causing malaria in humans ^[2]. Malaria is most common in tropical regions with high temperatures and humidity such as Africa, Central South America, the Dominican Republic, areas in the Caribbean, Eastern Europe, South Asia, and islands in the central and South Pacific Ocean ^[3]. However, other regions such as Southeast Asia, the Eastern Mediterranean, the Western Pacific and America also reported significant numbers of cases and deaths ^[4]. Information gathered by the World Health Organization (WHO) ^[4] indicated that there was an estimated 241 million cases of malaria in 2020 while deaths stood at 627,000, with most cases and deaths being recorded in Sub-Saharan Africa. The Sub-Saharan African region accounted for $95\%$ of cases and $96\%$ of malaria deaths with children under the age of five years, being the most affected, accounting for about $80\%$ of malaria deaths ^[4], which has motivated us to carryout this study. People who live in poverty and do not have access to healthcare are more likely to have complications from the disease ^[3]. It is less common for pregnant women with malaria to transfer the parasite to their children before or during birth ^[3]. After a person has been infected, malaria symptoms usually appear within $10$ to $30$ days ^[3]. However, some people do not feel sick for up to a year after infection. Some studies suggest that parasites can live in the host body for several years without causing any symptoms ^[3]. Malaria often has mild flu-like symptoms such as headache, cough, and nausea, which may be difficult to recognize as malaria ^[4]. If not attended to in 24 hours, susceptible specimens like P. falciparum can progress to severe illness with symptoms such as high fever, fatigue, chest pains, breathing problems, diarrhea and vomiting ^[5], and when not treated, P. falciparum can cause anemia and jaundice ^[3], which can lead to death ^[4]. Children with severe malaria frequently develop severe anemia, respiratory distress to metabolic acidosis or cerebral malaria ^[4]. People residing in malaria-endemic areas may develop partial immunity, allowing asymptomatic infections to develop ^[4]. Some population groups such as infants, children under the age of five, patients living with HIV/AIDS, pregnant women and nonimmune migrants are considered to have a higher risk of contracting malaria and having a severe disease ^[6]. The first ordinary differential equation (ODE) model of malaria transmission was proposed and analyzed by Ross ^[7] and refinements of this model were produced by MacDonald ^[8]. The result became known as the Ross-MacDonald model, which has become a basis for modeling of mosquito-borne diseases. Several models of malaria have been studied and improvements on the formulations and analysis of the models keep on evolving. Models of malaria disease with vaccination have been studied by several authors; for instance, Mohammed ^[9], Tchoumi ^[10] and models with infected immigrants presented by, for instance; Witbooi ^[11] and Traoré ^[12]. De la Sen et al. (2015) ^[14] studied a susceptible, exposed, infectious and recovered (SEIR) epidemic model with finitely distributed delays and eventual vaccination with a pseudo-mass action type, where the infective transmission rate does not depend directly on the total population. The proposed model was analyzed qualitatively, and they further proposed general time varying feedback vaccination rules. In this paper, we present a new deterministic model with an influx of malaria parasites via infected immigrants as well as vaccination. We investigate the effectiveness of the malaria vaccine called RTS, S/AS01 or Mosquirix, in the presence of infected immigrants. The name RTS is due to the fact that the vaccine was manufactured using genes from the repeat ("R") and T-cell epitope ("T") of the pre-erythrocytic circumsporozoite protein (CSP) of the P. falciparum, together with a surface antigen ("S") of the hepatitis B virus, and AS01 is a chemical adjuvant added to it to increase the immune system response ^[13]. Our proposed model consists of a system of first-order ordinary differential equations. When there is zero influx of infectives and ${R_0 < 1},$ the model was shown to have a disease-free equilibrium point that is globally asymptotically stable. The existence of endemic equilibria for the system is deduced and the local dynamics are proved using the center manifold theorem. Numerical simulations are also performed to investigate the effect of immigration in controlling the malaria disease.

1.1.   Malaria vaccination program in Sub-Saharan Africa

Malaria interventions have been of great use over the years but have mostly not succeeded in eliminating the disease, and these efforts are being challenged by insecticide and drug resistance ^[15]. In complementing anti-malarial methods, the European Medicines Agency (EMA) ^[16] approved the use of a malaria vaccine known as RTS, S/AS01, or Mosquirix. This vaccine candidate is intended for use against malaria caused by the P. falciparum parasite to reduce severe cases of the disease in infants aged from six weeks to 17 months in areas where malaria is dominant ^[16]. Formal recommendations were issued by the WHO ^[17] on the pilot implementations of the four-dose schedule of Mosquirix in children below the age of five years. This pilot project was carried out in three Sub-Saharan African countries, namely Ghana, Kenya and Malawi, where three doses were given to children between one and 59 months old, at least four weeks apart from the first dose, with the fourth dose being issued 15–18 months after the third dose ^[18].

1.2.   Model formulation

At a particular time $t,$ the total human population under study is of size $N(t).$ The human population is divided into susceptible individuals $(S)$ who can be infected; fully vaccinated individuals $(V)$ who have a reduced chance of being infected; infectious individuals $(I)$ who have malaria and can infect mosquitoes and recovered individuals $(R)$ who are free of malaria and have gained a partial immunity, such that $N(t) = S(t)+V(t)+I(t)+R(t).$ The total vector population is constant and is given by $L(t).$ It is divided into susceptible mosquitoes $(M)$ and infectious mosquitoes $(J)$ such that $L(t) = M(t)+J(t).$ When susceptible and vaccinated individuals interact with infectious mosquitoes, they get infected. The rate at which vaccinated individuals become infected is much less than the rate at which unvaccinated susceptible individuals may be infected. The movement of individuals and mosquitoes between compartments and within the two sub-populations is given in Figure 1. Our model differs from existing models in the literature, such as, ^{[10,11,19,20,21]} in that it studies the dynamics of the malaria disease in the presence of immigration and vaccination simultaneously. This makes our model more realistic with respect to the Sub-Saharan African region, as the region is characterized by movement of both infected and uninfected individuals across different regions; at the same time, authorities are working on curbing malaria infections by using different control measures including vaccination.

1.3.   Model equations

The population of susceptible individuals is increased by a constant recruitment $\mu H$ and the natural recovery of infected individuals at a rate $h.$ We assume that susceptible individuals can contract the malaria parasite through a successful bite by an infectious female Anopheles mosquito at a rate $a,$ where each bite has a probability $b$ of bringing an infection. These individuals will die a natural death at a rate of $\mu.$ The class of vaccinated children is increased by vaccination of susceptible humans at a rate $\nu.$ Vaccinated children contract the disease when they are bitten by an infectious mosquito at a rate $a,$ where each bite has a probability $f$ of bringing an infection, where $f < < b.$ These individuals are decreased by a proportion progressing back to being susceptible when the antibody titers decrease at a rate $g$ and a natural mortality, at a rate $\mu.$ The compartment of infectious individuals is increased by both susceptible and vaccinated individuals who have been bitten by infectious mosquitoes $J$ at a rate $a$. Infected immigrants increase the population of humans in this class at a rate $\mu G.$ This class is reduced by a natural death at a rate $\mu,$ death due to the disease at a rate $\delta$ and recovery at a rate $k.$ The class of recovered individuals is increased by the recovery of infectious individuals at a rate of $k.$ We assume a recovery with temporary immunity, meaning recovered individuals have a chance of being susceptible again at a rate $h.$ Recovered individuals may die a natural death at a rate $\mu.$ On the other hand, the population of susceptible mosquitoes is increased by a constant birth rate at a rate $\theta L.$ When a susceptible mosquito bites an infectious individual $(I)$ at a rate $a$, with a probability of bringing an infection $c,$ it eventually become infectious and will progress to the class of infectious mosquitoes $(J).$ Susceptible mosquitoes will die a natural death at a rate $\theta$ or will become infectious through a force of infection $ac.$ The population of infectious mosquitoes increases when susceptible mosquitoes interact with infectious humans. These mosquitoes acquire the infection through a successful bite at a rate $a,$ with a probability of infection $c.$ The model equations describing the dynamics of malaria are given by system (1.1):

with initial conditions $S(0) = S_{0},$ $V(0) = V_{0},$ $I(0) = I_{0},$ $R(0) = R_{0},$ $M(0) = M_0,$ $J(0) = J_0,$ where

2.   Model parameters and variables

In this section we give a summary of the descriptions of the parameters and variables used in system (1.1). We also give the values of the parameters. The natural death rate of humans is calculated as

where the life expectancy is 66.2 years, with Kenya as an example ^[22].

3.   Positivity of solutions

System (1.1) describes the population of living organisms; therefore, it is important to show that all state variables ${S(t), V(t), I(t), R(t), M(t), J(t)}$ are always positive for every time ${t\geq 0}.$

Theorem 3.1. The solutions $S(t), V(t), I(t), R(t), M(t), J(t)$ of system (1.1) exist for all time and are always positive, smooth and unique for all $t > 0,$ subject to v,G >0 and the initial conditions $S_0 > 0, V_0 > 0, I_0 > 0, R_0 > 0, M_0 > 0, J_0 > 0.$

Lemma 3.2. Let $U \subseteq \mathbb{R}^n$ be open and consider a solution ${\boldsymbol{x}} = (x_1, \dots, x_n) \in C^1((0, T);U) \cap C([0, T);U)$ to the ODE

where ${\boldsymbol{x}}_0 = (x_{1, 0}, \dots, x_{n, 0}) \in U,$ with $x_{i, 0} > 0$ for all $i \in \{1, \dots, n\}$ and ${\boldsymbol{f}} = (f_1, \dots, f_n):(0, T) \times U \to \mathbb{R}^n$. With the property that $f_i(t, {\boldsymbol{y}}) > 0$ for all y∈U such that $\boldsymbol{y} \in \mathbb{R}_{>0}^{i-1} \times\{0\} \times \mathbb{R}_{\geq 0}^{n-i}$ , $t \in(0, T)$ and $i \in \{1, \dots, n\},$ then ${\boldsymbol{x}}$ is positive for all time, i.e., $x_i(t) > 0$ for all $t \in [0, T)$ and $i \in \{1, \dots, n\}$.

Proof. of Lemma 3.2.

Suppose ${\boldsymbol{x}}$ is not positive for all time, then by continuity there exists $s \in (0, T)$ and $j \in \{1, \dots, n\}$ such that $x_j(s) = 0$ and, thus, the set of all the zeros of all the $x_i;$ that is, $\cup_{i = 1}^n x_i^{-1}(\{0\})$, is closed, nonempty and bounded below and, thus, has a minimum zero, for $\tau \in (0, T)$. Therefore, we have that $x_i(\tau) \geq 0$, $x_i(t) > 0$ for all $t \in (0, \tau)$ and $i \in \{1, \dots, n\}$.

Thus, we have that $\left(x_2(\tau), \ldots, x_n(\tau)\right) \in \mathbb{R}_{\geq 0}^{n-1}$ and, if $x_1(\tau) = 0$, then we have satisfied the conditions to apply the property of ${\boldsymbol{f}}$ to get $f_1(\tau,{\boldsymbol{x}}(\tau)) > 0$. Thus,

By continuity, there exists $0 < \varepsilon < \tau$ such that $\dot{x}_1(t) > 0$, for all $t \in (\tau - \varepsilon,\tau + \varepsilon)$. Thus, by the fundamental theorem of calculus, we have that

This is a contradiction, as ${x_1(t) \geq 0}$ for all ${t \in (0, \tau)}$ and, thus, ${x_1(\tau) > 0}.$

Suppose $\left(x_1(\tau), \ldots, x_{i-1}(\tau)\right) \in \mathbb{R}_{>0}^{i-1}$ and $\left(x_i(\tau), \ldots, x_n(\tau)\right) \in \mathbb{R}_{\geq 0}^{n-i+1}$ for $i \in \{2, \dots, n\}$. If $x_i(\tau) = 0$, then similarly to before, we satisfy the conditions to apply the property of ${\boldsymbol{f}}$ to get $f_i(\tau, {\boldsymbol{x}}(\tau)) > 0$. The same as before, we derive a contradiction and conclude $x_i(\tau) > 0$ and that $\left(x_1(\tau), \ldots, x_i(\tau)\right) \in \mathbb{R}_{>0}^i$ and $\left(x_{i+1}(\tau), \ldots, x_n(\tau)\right) \in \mathbb{R}_{\geq 0}^{n-i}$. Thus, by induction, we obtain that $x_1(\tau), \dots, x_{n}(\tau) > 0$. This contradicts that $\tau \in \cup_{i = 1}^n x_i^{-1}(\{0\})$; thus, there must be no zeros, so ${\boldsymbol{x}}$ is positive for all time.

□

Next, we prove Theorem 3.1.

Proof. of Theorem 3.1

We consider the ODE system

We define the set $U = B_{N_0 + L + 1}$ as the open ball of radius $N + L + 1$, where $N_0$ is the human population and $L$ is the mosquito population; we can see from the ODE that these remain constant. We see that $x_0 \in U,$ and that on $U$, $f$ is a bounded function, Lipschitz continuous and is smooth. By the Picard-Lindelöf theorem, there exists a $T > 0$ and unique solution ${\boldsymbol{x}} = (x_1, \dots, x_n) \in C^1((0, T);U) \cap C([0, T);U)$ to the system of ODEs with given initial conditions. Since the system of ODEs is in the format of Lemma 3.2 with $n = 6$, we have the positivity of the solutions of system (1.1). Since we have positivity of solutions and $N(t) = S(t)+V(t)+I(t)+R(t)$ and $L = M(t) + J(t)$, then ${\boldsymbol{x}}(t) \in \overline{B_{N_0+L+H+G}} \Subset U$ bu Eq (4.2), so by the Corollary in section 17.4 of the book ^[24], we have that the solution is defined for all time and we have a global positive unique solution. By applying a bootstrap argument, we obtain that the solution is smooth in time. That is, since ${\boldsymbol{f}}$ is smooth and

if ${\boldsymbol{x}}(t)$ is $C^n$, then from (3.1) we have that ${\boldsymbol{x}}(t)$ is $C^{n+1}$. Thus, we obtain that ${\boldsymbol{x}}(t)$ is $C^\infty$. □

4.   Well-Posedness of the model in a feasible region

We aim to analyze system (1.1) in a suitable feasible region $\mathcal{D}$ of biological interest.

Theorem 4.1. The feasible region $\mathcal{D} = \mathcal{D}_H\times \mathcal{D}_V,$ where

is attracting and positively invariant with respect to system (1.1) for all values of $t\geq 0.$

Proof.

Substituting values of the derivatives of state variables of system (1.1), we get

Solving Eq (4.1) by using the integrating factor ${I_f = e^{\int_0^t \mu d\tau} = e^{\mu t}},$ we have:

Evaluating the integrals on both sides, we get:

which reduces to:

For all $t\geq 0,$ the total population $N(t) \leq (H+G)$ if $N(0) \leq (H+G).$ Thus, the region $\mathcal{D}_H$ is positively invariant. If $N(0) > (H+G),$ then $N(t)$ approaches $\mathcal{D}_H$ asymptotically. Thus, the region $\mathcal{D}_H$ attracts all solutions in $\mathbb{R}^4_+.$

Similarly, for the vector population,

The mosquito population will remain constant and will stay in $\mathcal{D}_V$ for all $t\geq 0.$

Thus, $\mathcal{D}$ is positively invariant and system (1.1) is well-posed mathematically and biologically; hence, it is sufficient to study the dynamics of system (1.1) in $\mathcal{D}.$ □

5.   Disease-free equilibrium

It is very important to check if the disease under study has steady states. These are points at which the disease is neither increasing nor decreasing. We equate the changes in the state variables to zero and we solve for the corresponding state variables such that:

Observing that ${M+J = L\implies M = L-J},$ system (1.1) can be reduced to a system of five equations by eliminating the equation for ${\frac{dM}{dt}}$ so that we have

The solutions ${(S^*, V^*, I^*, R^*, M^*, J^*)}$ of equations described in Eq (5.2) are the steady states of system (1.1).

NB: In the rest of this section we restrict to the case when there is no influx of infected immigrants, i.e., $(G = 0).$

In this case, we have a disease-free equilibrium point

The reproductive number can be defined as the average number of new infections of the disease caused by one typically infected individual or mosquito in a completely susceptible population. We make use of the next-generation matrix method by Van den Driessche and Watmough ^[25]. This is done by considering transmissions and transitions between and within classes. Considering infectious classes only, we have

Taking partial derivatives, we define the matrix $F$ as the matrix of new infections and matrix $V$ as the matrix of transitions, such that

$R_0 = \rho(FV^{-1}),$ where $\rho$ is the spectral radius or largest eigenvalue. Thus,

● The term ${R_H = \sqrt{\frac{aH(bg+b\mu+f\nu)}{(k+\delta+\mu)(g+\mu+\nu)}}}$ is the average number of new infections caused by one infectious human in a completely susceptible mosquito population.

● The term ${R_V = \sqrt{\frac{acL}{\theta}}}$ is the average number of new infections caused by one infectious mosquito in a completely susceptible human population.

Theorem 5.1. The disease-free equilibrium point is locally asymptotically stable when $R_0 < 1,$ otherwise it is unstable.

Proof. To prove the local dynamics of our proposed model, we compute the Jacobian matrix $J$ evaluated at the disease-free equilibrium point and we obtain:

The eigenvalues of $J$ are:

Following Allen ^[26], all eigenvalues are negative if they are real, or have negative real parts if they are complex, except for $\lambda_6,$ which is negative or has a negative real part whenever the term

Thus the disease-free equilibrium point is locally asymptotically stable when $R_0 < 1.$ □

Theorem 5.2. The disease-free equilibrium point is globally asymptotically stable when $R_0 < 1.$

Proof. We use a function of the form ${P(I, J) = \alpha_1 I+\alpha_2 J},$ where $\alpha_1, \alpha_2 > 0$ are constants to be determined. Since we proved in Lemma (3.1) that $I, J > 0,$ then $P$ is positive definite. Taking first derivatives on $P,$ we get

Substituting the right hand side of the third and sixth equations of system (1.1), we get;

Evaluating $\dot{P}$ at the disease-free equilibrium point, we get

Equating the coefficients of ${I}$ and ${J}$ to zero and solving for the nontrivial values of ${\alpha_1}$ and ${\alpha_2}$ by setting $\alpha_1 = 1,$ we have ${\alpha_2 = \frac{aH[bg+b\mu+f\nu]}{\theta(g+\mu+\nu)}}.$

Simplifying $\dot{P},$ we have

We note that ${\dot{P} < 0}$ if, and only if, ${R^2_0 < 1}\implies {R_0 < 1}.$ When $R_0 = 1,$ then $P = 0.$ Thus, by ^[27], the disease free equilibrium point is globally asymptotically stable when ${R_0 < 1}.$ □

6.   Endemic equilibrium points

In our study of the existence of endemic equilibrium points, we encounter the number $R_i^2$, which is:

This invariant is the indicator that determines the existence, or not, of a unique endemic equilibrium point. Note that if $G = 0$, then $R_i = R_0$.

Theorem 6.1. (a) Suppose that $G = 0$. There exists a unique positive equilibrium point if $R_0 > 1$.

(b) Suppose that $G > 0$. There exists a unique positive equilibrium point if, and only if, $R_i > 1$.

Proof. At the steady states, we equate the derivatives of the state variables of system (1.1) to zero, with $G > 0$. Since $M = L-J,$ we can eliminate the fifth equation of system (1.1) so that we have

From the last equation of system (6.2), we have

By back substitution, we get

${V^* = \frac{\nu(acI^*+\theta)^2[\mu H(\mu+h)+hkI^*]}{(\mu+h)[(a^2bcLI^*+(\mu+\nu)(acI^*+\theta))(a^2fcLI^*+(\mu+g)(acI^*+\theta))-\nu g(acI^*+\theta)^2]}},$

${S^* = \frac{[a^2fcLI^*+(\mu+g)(acI+\theta)][\mu H(\mu+h)+hkI^*](acI^*+\theta)}{(\mu+h)[(a^2bcLI^*+(\mu+\nu)(acI^*+\theta))(a^2fcLI^*+(\mu+g)(acI^*+\theta))-\nu g(acI^*+\theta)^2]}}.$

Substituting the values of ${S^*, V^*, R^*, J^*}$ into the equation for infective humans, we get the polynomial $P(I^*)$, with

$P(I^*) \ = \ G \mu (h+\mu) \{[a^2 b c L I^*+(\mu +\nu) (a c I^*+\theta)][a^2 c f L I^*+(g+\mu) (a c I^*+\theta)]$

$-g \nu (a c I^*+\theta)^2\}+a^2 b c L I^* \{[H \mu (h+\mu)+h k I^*] ([a^2 c f L I^*+(g+\mu) (a c I^*+\theta)]\}$

$-(\delta+k+\mu) I^* (h+\mu) \{[a^2 b c L I^*+(\mu +\nu) (a cI^*+\theta)] [a^2 c f L I^*+(g+\mu) (a c I^*+\theta)]$

$-g \nu (a c I^*+\theta)^2\}+a^2 c f L \nu I^* (a c I^*+\theta) [H \mu (h+\mu)+h k I^*] = 0.\quad$ This simplifies to,

where

$\mathcal{A}_3\ = \ -a^2c^2\{[(\mu+\delta)(\mu +h) +\mu k] [a^2bfL^2+abL(\mu+g)+afL\nu]+ \mu(\mu+h)(\mu+\delta+k)[afL+\mu+g +\nu] < 0,$

$\mathcal{A}_2 = a c (a (b L (h \mu (a c (G+H) (a f L+g)+\theta k)+a c \mu ^2 (G+H) (a f L+g+h)$

$+a c \mu ^3 (G+H)+g h \theta k)+c \mu (h+\mu) (G \mu (a f L+g+\mu)+a f L \nu (G+H)$

$+G \mu \nu)+f h \theta k L \nu)-\theta (\delta+k+\mu) (h+\mu) (a L (b (g+\mu)+f (\mu +\nu))+2 \mu (g+\mu +\nu))),$

${\mathcal{A}_1 = ac\theta^2\mu(h+\mu)(g+\mu+\nu)(\delta+k+\mu)\left[R_i^2-1\right]},$

${\mathcal{A}_0 = G \theta ^2 \mu ^2 (h+\mu) (g+\mu +\nu) > 0}.$

For case (b), the number of positive solutions of Eq (6.3) depends on the signs of $\mathcal{A}_3, \mathcal{A}_2, \mathcal{A}_1$ and $\mathcal{A}_0.$ We apply Descartes rule of signs by ^[28] to investigate the number of sign changes of Eq (6.3), which will be the possible solutions of the polynomial. Possible signs of the coefficients of Eq (6.3) are shown in Table 2. It can be observed that $\mathcal{A}_0 > 0$ and $\mathcal{A}_3 < 0$. With $R_i > 1$, we have $\mathcal{A}_1 > 0$. By the Descartes rule of signs, it follows that $P(I^{*})$ has a unique positive root. Thus, the statement (b) is proved.

Now, for case (a) with $G = 0$, the polynomial of interest becomes:

It can be observed that $\mathcal{A}_3 < 0$ with $R_0 > 1,$ thus $\mathcal{A}_1 > 0$. By the Descartes rule of signs, it follows that $Q(I^{*})$ has a unique positive root. Thus, the claim (a) is also true. □

7.   Bifurcation analysis

We first present the case when there is no immigration of infectives, i.e. $G = 0$. This is done by applying a change of variable to system (1.1):

We choose the bifurcation point $a = \phi^*\implies a^2 = \phi^{*2}$ and $R_0 = 1.$

System (1.1) becomes

The Jacobian matrix of system (7.1) evaluated at the disease-free equilibrium point is

The eigenvalues of the system at $R_0 = 1$ are

We note that the system has a simple eigenvalue, while the other eigenvalues have negative real parts. The right eigenvector associated with the simple eigenvalue is given by $w = (w_1, w_2, w_3, w_4, w_5, w_6),$ where

Similarly, the left eigenvector corresponding the simple eigenvalue is given by $v = (v_1, v_2, v_3, v_4, v_5, v_6),$ where

The nonzero second order partial derivatives are:

Substituting the values of ${w_1, w_2, w_3, w_5, w_6, v_3, v_6}$ and evaluating $b$ at the disease-free equilibrium point, we get,

We note that $b$ is always positive. The type of bifurcation of system (1.1) depends on the sign of $a.$

Theorem 7.1. Following ^[29],

(i) When $a > 0,$ meaning,

system (1.1) undergoes a backward bifurcation when $R_0 = 1,$ $a = \phi^*.$

(ii) When $a < 0,$ meaning,

system (1.1) undergoes a transcritical bifurcation when $R_0 = 1,$ $a = \phi^*.$

Therefore the endemic equilibrium point is locally asymptotically stable when $R_{0} > 1$ but sufficiently close to one.

Now we consider the case $G > 0$.

Since $\mu G$ is a constant, its partial derivatives will be zero. We inherit the proof of the local dynamics of the endemic equilibrium points when $G > 0$ and the bifurcation analysis of system (1.1) from the center manifold studied for the case when $G = 0.$ Thus, the endemic equilibrium point for system (1.1) when $G > 0$ is locally asymptotically stable.

8.   Sensitivity analysis

In order to implement control strategies, it is important to know the parameters that influence the basic reproductive number. This will help in providing a clear picture regarding the factors that need to be addressed in order to reduce malaria incidents. In the event of a low budget, sensitivity analysis helps us identify key factors that can be prioritized in the control of a disease. In this section we carry out sensitivity analysis of system (1.1) using the Latin hypercube sampling (LHS) method described in ^{[30,31,32,33]}, with 1000 simulations per run. The partial rank correlation coefficient (PRCC) is an efficient and reliable sampling-dependent method that performs global sensitivity analysis, which allows us to investigate the monotonicity between parameters of the model and the model output, or in this case, the reproductive number when all parameters are varied ^[33,34]. In general, a standard correlation coefficient $\rho$ for any two variables, $x$ and $y,$ is calculated using the formula:

where $\{(x_i, y_i):x_i\in x, \; y_i\in y\}$ are a set of paired sampled data ^[33]. The method of PRCCs is explained in detail in ^[34]. Figure 2 shows the tornado plot of the PRCCs of the parameters of our proposed model.

Parameters with positive PRCCs increase the basic reproductive number $R_0$ when they are increased, while parameters with negative PRCCs decrease the reproductive number when they are increased. Parameter values used are displayed in Table 2. It can be observed that factors such as the mosquito biting rate, probabilities of successful bites leading to both new mosquito infections and human infections are the most sensitive factors. Interventions should be introduced to reduce these factors. On the other hand, factors such as the mosquito mortality rate and vaccination rate are also sensitive to the reproductive number. Increasing these factors will decrease the reproductive number and, hence reduce malaria, provided there is no immigration of infected humans.

The scatter plots shown in Figure 3 show the relationship between the most sensitive parameters and the $\log R_0.$ It can been seen that the mosquito biting rate and the probability of new mosquito infections are monotonically increasing, meaning if we increase these parameters, the disease will also increase. Similarly, the mosquito mortality rate and vaccination rate monotonically decreasing and increasing their values will lead to a decrease in the number of new cases.

9.   Numerical simulations

In this section we carry out numerical simulations for our system (1.1). Parameters used in the simulations are obtained from literature, with some being reasonable estimates. The values are listed in Table 2. We start by investigating the impact of immigration of infectives in the population.

10.   Effect of controlling the influx of infected immigrants

One of our objectives is to investigate the impact of influx of infected immigrants into the population. Figure 4 below shows the population of infectious individuals when there is a constant influx of infected immigrants and the case when there is no influx of infected immigrants. The Python notebook for the codes used to produce the plot below can be obtained here.

It can be observed that the population of infectious human beings can be reduced when there is no influx of infected immigrants. In the presence of a constant flow of infected immigrants, the population of infectious humans will increase and malaria will remain a pandemic and a huge threat to human lives. The shaded region between the two graphs indicates the number of infections that can be prevented by curbing or reducing the influx of infected humans.

10.1.   Vaccine coverage when there is no immigration

We investigate the impact of introducing the mosquirix malaria vaccine when there is no influx of infected immigrants. This is done by varying values of $\nu,$ the rate at which susceptible humans are vaccinated. Figure 5 below indicates the population of infectious humans in the absence of the influx infected immigrants. The Python notebook for the code used to produce this plot can be obtained here.

10.2.   Vaccine coverage when there is immigration

Immigration is part of our daily lives. It is difficult for people not to migrate because of other factors that can include poverty and search for greener pastures. We investigate the effectiveness of the malaria vaccine when there is an influx of infected immigrants.

It can be observed that in the presence of infected immigrants, a higher vaccine coverage also helps in reducing malaria infections. This is evidenced by a decrease in the number of cases, as shown in Figure 6. The disease will remain endemic because of the influx of infected humans, but vaccination reduces the rate at which new infections emerge in the population. The code used to produce this plot can be obtained on this link.

11.   Optimal control

In this section we carry out an optimal control analysis of different strategies to control malaria. We wish to choose the best scenario that, if well implemented, could help control the spread of malaria. To investigate the potential impact of the implemented intervention measures, the following control variables are incorporated into our model system:

● $p_1(t):$ Represent the use of personal protection measures to prevent mosquito bites during the day and night, such as the use of insecticide-treated nets, application of repellents to skin or spraying of insecticides,

● $p_2(t):$ Represent the treatment and

● $p_3(t):$ Represent the use of vaccination to prevent malaria.

The model with optimal control is shown in Figure 7 and system (11.1):

We want to investigate the impact of the control strategies at different levels of immigration. We consider the following quadratic objective functional:

The objective functional measures the cost of our proposed strategies from time $t = 0$ to the terminal time $t = T.$ For realism, we need to know the unit cost of the intervention, but for now, the cost is not available and, hence, we base our experiments on nominal values. The constants $M_1, M_2, M_3$ are the weight constants for the state variables and are positive for all time $t,$ while $\omega_1, \omega_2, \omega_3$ are the constants for the three strategies. We derive the necessary conditions that must be satisfied by an optimal control problem from the Pontryagin's maximum principle ^[35,36], which converts system (11.1) and the objective functional defined in Eq (11.2) into the problem of minimizing pointwise Hamiltonian $\mathcal{H},$ with respect to the control variables ${\left(p_1, \; p_2, \; p_3\right)\in \mathcal{U}}$ at a minimum cost:

where ${\mathcal{U}}$ is the set of admissible controls. The Hamiltonian $(\mathcal{H})$ is given by

where the variables $\lambda_1, \lambda_2, \lambda_3, \lambda_4, \lambda_5, \lambda_6$ are the adjoints to the state variables $S, V, I, R, M, J,$ respectively. The adjoint variables should satisfy

implying that:

with transversality conditions

Moreover, $p_{1}^{*}, p_{2}^{*}, p_{3}^{*}$ satisfy the conditions

implying that;

in the interior of the control set $\mathcal{U}.$ Solving Eq (11.4) for ${(p_{1}, p_{2}, p_{3})}$ in ${\mathcal{U}},$ we get the optimal solution as:

11.1.   Numerical simulations of the optimal control

Using parameters described in Table 1, we obtain the following simulations for the different strategies under consideration. The code used in these plots can be obtained here.

It can be observed that even though the use of personal protective equipment, treatment and vaccination of young children can help reduce the burden of malaria, the combined use of these strategies is the best way to achieve better results.

When there is immigration of infected humans, the three strategies do reduce malaria incidence, but malaria remains endemic in the population. Even the combined use of all strategies at the same time is not enough to completely eliminate malaria because all these available strategies are not $100\%$ efficient. Thus, the best strategy is to use the three strategies at the same time, but making sure that immigration of infected humans is stopped.

12.   Conclusions and discussions

We have presented a new model for which we have proved that the solutions are well-behaved and the equilibria has good stability properties. Existence of endemic equilibria for system (1.1) was deduced and the local dynamics were proven using the center manifold theorem. When there is immigration, the model had no disease free equilibrium point, indicating that malaria would remain endemic in the population. Numerical simulations were also performed to investigate the effect of immigration in controlling the malaria disease. The results showed that completely stopping the influx of infected immigrants will be very useful in reducing the burden of malaria in malaria endemic regions. We performed a sensitivity analysis to $R_0$ and our results showed that interventions must be introduced with an aim of reducing contacts between mosquitoes and humans. An investigation of the vaccine coverage indicated that when there is no influx of infected immigrants, vaccination will be essential in eliminating malaria. We then analyzed the effectiveness of the vaccine by considering different rates of coverage of the vaccine when there is an influx of immigrants. Results from this analysis showed that a wide coverage rate will reduce infections, although malaria will remain endemic. Following the results from the sensitivity analysis, we introduced optimal control strategies to curb malaria incidence. The strategies proposed are the use of personal protection measures to prevent mosquito bites during the day and night, such as the use of insecticide-treated nets, application of repellents to skin or spraying of insecticides, treatment of infected humans and vaccination of young children below the age of five years. Results from this optimal control showed that the combined use of the three strategies is the best way to reduce malaria incidence, provided the influx of infected humans is completely stopped. In conclusion, it is evident that a malaria free population can be achieved when more children get vaccinated, infected people get treatment, people living in malaria endemic regions continue using personal protective measures that prevent mosquito bites during the day and night and when immigration of infected humans is completely stopped. The modeling can be improved by also considering the effect of climatic factors on malaria dynamics. Several such models, but without vaccination or immigration of infected humans, are already in the literature, e.g., ^[19,20,37]. Also, when malaria spreads within a population, individuals gain knowledge and a better understanding of the disease. Thus, this work can be extended by means of investigating the memory effect on the dynamics of malaria, using fractional order differential equations described in ^{[38,39,40,41]}. It is also interesting to consider stochastic models and their optimal control, as described in ^[42].

Use of AI tools declaration

The authors declare they have not used any Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

Pride Duve would like to acknowledge the African Institute for Mathematical Sciences (AIMS) South Africa for the scholarship he obtained during his Masters degree studies, which made it possible for him to carry out this work. Furthermore, the project is funded by the Federal Ministry of Education and Research of Germany (BMBF) under the project NEED (Grant Number 01Kl2022) and the German Research Foundation (JO 1276/5-1).

Conflict of interest

All authors declare no conflicts of interest in this paper.