Dynamical analysis of a nonlocal delays spatial malaria model with Wolbachia-infected male mosquitoes release

1.   Introduction

Malaria remains a major global health challenge, primarily affecting tropical and subtropical regions. The disease is caused by Plasmodium parasites, transmitted through bites of infected Anopheles mosquitoes. Symptoms include fever, chills, and flu-like illness, and without treatment, malaria can lead to severe complications and death ^[1]. According to the World Health Organization (WHO), in 2023 there were an estimated 263 million malaria cases and 597,000 malaria deaths worldwide. The African region carries a disproportionately high share of the global malaria burden, accounting for 94% of malaria cases and 95% of malaria deaths ^[2]. Despite ongoing efforts to control malaria, it continues to be a major global health issue that is especially prevalent in the sub-Saharan African region ^[3].

Mathematical models are essential for depicting the dynamics of disease transmission and developing effective control strategies ^[4,5,6]. Early models often assumed a homogeneous environment, ignoring spatial variations such as population density, mosquito distribution, or local ecological conditions. However, newer studies have focused on the importance of spatial heterogeneity in malaria transmission and the need for reaction-diffusion models that incorporate these variations ^[7]. Spatial heterogeneity significantly impacts malaria transmission, especially in regions with diverse ecological factors. Moreover, delays, such as the incubation periods of malaria in humans and mosquitoes, are crucial components of modern models. These time lags affect the stability and transmission dynamics of malaria. By including these delays in reaction-diffusion malaria models, researchers can better simulate real-world scenarios and assess the effectiveness of control measures over time ^[8]. For articles on the variable malaria model, please refer to ^[9,10].

Recent studies have shown promising results with the use of Wolbachia in controlling mosquito-borne diseases, including malaria and dengue. Field trials conducted in various parts of the world have demonstrated that releasing Wolbachia-infected Aedes aegypti mosquitoes can significantly reduce dengue transmission ^[11]. Similarly, the potential of Wolbachia to control malaria transmission has been explored through laboratory and field studies involving Anopheles mosquitoes, the primary malaria vectors ^[12]. Moreover, studies in malaria-endemic areas, including regions in Indonesia and Australia, have demonstrated that Wolbachia-infected mosquitoes can successfully integrate into local mosquito populations. This indicates a viable and potentially large-scale strategy for managing malaria ^[13,14]. Consequently, it is of great importance to investigate the effects of Wolbachia-based methods that disrupt insect reproduction on malaria control efforts.

This paper aims to study the effect of releasing Wolbachia-infected male mosquitoes on malaria control. The paper proceeds with the following structure. The next section details the development of our mathematical model. The well-posedness of the model is studied in Section 3. Section 4 delves into the introduction of the basic reproduction number, denoted as $R_0$. Section 5 focuses on analyzing the threshold dynamics of our model, which are contingent upon the value of $R_0$. Section 6 presents numerical simulations to elucidate the theoretical results, explore the impact of the release ratio of Wolbachia-infected males on transmission and control of malaria, and perform a sensitivity analysis. The paper concludes with a summary in the final section.

2.   Model formulation

2.1.   Age-structured model with mosquitoes

The mosquito population has two subclasses: the aquatic population and the winged population. Let $L_1(x, a_1, t)$ represent wild mosquitoes at time $t$, position $x$, and chronological age $a_1. \; A(x, t)$ aquatic mosquitoes at time $t$ and position $x$. $W(x, t)$ represents female winged mosquitoes at time $t$ and position $x$. $\Omega$ is a bounded region with a smooth boundary $\partial\Omega.$ In the aquatic stage, limited habitat space and food resources lead to intense intraspecific competition among the larvae. Consistent with the approach proposed in ^[15], we utilize the following classic age-structured equations and introduce an additional nonlinear factor to characterize the evolutionary dynamics of mosquito populations under such competition. For $t\geq 0, \; a_1 > 0$

where $\nabla$ indicates the gradient, $\mathbf{n}$ is the outward unit normal vector on $\partial\Omega$, $d(x, a_1)$ denotes the diffusion coefficient, and $g(L_1)$ represents a continuous function of $L$. The expressions for $d(x, a_1)$ and $g(L_1)$ are as follows:

Here, $d_w$ signifies the diffusion coefficient for adult mosquitoes. Meanwhile, $c$ indicates the intraspecific competition rate among the larvae, and $\mu_a$ and $\mu_w$ correspond to the natural mortality rates of larvae and adult mosquitoes, respectively. Then $A(x, t)$ and $W(x, t)$ can be determined by integrating the population density over the specified age ranges, as detailed below

