Optimal control of a two-group malaria transmission model with vaccination

Abstract

Malaria is a vector-borne disease that poses major health challenges globally, with the highest burden in children less than 5 years old. Prevention and treatment have been the main interventions measures until the recent groundbreaking highly recommended malaria vaccine by WHO for children below five. A two-group malaria model structured by age with vaccination of individuals aged below 5 years old is formulated and theoretically analyzed. The disease-free equilibrium is globally asymptotically stable when the disease-induced death rate in both human groups is zero. Descarte’s rule of signs is used to discuss the possible existence of multiple endemic equilibria. By construction, mathematical models inherit the loss of information that could make prediction of model outcomes imprecise. Thus, a global sensitivity analysis of the basic reproduction number and the vaccination class as response functions using Latin-Hypercube Sampling in combination with partial rank correlation coefficient are graphically depicted. As expected, the most sensitive parameters are related to children under 5 years old. Through the application of optimal control theory, the best combination of interventions measures to mitigate the spread of malaria is investigated. Simulations results show that concurrently applying the three intervention measures, namely: personal protection, treatment, and vaccination of childreen under-five is the best strategy for fighting against malaria epidemic in a community, relative to using either single or any dual combination of intervention(s) at a time.

1 Introduction

Malaria is a vector-borne disease transmitted during blood meal by infectious Anopheles mosquitoes via insertion of sporozoites in the blood of susceptible humans (Niare et al. 2016). The worldwide casualties of malaria are huge, particularly in sub-Sahara Africa, and especially in the children population (WHO 2022c). The World Health Organization (WHO) estimated about 241 million malaria cases globally in 2020 that resulted in 627, 000 deaths. The tropical and sub-tropical regions (mainly Africa) carry a disproportionately high share of the global malaria burden with as much as $95\%$ of malaria cases, $96\%$ of malaria deaths with about $80\%$ of these deaths are children less than 5 years of age (WHO 2022b).

Malaria is a severe life-threatening vector-borne disease, and with the ongoing COVID-19 pandemic, malaria morbidity and mortality could increase (Tchoumi et al. 2022). While several prevention and therapeutic measures have been implemented to fight against this deadly disease, the recent groundbreaking malaria vaccine for children now recommended by the WHO (2022c) after successful pilots in Ghana, Kenya, and Malawi, is expected to help strengthen the fight against malaria infection (WHO 2022a). Also, the development of Transmission Blocking Drugs (TBDs) conferring protection against malaria will play a significant role in mitigating morbidity and mortality in malaria-prone regions (Woldegerima et al. 2021). It is expected that the widespread use of the RTS malaria vaccine (trade name Mosquirix), a recombinant protein-based malaria vaccine in children in sub-Saharan Africa and other areas with moderate or high transmission of Plasmodium falciparum malaria along with other preventive measures could help mitigate the spread and eventually the eradication of the disease (WHO 2022a).

Since the eighteen century, the development of mathematical models has been critical to provide framework and understanding of the dynamics of infectious diseases (Bernoulli 1760; Kermack and McKendrick 1927; Dietz and Heesterbeek 2002). Several mathematical models of malaria dynamics investigating various aspects of the disease have flourished in the literature Herdicho et al. (2021) . While mosquito’s population fluctuates between climatic seasons, seasonal factor tends to impact on the dynamics of infected mosquitoes and human populations in regions with hot climate (Herdicho et al. 2021; Tasman 2021). Because of mathematical tractability and convenience, though important, we will not account for seasonality in the birth rate of mosquitoes. For several decades, concerted global stringent efforts have been underway to develop effective and safe vaccine for use against malaria in humans (Forouzannia and Gumel 2014; Karunamoorthi 2014), with several candidate vaccines targeting different stages of the malaria parasite’s lifecycle (White et al. 2015). An overview of integrated mathematical models for predicting the epidemiologic and economic effects of malaria vaccines on the clinical epidemiology and natural history of Plasmodium falciparum malaria both at the individual and population level has been reported in Smith et al. (2006) and Atcheson et al. (2019). It is noted that these models provide a unique platform for predicting both the short- and long-term effects of malaria vaccines on the burden of disease, allowing for the temporal dynamics of effects on immunity and transmission. With mathematical models being increasingly used to inform decisions throughout product development pathways from pre-clinical studies to country implementation of novel health interventions, Galactionova et al. (2021) illustrate the utility of simulation approaches by reviewing malaria vaccine modelling studies. A mathematical model of vaccine combination adapted to murine malaria studies based on simple probabilistic assumptions was developed in Atcheson et al. (2019).

Transmission-blocking vaccines of malaria have been investigated in Zhao et al. (2022) and Takashima et al. (2021) where as expected, vaccination has a positive impact on reducing the disease burden, while malaria could be controlled if the duration of efficacy is in the order of a human life-span (Koella 1991) . Public health impact of dynamics model of a transmission-blocking vaccine alongside existing interventions suggests that school-aged children are an attractive population to target for vaccination (Challenger et al. 2021). That is, benefit of vaccination distributed across the population averts the greatest number of cases in younger children. Even an imperfect anti-malaria vaccine (with a modest efficacy and coverage rate) can lead to effective disease control (Teboh-Ewungkem et al. 2010). Handari et al. (2022) analyze an optimal control model of malaria incorporating pre-erythrocytic vaccine and transmission-blocking treatment.

Previous studies have been hypothetical because there was no approved/licensed malaria vaccine, and also vaccination was applied to the entire population. As predicted in Challenger et al. (2021) and Hill (2011), the first approved malaria vaccine is for children age five and below, and consequently, we propose a mathematical model of malaria transmission dynamic by extending our previous work Tchoumi et al. (2022), incorporating vaccination in the class of children less than 5 years old, but with no seasonal birth rate (Herdicho et al. 2021). The population of children 5 months (minimum age to take the vaccine) to 5 years are vaccinated at a certain rate. The vaccinated class comprises fully vaccinated children (i.e., those who have taken all the 4 doses of the S/AS01 (RTS,S) vaccine). It is important to note that in Tchoumi et al. (2022), the authors studied a two-group malaria model structured by age with no vaccination, while herein, to realistically capture what is currently known about the state of the vaccine development, we incorporate vaccination for children less than 5 years old only to capture the recently approved children vaccine. To the best of our knowledge, this is the first mathematical model investigating the impact of the newly approved children malaria vaccine on the dynamics of the disease.

Here is the outline of the rest of the paper. The mathematical model of our proposed malaria model with vaccination of children under 5 years old with no seasonality is formulated in Sect. 2.1. Theoretical analyses of the model using the fundamental theory of dynamical systems is carried out in Sect. 3. In Sect. 4, we formulate an optimal control problem to investigate the impact of the optimal control strategy on mitigating the spread of the disease. Conditions for the existence of optimal control and the optimality system are established using the Pontryagin’s Maximum Principle. Numerical simulations along with global sensitivity analysis using the vaccinated class as our response variable are carried out in Sect. 5. while Sect. 6 is the conclusion.

2 Malaria model without control

2.1 Model formulation

A deterministic compartmental modeling approach is used to describe the disease transmission dynamics. The model flow is a susceptible-vaccinated-exposed-infected-recovered-susceptible SVEI(R)S malaria model in the human population, and SEI in the mosquitoe population. Susceptibles ($S_e$) under 5 years old are recruited at the constant rate $\Lambda _e$ and can die naturally at the rate $\mu _h$, or grow to become susceptible over 5 years old at the rate $\xi$, becoming vaccinated ($V_e$) at the rate $\vartheta _e$, or infected ($E_e$) after a bite from an infectious mosquito with a strength of infection $\lambda _e$. After a latency period $\dfrac{1}{\sigma _e}$, infected individuals become infectious ($I_e$), they can recover and become susceptible again at the rate $\omega _e$, or die naturally at the rate $\mu _h$ or as a result of an illness at rate $d_1$. Entry into susceptible humans over 5 years old ($S_a$) comes from the growth of susceptible humans under 5 years old, who can die naturally at the rate $\mu _h$, become infected ($E_a$) with the infection strength equal to $\lambda _a$. Infected individuals become infectious ($I_a$) at the rate $\sigma _a$ or die naturally at the rate $\mu _h$. An infectious over 5 years recovers at the rate $\omega _a$ and becomes susceptible ($S_a$), or recovered ($R_a$) with probability $1-p$ and p respectively, or dies naturally at the rate $\mu _h$ or due to illness at rate $d_2$. A treated human over 5 years old can die at rate $\mu _h$, or lose immunity at rate $\delta _a$ to become susceptible. Similarly, following an infectious bite, susceptible mosquitoes ($S_v$) can become infected ($E_v$) with the strength of the infection $\lambda _v$ before becoming infectious themselves ($I_v$) at rate $\sigma _v$. Mosquitoes do not recover from malaria, they all die naturally at the rate $\mu _v$. The schematic flow diagrams of human and mosquito components of the model are depicted in Figs. 1 and 2.

Fig. 1

Model flow diagram of the human component of the model

Fig. 2

Model flow diagram of the mosquito component of the model

2.2 The model

The variables and parameters values of the model are presented in the following Tables 1 and 2, respectively.

Table 1 Variable of model

Table 2 Model parameter

Based on our model description and assumptions, we establish the following system of non-linear ordinary differential equations.

