This section presents a compartmental epidemiological model to capture the interdependent dynamics of malaria transmission between human and mosquito populations, a necessity given the disease’s reliance on cross-species interaction. The human population is divided into three compartments: susceptible (\(S_h\)), infected (\(I_h\)), and recovered (\(R_h\)). In parallel, the mosquito population is categorized into susceptible (\(S_m\)) and infected (\(I_m\)) compartments. The model accounts for bidirectional transmission mechanisms: infected mosquitoes (\(I_m\)) transmit the parasite to susceptible humans (\(S_h\)), increasing the \(I_h\) population, while infected humans (\(I_h\)) subsequently infect susceptible mosquitoes (\(S_m\)), driving the rise of \(I_m\). These interactions are governed by density-dependent transmission rates, reflecting real-world contact patterns. The model assumes constant birth and death rates, homogeneous mixing of the population, and no demographic or spatial heterogeneity. Figure 2 schematically represents the coupled transmission pathways, emphasizing the feedback loop central to malaria’s persistence.
Schematic diagram of SIR-SI system.
For the human population, we consider the standard SIR model, and for the mosquito population, we consider the SI model. The following is the system of differential equations.
$$\begin{aligned} \frac{dS_{h}}{dt}= & \Gamma _{h}N_{h}-\frac{\beta _{h}S_{h}I_{m}}{N_{h}}- \mu _{h}S_{h} \nonumber \\ \frac{dI_{h}}{dt}= & \frac{\beta _{h}S_{h}I_{m}}{N_{h}} – (\gamma _{h} + \mu _{h})I_{h}\nonumber \\ \frac{dR_{h}}{dt}= & \gamma _{h} I_{h} -\mu _{h} R_{h}\nonumber \\ \frac{dS_{m}}{dt}= & \Gamma _{m}N_{m}-\frac{\beta _{m}S_{m}I_{h}}{N_{m}}- \mu _{m}S_{m} \nonumber \\ \frac{dI_{m}}{dt}= & \frac{\beta _{m}S_{m}I_{h}}{N_{m}} – \mu _{m}I_{m}. \end{aligned}$$
(1)
The model used in this work is the SIR-SI model. The SIR model is used for the human population, and the SI model is used for the mosquito population. Table 1 and 2 explain the compartments and the parameters used.
Disease-free steady state analysis
In this section, we perform the steady-state analysis. Steady-state solutions play an important role when the analytical solution is not known, and we want to study the qualitative properties of solutions.
We define the basic reproduction number for the model (1) as
$$R_{0} = \frac{\beta _{h}\beta _{m}\Gamma _{h}\Gamma _{m}}{\mu _{h}\mu _{m}^{2}(\mu _{h}+\gamma _{h})}.$$
Observe that the basic reproduction number depends on recruitment rates, infection rates, recovery rates, and mortality rates.
Theorem 1
If \(R_{0}<1\), the disease-free steady state is locally stable.
Proof
For the analysis of the disease-free steady state, we need to equate the infected and recovered populations of both species to zero. Therefore, we obtain \(S_{h} = \frac{\Gamma _{h}}{\mu _{h}}\), \(I_{h} = 0\), \(R_{h} = 0\), \(S_{m} = \frac{\Gamma _{m}}{\mu _{m}}\), \(I_{m} = 0.\)
Dividing the first three equations \(N_{h}\) and the last two equations with \(N_{m}\) in (1), we get
$$\begin{aligned} \frac{dS_{h}}{dt}= & \Gamma _{h} – m\beta _{h}S_{h}I_{m} – \mu _{h}S_{h} \nonumber \\ \frac{dI_{h}}{dt}= & m\beta _{h}S_{h}I_{m} -(\gamma _{h} – \mu _{h})I_{h} \nonumber \\ \frac{dR_{h}}{dt}= & \gamma _{h}I_{h} – \mu _{h}R_{h} \nonumber \\ \frac{dS_{m}}{dt}= & \Gamma _{m} – \frac{\beta _{m}S_{m}I_{h}}{m} – \mu _{m}S_{m} \nonumber \\ \frac{dI_{m}}{dt}= & \frac{\beta _{m}S_{m}I_{h}}{m} – \mu _{m}I_{m}. \end{aligned}$$
(2)
Now our task is to linearize the system (2) around the disease-free steady state. After linearization, we obtain the following system
$$\begin{aligned} \frac{dS_{h}}{dt}= & \Gamma _{h} – m\beta _{h}\frac{\Gamma _{h}}{\mu _{h}}I_{m} – \mu _{h}S_{h} \\ \frac{dI_{h}}{dt}= & m\beta _{h}\frac{\Gamma _{h}}{\mu _{h}}I_{m} – (\gamma _{h} + \mu _{h})I_{h}\\ \frac{dR_{h}}{dt}= & \gamma _{h}I_{h} – \mu _{h}R_{h}\\ \frac{dS_{m}}{dt}= & \Gamma _{m} – \beta _{m}\frac{\Gamma _{m}}{m\mu _{m}}I_{h} – \mu _{m}S_{m} \\ \frac{dI_{m}}{dt}= & \frac{\Gamma _{m}\beta _{m}I_{h}}{m\mu _{m}} – \mu _{m}I_{m}. \end{aligned}$$
Constructing the Jacobian matrix, we get
$$J = \begin{bmatrix} -\mu _{h} & 0 & 0 & 0 & -\frac{m\beta _{h}\Gamma _{h}}{\mu _{h}} \\ 0 & -(\gamma _{h}+ \mu _{h}) & 0 & 0 & \frac{m\beta _{h}\Gamma _{h}}{\mu _{h}}\\ 0 & \gamma _{h} & -\mu _{h} & 0 & 0 \\ 0 & \frac{-\beta _{m}\Gamma _{m}}{m\mu _{m}} & 0 & -\mu _{m} & 0 \\ 0 & \frac{\beta _{m}\Gamma _{m}}{m\mu _{m}}& 0 & 0 & -\mu _{m} \\ \end{bmatrix}.$$
The characteristic polynomial of the above matrix is
$$\begin{aligned} f(\lambda ) = -(\lambda + \mu _{h})(\lambda + \mu _{h})(\lambda + \mu _{m}) \left[ (\lambda + \mu _{m})(\lambda + \mu _{h} + \gamma _{h}) – \frac{\beta _{h}\beta _{m}\Gamma _{h}\Gamma _{m}}{\mu _{m}\mu _{h}}\right] . \end{aligned}$$
For the system to be stable, all the eigenvalues must have negative real parts. It is clear that three of the eigenvalues are negative; we need to check the roots of the quadratic polynomial for the remaining two eigenvalues. By solving the quadratic equation, we get the roots as
$$\begin{aligned} & \frac{-(\gamma + \mu _{h}+\mu _{m}) + \sqrt{(\mu _{m}+ \mu _{h} + \gamma _{h})^{2} – 4(\mu _{m}(\mu _{h}+ \gamma _{h}) – \frac{\beta _{h}\beta _{m}\Gamma _{h}\Gamma _{m}}{\mu _{h}\mu _{m}}) }}{2}\\ & \frac{-(\gamma + \mu _{h}+\mu _{m}) – \sqrt{(\mu _{m}+ \mu _{h} + \gamma _{h})^{2} – 4(\mu _{m}(\mu _{h}+ \gamma _{h}) – \frac{\beta _{h}\beta _{m}\Gamma _{h}\Gamma _{m}}{\mu _{h}\mu _{m}}) }}{2}. \end{aligned}$$
If we observe, the second root is always negative, and thus we need to find out the condition for which the first root is negative, and that is
$$\mu _{m}( \mu _{h} + \gamma _{h}) – \frac{\beta _{m}\beta _{h}\Gamma _{m}\Gamma _{h}}{\mu _{h}\mu _{m}}> 0.$$
Rearranging the above inequality, we get the expression
$$\frac{\beta _{h}\beta _{m}\Gamma _{h}\Gamma _{m}}{\mu _{h}\mu _{m}^{2}(\mu _{h}+\gamma _{h})}< 1.$$
So, if \(R_{0}<1,\) it will ensure us that the disease-free steady state is stable. \(\square\)
Endemic steady state analysis
In this subsection, we aim to study the stability of the endemic steady-state.
The given system of equations is
$$\begin{aligned} \frac{dS_{h}}{dt}= & \Gamma _{h} – m\beta _{h}\frac{\Gamma _{h}}{\mu _{h}}I_{m} – \mu _{h}S_{h}\\ \frac{dI_{h}}{dt}= & m\beta _{h}\frac{\Gamma _{h}}{\mu _{h}}I_{m} – (\gamma _{h} + \mu _{h})I_{h}\\ \frac{dR_{h}}{dt}= & \gamma _{h}I_{h} – \mu _{h}R_{h}\\ \frac{dS_{m}}{dt}= & \Gamma _{m} – \beta _{m}\frac{\Gamma _{m}}{m\mu _{m}}I_{h} – \mu _{m}S_{m} \\ \frac{dI_{m}}{dt}= & \frac{\Gamma _{m}\beta _{m}I_{h}}{m\mu _{m}} – \mu _{m}I_{m}. \end{aligned}$$
Before proceeding with the next theorem, it is essential to define the following quantities
$$\begin{aligned} b= & \mu _{h} + \frac{\beta _{h}I_{m}^{*}}{N_{h}}+ \mu _{h} + \gamma + \frac{\beta _{h}I_{h}^{*}}{N_{m}} + \mu _{m}\\ c= & (\mu _{h} + \frac{\beta _{h}I_{m}^{*}}{N_{h}})(\mu _{h} + \gamma ) + (\mu _{h} + \gamma )(\frac{\beta _{h}I_{h}^{*}}{N_{m}} + \mu _{m})+(\frac{\beta _{h}I_{h}^{*}}{N_{m}} + \mu _{m})(\mu _{h} + \frac{\beta _{h}I_{m}^{*}}{N_{h}})- \frac{\beta _{m}\beta _{h}S_{m}^{*}S_{h}^{*}}{N_{m}N_{h}}\\ d= & \left( \mu _{h} + \frac{\beta _{h}I_{m}^{*}}{N_{h}}\right) \left[ (\mu _{h} + \gamma )(\frac{\beta _{h}I_{h}^{*}}{N_{m}} + \mu _{m})-\frac{\beta _{m}\beta _{h}S_{m}^{*}S_{h}^{*}}{N_{m}N_{h}} + \frac{\beta _{m}\beta _{h}^{2}S_{m}^{*}S_{h}^{*}I_{m}^{*}}{N_{m}N_{h}^{*}}\right] , \end{aligned}$$
where \(S^{*}_{m},S^{*}_{h},I^{*}_{m},I^{*}_{h}\) are the non trivial equilibrium solutions.
Theorem 2
If \(R_{0}>1\) and \(b, c, d, bc-d >0\), the endemic steady state is locally stable.
Proof
The equilibrium points will be obtained by equating all of the above time derivatives to zero, and solving them for non-trivial solutions, we obtain the solutions as
$$\begin{aligned} S_{h}^{*}= & \frac{N_{h}(\gamma + \mu _{h})(\mu _{m}N_{m}+ \frac{\Gamma _{h}N_{h}\beta _{m}}{\gamma + \mu _{h}})}{\beta _{m}(\mu _{h}N_{h}+ \frac{\Gamma _{m}N_{m}\beta _{h}}{\mu _{m}})} \\ I_{h}^{*}= & \frac{N_{h}N_{m}\mu _{m}\mu _{h}(R_{0}-1)}{\beta _{m}(\mu _{h}N_{h}+ \frac{\Gamma _{m}N_{m}\beta _{h}}{\mu _{m}})} \\ R_{h}^{*}= & \frac{N_{h}N_{m}\mu _{m}\gamma _{h}(R_{0}-1)}{\beta _{m}(\mu _{h}N_{h}+ \frac{\Gamma _{m}N_{m}\beta _{h}}{\mu _{m}})} \\ S_{m}^{*}= & \frac{N_{m}\mu _{m}(\mu _{h}N_{h}+ \frac{\Gamma _{m}N_{m}\beta _{h}}{\mu _{m}})}{\beta _{h}(\mu _{m}N_{m}+ \frac{\Gamma _{h}N_{h}\beta _{m}}{\gamma + \mu _{h}})}\\ I_{m}^{*}= & \frac{N_{h}N_{m}\mu _{m}\mu _{h}(R_{0}-1)}{\beta _{h}(\mu _{m}N_{m}+ \frac{\Gamma _{h}N_{h}\beta _{m}}{\gamma + \mu _{h}})}. \end{aligned}$$
Linearizing the model (2) around the endemic equilibrium point, we get the following system
$$\begin{aligned} \frac{dS_{h}}{dt}= & \Gamma _{h}N_{h}-\frac{\beta _{h}}{N_{h}}(S_{h}^{*}I_{m} + I_{m}^{*}S_{h}- S_{h}^{*}I_{m}^{*})- \mu _{h}S_{h}\\ \frac{dI_{h}}{dt}= & \frac{\beta _{h}}{N_{h}}(S_{h}^{*}I_{m} + I_{m}^{*}S_{h}-S_{h}^{*}I_{m}^{*})- \gamma I_{h} -\mu _{h}I_{h} \\ \frac{dR_{h}}{dt}= & \gamma I_{h} – \mu _{h}R_{h} \\ \frac{dS_{m}}{dt}= & \gamma _{m}N_{m}-\frac{\beta _{m}}{N_{m}}(S_{m}^{*}I_{h}+ S_{m} I_{h}*) – \mu _{m}S_{m} \\ \frac{dI_{m}}{dt}= & \frac{\beta _{m}}{N_{m}}(S_{m}^{*}I_{h}+ S_{m} I_{h}*) – \mu _{m}I_{m}. \end{aligned}$$
Writing the Jacobian matrix, we get
$$J = \begin{bmatrix} -\frac{\beta _{h}I_{m}^{*}}{N_{h}}-\mu _{h} & 0 & 0 & 0 & -\frac{\beta _{h}S_{h}^{*}}{N_{h}}\\ \frac{\beta _{h}I_{m}^{*}}{N_{h}} & -(\gamma _{h}+ \mu _{h}) & 0 & 0 & \frac{\beta _{h}S_{h}^{*}}{N_{h}} \\ 0 & \gamma _{h} & -\mu _{h} & 0 & 0 \\ 0 & -\frac{\beta _{h}S_{m}^{*}}{N_{m}} & 0 & -\frac{\beta _{h}S_{m}^{*}}{N_{m}}-\mu _{m} & 0 \\ 0 & \frac{\beta _{h}S_{m}^{*}}{N_{m}} & 0 & \frac{\beta _{h}S_{m}^{*}}{N_{m}} & -\mu _{m} \\ \end{bmatrix}$$
The characteristic polynomial of the above matrix is
$$f(\lambda ) = -\left( \lambda + \mu _{h}\right) (\lambda + \mu _{m}) \left( \left[ \lambda + \mu _{h} + \frac{\beta _{h}I_{m}^{*}}{N_{h}} \right] \left[ \left( \lambda + \mu _{h}+\gamma \right) \left( \lambda + \frac{\beta _{m}I_{h}^{*}}{N_{m}}+\mu _{m} \right) -\frac{\beta _{m}\beta _{h}S_{m}^{*}S_{h}^{*}}{N_{h}N_{m}} \right] + \frac{\beta _{h}^{2}\beta _{m}S_{h}^{*}I_{m}^{*}S_{m}^{*}}{N_{h}^{2}N_{m}} \right) .$$
From the above equation, we can see that we have two linear factors and a cubic factor. To analyze the cubic factor, let us state the following lemma:
Lemma 1
Let \(f(x) = ax^3 + bx^2 + cx + d\) be a cubic polynomial. For f(x) to have all negative roots or complex roots with negative real parts, the following conditions are necessary:
$$a> 0, b> 0, c> 0, d>0, bc – ad >0.$$
Now, by using the above-stated lemma, we can arrive at our required result. \(\square\)
Remark 1
Equilibrium points are crucial for evaluating malaria control efforts. A disease-free equilibrium indicates that transmission can be halted through vector control and drug administration, while an endemic equilibrium suggests sustained transmission, requiring ongoing interventions. Stability analysis helps identify critical thresholds, such as intervention coverage or mosquito density, to guide effective control strategies.
Numerical validation of the stability theorems
In the previous section, we established that the disease-free steady state is attained when \(R_0 < 1\). Accordingly, we selected \(R_0 = 0.04, 0.16,\) and \(0.51\), ensuring \(R_0 < 1\). As shown in Figs. 3, 4, and 5, the infected populations of humans and mosquitoes decline over time.
For very small value of \(R_0 = 0.04\), infected population converges to zero.
For \(R_0 = 0.16\), firstly infected population increases and then diminishes to zero.
Behaviour of infected population for \(R_0 = 0.51\). It reveals that the number of infected mosquitoes initially increases for a brief period before diminishing and eventually converging to zero.
From Figs. 3, 4, and 5, we observe that the proposed model for malaria spread exhibits stable dynamics when the reproduction number remains below 1. Furthermore, as the reproduction number increases, the peak infection level intensifies, indicating a heightened disease burden within the population.
Importance of steady-state analysis and equilibrium points
Steady-state analysis in malaria modeling provides insights into the long-term behavior of transmission dynamics, focusing on equilibrium points where the rate of new infections and recoveries stabilizes. These equilibrium points are critical in understanding malaria persistence and eradication potential in specific regions. Identifying whether the system reaches a disease-free equilibrium or a sustained endemic state allows for targeted intervention strategies.
For malaria control, equilibrium points are essential for assessing the effectiveness of interventions. A disease-free equilibrium indicates that malaria transmission can be halted, typically through a combination of vector control measures such as insecticide-treated nets, indoor spraying, and mass drug administration. Conversely, an endemic equilibrium suggests that transmission is sustained within the population despite interventions, indicating that continuous, long-term strategies such as routine treatment, surveillance, and seasonal control programs are needed to keep malaria prevalence low. Stability analysis of these equilibrium points helps to identify critical thresholds, such as coverage levels for interventions or mosquito density, that determine whether malaria transmission will be suppressed or continue to persist, guiding the implementation of more effective control measures.
$$\frac{\beta _{h}\beta _{m}\Gamma _{h}\Gamma _{m}}{\mu _{h}\mu _{m}^{2}(\mu _{h}+\gamma _{h})}< 1.$$
This inequality indicates that if the basic reproduction number \(R_0\) is less than 1, it ensures that the disease-free steady state is stable, further supporting the potential for malaria eradication through effective control strategies.
Temperature and altitude dependence of the transmission rate
Several factors influence malaria transmission; however, Patz et al.44 identify temperature and altitude as among the most significant. Therefore, we model transmission as a function of temperature and altitude as follows:
$$\begin{aligned} \beta (T,h) = \beta _{0}e^{-\frac{(T-25)^{2}}{\eta ^{2}}}e^{-\frac{h^{2}}{\xi ^{2}}}(1-e^{-\frac{h^{2}}{\xi ^{2}}}), \end{aligned}$$
(3)
where \(\beta _{0}\) is the transmission constant of the region, T is the temperature of the region, h is the altitude, \(\eta\) and \(\xi\) are the constants associated with the region’s temperature and height, respectively. The Gaussian function is centered at \(25^{\circ }C\), as this temperature is biologically optimal for mosquito survival, leading to maximum transmission. This aligns with the findings of Shapiro et al.45, which also identify a similar optimal temperature for malaria mosquito breeding.
Here \(e^{-\frac{(T-25)^{2}}{\eta ^{2}}}\) is used to model the temperature and \(e^{-\frac{h^{2}}{\xi ^{2}}}(1-e^{-\frac{h^{2}}{\xi ^{2}}})\) is used to model the height. Malaria transmission will be very minimal when the temperature is either extremely high or it is extremely low and thus to model this variation, the Gaussian function is used and the reason for shifting it by 25 is because the optimum temperature for mosquito’s existence and malaria transmission is \(25\,^\circ \textrm{C}\). Also, malaria transmission is completely zero when the altitude is zero since there won’t be any mosquitoes in the sea and in the same way when the altitude is extremely high again the transmission is completely zero since there are no mosquitoes in the space and thus to model both of these conditions the negative exponential function is used in this manner.
At a constant altitude of 75 meters, temperature variations above and below \(25^\circ \textrm{C}\) were rigorously assessed to determine their influence on transmission dynamics. Our analysis reveals that \(25^\circ \textrm{C}\) serves as the optimal temperature for maximizing transmission rates. Analysis demonstrates a consistent decline in both infected mosquito and human populations as temperatures rise beyond this threshold. Furthermore, population trends for infected mosquitoes and humans align with a Gaussian distribution, reflecting a symmetrical, bell-shaped relationship with temperature. These findings underscore the critical role of temperature regulation in transmission mitigation strategies.
Assuming the temperature is from \(T_1\) to \(T_2\) and the height is from \(h_1\) to \(h_2\), we can write
$$\begin{aligned} \beta _{avg} = \int _{T_{1}}^{T_{2}} \int _{h_{1}}^{h_{2}}e^{-\frac{(T-25)^{2}}{\eta ^{2}}}e^{-\frac{h^{2}}{\xi ^{2}}}(1-e^{-\frac{h^{2}}{\xi ^{2}}})dh dT. \end{aligned}$$
(4)
Here, the effect of transmission rate is studied (Figs. 6 – 8) by changing the temperature values for a fixed height. The parameters considered are the following: \(\beta _{0} = 10\), and for human we have taken \(\eta =200 ,\xi = 20000\) and for mosquitoes, \(\eta = 400 , \xi = 40000\). The different temperature values which we used are \(25^\circ\), \(30^\circ\), \(35^\circ\), \(40^\circ\), \(45 ^\circ\).
When the altitude was maintained at \(100\; m\), temperature variations were analyzed to assess their impact on transmission dynamics. The findings confirm that \(25^\circ \textrm{C}\) is the optimal temperature for maximizing the transmission rate. Additionally, the population density of infected mosquitoes was higher compared to an altitude of \(75\; m\), though the variation was marginal.
At a fixed altitude of \(125\; m\), the impact of temperature variations on transmission dynamics was analyzed. The findings consistently indicate that \(25^\circ \textrm{C}\) remains the optimal temperature for maximizing the transmission rate. A similar trend was observed, with a slight variation in the population of infected mosquitoes, reaching its peak at a higher point.
At a constant temperature of \(28^\circ \textrm{C}\), altitude variations were studied. Population density trends for infected mosquitoes and humans exhibit a Gaussian distribution, reflecting a symmetrical, bell-shaped relationship with altitude. These findings highlight the significant impact of altitude on disease transmission and emphasize the importance of considering elevation in vector control strategies.
With the temperature maintained at \(35^\circ \textrm{C}\), the impact of altitude fluctuations on transmission dynamics was further examined. A similar trend was observed; however, at \(35^\circ \textrm{C}\), the peak population of infected mosquitoes occurred at lower altitudes.
The impact of altitude variations on transmission dynamics was further analyzed at a fixed temperature of \(42^\circ \textrm{C}\). In this case, a more significant decline in the population of infected mosquitoes was observed. These findings highlight the critical role of altitude in shaping transmission patterns and optimizing intervention strategies.
From Figs. 9 to 11, we observe that the transmission rate is highest when the altitude is 150. This represents the maximum value among the altitudes considered. Therefore, we conclude that \(h = 150\) may serve as a threshold value for altitude in the model. The function \(f(h) = e^{-\frac{h^{2}}{\xi ^{2}}}(1-e^{-\frac{h^{2}}{\xi ^{2}}})\) achieves its highest value when \(h = \xi {\sqrt{\ln (2)}}\), and thus the transmission rate increases as the height value approaches \(h = {\xi \sqrt{\ln (2)}}\). By substituting both values of \(\xi\) and taking the average, we get approximately 142.1; thus, the optimal value for altitude can be concluded to be in the range \(140 – 150\), and since it can be observed that the trajectory for all of these values is nearly the same.