2.2.   Releasing Wolbachia-infected male mosquitoes

The female male proportion is assumed to be $1 :1$. Assuming that Wolbachia-infected male mosquitoes are deployed under a proportional release strategy, that is, the released amount of Wolbachia-infected male mosquitoes is proportional to the current wild male mosquitoes. Here, $q(x)$ represents the ratio of released Wolbachia-infected males to wild males. Let $W_c(x, t)$ be the Wolbachia-infected male mosquitoes. Then we have $W_c(x, t) = q(x)W(x, t)$.

We assume that Wolbachia-infected male mosquitoes have the same mating ability as wild male mosquitoes. Thus, $\frac{W_c}{W+W_c} = \frac{q(x)}{1+q(x)}$ describes the probability of a wild female mosquito mating with Wolbachia-infected male mosquitoes. Then $1-\frac{q(x)}{1+q(x)} = \frac{1}{1+q(x)}$ is the probability of a wild female mosquito mating with wild male mosquitoes. Denote $\rho(x)$ as the egg-laying rate of each wild females mosquito. We have

From (2.1), we can get

It is natural to assume that $L_1(x, +\infty, t) = 0$. By the characteristic line method, we can obtain the expression of $L_1(x, \tau_1, t).$ If we set $w_1(x, a_1, s) = L_1(x, a_1, a_1+s), \; a_1\in(0, \tau_1],$ we have

We can obtain

If we let $a_1 = \tau_1$, then $s = t-\tau_1$, and we have

Thus, one has

2.2.1.   Human mosquito model

Let $S_w(x, t)$, $E_w(x, t)$, and $I_w(x, t)$ represent susceptible, exposed, and infectious (female) winged mosquitoes. Thus we have $W(x, t) = S_w(x, t)+E_w(x, t)+I_w(x, t)$. Let $S_h(x, t)$, $E_h(x, t)$, $I_h(x, t)$ and $R_h(x, t)$ represent susceptible, exposed, infectious and recovered humans, respectively. The total density of human population can be expressed by $N_h(x, t) = S_h(x, t)+E_h(x, t)+I_h(x, t)+R_h(x, t)$. Assume that $N_{h}(x, t)$ satisfies the following equation:

where $\Lambda$ and $\mu_h$ correspond to the influx rate and natural mortality rate of humans, respectively. According to Lemma 1 in ^[7], system (2.5) has a globally attractive positive steady state, denoted by $N_h^{*}(x)$. In this study, we assume $N_{h}(x, t) \equiv N_h^{*}(x)$ for $\forall\; t \geq0$ and $x\in \Omega$.

Let $a_w$ be infection age of winged mosquitoes and $L_w(x, a_w, t)$ be the density of infected winged mosquitoes. Then

So

From (2.6), we can get

It is natural to assume that $L_w(x, +\infty, t) = 0$. By the characteristic line method, we can obtain the expression of $L_w(x, \tau_w, t).$ Set $w_2(x, a_w, s) = L_w(x, a_w, a_w+s), \; a_w\in(0, \tau_w].$ We have

We can obtain

where $\Gamma_w$ is the Green function of the operator $\nabla\cdot(d_w(\cdot)\nabla)-\mu_w(\cdot)$ associated with the Neumann boundary condition, $b(x)$ represents bite rate of mosquitoes, $\beta_{hw}(x)$ represents the transmission probability from infectious humans to adult mosquitoes.

Let $a_w = \tau_w$, then $s = t-\tau_w$, and

We then have

Let $a_h$ be the infection age of humans. Similar to the derivation of $\frac{\partial E_w(x, t)}{\partial t}$ and $\frac{\partial I_w(x, t)}{\partial t}$, we can obtain the following expressions of $\frac{\partial E_h(x, t)}{\partial t}$ and $\frac{\partial I_h(x, t)}{\partial t}$:

where $\Gamma_h$ is the Green function of the operator $\nabla\cdot(d_h(\cdot)\nabla)-\mu_h(\cdot)$ associated with the Neumann boundary condition, $d_h$ represents the diffusion coefficient for humans, $\beta_{wh}$ represents the transmission probability from infectious mosquitoes to humans, and $r_h$ and $\mu_h$ represent recovery rate and death rate of infectious humans, respectively.

In a word, we can obtain the following nonlocal delays reaction-diffusion system:

We know that $A(x, t), E_h(x, t)$, and $R_h(x, t)$ are decoupled from the other equations. Then it can be applied to the following system:

where $p(x) = \frac{1}{2(1+q(x))}, \; \iota(x) = \frac{c(x)}{\mu_a(x)}$, $\beta_{w}(x) = \frac{b(x)\beta_{hw}(x)}{N^*_h(x)}, \; \beta_{h}(x) = \frac{b(x)\beta_{wh}(x)}{N^*_h(x)}.$

3.   The well-posedness

Let $\mathcal{X} = C(\overline{\Omega}, \mathbb{R}^{5})$ be the Banach space with the supremum norm $\|\cdot\|_{\mathcal{X}}$, and $\mathcal{X}^{+} = C(\overline{\Omega}, \mathbb{R}^{5}_{+})$. Let $\tau = \max\{\tau_1, \tau_w, \tau_h\}.$ Define $\mathcal{C} = C([-\tau, 0], \mathcal{X})$ as the Banach space with the norm $\|\varphi\| = \max\limits_{\theta\in [-\tau, 0]}\|\varphi(\theta)\|_{\mathcal{X}}$ for $\varphi\in \mathcal{C}$, and $\mathcal{C}^{+} = C([-\tau, 0], \mathcal{X}^{+})$. Then $(\mathcal{X}, \mathcal{X}^{+})$ and $(\mathcal{C}, \mathcal{C}^{+})$ are strongly ordered Banach spaces. For $\varsigma > 0$ and a function $z:[-\tau, \varsigma)\rightarrow \mathcal{X},$ we define $z_t\in\mathcal{C}$ by $z_t(\theta) = z(t+\theta), \; \theta\in[-\tau, 0].$ Denote $\mathcal{T}_{w}(t), \; \mathcal{T}_{s}(t)$, and $\mathcal{T}_{h}(t)$ : $C(\overline{\Omega}, \mathbb{R})\rightarrow C(\overline{\Omega}, \mathbb{R})$ are the evolution operators associated with

Then, for $\vartheta\in C(\overline{\Omega}, \mathbb{R}),$ one has

where $\Gamma_s$ is the Green function of the operator $\nabla\cdot(d_h(\cdot)\nabla)-\mu_h(\cdot)$ associated with the Neumann boundary condition. We can see that $\mathcal{T}_{j}(t)$ is strongly positive and compact, where $j = w, s, h$. Set $\mathcal{T} =$ diag$\{\mathcal{T}_{w}(t), \mathcal{T}_{w}(t), \mathcal{T}_{w}(t), \mathcal{T}_{s}(t), \mathcal{T}_{h}(t) \}$. Then $\mathcal{T} :\mathcal{X}\rightarrow \mathcal{X}$ is a semigroup generated by the operator $\mathcal{A} =$ diag $\{\mathcal{A}_1, \mathcal{A}_1, \mathcal{A}_1, \mathcal{A}_2, \mathcal{A}_3\}$ defined on $D(\mathcal{A}) = D(\mathcal{A}_1)\times D(\mathcal{A}_1) \times D(\mathcal{A}_1) \times D(\mathcal{A}_2) \times D(\mathcal{A}_3),$ in which

For $\varphi = (\varphi_1, \varphi_2, \varphi_3, \varphi_4, \varphi_5)\in \mathcal{C}^{+}$, define $\mathcal{F} = (\mathcal{F}_1, \mathcal{F}_2, \mathcal{F}_3, \mathcal{F}_4, \mathcal{F}_5):\mathcal{C}^{+}\rightarrow \mathcal{X}$ as

Then system (2.11) can be rewritten as an abstract functional equation as follows:

where $u: = (S_w, E_w, I_w, S_h, I_h)$.

According to Corollary 8.1.3 in ^[16] and Corollary 4 in ^[17], we can get the following lemma.

Lemma 3.1. System (2.11) with an initial value function $\varphi\in \mathcal{C}^+$ has a unique mild solution $u(\cdot, t, \varphi)$ on its maximal interval of existence $[0, t_\varphi)$. Moreover, $u(\cdot, t, \varphi)\in \mathcal{C}^+$ for any $t\in [0, t_\varphi)$, and $u(\cdot, t, \varphi)$ is a classical solution of system (2.11) for $t > \tau.$

Lemma 3.2. System (2.11) with an initial value function $\varphi\in \mathcal{C}^+$ has a unique global classical solution $u(\cdot, t, \varphi)$ with $u^{0} = \varphi$ for $t\in [0, \; \infty)$. Moreover, the solution semiflow $\Upsilon(t) = u_t(\cdot):\mathcal{C}^+\rightarrow \mathcal{C}^+$ has a compact global attractor in $\mathcal{C}^+$.

Proof. From the first and fourth equations of (2.11), we have