$$\begin{aligned} \left\{ \begin{array}{l} S_{e}' = \Lambda _{e}+\omega _{e}I_{e}-(\lambda _{e}+\vartheta _e+\xi +\mu _{h})S_{e},\\ V_{e}' = \vartheta _e S_{e} -((1-\varepsilon ) \lambda _e +\mu _h)V_e,\\ S_{a}' = \xi S_{e}+u\omega _{a}I_{a}+\delta _a R_{a}-(\lambda _{a}+\mu _{h})S_{a},\\ E_{e}' = \lambda _{e}(S_{e}+(1-\varepsilon ) V_e)-(\sigma _{e}+\mu _{h}) E_{e},\\ E_{a}' = \lambda _{a}S_{a}-(\sigma _{a}+\mu _{h}) E_{a},\\ I_{e}' = \sigma _{e}E_{e}-(\omega _{e}+\mu _{h}+d_1) I_{e},\\ I_{a}' = \sigma _{a}E_{a} -( \omega _{a}+\mu _{h}+d_2) I_{a},\\ R_{a}' = (1-u)\omega _{a}I_{a}-(\delta _a+\mu _{h})R_{a},\\ S_{v}' = \Lambda _{v}-(\lambda _{v}+\mu _{v})S_{v},\\ E_{v}' = \lambda _{v}S_{v}-(\sigma _{v}+\mu _{v})E_{v},\\ I_{v}' = \sigma _{v}E_{v}-\mu _{v} I_{v}, \end{array} \right. \end{aligned}$$

(1)

with initial conditions

$$\begin{aligned} (S_{e}(0),V_{e}(0),S_{a}(0),S_{v}(0),E_{e}(0),E_{a}(0),E_{v}(0),I_{e}(0),I_{a}(0),I_{v}(0),R_{a}(0))\in {\mathbb {R}}_{+}^{11} \end{aligned}.$$

The forces of infections are given by

$$\begin{aligned} \lambda _{e}=\dfrac{b_{2}\beta _{e}}{N_{h}}I_{v},\,\lambda _{a}=\dfrac{b_{1}\beta _{a}}{N_{h}}I_{v},\, \lambda _{v}=\dfrac{\beta _{v}[b_{2}I_{e}+b_{1}I_{a}]}{N_{h}}. \end{aligned}$$

3 Model analysis

3.1 Existence, uniqueness and positivity of solutions

The functions on the right-hand side of the system (1) are Lipschitz continuous, therefore by Picard’s existence Theorem, the system (1) has a solution.

Theorem 3.1

The solutions of the model system (1) with non-negative initial conditions are all non-negative.

Proof

Assume that there exists a time $\tilde{t}$ such that $S_e(\tilde{t})=0,$ $S_e'(\tilde{t})<0,$ $S_e(t)>0,$ $S_a(t)>0,$ $I_e(t)>0,$ $I_a(t)>0,$ $V_e(t)>0,$ $E_e(t)>0,$ $E_a(t)>0,$ $E_v(t)>0,$ $R_a(t)>0,$ $S_v(t)>0,$ $I_v(t)>0,$ for $0<t<\tilde{t}.$ From the first equation of (1), we have

$$\begin{aligned} \dfrac{dS_e(\tilde{t})}{dt} =\Lambda _e+w_e I >0. \end{aligned}$$

This contradicts the assumption that $S_e'(\tilde{t})<0.$ Therefore S(t) is positive.

Similarly, $S_a(t),$ $I_e(t),$ $I_a(t),$ $V_e(t),$ $E_e(t),$ $E_a(t),$ $E_v(t),$ $R_a(t),$ $S_v(t),$ $I_v(t)$ are all positive. $\square$

Theorem 3.2

Solutions the model system (1) are bounded in the invariant region

$$\begin{aligned} \Omega =\left\{ (S_e, S_a, I_e, I_a, V_e, E_e, E_a, E_v, R_a, S_v,I_v,) \in {\mathbb {R}}^{11}_{+}:N_h(t)\le \displaystyle \frac{\Lambda _e}{\mu _h},\, N_v(t)\le \dfrac{\Lambda _v}{\mu _v} \right\} . \end{aligned}$$

Proof

Given $N_h(t)=S_a(t)+S_e(t)+I_e(t)+I_a(t)+V_e(t)+E_e(t)+E_a(t)+R_a(t)$ and $N_v(t)= +E_v(t)+S_v(t)+I_v(t).$

$$\begin{aligned} \dfrac{dN_h}{dt}= & {} (\Lambda _e-\mu _h N_h)-(d_1I_e+d_2I_a), \\\le & {} \Lambda _e-\mu _h N_h. \end{aligned}$$

Solving the differential inequality, we obtain

$$\begin{aligned} N_h(t)\le \dfrac{\Lambda _e}{\mu _h}+\left( \dfrac{\Lambda _e}{\mu _h}-N_{h}(0)\right) \exp {(-\mu _h t)}. \end{aligned}$$

Therefore,

$$\begin{aligned} \limsup _{t\longrightarrow \infty }N_h(t)=\dfrac{\Lambda _e}{\mu _h}. \end{aligned}$$

Since $N_h(t)=S_a(t)+S_e(t)+I_e(t)+I_a(t)+V_e(t)+E_e(t)+E_a(t)+R_a(t),$ it follows that $S_e(t)\le \dfrac{\Lambda _e}{\mu _h},$   $S_a(t)\le \dfrac{\Lambda _e}{\mu _h},\,\,$ $E_e(t)\le \dfrac{\Lambda _e}{\mu _h},$   $E_a(t)\le \dfrac{\Lambda _e}{\mu _h},$   $I_e(t)\le \dfrac{\Lambda _e}{\mu _h},$   $I_a(t)\le \dfrac{\Lambda _e}{\mu _h},$   $V_e(t)\le \dfrac{\Lambda _e}{\mu _h},$   $R_a(t)\le \dfrac{\Lambda _e}{\mu _h}.$

$$\begin{aligned} \dfrac{dN_v}{dt}= & {} \Lambda _v-\mu _v N_v. \end{aligned}$$

Thus, the solution of this differential equation is

$$\begin{aligned} N_v(t) = \dfrac{\Lambda _v}{\mu _v}-\left( \dfrac{\Lambda _v}{\mu _v}-N_{v}(0)\right) \exp {(-\mu _v t)}. \end{aligned}$$

Therefore,

$$\begin{aligned} \limsup _{t\longrightarrow \infty }N_v(t)=\dfrac{\Lambda }{\mu _v}, \end{aligned}$$

and it follows that $S_v(t)\le \dfrac{\Lambda _v}{\mu _v},$  $E_v(t)\le \dfrac{\Lambda _v}{\mu _v}$  and  $I_v(t)\le \dfrac{\Lambda _v}{\mu _v}.$ $\square$

In what follows, for simplicity, let $k=(\vartheta _e+\xi +\mu _{h}),$ $k_{0}=(\sigma _{e}+\mu _{h}),$ $k_{1}=(\sigma _{a}+\mu _{h}) ,$ $k_{2}=(\omega _{e}+\mu _{h}+d_1),$ $k_{3}=( \omega _{a}+\mu _{h}+d_2) ,$ $k_{4}=(\mu _{h}+\delta _a),$ and  $k_{5}=(\sigma _{v}+\mu _{v})$.

3.2 Disease-free equilibrium and basic reproduction number

We first show the existence of a trivial disease-free equilibrium (DFE) for our malaria model without control (MMWC) which is used in computing the basic reproduction number. The MMWC DFE is

$$\begin{aligned} \displaystyle {\mathcal {M}}^{0}=\displaystyle \left( S_{e}^0,V_{e}^0,S_{a}^0,S_{v}^0,E_{e}^0,E_{a}^0,E_{v}^0,I_{e}^0,I_{a}^0,I_{v}^0,R_{a}^0\right) =\left( \dfrac{\Lambda _{e}}{k},\dfrac{\vartheta _e\Lambda _{e}}{\mu _{h}k},\dfrac{\xi \Lambda _{e}}{\mu _{h}k},\dfrac{\Lambda _{v}}{\mu _{v}},0,0,0,0,0,0,0\right) . \end{aligned}$$

To calculate the basic reproduction number, we apply the next generation method in van den Driessche and Watmough (2002). We now have that

$$\begin{aligned} F= & {} \left( \begin{array}{cccccccccc} \displaystyle \lambda _{e}(S_{e}+(1-\varepsilon ) V_e)\\ \displaystyle \lambda _{a}S_{a}\\ \displaystyle \lambda _{v}S_{v}\\ \displaystyle 0 \\ \displaystyle 0\\ \displaystyle 0\\ \end{array} \right) ,\,\,\text {and}\,\, V= \left( \begin{array}{cccccccccc} \displaystyle k_{0}E_{e}\\ \displaystyle k_{1}E_{a}\\ \displaystyle k_{5}E_{v}\\ \displaystyle k_{2}I_{e}-\sigma _{e}E_{e}\\ \displaystyle k_{3}I_{a}-\sigma _{a}E_{a}\\ \displaystyle \mu _{v} I_{v}- \sigma _{v}E_{v}\\ \end{array} \right) \\ {\mathcal {F}}= & {} \left( \begin{array}{cccccccccc} \displaystyle 0&{}0&{}0&{}0&{}0&{}\dfrac{\beta _{e}b_{2}\left( k+(1-\varepsilon )\vartheta _e\right) }{k}\\ \displaystyle 0&{}0&{}0&{}0&{}0&{}\dfrac{\beta _{a}b_{1}\xi }{k}\\ \displaystyle 0&{}0&{}0&{}\dfrac{\beta _{v}b_{2}\Lambda _{v}\mu _{h}}{\mu _{v}\Lambda _{e}}&{}\dfrac{\beta _{v}b_{1}\Lambda _{v}\mu _{h}}{\mu _{v}\Lambda _{e}}&{}0\\ \displaystyle 0&{}0&{}0&{}0&{}0 &{}0\\ \displaystyle 0&{}0&{}0&{}0&{}0&{}0 \\ \displaystyle 0&{}0&{}0&{}0&{}0&{}0 \\ \end{array} \right), \,\,\text {and}\,\, {\mathcal {V}}= \left( \begin{array}{cccccccccc} \displaystyle k_{0}&{}0&{}0&{}0&{}0&{}0\\ \displaystyle 0&{}k_{1}&{}0&{}0&{}0&{}0\\ \displaystyle 0&{}0&{}k_{5}&{}0&{}0&{}0\\ \displaystyle -\sigma _{e}&{}0&{}0&{}k_{2}&{}0&{}0\\ \displaystyle 0&{}-\sigma _{a}&{}0&{}0&{}k_{3} &{}0 \\ \displaystyle 0&{}0&{}-\sigma _{v}&{}0&{}0&{}\mu _{v}\\ \end{array} \right) . \end{aligned}$$

Hence,

$$\begin{aligned} \displaystyle {\mathcal {R}}_{C}^{m}=\displaystyle \sqrt{\dfrac{\Lambda _v \beta _v\mu _h\sigma _v \left[ b_2^2 \beta _ek_1 k_3 \sigma _e(\mu _h+(1-\varepsilon )\vartheta _e)+b_1^2\beta _ak_0k_2\sigma _a\xi \right] }{\Lambda _e \mu _v^2 k k_0 k_1 k_2 k_5} }\cdot \end{aligned}$$

3.3 Global stability of the DFE

To prove the global asymptotically stability (GAS) of the DFE, we use the approach described in Castillo-Chavez et al. (2002). We then re-write the malaria model (1) as follows

$$\begin{aligned} \left\{ \begin{array}{lll} \dfrac{dX}{dt}={F}\left( X,I\right) ,\\ \\ \dfrac{d I}{dt}={\mathcal {G}}\left( X,I\right) ,\,\,{\mathcal {G}}\left( X,0\right) =0,\\ \end{array} \right. \end{aligned}$$

(2)

in which $X=\left( S_{e},V_{e},S_{a},R_{a},S_{v}\right) \in {\mathbb {R}}^{5}$ and $I=\left( E_{e},E_{a},E_{v},I_{e},I_{a},I_{v}\right) \in {\mathbb {R}}^{6}.$ We note here that X and I represents the classes of the un-infectious and infectious individuals, respectively. Let the DFE from Sect. 3.2 be

$$\begin{aligned} \displaystyle {\mathcal {M}}^{0}=(X_{0},0)=\displaystyle \left( S_{e}^0,V_{e}^0,S_{a}^0,R_{a}^0,S_{v}^0,E_{e}^0,E_{a}^0,E_{v}^0,I_{e}^0,I_{a}^0,I_{v}^0\right) =\left( \dfrac{\Lambda _{e}}{k},\dfrac{\vartheta _e\Lambda _{e}}{\mu _{h}k},\dfrac{\xi \Lambda _{e}}{\mu _{h}k},0,\dfrac{\Lambda _{v}}{\mu _{v}},0,0,0,0,0,0\right) . \end{aligned}$$

For the model to be GAS at ${\mathcal {M}}^{0},$ it needs to satisfy the following conditions from Castillo-Chavez et al. (2002), which are

$({\mathcal {C}}_{1})$:

Local stability is guaranteed at ${\mathcal {M}}^{0}$ whenever $({\mathcal {R}}_{C}^{m}<1.$

${\mathcal {C}}_{1})$:

At $\dfrac{dX}{dt}=F(X_{0},0)$ the DFE is globally asymptotically stable.

$({\mathcal {C}}_{3})$:

${\mathcal {G}}(X,I)={\mathcal {A}}I-\hat{{\mathcal {G}}}(X,I),\,\hat{{\mathcal {G}}}(X,I)\ge 0$ for $(X,I)\in \Omega ,$ where ${\mathcal {A}}={\mathcal {D}}_{I}{\mathcal {G}}({\mathcal {M}}^{0})$ is an Metzler matrix, and $\Omega$ is the model biologically feasible region defined earlier.

Theorem 3.3

If the disease-induced death rate is zero $(d_1=d_2=0)$, then the disease-free equilibrium ${\mathcal {M}}^{0}$ is globally asymptotically (GAS) stable when ${\mathcal {R}}_{C}^{m}<1$.

Proof

To prove that the DFE is GAS when ${\mathcal {R}}_{C}^{m}<1$, we have to verify the conditions ${\mathcal {C}}_{1}$ to ${\mathcal {C}}_{3}.$

Using the approach in van den Driessche and Watmough (2002), we obtain that the DFE ${\mathcal {M}}^{0}$ is LAS when ${\mathcal {R}}_{C}^{m}<1$, so the condition ${\mathcal {C}}_{1}$ is verified.

Next, we re-write the model system (1) in the form given in (2) as

$$\begin{aligned} \dfrac{dX}{dt}=F(X,I)= \begin{pmatrix} \displaystyle \Lambda _{e}+\omega _{e}I_{e}-(\lambda _{e}+k)S_{e}\\ \displaystyle \vartheta _e S_{e} -((1-\varepsilon ) \lambda _e +\mu _h)V_e\\ \displaystyle \xi S_{e}+u\omega _{a}I_{a}+\delta _aR_{a}-(\lambda _{a}+\mu _{h})S_{a}\\ \displaystyle \Lambda _{v}-(\lambda _{v}+\mu _{v})S_{v},\\ \displaystyle (1-u)\omega _{a}I_{a}-k_{4}R_{a}\\ \end{pmatrix}, \,\text {and}\, \dfrac{dI}{dt}=G(X,I)= \begin{pmatrix} \displaystyle \lambda _{e}(S_{e}+(1-\varepsilon ) V_e)-k_{0} E_{e}\\ \displaystyle \lambda _{a}S_{a}-k_{1} E_{a}\\ \displaystyle \lambda _{v}S_{v}-k_{5}E_{v}\\ \displaystyle \sigma _{e}E_{e}-k_{2} I_{e}\\ \displaystyle \sigma _{a}E_{a} -k_{3} I_{a}\\ \displaystyle \sigma _{v}E_{v}-\mu _{v} I_{v} \end{pmatrix}. \end{aligned}$$

(3)

$$\begin{aligned} \dfrac{dX}{dt}=F(X_{0},0) \Leftrightarrow \left\{ \begin{array}{l} {\dot{S}}_e=\Lambda _{e}-kS_{e}, \\ {\dot{V}}_e=\vartheta _e S_{e} -\mu _hV_e, \\ {\dot{S}}_a =\xi S_{e}+\delta _aR_{a}-\mu _{h}S_{a},\\ {\dot{R}}_a=-k_{4}R_{a},\\ {\dot{S}}_v=-\mu _v S_v. \end{array} \right. \end{aligned}$$

(4)

This equation has a unique equilibrium point $\left( \dfrac{\Lambda _{e}}{k},\dfrac{\vartheta _e\Lambda _{e}}{\mu _{h}k},\dfrac{\xi \Lambda _{e}}{\mu _{h}k},0,\dfrac{\Lambda _{v}}{\mu _{v}}\right)$ which is globally asymptotically stable. Therefore, the condition ${\mathcal {C}}_{2}$ is satisfied.

Linearizing the second matrix in equation (3) gives the Metzler Matrix

$$\begin{aligned} {\mathcal {A}}={\mathcal {D}}_{Z}{\mathcal {G}}({\mathcal {M}}^{0})= \begin{pmatrix} \displaystyle -k_{0}&{}0&{}0&{}0&{}0&{}\dfrac{\beta _{e}b_{2}}{N_{h}^0}(S_{e}^0+(1-\varepsilon ) V_e^0)\\ \displaystyle 0&{}-k_{1}&{}0&{}0&{}0&{}\dfrac{\beta _{a}b_{1}}{N_{h}^0}S_a^0\\ \displaystyle 0&{}0&{}-k_{5}&{}\dfrac{\beta _{v}b_{1}}{N_{h}^0}S_v^0&{}\dfrac{\beta _{v}b_{2}}{N_{h}^0}S_v^0&{}0\\ \displaystyle \sigma _{e}&{}0&{}0&{}-k_{2}&{}0&{}0\\ \displaystyle 0&{}\sigma _{a}&{}0&{}0&{}-k_{3}&{}0\\ \displaystyle 0&{}0&{}\sigma _{v}&{}0&{}0&{}-\mu _{v}\\ \end{pmatrix}. \end{aligned}$$

Computing $\hat{{\mathcal {G}}}(X,Z)$ and after some algebraic simplifications, we have

$$\begin{aligned} \hat{{\mathcal {G}}}(X,I)={\mathcal {A}}I- {\mathcal {G}}(X,I)= \begin{pmatrix} \displaystyle \beta _{e}b_{2}I_{v} \left[ \dfrac{S_e^0}{N_h^0}-\dfrac{S_e}{N_h}+(1-\varepsilon ) \left( \dfrac{V_e^0}{N_h^0}-\dfrac{V_e}{N_h} \right) \right] \\ \displaystyle \beta _{a}b_{1}I_{v} \left[ \dfrac{S_a^0}{N_h^0}-\dfrac{S_a}{N_h} \right] \\ \displaystyle \beta _{v}\left( b_{2}I_{e} + b_{1}I_{a}\right) \left[ \dfrac{S_v^0}{N_h^0}-\dfrac{S_v}{N_h} \right] \\ \displaystyle 0\\ \displaystyle 0\\ \displaystyle 0\\ \end{pmatrix}. \end{aligned}$$

Thus,

$$\begin{aligned} \hat{{\mathcal {G}}}(X,I)\ge \begin{pmatrix} \displaystyle \beta _{e}b_{2}I_{v}\left[ S_e^0+(1-\varepsilon )V_e^0\right] \left( \dfrac{1}{N_h^0}-\dfrac{1}{N_h} \right) \\ \\ \displaystyle \beta _{a}b_{1}I_{v}S_a^0\left( \dfrac{1}{N_h^0}-\dfrac{1}{N_h} \right) \\ \displaystyle \beta _{v}\left( b_{2}I_{e} + b_{1}I_{a}\right) S_v^0\left( \dfrac{1}{N_h^0}-\dfrac{1}{N_h} \right) \\ \displaystyle 0\\ \displaystyle 0\\ \displaystyle 0\\ \end{pmatrix}. \end{aligned}$$

Since $\left( \dfrac{1}{N_h^0}-\dfrac{1}{N_h} \right) =\dfrac{N_h-N_h^0}{N_hN_h^0}$, when $d_1=d_2=0$, $N_h(t)-N_h^0=\left( \dfrac{\Lambda _e}{\mu _h}-N_{h}(0)\right) \exp {(-\mu _h t)}$ is positive, and we obtain $\hat{{\mathcal {G}}}(X,I) \ge 0$. The condition ${\mathcal {C}}_3$ is satisfied. We can conclude that if $d_1=d_2=0$, then, the DFE is GAS when ${\mathcal {R}}_{C}^{m}<1$.

$\square$

Remark 3.1

When the disease-induced death rate is not zero, the condition ${\mathcal {R}}_{C}^{m}<1$ is not sufficient for the global stability of disease-free equilibrium and the phenomenon of backward bifurcation may occur. In this case, to mitigate the spread of the disease, it is necessary to reduce ${\mathcal {R}}_{C}^{m}$ to less than another threshold, say ${\mathcal {R}}_{C}^{\#}<{\mathcal {R}}_{C}^{m}<1$.

3.4 Existence of the endemic equilibrium

By setting the left-hand side of the model system (1) to zero and solving it, we obtain the endemic equilibrium point as follows:

$$\begin{aligned} S_e^*=\dfrac{\left( \Lambda _{v}^* b_{2} \beta _{e} (1-\varepsilon ) \lambda _{v}^* \sigma _{v} + N_{h}^* k_{5} \lambda _{v}^* \mu _{h} \mu _{v} + N_{h}^* k_{5} \mu _{h} \mu _{v}^{2}\right) \Lambda _{e} N_{h}^* k_{0} k_{2} k_{5} {\left( \lambda _{v}^* + \mu _{v}\right) } \mu _{v}}{\Lambda _e N_h^*b_2\beta _e\sigma _v \lambda _v^* \left[ k_5\mu _v(\lambda _v^*+\mu _v)((1-\varepsilon )(kk_0k_2-\omega _e\sigma _e\vartheta _e)+\mu _h(k_0k_2-\omega _e\sigma _e)) \right] +k_6}, \end{aligned}$$

(5)

where   $k_6=(\Lambda _vb_2\beta _e \sigma _v \lambda _v)^2(1-\varepsilon )(k_0k_2-\omega _e\sigma _e)+N_h^*kk_0k_2k_5\mu _h\mu _v^2(\lambda _v+\mu _v)^2$, and $N_h^*=\dfrac{\Lambda _e}{\mu _h}.$

We also have

$$\begin{aligned} \begin{array}{cccc} S_v^*=\dfrac{\Lambda _v}{\lambda _v+\mu _v}, &{} E_v^*=\dfrac{\lambda _v S_v^*}{k_5}, &{} I_v^*=\dfrac{\sigma _v E_v^*}{\mu _v}, &{} \\ \\ \lambda _e^*=\dfrac{b_2 \beta _e I_v^*}{N_h^*}, &{} V_e^*=\dfrac{\vartheta S_e^*}{(1-\varepsilon )\lambda _e^*+mu_h)} ,&{} E_e^*=\dfrac{\lambda _e^*(S_e^*+(1-\varepsilon )V_e^*)}{k_0}, &{} I_e=\dfrac{\sigma _e E_e^*}{k_2} \\ \end{array}, \end{aligned}$$

(6)

$$\begin{aligned} S_a^*=\dfrac{\Lambda _eS_e^*k_1k_3k_4k_5\xi (\lambda _v^*+\mu _v)}{\mu _v\left[ \Lambda _ek_1k_3k_4k_5\mu _v(\lambda _v^*+\mu _v)+\Lambda _vb_1\beta _a\sigma _v\lambda _v^*(k_1k_3k_4-\delta _a\omega _a\sigma _a-u\omega _a\sigma _a\mu _h) \right] }, \end{aligned}$$

(7)

$$\begin{aligned} \begin{array}{cccc} \lambda _a^*=\dfrac{b_1 \beta _a I_v^*}{N_h^*},&E_a^*=\dfrac{\lambda _a^* S_a^*}{k_1},&I_a^*=\dfrac{\sigma _a E_a^*}{k_3},&R_a^* =\dfrac{(1-u)\omega _a I_a^*}{k_4}. \end{array} \end{aligned}$$

(8)

Substituting the expressions for $I_a^*$ and $I_e^*$ into the expression for $\lambda _v,$ we obtain the polynomial

$$\begin{aligned} a_3\lambda _v^{*3}+a_2\lambda _v^{*2}+a_1\lambda _v^*+a_0, \end{aligned}$$

(9)

where the expressions of $a_0, a_1, a_2$ and $a_3$ are given in the Appendix.

Because the model monitors human populations, all associated state variables should be non-negative for all time $t \ge 0$. We therefore use Descarte’s rule of signs to discuss the existence of possible positive roots of Eq. (9). The results are summarized in Table 3.

Table 3 Number of possible positive roots of equation (9) using Descarte’s rule of signs

The following remark summarizes the results in Table 3

Remark 3.2

• If all the coefficients of the polynomial (9) have the same signs, then the model system (1) has no endemic equilibrium point,

• If $\mathcal{R}_C^m \ne 1,$ then the system has zero, one, two or three endemic equilibrium points,

• If $\mathcal{R}_C^m=1,$ then the system has either zero, one or two endemic equilibrium points.

The last two conditions above imply that the model system (1) could exhibit the phenomenon of backward/subcritical bifurcation, when a stable DFE co-exists with a stable endemic equilibrium. This is an epidemiological situation in which the classical requirement of having the basic reproduction number less than unity although necessary is not sufficient to eliminate the disease.

4 Optimal control model

The application of optimal control enables us to forecast or choose the best scenario that if well implemented could help mitigate the spread of the disease. Thus, to investigate the potential impact of the implemented intervention measures, the following control variables are incorporated into the model system (1):

• $c_1 (t)$ representing the use of personal protection measures to prevent mosquitoes bites during the day and the night such as the use of insecticide-treated nets, application of repellents to skin or spraying of insecticides,

• $c_2(t)$ representing the treatment, and

• $c_3(t)$ representing the use of vaccination to prevent malaria,

as follows

$$\begin{aligned} \left\{ \begin{array}{lll} S_{e}'=\Lambda _{e}+c_2\omega _{e}I_{e}-((1-c_1)\lambda _{e}+c_3\vartheta _e+\xi +\mu _{h})S_{e},\\ \\ V_{e}'=c_3 \vartheta _e S_{e} -((1-c_1)(1-\varepsilon ) \lambda _e +\mu _h)V_e,\\ \\ S_{a}'=\xi S_{e}+(1-p)c_2\omega _{a}I_{a}+\delta _a R_{a}-((1-c_1)\lambda _{a}+\mu _{h})S_{a},\\ \\ E_{e}'=(1-c_1)\lambda _{e}(S_{e}+(1-\varepsilon ) V_e)-(\sigma _{e}+\mu _{h}) E_{e},\\ \\ E_{a}'= (1-c_1)\lambda _{a}S_{a}-(\sigma _{a}+\mu _{h}) E_{a},\\ \\ I_{e}'=\sigma _{e}E_{e}-(c_2\omega _{e}+\mu _{h}+d_1) I_{e},\\ \\ I_{a}'=\sigma _{a}E_{a} -(c_2 \omega _{a}+\mu _{h}+d_2) I_{a},\\ \\ R_{a}'=pc_2\omega _{a}I_{a}-(\mu _{h}+\delta _a)R_{a},\\ \\ S_{v}'=\Lambda _v-((1-c_1)\lambda _{v}+\mu _{v})S_{v},\\ \\ E_{v}'=(1-c_1)\lambda _{v}S_{v}-(\sigma _{v}+\mu _{v})E_{v},\\ \\ I_{v}'=\sigma _{v}E_{v}-\mu _{v} I_{v}.\\ \end{array} \right. \end{aligned}$$

(10)

Consider the following quadratic objective functional which measures the cost of the control. This cost includes the above interventions. The the nonlinear objective function is

$$\begin{aligned} J\left( c_1, c_2, c_3 \right) = \int _{0}^{T} \left[ A_1 E_{e}(t) + A_2 E_{a}(t) +A_3 I_{e}(t) +A_4 I_(t) + A_5 N_{v}(t) + \frac{w_1}{2}c_1^2 + \frac{w_2}{2}c_2^2+ \frac{w_3}{2}c_3^2 \right] dt, \end{aligned}$$

(11)

where T is the final time, $A_i, \ i=1,\cdots ,5$ are positive weight constants, $N_v = S_v+ + E_v + I_v$ and $w_i, \ i=1,\cdots ,3$ are weight constants for the strategies and treatments against proliferation of malaria. The fact that the controls are linearly in (10) and quadratic in the objective functional allows the Hamiltonian associated to the optimal control problem to be maximized. Therefore, we seek to find, using the Pontryagin Maximum Principle (Pontryagin et al. 1986), an optimal control $(c_1^*, c_2^*, c_3^*) \in U$ satisfying (10), such that

$$\begin{aligned} J\left( c_1^*, c_2^*, c_3^* \right) = \min \left\{ J\left( c_1, c_2, c_3\right) \ |\ (c_1, c_2, c_3) \in U\right\} . \end{aligned}$$

(12)

The associated pseudo-Hamiltonian is

$$\begin{aligned} \displaystyle H= & {} \,\displaystyle A_1 E_{e}(t) + A_2 E_{a}(t) +A_3 I_{e}(t) +A_4 I_{a}(t) + A_5 N_{v}(t) + \frac{w_1}{2}c_1^2 + \frac{w_2}{2}c_2^2+ \frac{w_3}{2}c_3^2 \nonumber \\{} & {} + \displaystyle \xi _1\left[ \Lambda _{e}+c_2\omega _{e}I_{e}-((1-c_1)\lambda _{e}+c_3\vartheta _e+\xi +\mu _{h})S_{e}\right] \nonumber \\{} & {} +\displaystyle \xi _2 \left[ c_3 \vartheta _e S_{e} -((1-c_1)(1-\varepsilon ) \lambda _e +\mu _h)V_e \right) \nonumber \\{} & {} +\displaystyle \xi _3\left[ \xi S_{e}+(1-u)c_2\omega _{a}I_{a}+\delta _a R_{a}-((1-c_1)\lambda _{a}+\mu _{h})S_{a}, \right] \nonumber \\{} & {} +\displaystyle \xi _4 \left[ (1-c_1)\lambda _{e}(S_{e}+(1-\varepsilon ) V_e)-(\sigma _{e}+\mu _{h}) E_{e} \right] \nonumber \\{} & {} +\displaystyle \xi _5 \left[ (1-c_1)\lambda _{a}S_{a}-(\sigma _{a}+\mu _{h}) E_{a} \right] \nonumber \\ \nonumber \\{} & {} +\displaystyle \xi _6 \left[ \sigma _{e}E_{e}-(c_2\omega _{e}+\mu _{h}+d_1) I_{e} \right] \nonumber \\ \nonumber \\{} & {} +\displaystyle \xi _7 \left[ \sigma _{a}E_{a} -(c_2 \omega _{a}+\mu _{h}+d_2) I_{a} \right] \nonumber \\ \nonumber \\{} & {} +\displaystyle \xi _8 \left[ uc_2\omega _{a}I_{a}-(\mu _{h}+\delta _a)R_{a} \right] \nonumber \\ \nonumber \\{} & {} +\displaystyle \xi _9 \left[ \Lambda _v-((1-c_1)\lambda _{v}+\mu _{v})S_{v} \right] \nonumber \\ \nonumber \\{} & {} +\displaystyle \xi _{10} \left[ (1-c_1)\lambda _{v}S_{v}-(\sigma _{v}+\mu _{v})E_{v}\right] \nonumber \\ \nonumber \\{} & {} +\displaystyle \xi _{11} \left[ \sigma _{v}E_{v}-\mu _{v} I_{v} \right] , \end{aligned}$$

(13)

where the $\, \xi _1 ,\, \xi _2 ,\,\xi _3 ,\, \xi _4, \, \xi _5 ,\, \xi _6,\,\xi _7 ,\, \xi _8 ,\,\xi _9 ,\, \xi _{10}, \, \xi _{11}$ are the associated adjoints for the states $S_{e},V_{e},S_{a},E_{e},E_{a},I_{e},I_{a},R_{a}, S_{v},E_{v},I_{v}$. The system of equations is found by taking the appropriate partial derivatives of the Hamiltonian (13) with respect to the associated state variable.

Where $\xi _i, \ i=1,\cdots ,11$ are the adjoint variables satisfying

$$\begin{aligned} \begin{aligned} \displaystyle \xi ^{'}_1&= - \displaystyle \frac{\partial {\mathbb {H}}}{\partial S_{e}}\qquad \qquad&\xi ^{'}_2 = - \displaystyle \frac{\partial {\mathbb {H}}}{\partial V_{e}}, \\ \displaystyle \xi ^{'}_3&= - \displaystyle \frac{\partial {\mathbb {H}}}{\partial S_{a}}\qquad \qquad&\xi ^{'}_4 = - \displaystyle \frac{\partial {\mathbb {H}}}{\partial E_{e}}, \\ \displaystyle \xi ^{'}_5&= - \displaystyle \frac{\partial {\mathbb {H}}}{\partial E_{a}}\qquad \qquad&\xi ^{'}_6 = - \displaystyle \frac{\partial {\mathbb {H}}}{\partial I_{e}}, \\ \displaystyle \xi ^{'}_7&= - \displaystyle \frac{\partial {\mathbb {H}}}{\partial I_{a}}\qquad \qquad&\xi ^{'}_8 = - \displaystyle \frac{\partial {\mathbb {H}}}{\partial R_{a}}, \\ \displaystyle \xi ^{'}_9&= - \displaystyle \frac{\partial {\mathbb {H}}}{\partial S{v}}\qquad \qquad&\xi ^{'}_{10} = - \displaystyle \frac{\partial {\mathbb {H}}}{\partial E_v}.\\ \displaystyle \xi ^{'}_7&= - \displaystyle \frac{\partial {\mathbb {H}}}{\partial I_{v}}\qquad \qquad .&\end{aligned} \end{aligned}$$

(14)

That is,

$$\begin{aligned} \begin{aligned} \displaystyle \xi ^{'}_1&= \displaystyle \lambda _{e}(1-c_1)(\xi _{1} -\xi _{4}) + c_3\vartheta _e(\xi _{1} -\xi _{2}) +\xi (\xi _{1} -\xi _{3})+\dfrac{\lambda _{e}S_{e}}{N_{h}}(1-c_1)(\xi _{4} -\xi _{1})+ \dfrac{\lambda _{a}S_{a}}{N_{h}}(1-c_1)(\xi _{5} -\xi _{3})\\&\quad + \dfrac{\lambda _{v}S_{v}}{N_{h}}(1-c_1)(\xi _{10} -\xi _{9})+ \dfrac{(1-\varepsilon )\lambda _{e}V_{e}}{N_{h}}(1-c_1)(\xi _{4} -\xi _{2}) +\mu _{h}\xi _{1},\\ \displaystyle \xi ^{'}_2&= \displaystyle \dfrac{\lambda _{e}S_{e}}{N_{h}}(1-c_1)(\xi _{4} -\xi _{1}) + \dfrac{\lambda _{e}V_e}{N_{h}}(1-c_1)(1-\varepsilon )(\xi _{4} -\xi _{2}) + \lambda _{e}(1-c_1)(1-\varepsilon ) (\xi _{2}-\xi _{4})\\&\quad + \dfrac{\lambda _{v}S_{v}}{N_{h}}(1-c_1)(\xi _{10} -\xi {_9})+ \dfrac{\lambda _{a}S_{a}}{N_{h}}(1-c_1)(\xi _{5} -\xi _{3}) +\mu _{h}\xi _{2}, \\ \displaystyle \xi ^{'}_3&= \displaystyle \dfrac{\lambda _{e}S_{e}}{N_{h}}(1-c_1)(\xi _{4} -\xi _{1}) + \dfrac{\lambda _{e}V_e}{N_{h}}(1-c_1)(1-\varepsilon )(\xi _{4} -\xi _{2}) + \dfrac{\lambda _{a}S_a}{N_{h}}(1-c_1)(\xi _5- \xi _{3}) +(1-c_1)\lambda _{a}(\xi _3-\xi _5) \\&\quad + \dfrac{\lambda _{v}S_{v}}{N_{h}}(1-c_1)(\xi _{10} -\xi _{9}) +\mu _{h}\xi _{3}, \\ \displaystyle \xi ^{'}_4&= \displaystyle \dfrac{\lambda _{e}S_{e}}{N_{h}}(1-c_1)(\xi _{4} -\xi _{1}) + \dfrac{\lambda _{e}V_e}{N_{h}}(1-c_1)(1-\varepsilon )(\xi _{4} -\xi _{2}) + \dfrac{\lambda _{a}S_{a}}{N_{h}}(1-c_1) (\xi _{5}-\xi _{3}) +\sigma _{e}(\xi _{4}-\xi _{6})\\&\quad + \dfrac{\lambda _{v}S_{v}}{N_{h}}(1-c_1)(\xi _{10} -\xi _{9}) +\mu _{h}\xi _{4}-A_1, \\ \displaystyle \xi ^{'}_5&= \displaystyle \dfrac{\lambda _{e}S_{e}}{N_{h}}(1-c_1)(\xi _{4} -\xi _{1}) + \dfrac{\lambda _{e}V_e}{N_{h}}(1-c_1)(1-\varepsilon )(\xi _{4} -\xi _{2}) + \dfrac{\lambda _{a}S_a}{N_{h}}(1-c_1) (\xi _{5}-\xi _{3}) +\sigma _{a}(\xi _{5}-\xi _{7})\\&\quad + \dfrac{\lambda _{v}S_{v}}{N_{h}}(1-c_1)(\xi _{10} -\xi _{9}) +\mu _{h}\xi _{5}-A_2, \\ \displaystyle \xi ^{'}_6&= \displaystyle \dfrac{\lambda _{e}S_{e}}{N_{h}}(1-c_1)(\xi _{4} -\xi _{1}) + \dfrac{\lambda _{e}V_e}{N_{h}}(1-c_1)(1-\varepsilon )(\xi _{4} -\xi _{2}) + \dfrac{\lambda _{a}S_a}{N_{h}}(1-c_1) (\xi _{5}-\xi _{3}) + c_2\omega _{e}(\xi _{6}-\xi _{1})\\&\dfrac{ \lambda _{v}S_v}{N_{h}}(1-c_1)(\xi _{10} -\xi _{9})+ \dfrac{b_2 \beta _{v}S_v}{N_{h}}(1-c_1)(\xi _{9} -\xi _{10}) +(\mu _{h}+d_{1})\xi _{6}-A_3, \\ \displaystyle \xi ^{'}_7&= \displaystyle \dfrac{\lambda _{e}S_{e}}{N_{h}}(1-c_1)(\xi _{4} -\xi _{1}) + \dfrac{\lambda _{e}V_{e}}{N_{h}}(1-c_1)(1-\varepsilon ) (\xi _{4} -\xi _{2}) + \dfrac{\lambda _{a}S_{a}}{N_{h}}(1-c_1)(\xi _{5}-\xi _{3}) +\dfrac{b_1 \beta _v S_v}{N_{h}}(1-c_1)(\xi _9-\xi _{10})\\&\quad -c_2\omega _{a}(u\xi _{8}+(1-u)\xi _{3})+ \dfrac{\lambda _{v}S_v}{N_{h}}(1-c_1)(\xi _{10} -\xi _{9}) +(c_2\omega _a+\mu _{h}+d_{2})\xi _{7}-A_4, \\ \displaystyle \xi ^{'}_8&= \displaystyle \dfrac{\lambda _{e}S_{e}}{N_{h}}(1-c_1)(\xi _{4} -\xi _{1}) + \dfrac{\lambda _{e}V_{e}}{N_{h}}(1-c_1)(1-\varepsilon ) (\xi _{4} -\xi _{2}) + \dfrac{\lambda _{a}S_{a}}{N_{h}}(1-c_1)(\xi _{5}-\xi _{3}) + \delta _{a}(\xi _{8}-\xi _{3})\\&\quad + \dfrac{\lambda _v S_{v}}{N_{h}}(1-c_1)(\xi _{10} -\xi _{9}) + \mu _{h}\xi _{8}, \\ \displaystyle \xi ^{'}_9&= \displaystyle \lambda _v (1-c_1) (\xi _{9} -\xi _{10}) + \mu _{v}\xi _{9}- A_5, \\ \displaystyle \xi ^{'}_{10}&= \displaystyle \sigma _{v}(\xi _{10} -\xi _{11}) + \mu _{v}\xi _{10}- A_5, \\ \displaystyle \xi ^{'}_{11}&= \displaystyle \lambda _{e}S_{e}(1-c_1)(\xi _{1} -\xi _{4}) +\lambda _{e}V_{e}(1-\varepsilon )(1-c_1)(\xi _{2} -\xi _{4})+\lambda _{a}S_{a}(1-c_1)(\xi _{3} -\xi _{5}) + \mu _{v}\xi _{11}-A_5, \\ \end{aligned} \end{aligned}$$

Considering the optimality conditions, the Hamiltonian function is differentiated with respect to the control variables resulting in

$$\begin{aligned} \begin{aligned} 0=\displaystyle \frac{\partial H}{\partial c_1}&=\displaystyle w_{1}c_1 +(\xi _{1}-\xi _{4})\lambda _{2}S_{e} +(\xi _{2}-\xi _{4})\varepsilon \lambda _{3} V_{e}+\xi _{4}\lambda _{e}V_{e}+(\xi _{3}-\xi _{5})\lambda _{a}V_{a}+(\xi _{9}-\xi _{10})\lambda _{v} S_{v} \\ 0=\displaystyle \frac{\partial H}{\partial c_2}&=\displaystyle w_{2}c_2+(\xi _{1}-\xi _{6})\omega _{e}I_{e}+(\xi _{8}-\xi _{3})p\omega _{a}I_{a} +(\xi _{3}-\xi _{7})\omega _{a}I_{a}. \\ 0=\displaystyle \frac{\partial H}{\partial c_3}&=\displaystyle w_{3}c_3+(\xi _{2}-\xi _{1})\vartheta _{e}S_{e}, \end{aligned} \end{aligned}$$

(15)

on the interior of the control set ${\mathcal {U}}$. Then, solving for $c_1^*$ (on the interior of the control set) gives

$$\begin{aligned} c_1^*= & {} \,\displaystyle \dfrac{\lambda _{2}S_{e} (\xi _{4}-\xi _{1}) +(1-\varepsilon ) \lambda _{e} V_{e}(\xi _{4}-\xi _{2})+\lambda _{a}S_{a}(\xi _{5}-\xi _{3})+\lambda _{v} S_{v} (\xi _{10}-\xi _{9})}{\displaystyle w_1}, \nonumber \\ c_2^*= & {} \, \displaystyle \dfrac{\displaystyle \omega _{e}I_{e}(\xi _{6}-\xi _{1})+\omega _{a}I_{a}(\xi _{7}-\xi _{3})+u\omega _{a}I_{a}(\xi _{3}-\xi _{8}) }{\displaystyle w_2 },\nonumber \\ c_3^*= & {} \, \displaystyle \dfrac{\displaystyle \vartheta _{e}S_{e}(\xi _{1}-\xi _{2}) }{\displaystyle w_3 } \end{aligned}$$

(16)

Using the bounds on the controls, we obtain the characterization, and hence

$$\begin{aligned} c_1^*= & {} min\left\{ b,max\left[ a,\,\displaystyle -\frac{ (\xi _{1}-\xi _{4})\lambda _{2}S_{e} +(\xi _{2}-\xi _{4})\varepsilon \lambda _{3} V_{e} +(\xi _{3}-\xi _{5})\lambda _{a}V_{a}+(\xi _{9}-\xi _{10})\lambda _{v} S_{v} }{\displaystyle w_1} \right\} ,\right. \\ c_2^*= & {} min\left\{ d,max\left[ c,\,\displaystyle -\frac{\displaystyle (\xi _{1}-\xi _{6})\omega _{e}I_{e}+(\xi _{8}-\xi _{3})p\omega _{a}I_{a} +(\xi _{3}-\xi _{7})\omega _{a}I_{a} }{\displaystyle w_2 } \right\} \right. \\ c_3^*= & {} min\left\{ d,max\left[ c,\,\displaystyle -\frac{\displaystyle (\xi _{2}-\xi _{1})\vartheta _{e}S_{e} }{\displaystyle w_3 } \right\} .\right. \end{aligned}$$

5 Numerical simulations

Graphical representations using the model parameter values in Table 4 are illustrated below. When $c_1(t) = 0$ there is no vaccination at all, we set the lower bound of the controls to 0 and the upper bound to 1, that is, $a = c = 0$, $b = d = 1$. Thus, $0 \le c_1(t), c_2(t) \le 1$. The model parameter values in Table Table 4 are taken from the literature, or assumed for illustrative purpose.

Table 4 Model parameter values

5.1 Global sensitivity analysis of $\mathcal{R}_C^m$ and $V_{e}$

This subsection is devoted to the global sensitivity analyses of the model reproduction number $\mathcal{R}_C^m$ and the vaccination compartment $V_{e}$. The threshold $\mathcal{R}_C^m$ is chosen because of its crucial role in forecasting the spread/persistence of an epidemic, while the vaccination class is chosen considering ’prevention is better than cure’. In general, mathematical models possess some uncertainties due to variations such as demography and geographical location incorporated during the model formulation. Due to these uncertainties, we have an inherent epistemic uncertainty in our model parameterization for those estimated or calculated (Marino et al. 2008). For this reason, we investigate these uncertainties in the model parameters by applying the method of Latin-Hyper-Cube Sampling (LHS) with a combination of partial rank correlation coefficient (PRCC) to bypass unbiased estimates of the parameters (Marino et al. 2008; Bauer et al. 2008). This method takes an input variable and generates an output containing the tornado and/or scatter plots (Bauer et al. 2008). More details on this can be found in Herdicho et al. (2021), Chukwu and Nyabadza (2020) and the references therein. Following Bauer et al. (2008), parameters with PRCC value output less than $-0.5$ or greater than 0.5 are assumed to be most sensitive while below these ranges are less significant.  The global sensitivity analysis is carried out using Matlab and R softwares.

Figure 3 is the tornado plot for the PRCC of the model parameters from the reproduction number $\mathcal{R}_C^m$, while Fig.  4 depicts the scatter plots obtained against the vaccination compartment.

Fig. 3

Tonardo plot showing all the model parameters against $\mathcal{R}_C^m.$ The longer the bar, the more sensitive if the corresponding parameter

Table 5 PRCC values and p-values with their significant impact

Fig. 4

Scatter plots simulation showing the PRCC values of parameters $\lambda _{a}$ $\vartheta _{e},\xi$ and $\mu _{h}$ against the vaccination class $V_{e}$

From Fig. 3, the parameter $b_{2}$, $\beta _{e}$ and $\beta _{v}$ are strongly positively correlated with positive PRCC, while the parameter $\mu _{v}$ has a strong negative PRCC value. The biological implication of the positive PRCC values implies their increase will certainly increase the numerical value of the $\mathcal{R}_C^m$, and the contrast is true for the negative PRCC values. Last, with a 1000 sample size and a unit 1, we performed a global sensitivity using the $V_{e}$ class as our response variable. The results are presented using a scatter plot for visualization of the sensitivity of the model parameters. It can be seen that $\lambda _{a}$ has a strong positive correlation, while $\vartheta _{e}, \xi$ and $\mu _{h}$ have a strong negative PRCC values. These suggest a need to decrease the transmission parameters and increase vaccination of children less than 5 years old, a strategy that could help contain the spread of malaria in the human population. Importantly, Table 5 gives the PRCC values of each model parameter in the $\mathcal{R}_C^m$ with the associated p-values showing the level statistically significance. Those parameters with p-values less than 0.05 are said to be significant and have more effect on the reproduction number when compared to other parameters.

5.2 Effect of implementing control measures

5.2.1 Personal protection $c_1\ne 0$, treatment $c_2\ne 0$, and vaccination $c_3\ne 0$

Figure 5 depicts the time series of the scenario when the three control measures, namely personal protection $c_1\ne 0$,  treatmentt $c_2\ne 0$, and vaccination of children under 5 years old $c_3\ne 0$ are implemented.

Fig. 5

Effect of implementing controls of the model state variables for (a) Susceptible humans under 5-years. (b) Susceptible humans over 5-years. (c) Infectious humans under 5-years. (d) Infectious humans over 5-years. (e) Vaccinated humans under 5-years

Results of the strategy that combines the three intervention measures indicate that the number of susceptible humans under over 5-years of age decreases, whereas the number of susceptible individuals over 5-years increase in time, as seen in Fig. 5a, b respectively. On the other hand, in the presence of control measures, we have a reduced number of infectious individuals in both of these sub-groups of less and greater than five as shown in Fig. 5c, d. This very interesting result suggests a need for simultaneous implementation of treatment, vaccination of less than 5-years old, and personal protective measures to help mitigate the spread of malaria outbreak(s). In addition, Fig. 5e shows that the use of vaccination does have a significant positive impact on individuals under the age of five, this is very important as this group bears the highest burden of malaria morbidity and mortality in the community. This results further asserts the WHO recommendation to vaccinate children under five (WHO 2022c), and thanks to the availability of the vaccine currently for this age group.

Fig. 6

Control profiles for the use of optimal control variables (a) $c_{1}$-personal protective measures. (b) $c_{2}$-Treatment effort (c) $c_{3}$-vaccination effort

Figure 6 shows the control profiles for three types of controls considered in this paper. It can be seen in Fig. 6a, b that it is important to keep the use of personal protective measures and the treatment effort at its maximum level throughout the modeling time to achieve the control of malaria. In contrast, vaccination is effective at the start of the simulation for about 200 days, but start to reduce, Fig. 6c. This may not be surprising as mass vaccination is expected because the vaccine has just been available for the first time recently, but it is expected that after a while, vaccination rate will decrease. This results can be ascertain as well based on the Covid19 vaccination campaign in the early 2021 which has now drastically fade down.

5.2.2 Personal protection $c_1\ne 0$, treatment $c_2 = 0$, and vaccination $c_3 = 0$

The implementation of personal protection alone has a positive impact on reducing the malaria spread in the community as depicted in Fig. 7. The control profile of $c_1$ is at its maximum value through the intervention period, Fig. 7f.

Fig. 7

Effect of implementing only control $c_1$ of the model state variables for (a) Susceptible humans under 5-years. (b) Susceptible humans over 5-years. (c) Infectious humans under 5-years. (d) Infectious humans over 5-years. (e) Vaccinated humans under 5-years

5.2.3 Personal protection $c_1 = 0$, treatment $c_2 \ne 0$, and vaccination $c_3 = 0$

The implementation of treatment as a sole control measure is depicted in Fig. 8. As in the case of singly implementing personal protection only, though treatment alone has a positive population level impact, it is not sufficient to eradicate the disease. The control profile of $c_2$ is at its maximum value through the intervention period, Fig. 8f.

Fig. 8

Effect of implementing only control $c_2$ of the model state variables for (a) Susceptible humans under 5-years. (b) Susceptible humans over 5-years. (c) Infectious humans under 5-years. (d) Infectious humans over 5-years. (e) Vaccinated humans under 5-years

5.2.4 Personal protection $c_1 = 0$, treatment $c_2 = 0$, and vaccination $c_3 \ne 0$

Vaccination of children under five also has a positive population level impact, see Fig. 9. But, in this case, the vaccination should be kept at its maximum for about 600 days, which is about 3 times more than when the three intervention strategies are concurrently implemented, Fig. 9f.

Fig. 9

Effect of implementing only control $c_3$ of the model state variables for (a) Susceptible humans under 5-years. (b) Susceptible humans over 5-years. (c) Infectious humans under 5-years. (d) Infectious humans over 5-years. (e) Vaccinated humans under 5-years

When treatment is the only implemented control measure, it takes more than 3 years of treatment of infected individuals at the maximum control level as depicted in Fig. 8f to have a meaningful population-level impact. On the other hand, in Fig. 9f, when vaccination alone is implemented for children less than 5 years old, it takes about a year of optimally administering vaccination for a meaningful population-level impact.

6 Conclusion

Malaria is an infectious vector-borne disease of global public health concern, with children less than 5 years old bearing the highest burden. While several prevention and therapeutic measures have been implemented to fight against malaria, the recent groundbreaking RTS malaria vaccine (trade name Mosquirix), a recombinant protein-based malaria vaccine for children age five and below, recommended by the WHO could be a game changer.

A two-group malaria model structured by age with individuals above and below five years of age is formulated and analyzed. The basic reproduction number defined as the expected number of secondary infections generated by one infected individual during its entire period of infectiousness in a naive population is computed using the next generation method. The disease-free equilibrium is shown to be globally asymptotically stable when the disease-induced death rate in both human groups is zero. Descarte’s rule of signs is used to discuss the possible existence of multiple endemic equilibria, in which case the model undergoes a backward or subcritical bifurcation. The epidemiological implication of this situation when a stable DFE co-exists with a stable endemic equilibrium is that having the basic reproduction number less than unity although necessary is not sufficient to eliminate the disease.

By construction, mathematical models inherit the loss of information that could make model outcomes imprecise. Therefore, global sensitivity analysis of the basic reproduction number and the vaccination class using Latin-Hyper Cube Sampling (LHS) in combination with partial rank correlation coefficient (PRCC) are investigated and results depicted graphically. The most sensitive parameters are related to children under five years old such as the average biting rate of mosquitoes on susceptible individuals under 5 years old $b_2$, the progression rate from exposed to infectious individuals under 5 years old $\sigma _e$, and the probability of infection of susceptible under 5 years old per mosquito bite $\beta _e$.

To mitigate the spread of infectious diseases, intervention measures (both therapeutic and non-therapeutic) are paramount. Consequently, we extended the proposed model by incorporating three time-dependent controls, namely personal protection, treatment, and vaccination of children under-five. Using Pontryagin’s maximum principle, we prove the existence of an optimal control problem and find the optimal control combination. Numerical simulations reveal that concurrently applying the three intervention measures is the best scenario for fighting against malaria epidemic in a community, compared to using either single or any dual combination of intervention(s) at a time. While singly or dual intervention measure implemented at a time still has a positive population effect, they are all less effective than concurrently implementing the triple intervention strategy at a time.

The proposed model has some limitations. We assumed that individuals in the human population mix homogeneously. The model could be extended by developing an agent-based two-group malarial model. Also, since the new malaria vaccine for children less than five years old is not yet widely available, a stochastic version of the model is another avenue that warrant further investigation. As countries data become available, fitting the model to real data could enable better estimation of some model parameters, which herein are mainly extracted from existing literature. In general, the density of mosquitoes fluctuates between climatic seasons. For this reason, accounting for seasonal factor (in the birth rate of mosquitoes) as well as the influence of climate (temperature-dependent model), two important drivers of the malaria dynamics are important.

Appendix

The coefficients $a_0, a_1, a_2$ and $a_3$ of the polynomial 9 are given below.

$$\begin{aligned} a_3= & {} (1-\varepsilon )(\Lambda _v b_2 \beta _e\sigma _v)^2(k_0k_2-\omega _e\sigma _e)+\Lambda _e\Lambda _vb_2\beta _ek_5\mu _v\sigma _v(1-\varepsilon )(kk_0k_2-\omega _e\sigma _e\vartheta _e)\nonumber \\{} & {} +\Lambda _e\Lambda _vb_2\beta _ek_0k_2k_5\mu _h\mu _v\sigma _v\left[ \Lambda _vb_1\beta _a\sigma _v(k_4\omega _a\sigma _au+\delta _a\omega _a(1-u)-k_1k_3k_4)-\Lambda _ek_1k_3k_4k_5\mu _v \right] \nonumber \\{} & {} -\Lambda _e\Lambda _vb_2\beta _ek_5\mu _h\mu _v\omega _e\sigma _e\sigma _v+\Lambda _e^2kk_0k_2k_5^2\mu _v^2 \nonumber \\ a_2= & {} -\Lambda _{v}^{3} b_{1} b_{2}^{3} \beta _{a} \beta _{e}^{2} \beta _{v} (1-\varepsilon ) k_{4} \mu _{h}^{2} \omega _{a} \sigma _{a} \sigma _{e} \sigma _{v}^{3} u - \Lambda _{v}^{3} b_{1} b_{2}^{3} \beta _{a} \beta _{e}^{2} \beta _{v} \delta _{a} (1-\varepsilon ) \mu _{h}^{2} \omega _{a} \sigma _{a} \sigma _{e} \sigma _{v}^{3} (1-u) \nonumber \\{} & {} + \Lambda _{v}^{3} b_{1} b_{2}^{3} \beta _{a} \beta _{e}^{2} \beta _{v} (1-\varepsilon ) k_{1} k_{3} k_{4} \mu _{h}^{2} \sigma _{e} \sigma _{v}^{3}\nonumber \\{} & {} - \Lambda _{e} \Lambda _{v}^{2} b_{1} b_{2}^{2} \beta _{a} \beta _{e} \beta _{v} (1-\varepsilon ) k_{4} k_{5} \mu _{h} \mu _{v} \omega _{a} \sigma _{a} \sigma _{e} \sigma _{v}^{2} u \vartheta _e - \Lambda _{e} \Lambda _{v}^{2} b_{1} b_{2}^{2} \beta _{a} \beta _{e} \beta _{v} \delta _{a} (1-\varepsilon ) k_{5} \mu _{h} \mu _{v} \omega _{a} \sigma _{a} \sigma _{e} \sigma _{v}^{2} (1-u) \vartheta _e\nonumber \\{} & {} + \Lambda _{e} \Lambda _{v}^{2} b_{2}^{3} \beta _{e}^{2} \beta _{v} (1-\varepsilon ) k_{1} k_{3} k_{4} k_{5} \mu _{h}^{2} \mu _{v} \sigma _{e} \sigma _{v}^{2} + \Lambda _{e} \Lambda _{v}^{2} b_{1} b_{2} \beta _{a} \beta _{e} (1-\varepsilon ) k k_{0} k_{2} k_{4} k_{5} \mu _{v}^{2} \omega _{a} \sigma _{a} \sigma _{v}^{2} u\nonumber \\{} & {} - \Lambda _{e} \Lambda _{v}^{2} b_{1} b_{2}^{2} \beta _{a} \beta _{e} \beta _{v} k_{4} k_{5} \mu _{h}^{2} \mu _{v} \omega _{a} \sigma _{a} \sigma _{e} \sigma _{v}^{2} u + \Lambda _{e} \Lambda _{v}^{2} b_{1} b_{2} \beta _{a} \beta _{e} \delta _{a} (1-\varepsilon ) k k_{0} k_{2} k_{5} \mu _{v}^{2} \omega _{a} \sigma _{a} \sigma _{v}^{2} (1-u)\nonumber \\{} & {} - \Lambda _{e} \Lambda _{v}^{2} b_{1} b_{2}^{2} \beta _{a} \beta _{e} \beta _{v} \delta _{a} k_{5} \mu _{h}^{2} \mu _{v} \omega _{a} \sigma _{a} \sigma _{e} \sigma _{v}^{2} (1-u) + \Lambda _{e} \Lambda _{v}^{2} b_{1} b_{2}^{2} \beta _{a} \beta _{e} \beta _{v} (1-\varepsilon ) k_{1} k_{3} k_{4} k_{5} \mu _{h} \mu _{v} \sigma _{e} \sigma _{v}^{2} v \nonumber \\{} & {} - \Lambda _{e} \Lambda _{v}^{2} b_{1} b_{2} \beta _{a} \beta _{e} (1-\varepsilon ) k_{4} k_{5} \mu _{v}^{2} \omega _{a} \omega _{e} \sigma _{a} \sigma _{e} \sigma _{v}^{2} u \vartheta _e - \Lambda _{e} \Lambda _{v}^{2} b_{1} b_{2} \beta _{a} \beta _{e} \delta _{a} (1-\varepsilon ) k_{5} \mu _{v}^{2} \omega _{a} \omega _{e} \sigma _{a} \sigma _{e} \sigma _{v}^{2} (1-u) \vartheta _e\nonumber \\{} & {} + \Lambda _{e} \Lambda _{v}^{2} b_{1}^{2} b_{2} \beta _{a} \beta _{e} \beta _{v} (1-\varepsilon ) k_{0} k_{2} k_{4} k_{5} \mu _{h} \mu _{v} \sigma _{a} \sigma _{v}^{2} \xi - \Lambda _{e} \Lambda _{v}^{2} b_{1} b_{2} \beta _{a} \beta _{e} (1-\varepsilon ) k k_{0} k_{1} k_{2} k_{3} k_{4} k_{5} \mu _{v}^{2} \sigma _{v}^{2} \nonumber \\{} & {} - \Lambda _{e} \Lambda _{v}^{2} b_{2}^{2} \beta _{e}^{2} (1-\varepsilon ) k_{0} k_{1} k_{2} k_{3} k_{4} k_{5} \mu _{h} \mu _{v}^{2} \sigma _{v}^{2} + \Lambda _{e} \Lambda _{v}^{2} b_{1} b_{2}^{2} \beta _{a} \beta _{e} \beta _{v} k_{1} k_{3} k_{4} k_{5} \mu _{h}^{2} \mu _{v} \sigma _{e} \sigma _{v}^{2} \nonumber \\{} & {} + \Lambda _{e} \Lambda _{v}^{2} b_{2}^{2} \beta _{e}^{2} (1-\varepsilon ) k_{1} k_{3} k_{4} k_{5} \mu _{h} \mu _{v}^{2} \omega _{e} \sigma _{e} \sigma _{v}^{2} + \Lambda _{e} \Lambda _{v}^{2} b_{1} b_{2} \beta _{a} \beta _{e} k_{0} k_{2} k_{4} k_{5} \mu _{h} \mu _{v}^{2} \omega _{a} \sigma _{a} \sigma _{v}^{2} u\nonumber \\{} & {} - \Lambda _{e} \Lambda _{v}^{2} b_{1} b_{2} \beta _{a} \beta _{e} k_{4} k_{5} \mu _{h} \mu _{v}^{2} \omega _{a} \omega _{e} \sigma _{a} \sigma _{e} \sigma _{v}^{2} u + \Lambda _{e} \Lambda _{v}^{2} b_{1} b_{2} \beta _{a} \beta _{e} \delta _{a} k_{0} k_{2} k_{5} \mu _{h} \mu _{v}^{2} \omega _{a} \sigma _{a} \sigma _{v}^{2} (1-u) \nonumber \\{} & {} - \Lambda _{e} \Lambda _{v}^{2} b_{1} b_{2} \beta _{a} \beta _{e} \delta _{a} k_{5} \mu _{h} \mu _{v}^{2} \omega _{a} \omega _{e} \sigma _{a} \sigma _{e} \sigma _{v}^{2} (1-u) + \Lambda _{e}^{2} \Lambda _{v} b_{2}^{2} \beta _{e} \beta _{v} (1-\varepsilon ) k_{1} k_{3} k_{4} k_{5}^{2} \mu _{h} \mu _{v}^{2} \sigma _{e} \sigma _{v} \vartheta _e\nonumber \\{} & {} + \Lambda _{e} \Lambda _{v}^{2} b_{1} b_{2} \beta _{a} \beta _{e} (1-\varepsilon ) k_{1} k_{3} k_{4} k_{5} \mu _{v}^{2} \omega _{e} \sigma _{e} \sigma _{v}^{2} v - 2 \, \Lambda _{e}^{2} \Lambda _{v} b_{2} \beta _{e} (1-\varepsilon ) k k_{0} k_{1} k_{2} k_{3} k_{4} k_{5}^{2} \mu _{v}^{3} \sigma _{v} \nonumber \\{} & {} + \Lambda _{e}^{2} \Lambda _{v} b_{2}^{2} \beta _{e} \beta _{v} k_{1} k_{3} k_{4} k_{5}^{2} \mu _{h}^{2} \mu _{v}^{2} \sigma _{e} \sigma _{v} - \Lambda _{e} \Lambda _{v}^{2} b_{1} b_{2} \beta _{a} \beta _{e} k_{0} k_{1} k_{2} k_{3} k_{4} k_{5} \mu _{h} \mu _{v}^{2} \sigma _{v}^{2} \nonumber \\{} & {} + \Lambda _{e} \Lambda _{v}^{2} b_{1} b_{2} \beta _{a} \beta _{e} k_{1} k_{3} k_{4} k_{5} \mu _{h} \mu _{v}^{2} \omega _{e} \sigma _{e} \sigma _{v}^{2} + 2 \, \Lambda _{e}^{2} \Lambda _{v} b_{1} \beta _{a} k k_{0} k_{2} k_{4} k_{5}^{2} \mu _{v}^{3} \omega _{a} \sigma _{a} \sigma _{v} u + 2 \, \Lambda _{e}^{2} \Lambda _{v} b_{1} \beta _{a} \delta _{a} k k_{0} k_{2} k_{5}^{2} \mu _{v}^{3} \omega _{a} \sigma _{a} \sigma _{v} (1-u) \nonumber \\{} & {} + 2 \, \Lambda _{e}^{2} \Lambda _{v} b_{2} \beta _{e} (1-\varepsilon ) k_{1} k_{3} k_{4} k_{5}^{2} \mu _{v}^{3} \omega _{e} \sigma _{e} \sigma _{v} v + \Lambda _{e}^{2} \Lambda _{v} b_{1}^{2} \beta _{a} \beta _{v} k_{0} k_{2} k_{4} k_{5}^{2} \mu _{h} \mu _{v}^{2} \sigma _{a} \sigma _{v} \xi - 2 \, \Lambda _{e}^{2} \Lambda _{v} b_{1} \beta _{a} k k_{0} k_{1} k_{2} k_{3} k_{4} k_{5}^{2} \mu _{v}^{3} \sigma _{v} \nonumber \\{} & {} - 2 \, \Lambda _{e}^{2} \Lambda _{v} b_{2} \beta _{e} k_{0} k_{1} k_{2} k_{3} k_{4} k_{5}^{2} \mu _{h} \mu _{v}^{3} \sigma _{v} + 2 \, \Lambda _{e}^{2} \Lambda _{v} b_{2} \beta _{e} k_{1} k_{3} k_{4} k_{5}^{2} \mu _{h} \mu _{v}^{3} \omega _{e} \sigma _{e} \sigma _{v} - 3 \, \Lambda _{e}^{3} k k_{0} k_{1} k_{2} k_{3} k_{4} k_{5}^{3} \mu _{v}^{4}\nonumber \\ a_1= & {} -\Lambda _{e} \Lambda _{v}^{2} b_{1} b_{2}^{2} \beta _{a} \beta _{e} \beta _{v} (1-\varepsilon ) k_{4} k_{5} \mu _{h} \mu _{v}^{2} \omega _{a} \sigma _{a} \sigma _{e} \sigma _{v}^{2} u \vartheta _e - \Lambda _{e} \Lambda _{v}^{2} b_{1} b_{2}^{2} \beta _{a} \beta _{e} \beta _{v} \delta _{a} (1-\varepsilon ) k_{5} \mu _{h} \mu _{v}^{2} \omega _{a} \sigma _{a} \sigma _{e} \sigma _{v}^{2} (1-u) \vartheta _e\nonumber \\{} & {} + \Lambda _{e} \Lambda _{v}^{2} b_{2}^{3} \beta _{e}^{2} \beta _{v} (1-\varepsilon ) k_{1} k_{3} k_{4} k_{5} \mu _{h}^{2} \mu _{v}^{2} \sigma _{e} \sigma _{v}^{2} - \Lambda _{e} \Lambda _{v}^{2} b_{1} b_{2}^{2} \beta _{a} \beta _{e} \beta _{v} k_{4} k_{5} \mu _{h}^{2} \mu _{v}^{2} \omega _{a} \sigma _{a} \sigma _{e} \sigma _{v}^{2} u\nonumber \\{} & {} - \Lambda _{e} \Lambda _{v}^{2} b_{1} b_{2}^{2} \beta _{a} \beta _{e} \beta _{v} \delta _{a} k_{5} \mu _{h}^{2} \mu _{v}^{2} \omega _{a} \sigma _{a} \sigma _{e} \sigma _{v}^{2} (1-u) + \Lambda _{e} \Lambda _{v}^{2} b_{1} b_{2}^{2} \beta _{a} \beta _{e} \beta _{v} (1-\varepsilon ) k_{1} k_{3} k_{4} k_{5} \mu _{h} \mu _{v}^{2} \sigma _{e} \sigma _{v}^{2} \vartheta _e \nonumber \\{} & {} + \Lambda _{e} \Lambda _{v}^{2} b_{1}^{2} b_{2} \beta _{a} \beta _{e} \beta _{v} (1-\varepsilon ) k_{0} k_{2} k_{4} k_{5} \mu _{h} \mu _{v}^{2} \sigma _{a} \sigma _{v}^{2} \xi + \Lambda _{e} \Lambda _{v}^{2} b_{1} b_{2}^{2} \beta _{a} \beta _{e} \beta _{v} k_{1} k_{3} k_{4} k_{5} \mu _{h}^{2} \mu _{v}^{2} \sigma _{e} \sigma _{v}^{2} \nonumber \\{} & {} + 2 \, \Lambda _{e}^{2} \Lambda _{v} b_{2}^{2} \beta _{e} \beta _{v} (1-\varepsilon ) k_{1} k_{3} k_{4} k_{5}^{2} \mu _{h} \mu _{v}^{3} \sigma _{e} \sigma _{v} v - \Lambda _{e}^{2} \Lambda _{v} b_{2} \beta _{e} (1-\varepsilon ) k k_{0} k_{1} k_{2} k_{3} k_{4} k_{5}^{2} \mu _{v}^{4} \sigma _{v}\nonumber \\{} & {} + 2 \, \Lambda _{e}^{2} \Lambda _{v} b_{2}^{2} \beta _{e} \beta _{v} k_{1} k_{3} k_{4} k_{5}^{2} \mu _{h}^{2} \mu _{v}^{3} \sigma _{e} \sigma _{v} + \Lambda _{e}^{2} \Lambda _{v} b_{1} \beta _{a} k k_{0} k_{2} k_{4} k_{5}^{2} \mu _{v}^{4} \omega _{a} \sigma _{a} \sigma _{v} u + \Lambda _{e}^{2} \Lambda _{v} b_{1} \beta _{a} \delta _{a} k k_{0} k_{2} k_{5}^{2} \mu _{v}^{4} \omega _{a} \sigma _{a} \sigma _{v} (1-u) \nonumber \\{} & {} + \Lambda _{e}^{2} \Lambda _{v} b_{2} \beta _{e} (1-\varepsilon ) k_{1} k_{3} k_{4} k_{5}^{2} \mu _{v}^{4} \omega _{e} \sigma _{e} \sigma _{v} \vartheta _e + 2 \, \Lambda _{e}^{2} \Lambda _{v} b_{1}^{2} \beta _{a} \beta _{v} k_{0} k_{2} k_{4} k_{5}^{2} \mu _{h} \mu _{v}^{3} \sigma _{a} \sigma _{v} \xi - \Lambda _{e}^{2} \Lambda _{v} b_{1} \beta _{a} k k_{0} k_{1} k_{2} k_{3} k_{4} k_{5}^{2} \mu _{v}^{4} \sigma _{v} \nonumber \\{} & {} - \Lambda _{e}^{2} \Lambda _{v} b_{2} \beta _{e} k_{0} k_{1} k_{2} k_{3} k_{4} k_{5}^{2} \mu _{h} \mu _{v}^{4} \sigma _{v} + \Lambda _{e}^{2} \Lambda _{v} b_{2} \beta _{e} k_{1} k_{3} k_{4} k_{5}^{2} \mu _{h} \mu _{v}^{4} \omega _{e} \sigma _{e} \sigma _{v} - 3 \, \Lambda _{e}^{3} k k_{0} k_{1} k_{2} k_{3} k_{4} k_{5}^{3} \mu _{v}^{5} \nonumber \\ a_0= & {} N_h^{*2}kk_0k_1k_2k_3k_4k_5^3\mu _h^2\mu _v^6 \left( {\mathcal {R}}_C^m-1 \right) . \end{aligned}$$

(17)