So, for any $\varphi\in \mathcal{C}^+$, $K_w, K_h > 0$ and $t_0 = t_0(\varphi) > 0$ exist such that $S_w(\cdot, t, \varphi)\leq K_w$ and $S_h(\cdot, t, \varphi)\leq K_h$ for $\forall t\geq t_0.$

Assume $M_1(t) = \int_\Omega S_w(x, t)\mathrm{d}x, \; M_2(t) = \int_\Omega E_w(x, t)\mathrm{d}x, \; M_3(t) = \int_\Omega I_w(x, t)\mathrm{d}x, \; M_4(t) = \int_\Omega S_h(x, t)\mathrm{d}x, \; M_5(t) = \int_\Omega I_h(x, t)\mathrm{d}x.$ Integrating $S_w(x, t)$ in system (2.11), we can obtain

where $K_1 = \frac{ e^{-\underline{\mu_a}\tau_1} }{ \underline{\iota}(1-e^{-\underline{\mu_a}\tau_1})}$. Thus

Integrating $E_w(x, t)$ in system (2.11), we can obtain

We can then get

It implies that

So, $M_2(t)$ is uniformly bounded.

Since $\Gamma_w(\cdot, \cdot, \cdot)$ is bounded, by integrating $I_w(x, t)$ in system (2.11), we can obtain

for some positive constant $K_2.$ One then has

where $\underline{\mu_1} = \min\{\underline{\mu_w}, \underline{\mu_a}\}$. It imply that there are a $\varphi-$dependent positive constant $K_3$ and a $\varphi-$independent positive constant $K_4$ such that

So, $M_3(t)$ is uniformly bounded.

Similarly, we can obtain

Thus

Since $\Gamma_h(\cdot, \cdot, \cdot)$ is bounded, by integrating $I_h(x, t)$ in system (2.11), we can obtain

for some positive constant $K_5.$ One then has

It imply that there are a $\varphi-$dependent positive constant $K_6$ and a $\varphi-$independent positive constant $K_7$

So, $M_5(t)$ is uniformly bounded. By the comparison theorem for delayed parabolic equations ^[16], we find that there exist $\varphi-$independent positive number $K_8, K_9, K_{10}$ and $t_1 = t_1(\varphi) > t_0(\varphi)+\tau$ such that $E_w(\cdot, t, \varphi)\leq K_8$, $I_w (\cdot, t, \varphi)\leq K_9$ and $I_h (\cdot, t, \varphi)\leq K_{10}$ for any $t\geq t_1$. So, the solution of system (2.11) is global and the solution semiflow $\Upsilon(t) = u_t(\cdot):\mathcal{C}^+\rightarrow \mathcal{C}^+$ is point dissipative by Theorem 2.1.8 in ^[16]. Based on Theorem 3.4.8 in ^[18], $\Upsilon(t)$ has a compact global attractor in $\mathcal{C}^+$ for $t\geq 0$. □

4.   The basic reproduction number

Set $E_w = I_w = I_h = 0$ in system (2.11). Then, $S_w(x, t)$ satisfies the following system:

By Theorem 11 in ^[19] and Lemma 2.1 in ^[20], system (4.1) admits a globally attractive positive steady state $W^*(x).$ It is easy to know that system (2.11) has an infection-free steady state $E_0(x) = (W^*(x), 0, 0, N_h^*(x), 0)$. Linearizing system (2.11) at $E_0(x)$, we can get

Consider the following linear nonlocal and cooperative reaction-diffusion system:

Substituting $I_w(x, t) = e^{\lambda t}\psi_1(x)$ and $I_h(x, t) = e^{\lambda t}\psi_2(x)$ into (4.2) and (4.3), we have

and

By similar methods to Theorem 7.6.1 in ^[21], we can see that (4.5) has a principal eigenvalue $\lambda_0(W^*)$ with strongly positive eigenfunctions. $\lambda_0(W^*)$ varies continuously under small perturbations $W^*$. Similar to Lemma 2.1 in ^[22], we can obtain the next lemma.

Lemma 4.1. The eigenvalue problem (4.4) admits principal eigenvalues associated with a strongly positive eigenfunction, denoted $\tilde{\lambda}_0(W^*, \tau)$. Moreover, $\tilde{\lambda}_0(W^*, \tau)$ and $\lambda_0(W^*)$ have the same sign.

Assume $\phi(x): = (\phi_w(x), \phi_h (x))$. Assume that $\phi(x)$ is the spatial distribution of the initial infective female winged mosquitos and humans. It is easy to see that the distribution of those infective individuals at time $t > 0$ is described by $\mathbb{T}(t)\phi(x): = (\mathcal{T}_w(t)\phi_w(x), \mathcal{T}_h(t)\phi_h (x))$. Then the distribution of new infectious female winged mosquitos can be expressed as

Therefore, the total distribution of infectious female winged mosquitos during the infective period is

Similarly, the total distribution of infectious humans during the infective period is

Define the operator $\mathbb{F}$ on $C(\overline{\Omega}, \mathbb{R})\times C(\overline{\Omega}, \mathbb{R})$ as

where

and

Then, we can define the next reproduction operator as

So, Based on the result in ^[23], the basic reproduction number $R_{0}$ can be defined by

Additionally, we can get the following lemma.

Lemma 4.2. $R_{0}-1$ has the same sign as $\lambda_0(W^*)$.

5.   Dynamic analysis

Theorem 5.1. If $R_{0} < 1$, then the infection-free steady state $E_{0}(x)$ of system (2.11) is globally attractive in $\mathcal{C}^+$.

Proof. When $R_{0} < 1$, according to Lemmas 4.1 and 4.2, we have $\tilde{\lambda}_0(W^*, \tau) < 0,$ Since $\lim\limits_{\zeta\rightarrow 0}\tilde{\lambda}_0(W^*+\zeta, \tau) = \tilde{\lambda}_0(W^*, \tau),$ a small enough $\zeta_0 > 0$ exists such that $\tilde{\lambda}_0(W^*+\zeta_0, \tau) < 0.$

Recall that $W^*(\cdot)$ is globally attractive for system (4.1), we have $\limsup\limits_{t\rightarrow \infty }S_w(\cdot, t)\leq W^*(\cdot).$ So, for some positive constant $t_0$, we have $S_w(\cdot, t)\leq W^*(\cdot)+\zeta_0$ for $t\geq t_0.$ Then, for all $t\geq t_0$, we can get

Set $\tilde{\psi}$ be the strongly positive eigenfunction corresponding to $\tilde{\lambda}_0(W^*+\zeta_0, \tau)$. Then, for $t > 0$, the following linear system

admits a solution $\tilde{u}(x, t) = e^{\tilde{\lambda}_0(W^*+\zeta_0, \tau)t}\tilde{\psi}(x).$ Then, for any given $\varphi\in \mathcal{C}^{+}$, there is some positive constant $\tilde{m}$, and one has $(I_w(\cdot, t, \varphi), I_h(\cdot, t, \varphi))\leq \tilde{m} \tilde{u}(\cdot, t)$ for all $t\in[t_0-\tau, t_0].$ According to the comparison principle, we can obtain

Therefore, $\lim\limits_{t\rightarrow \infty }(I_w(x, t, \varphi), I_h(x, t, \varphi)) = (0, 0)$ uniformly for $x\in\overline{\Omega}$. Then $S_w(x, t), E_w(x, t), S_h(x, t)$ in system (2.11) are asymptotic to the following equations

Therefore, $\lim\limits_{t\rightarrow \infty } E_w(x, t, \varphi) = 0$, $\lim\limits_{t\rightarrow \infty } S_w(x, t, \varphi) = W^*(x)$, $\lim\limits_{t\rightarrow \infty } S_h(x, t, \varphi) = N_h^*(x)$ uniformly for $x\in\overline{\Omega}$. This completes the proof. □

Lemma 5.1. For the solution $(S_{w}(x, t, \varphi), E_{w}(x, t, \varphi), I_{w}(x, t, \varphi), S_{h}(x, t, \varphi), I_{h}(x, t, \varphi))$ of system (2.11) with an initial value function $\varphi\in \mathcal{C}^{+}$.

(i) For any $t > 0$, $x\in\bar{\Omega}$, one has $S_{w}(x, t, \varphi) > 0$ and $S_{h}(x, t, \varphi) > 0$. Furthermore, for $\forall \; \varphi\in \mathcal{C}^{+}$, there is $\varphi-$independent positive constant $\zeta$ such that

(ii) Assume that there exists some $t^{*}\geq 0$ such that $I_{w}(\cdot, t^{*}, \varphi)\not\equiv 0$ and $I_{h}(\cdot, t^{*}, \varphi)\not\equiv 0$, then the solution satisfies

Proof. (i) From Lemma 4.1, for any $t > 0$, $x\in\bar{\Omega}$ and an initial value function $\varphi\in \mathcal{C}^{+}$, there is an $\breve{N}_1 > 0$ such that $S_{h}(x, t, \varphi)\leq \breve{N}_1$ and $I_{h}(x, t, \varphi)\leq \breve{N}_1$. Let $v_w(\cdot, t, \varphi)$ satisfy

It follows from the comparison principle that

for any $t > 0$ and $\varphi\in \mathcal{C}^{+}$. By Lemma 2.1 in ^[20], system (5.4) admits a globally attractive positive steady state $v_w^*(x).$ Set $\zeta_w = \min\limits_{ x\in\bar{\Omega} }v_w^*(x)$, then we have

Similarly, from Lemma 4.1, for any $t > 0$, $x\in\bar{\Omega}$, $\varphi\in \mathcal{C}^{+}$, there is an $\breve{N}_2 > 0$ such that $I_{w}(x, t, \varphi)\leq \breve{N}_2$.

It follows from the comparison principle that $S_h(\cdot, t, \varphi) > 0,$ and there is an $\zeta < \zeta_w$ such that $\liminf\limits_{t\rightarrow \infty } S_{h}(x, t, \varphi)\geq \zeta$ uniformly for $x\in\bar{\Omega}$. This implies that (i) holds.

(ii) From system (2.11), we can obtain $I_w$ and $I_h$, which satisfy

According to the maximum principle, assume some $t^{*}\geq 0$ exists such that $I_{w}(\cdot, t^{*}, \varphi)\not\equiv 0$ and $I_{h}(\cdot, t^{*}, \varphi)\not\equiv 0$, and then we have

This implies that (ii) holds. □

Theorem 5.2. The disease is uniformly persistent when $R_{0} > 1$; that is, $\varrho > 0$ exists such that for any initial value function $\varphi\in \mathcal{C}^{+}$ with $\varphi_3(\cdot, 0)\not\equiv 0$ and $\varphi_5(\cdot, 0)\not\equiv 0$, we have

uniformly for $x\in\bar{\Omega}$.

Proof. Let $(S_{w}(x, t, \varphi), E_{w}(x, t, \varphi), I_{w}(x, t, \varphi), S_{h}(x, t, \varphi), I_{h}(x, t, \varphi))$ be the solution of system (2.11) with an initial value function $\varphi\in \mathcal{C}^{+}$. Set

From Lemma 5.1, if $\varphi \in\partial\mathbb{S}_0$, then $I_{w}(x, t, \varphi) > 0$ and $I_{h}(x, t, \varphi) > 0$ for $x\in\bar{\Omega}$, $t > 0$. So, we have $\Upsilon(t)\mathbb{S}_0\subset \mathbb{S}_0.$ Denote $S_{\partial} = \{ \varphi \in\partial\mathbb{S}_0: \Upsilon(t)\varphi \in\partial\mathbb{S}_0, t\geq 0\}. \; \omega(\varphi)$ as the omega limit set of the forward orbit $\gamma^+(\varphi): = \{\Upsilon(t)\varphi:t\geq 0\}.$ We divide the proof into the following two claims.

Claim 1. For any $\varphi\in S_{\partial}$, the omega limit set $\omega(\varphi) = E_0(x)$.

For any $\varphi\in S_{\partial}$, we have $\Upsilon(t)\varphi \in\partial\mathbb{S}_0$. Hence, either $I_{w}(\cdot, t, \varphi)\equiv0$ or $I_{h}(\cdot, t, \varphi)\equiv0$ for all $t\geq 0.$ If $I_{w}(\cdot, t, \varphi)\equiv0$ for all $t\geq 0$, then, from the third equation of (2.11), we have $I_{h}(\cdot, t, \varphi)\equiv0$ for all $t\geq 0$. Then $E_{w}(\cdot, t, \varphi)$ satisfies

Thus, $\lim\limits_{t\rightarrow \infty }E_{w}(x, t, \varphi) = 0$, uniformly for $x\in \bar{\Omega}$. Then, $S_{w}(\cdot, t, \varphi)$ is asymptotic to the linear system (4.1). So, $\lim\limits_{t\rightarrow \infty }S_{w}(x, t, \varphi) = W^*(x)$ uniformly for $x\in \bar{\Omega}$. It is easy to know that $S_{h}(\cdot, t, \varphi)$ satisfies system (2.5). Then, $\lim\limits_{t\rightarrow \infty }S_{h}(x, t, \varphi) = N_h^*(x)$ uniformly for $x\in \bar{\Omega}$. In words,

Assume $I_{w}(\cdot, t_3, \varphi)\not \equiv 0$ for some $t_3 > 0$. Then $I_{h}(\cdot, t, \varphi)\equiv0$ for $t > t_3.$ From Lemma 5.1, we can obtain $I_{w}(\cdot, t, \varphi) > 0$ for $t > t_3$. Since $I_{h}(\cdot, t, \varphi)\equiv0$ for $t > t_3$, we have $I_{w}(\cdot, t, \varphi) \equiv 0$ for $t > t_3$ from the fifth equation of (2.11). This contradicts our assumption. Therefore, $\omega(\varphi) = E_0(x)$ for any $\varphi\in S_{\partial}$.

Claim 2. $E_{0}(x)$ is a uniform weak repeller for $\mathbb{S}_0$, i.e.,

Since $R_0 > 1$, from Lemma 4.2, $\lambda_0(W^*) > 0$. For any given $\sigma\in (0, \min\{\underline{W}^*, \underline{N_h}^*\}]$, set $\lambda_0(W^*, \sigma)$ be the principal eigenvalue of the nonlocal elliptic eigenvalue problem as follows:

We can then obtain $\lim\limits_{\sigma\rightarrow 0^+}\lambda_0(W^*, \sigma) = \lambda_0(W^*).$ So, we can find some sufficiently small constant $\sigma^*\in (0, \min\{\underline{W}^*, \underline{N_h}^*\}]$ such that $\lambda_0(W^*, \sigma^*) > 0.$

If (5.8) is not true, then $\varphi^* \in \mathbb{S}_0$ exists such that $\limsup\limits_{t\rightarrow +\infty}\|\Upsilon(t)\varphi- E_{0}(\cdot)\| < \sigma^*.$ Then some $t_4 > 0$ exists such that $\|\Upsilon(t)\varphi- E_{0}(\cdot)\| < \sigma^*$ for any $t\geq t_4.$ So, $S_{w}(\cdot, t, \varphi^*) > W^*(x)-\sigma^*$, $S_{h}(\cdot, t, \varphi^*) > N_h^*(x)-\sigma^*$, $0 < E_{w}(\cdot, t, \varphi^*) < \sigma^*$, $0 < I_{w}(\cdot, t, \varphi^*) < \sigma^*$, and $0 < I_{h}(\cdot, t, \varphi^*) < \sigma^*$ for any $t\geq t_4.$ Then we have

for all $t \geq t_5: = t_4+\tau.$

Set $\psi = (\psi_1, \psi_2)$ be the positive eigenfunction associated with the principal eigenvalue $\tilde{\lambda}_0(W^*, \tau, \sigma^*)$ of the following eigenvalue problem:

From Lemma 4.1, $\tilde{\lambda}_0(W^*, \tau, \sigma^*)$ and $\lambda_0(W^*, \sigma^*)$ have the same sign. So, $\tilde{\lambda}_0(W^*, \tau, \sigma^*) > 0.$ Consider the following the linear system:

for $t\geq t_5.$ It is easy to see that system (5.12) has a solution $(\nu_1(x, t), \nu_2(x, t)) = e^{\tilde{\lambda}_0(W^*, \tau, \sigma^*)t}((\psi_1(x), \psi_2(x))).$ Together (5.10) with the comparison principle, we can see that a sufficiently small positive number $\mathcal{L}$ exists such that

Since $\tilde{\lambda}_0(W^*, \tau, \sigma^*) > 0,$ we get

which is a contradiction. Therefore, (5.8) is true.

So, $E_0(\cdot)$ is an isolated invariant set $\mathcal{C}^+$, and $W^S(E_0(\cdot))\cap \mathbb{S}_0 = \emptyset,$ in which $W^S(E_0(\cdot))$ is a stable set of $E_0(\cdot)$. Define a continuous function $\mathfrak{f}: \mathcal{C}^{+}\rightarrow \mathbb{R}^+$ with the following form:

Clearly, $\mathfrak{f}^{-1}(0, +\infty) \subset \mathbb{S}_0$. From Lemma 5.1 (ii), we can see that if $\mathfrak{f}(\varphi) = 0$ for $\varphi \in \mathbb{S}_0$, or $\mathfrak{f}(\varphi) > 0$, then $\mathfrak{f}(\Upsilon(t)(\varphi)) > 0$ for any $t > 0$. So, $\mathfrak{f}$ is a generalized distance function for the semiflow $\Upsilon(t) :\mathcal{C}^+\rightarrow \mathcal{C}^+$. By Theorem 3 in ^[24], we can see that $\varrho_0 > 0$ exists such that $\min \limits_{\vartheta\in \omega(\varphi)} \mathfrak{f}(\vartheta)\geq \varrho_0$ for any $\varphi \in \mathbb{S}_0$. This implies uniform persistence. □

6.   Numerical simulations

In this section, numerical simulations are conducted to validate the analytical outcomes. By ^[25], we choose $H_*(x) = 53 (\text{km}^2)^{-1}$, the periodicity is set to $T = 12$ months, and $\Omega = (0, \pi).$ The parameter values are as follows: $N_h^* = 50, \; \mu_h = 0.00157$, $\Lambda = N_h^*\times\mu_h$, $r_h = 0.08, \; c = 0.000001, \; \rho = 50, \; \mu_a = 9, \; \mu_w = 2, \; \tau_1 = 0.5605, \; \tau_w = 17.25/30.4, \; \tau_h = 15/30.4, \; d_w = 0.0125, \; d_h = 0.1,$

6.1.   Dynamic behavior of model (2.11)

We set the initial functions as follows:

Choose $b = 5.$ Numerical simulation determines that $R_0\approx0.2844 < 1.$ Based on Theorem 5.1, $E_{0}(x)$ is globally attractive. As shown in Figure 1, the disease will eventually vanish. In other hand, we set $b = 18.$ The basic reproduction number can be calculated to be $R_0\approx1.024 > 1.$ The dynamic behavior of model (2.11) is illustrated in Figure 2, corroborating the theoretical result of Theorem 5.2. This indicates that an outbreak of the disease is likely to occur.

6.2.   The effect of the release ratio of Wolbachia-infected males on $R_0$

Let $q(x) = q^*(1-k\cos(2x)),$ and set $b = 18, \; \beta_{hw}(x) = 0.6, \; \beta_{wh}(x) = 0.2.$ Firstly, fix $k = 0.01$. We investigate the correlation between $q^*$ and $R_0.$ From Figure 3(a), as $q^*$ increases, $R_0$ gradually decreases and crosses the threshold of 1. Next, fix $q^* = 1.7.$ We analyze the influence of spatial heterogeneity, represented by $k$, on $R_0$. Figure 3(b) demonstrates that $R_0$ is an increasing function with respect to $k$. It should be noted that $q(x)$ represents a homogeneous distribution when $k = 0.$ Moreover, set $q = 1.5$, we explore the effects of the diffusion rates $d_w$ and $d_h$ on $R_0$ as shown in Figure 4.

6.3.   Sensitivity analysis

In this part, we explore the situation where the parameters of model (2.11) are constant. By ^[26], the basic reproduction number $R_0$ can be calculated as

Figure 5 displays scatter plots that illustrate the correlation between the parameters and the basic reproduction number $R_0$. Figure 6 displays a bar chart illustrating the partial rank correlation coefficient (PRCC) values for different parameters in relation to $R_0$. The bars represent the correlation coefficients for each parameter, with some bars above and some below the zero line, indicating positive and negative correlations, respectively. The parameters $\beta_{hw}$ and $\beta_{wh}$ show the strong positive correlations, while $q$, $r_h$ and $\mu_ w$ show the most negative correlations. The parameter $b$ shows a moderate positive correlation with $R_0$. Additionally, some parameters $\tau_1$, $\tau_w$, and $\mu_h$ show moderate negative correlations. The variations in the remaining parameters have a moderate impact on $R_0.$ The contour diagrams in Figure 7 illustrate the joint effect of two parameters on the basic reproduction number $R_0$.

7.   Conclusions

This paper has presented a comprehensive mathematical model for the spatial spread of malaria, incorporating the effects of nonlocal delays and the release of Wolbachia-infected male mosquitoes as a control strategy. The model's dynamics were primarily governed by the basic reproduction number $R_0$, which served as a threshold parameter. Through rigorous theoretical analysis, we established that when $R_0 < 1$, the infection-free steady state is globally attractive, indicating that malaria can be effectively controlled under this condition. Conversely, when $R_0 > 1$, the disease persists uniformly, suggesting that malaria will persist in the population. The simulations confirmed the theoretical predictions and provided insights into the impact of the release ratio of Wolbachia-infected males on the transmission and control of malaria. Future research could further explore the integration of Wolbachia release with other vector control methods to optimize malaria elimination efforts.

Use of AI tools declaration

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

Acknowledgements

The work is partially supported by the National Natural Science Foundation of China (No. 12201007), the Natural Science Foundation of Shandong Province (ZR2024QA021), the Natural Science Research Project of Anhui Educational Committee (No. 2022AH050961, 2023AH030021), the Teaching Research Project of Anhui Polytechnic University (No. 2022jyxm67), the Science and Technology Plan Program of Wuhu Municipal (No. 2024kj016).

Conflict of interest

The authors declare that there are no conflicts of interest.