The effect of migration on transmission of Wolbachia in Nilaparvata lugens

1.   Introduction

Brown planthopper Nilaparvata lugens (Stål) (Hereinafter referred to as N.lugens) is a monophagous insect species, which causes irreversible damage to rice mainly by either piercing or sucking the phloem juice of rice and transmitting Rice Ragged Stunt Viruses (RRSV) ^[1]; the frequent occurrence of N.lugens seriously affects food security worldwide. A number of studies have been performed to develop technology that controls agricultural pests by utilizing Wolbachia ^{[2,3,4,5,6,7,8,9,10]}. Researchers have identified an endosymbiotic Wolbachia strain, wStri, from a small planthopper Laodelphax striatellus (Fallén) that induces cytoplasmic incompatible (CI) phenotypes in its host ^[3,4,5]. CI is a phenomenon of either embryo death or offspring sterility when the male insects infected with Wolbachia mate with females that are either uninfected or infected with a different strain of Wolbachia ^{[5,6,7,8,9,10]}. A study from 2020 reported that Wolbachia wStri can cause cytoplasmic incompatibility in N.lugens and can significantly inhibit the infection and transmission of rice ragged stunt viruses (RRSV) in a laboratory environment ^[2]. This is a breakthrough achievement in the development of a control mechanism of rice planthopper by applying endosymbiotic bacteria. The artificial cultivation of the N.lugens line showed the required characteristics of either population replacement or population suppression in the laboratory.

Because N.lugens is a rice pest with the habit of long-distance migration ^[11], and the total number of migrating N.lugens is often in the hundreds of millions, they have a huge impact on the agricultural ecosystem. Over the years, the migration rule of N.lugens in China has been discovered mainly through modern monitoring technology and international cooperation research: under suitable conditions, N.lugens migrate from the southern rice area in spring, then fly northward alongside monsoons to the Huaibei rice area, and migrate back to the south in autumn ^{[12,13,14,15,16,17]}. The whole migration process spans across a large area in a large amount of time, which is one of the main causes for the large-scale N.lugens disaster. The existing dynamics models of N.lugens generally based the population migration on the diffusion mode between two patches ^{[18,19,20,21]}. In 2019, Vattikuti et al. established linear and nonlinear differential equation models in which the growth rate of N.lugens is described as a function of temperature, and further estimated the bioclimatological threshold under the influence of temperature ^[19]. Additionally, Thuy et al. proposed an ordinary differential equation model of N.lugens where they rapidly disperse in two patches with a stage structure ^[20]. In 2022, Nguyen et al. constructed a herbivore (N.lugens) -host plant (rice) model of N.lugens dispersal between two adjacent patches, and studied the effect of monsoons on the interaction between N.lugens and rice ^[21]. However, the seasonality of N.lugens migration was not fully reflected in the aforementioned models, and more importantly, these models did not consider the spread of Wolbachia in N.lugens. Although we previously studied the transmission of Wolbachia in N.lugens ^[5,6], we did not consider their migration.

In addition, there are some literatures that use the reaction diffusion equation model to study the transmission dynamics of Wolbachia in a mosquito population ^{[22,23,24,25]}. Ref.^[22,23] studied the spreading dynamics of Wolbachia in mosquitoes that activate and diffuse within a relatively isolated and limited region. The diffusion problem of mosquitoes within a free boundary region was also discussed in Ref.^[24]. Liu et al. proposed a coupled model to evaluate the efficiency of Wolbachia on a mosquito control ^[25]. However, the conception of the migration of N.lugens is different from that of the mosquitoes diffusion. The main characteristics of N.lugens migration are seasonality and the simultaneously uniform tropism (i.e., N.lugens migration in one direction is influenced by the unidirectional carrier air flow. In this case, the small random diffusion effect in other directions is negligible, where the carrier air flow is a southeast wind in summer and a northeast wind in autumn ^[26]). Therefore, N.lugens migration can not be studied simply by a research method based on diffusion behavior. In this article, we define a migration coefficient based on the population proportion to construct a Wolbachia transmission model that reflects the N.lugens migration. Through in-depth research on the basic dynamic characteristics of the model, we will analyze how to limit the control effect of Wolbachia during their migration on brown planthoppers.

This paper is organized as follows. In Section 2, the Wolbachia spreading dynamic model with the migration of N.lugens is introduced. In Sections 3 and 4, we analyze the existence and stability of equilibria for the main model. The global stabilities of three subsystems are studied in Section 5. Then, we provide numerical simulations for the stability of equilibria in Section 6. A discussion for the impact of migration on the dynamic properties of the system and the control effects of N.lugens are included in section 7. Finally, we give our conclusions in section 8.

2.   Model

2.1.   A general diffusion model with two patches

In order to establish the migration model of N.lugens, we give the following assumptions:

● Based on the data of the insect trap catches under lights and the data of field surveys of insects in Ref.^[11], we assume that N.lugens migrates between two areas, where area 1 is referred to as the main source area of the northward migrants of N.lugens, and area 2 is referred to as the main destination area of the northward migrants of N.lugens.

● N.lugens infected with wStri is only released in area 1.

● N.lugens mainly migrates northward from area 1 to area 2 in the spring, and it mainly migrates southward from area 2 to area 1 in the autumn.

Let $I^{(i)}(t)$ and $U^{(i)}(t)$ represent the number of the wStri-infected N.lugens and the uninfected N.lugens in time $t$ in the area $i$, respectively. Moreover, the total number of N.lugens within area 1 and area 2 are denoted as $N^{(i)} = I^{(i)}(t)+U^{(i)}(t)$, respectively, $i = 1, 2$.

When population migration is not considered, the N.lugens population only interacts within each patch. Thus, based on the modeling method in ^[5,6], it is easy to obtain the following interaction dynamic model of the N.lugens population in each patch:

where N.lugens infected with wStri has a perfect maternal transmission characteristic, and wStri can induce higher CI levels for the uninfected N.lugens female ^[2,5] and the CI intensity of the male N.lugens infected with wStri $I_M$ against the uninfected female N.lugens $U_F$ be $\xi\in (0, 1)$. According to a two-sex model (1) in Ref^[6], $b$ denotes the half natural birth rate of N.lugens. Note that N.lugens infected with wLug is not considered here. $d_1$ and $d_2$ denote the decay rate constants of the individuals infected with wStri and the uninfected, respectively.

When population migration is taken into the aforementioned system (2.1), it is transfered to the following general discrete diffusion system (2.2):

where $n_{ij}$ and $m_{ij}$ stand for the average migration rate of the wStri-infected N.lugens and the uninfected N.lugens from area $i$ to area $j$, respectively, where $i, j = 1, 2$. It must be specified that all parameters of model (2.2) have positive values. Notice that the average migration rates of N.lugens are all constant, and system (2.2) can not fully reflect the characteristics of N.lugens migration affected by a monsoon. Thus, we need to develop a N.lugens migration model that better reflects the effects of the season.

2.2.   Migration model

The diffusion coefficient is generally related to population density of the two patches. Many predator-prey models with diffusion assume that the diffusion coefficient is a linear function of the diffused species density; some of these models assume that only one species has this dispersal property ^[27], while others assume that there is a dispersal effect in only one monsoon direction ^[21]. The system we consider is markedly different from these models. First, two N.lugens species do not display a predator-prey relationship. Second, N.lugens migration is different from population diffusion, in which the directions of migration are different at different seasons. Simply considering N.lugens migration in one direction does not reflect the bidirectional characteristics of migration. Moreover, the migration coefficient is not a linear function of the population density of the initial source area because the linear function can not reflect the phased one-way migration phenomenon. Therefore, it is key to construct a migration coefficient of N.lugens which can reflect the influence of the season.

Since the migration of N.lugens is influenced by the season, its migration coefficient should reflect the seasonal factor. In the spring, N.lugens migrates northward due to a southeast wind from area 1 to area 2. At this time, the number of N.lugens in area 1 is much greater than in area 2. When N.lugens begins to migrate back due to a northeast wind from the northern rice area 2 to the southern rice area 1 in autumn, the number of N.lugens in area 2 is much greater than in area 1. We use the ratio $\frac{N^{(i)}}{N^{(j)}}$, which represent the number of N.lugens in two areas as a diffusion coefficient, to accurately reflect the seasonal migration characteristics of N.lugens ($i, j = 1, 2$). We call this coefficient the ratio dependent migration coefficient. Substituting $n_{ij}$ and $m_{ij}$ with $\frac{N^{(i)}}{N^{(j)}}$ in system (2.2), we obtain the following migration system (2.3):

with the initial value

where $a_{12}$ ($a_{21}$) and $\bar{a}_{12}$ ($\bar{a}_{21}$) denote the successful landing rate of the wStri-infected N.lugens and the uninfected N.lugens after a long-distance northward (or southward) migration, respectively, and $b_{12}$ ($b_{21}$) and $\bar{b}_{12}$ ($\bar{b}_{21}$) denote the successful take-off rate of the wStri-infected N.lugens and the uninfected N.lugens as northward (or southward) migration, respectively. $a_{12}$, $\bar{a}_{12}$, $a_{21}$, $\bar{a}_{21}$, $b_{12}$, $\bar{b}_{12}$, $b_{21}$, $\bar{b}_{21}\in [0, 1]$. Clearly, the successful take-off rate is more than the successful landing rate (i.e., $b_{12} > a_{12}, b_{21} > a_{21}, \bar{b}_{12} > \bar{a}_{12}$ and $\bar{b}_{21} > \bar{a}_{21}$). We will analyze the dynamics of system (2.3) with the initial value condition (2.4) under the aforementioned basic restriction conditions.

System (2.3) can be regarded as the phased migration model in the first and second half of the year. In the spring, the number of N.lugens in area 1 is much greater than in area 2. The major part of the last two expressions reflecting the migration in the first two equations of system (2.3) is the minus term, and the major part of the last two expressions reflecting the migration in the last two equations of system (2.3) is the plus term, which reveals the reality that there is a significant northward migration of N.lugens in spring, but not an obvious return to the south. On the contrary, the number of N.lugens in area 2 is much greater than in area 1 in autumn; this shows that N.lugens has a significant southward migration, but not an obvious northward migration. It essentially reflects the bidirectional property of N.lugens migration. Additionally, system (2.3) also reflects a simultaneously uniform tropism. Therefore, except for the migration coefficient being a function of population density, it is multiplied by the population density of the corresponding type of N.lugens, which reflects that the simultaneous uniform tropism effect increases when the population density increases during the same season.

Denote the right function of system (2.3) by $\textbf{F}(I^{(1)}, U^{(1)}, I^{(2)}, U^{(2)}) = (f_1, f_2, f_3, f_4)^T$. Since

the right-hand side function $f_i$ of system (2.3) can be continuously extended to the origin $(0, 0, 0, 0)$ by $f_i(0, 0, 0, 0) = 0$ ($i = 1, 2, 3, 4$). Thus, system (2.3) is well-defined in the following region:

3.   The nonnegativity and boundedness of solutions for system (2.3)

In this section, we will analyze the nonnegativity and boundedness of solutions for system (2.3).

Theorem 1. Suppose that $(I^{(1)}(t), U^{(1)}(t), I^{(2)}(t), U^{(2)}(t))$ is any solution of system (2.3) with the initial condition (2.4). Then, it is nonnegative and bounded.

Proof. Nonnegativity. When $I^{(1)} = 0$, $U^{(1)}\geq 0$, $I^{(2)}\geq 0$ and $U^{(2)}\geq 0$, we have the following:

Moreover, when $I^{(1)}\geq 0$, $U^{(1)} = 0$, $I^{(2)}\geq 0$ and $U^{(2)}\geq 0$, we have the following:

Similarly, when $I^{(1)}\geq 0$, $U^{(1)}\geq 0$, $I^{(2)} = 0$ and $U^{(2)}\geq 0$, we have the following:

when $I^{(1)}\geq 0$, $U^{(1)}\geq 0$, $I^{(2)}\geq 0$ and $U^{(2)} = 0$, we have

According to proposition B.7 in ^[28], it follows that

Thus, the solution of system (2.3) with the initial condition (2.4) is nonnegative.

Boundedness. Let $W(t) = I^{(1)}(t)+U^{(1)}(t)+I^{(2)}(t)+U^{(2)}(t)$. From the nonnegativity of solution which has been proven, and the basic constraints of parameters $b_{12} > a_{12}, b_{21} > a_{21}, \bar{b}_{12} > \bar{a}_{12}$ and $\bar{b}_{21} > \bar{a}_{21}$, we have the following:

Then, we have the following:

Therefore, we obtain the following:

namely, $\frac{{\rm d}W}{{\rm d}t}+bW\leq L$.

It follows from the differential inequality theory that $0\leq W(t) \leq \frac{L}{b}(1-e^{-bt})+W(0)e^{-bt}$. Therefore, we get $0\leq W(t)\leq L/b$ as $t\to +\infty$.

Consequently the solution $(I^{(1)}(t), U^{(1)}(t), I^{(2)}(t), U^{(2)}(t))$ starting from region $R_4^+$ is restricted to the following region:

where $L = \frac{2b^2}{d_1}+\frac{b^2(3-\xi)^2}{2d_2^2}$. □

The aforementioned region $\Lambda$ is the positive invariant set with respect to system (2.3). The nonnegative solution for system (2.3) shows that the population density of N.lugens can not be negative, and the boundedness of the solution for system (2.3) indicates that the population density of N.lugens in the wild can not increase indefinitely, which are consistent with reality.

4.   The existence and the stability of equilibria

4.1.   The existence of equilibria

Now we discuss the existence of equilibria of system (2.3). It is easy to see that system (2.3) always has a trivial equilibrium point $E_0(0, 0, 0, 0)$. In addition, system (2.3) may exist as equilibrium infected with wStri $E_1(\bar{I}^{(1)}, 0, \bar{I}^{(2)}, 0)$, the uninfected equilibrium $E_2(0, \hat{U}^{(1)}, 0, \hat{U}^{(2)})$ and the positive equilibrium $E_3(\tilde{I}^{(1)}, \tilde{U}^{(1)}, \tilde{I}^{(2)}, \tilde{U}^{(2)})$.

1). The existence analysis of the equilibrium infected with wStri $E_1$. Equilibrium $E_1(\bar{I}^{(1)}, 0, \bar{I}^{(2)}, 0)$ should satisfy the following algebraic equations:

where $\bar{I}^{(i)} > 0$, $i = 1, 2$. Simplifying the above equations gives two equivalent equations:

Transferring the terms $-d_1\bar{I}^{(1)}$ and $-d_1\bar{I}^{(2)}$ to the right-hand side of (4.1) and dividing the first equation with the second equation, we obtain the following:

Let $\frac{\bar{I}^{(1)}}{\bar{I}^{(2)}} = p$. By substituting it into formula (4.2), we obtain the following:

Simplifying and rearranging formula (4.3) produces a polynomial equation with respect to $p$:

Since $a_{12}$, $b+b_{12}$, $b+b_{21}$ and $a_{21}$ are greater than $0$, by ignoring those terms with the zero coefficient, the number of times of the coefficient signs changes for Eq (4.4) is $1$. According to Descartes' Rule of Signs, we assert that Eq (4.4) has only one positive real root, which is denoted as $c_1$. As for either the remaining complex roots or the non-positive real roots, if they exist, they do not meet the requirements.

Thus, by substituting $\frac{\bar{I}^{(1)}}{\bar{I}^{(2)}} = c_1$ into Eq (4.1), and noticing that $\bar{I}^{(1)} = c_1\bar{I}^{(2)}$, we have

or

Note that $\bar{I}^{(1)}$ in formula (4.5) is equivalent to that of formula (4.6), and so is also $\bar{I}^{(2)}$. Therefore, we obtain the existence of equilibrium $E_1$ as follows.

Theorem 2. If $a_{21}+bc_1^2-b_{12}c_1^3 > 0$ or $a_{12}c_1^3+bc_1-b_{21} > 0$, then system (2.3) exists as one equilibrium infected with $wStri \; E_1(\bar{I}^{(1)}, 0, \bar{I}^{(2)}, 0)$, where $\bar{I}^{(1)}$ and $\bar{I}^{(2)}$ are shown in formula (4.5) or (4.6), and $c_1$ is the positive real root of Eq (4.4).

2). The existence analysis of equilibrium $E_2(0, \hat{U}^{(1)}, 0, \hat{U}^{(2)})$. Similarly, equilibrium $E_2$ should satisfy the following equations:

where $\hat{U}^{(i)} > 0$, $i = 1, 2$.

Simplifying the equation above, we obtain that

Thus, we get

Let $\frac{{\hat{U}}^{(1)}}{{\hat{U}}^{(2)}} = q$. By substituting it into the formula above, the following polynomial equation can be easily obtained via simple algebraic operation:

Because of $\bar{a}_{12}$, $b+\bar{b}_{12}$, $b+\bar{b}_{21}$ and $\bar{a}_{21}$ are all greater than $0$, it follows from Descartes' Rule of Signs that Eq (4.8) has only a positive real root, which is denoted as $c_2$. As for either the remaining complex roots or non-positive real roots, if they exist, they do not meet the requirements.

Substituting $\frac{{\hat{U}}^{(1)}}{{\hat{U}}^{(2)}}$ with $c_2$ in formula (4.7), we have

or

It is easy to see that formulas (4.9) and (4.10) are actually equivalent from the calculation process above. Thus, we obtain the existence result of equilibrium $E_2$ as follows.

Theorem 3. If $\bar{a}_{21}+bc_2^2-\bar{b}_{12}c_2^3 > 0$ or $\bar{a}_{12}c_2^3+bc_2-\bar{b}_{21} > 0$, then system (2.3) has an uninfected equilibrium $E_2(0, \hat{U}^{(1)}, 0, \hat{U}^{(2)})$, where $\hat{U}^{(1)}$ and $\hat{U}^{(2)}$ are shown in formula (4.9) or (4.10), and $c_2$ satisfies Eq (4.8).

4.2.   The stability of equilibrium

In this subsection, we discuss the stability of equilibrium for system (2.3). For the sake of convenience, let symbols $E(I^{(1)}, U^{(1)}, I^{(2)}, U^{(2)})$ be any equilibrium of system (2.3), and let $J_E = (f_{ij})_{4\times 4}$ be the Jacobi matrix corresponding system (2.3) at equilibrium $E$, where

First, we consider the stability of system (2.3) at equilibrium $E_1$. The Jacobi matrix $J_{E_{1}}$ of the corresponding linearized system of system (2.3) at equilibrium $E_1$ is as follows:

Thus, the characteristic equation of $J_{E_{1}}$ is as follows:

where

Since $A_3 > 0$, to make all the roots of the characteristic equation (4.11) have the negative real part, the following condition (H1) must be satisfied:

(H1) $A_1 > 0, A_2 > 0$ and $A_4 > 0$.

Therefore, we have the stability result of system (2.3) at equilibrium $E_1$ as follows:

Theorem 4. Suppose that the condition (H1) holds; then, the equilibrium infected with wStri $E_1(\bar{I}^{(1)}, 0, \bar{I}^{(2)}, 0)$ of system (2.3) is locally asymptotically stable.

Then, we consider the stability of equilibrium $E_2(0, \hat{U}^{(1)}, 0, \hat{U}^{(2)})$ of system (2.3). The Jacobi matrix $J_{E_{2}}$ of the corresponding linear system of system (2.3) at $E_2$ is as follows:

Therefore, the characteristic equation of $J_{E_{2}}$ is as follows:

where

Due to $B_3 > 0$, the following condition (H2) must be satisfied to make all the roots of the characteristic equation of $J_{E_{2}}$ have a negative real part:

(H2) $B_1 > 0, B_2 > 0$ and $B_4 > 0$.

Therefore, we have the following the stability result of equilibrium $E_2$.

Theorem 5. Suppose that the condition (H2) holds; then, the uninfected equilibrium $E_2(0, \hat{U}^{(1)}, 0, \hat{U}^{(2)})$ of system (2.3) is locally asymptotically stable.

The existence of system (2.3) at the positive equilibrium will be analyzed by a numerical simulation method in Section 6. The existence and stability of $E_1$ indicates that wStri can fully invade the wild population even though there exists migration, while the existence and stability of the positive equilibrium $E_3$ suggest that migration may lead to an incomplete invasion of wStri; however, the existence and stability of $E_2$ suggest that migration may lead to a failure of the control method based on Wolbachia.

5.   The global stability of the degenerated system for system (2.3)

5.1.   The subsystem without N.lugens migration

N.lugens is a wing dimorphic insect with long and short wings. The long-winged morph has a migration habit, which allows them to escape adverse habitats and track changing resources ^[29]. The short-winged morph lacks a long distance migration ability, though its strong fecundity is conducive to the rapid proliferation of the population. Therefore, in order to theoretically reveal the interaction of the unmigrated N.lugens population, we need to study the interaction of the subsystem without migration.

When N.lugens migration is not considered, system (2.3) is degenerated into the following system (5.1):

Since the first two equations are independent of the second two equations in system (5.1), we study the following system of the interaction between Wolbachia and N.lugens when N.lugens infected with wStri is only released in area 1:

Obviously, system (5.2) always has equilibria $E_{00}(0, 0)$, $E_{10}(\frac{b}{d_1}, 0)$ and $E_{01}(0, \frac{b}{d_2})$. A straightforward calculation shows that system (5.2) has a unique positive equilibrium $E_{11}(\check{I}^{(1)}, \check{U}^{(1)})$ if and only if $0 < 1-\frac{d_2}{d_1} < \xi$ is true, where

By a direct calculation, the Jacobi matrix $J_{E_{00}}$ corresponding to system (5.2) at the $E_{00}$ is as follows:

The Jacobi matrix $J_{E_{10}}$ and $J_{E_{01}}$ corresponding to system (5.2) at the $E_{10}$ and $E_{01}$ are

and

respectively. Moreover, the Jacobi matrix of system (5.2) corresponding to $E_{11}$ is as follows:

and its determinant ${\rm det} J_{E_{11}} < 0$ from the existence condition of $E_{11}$. Therefore, the stability results of those equilibria for system (5.2) are obtained as follows:

Lemma 1. $E_{00}(0, 0)$ is always unstable. $E_{10}(\frac{b}{d_1}, 0)$ is locally asymptotically stable if $1-\xi < \frac{d_2}{d_1}$. $E_{01}(0, \frac{b}{d_2})$ is locally asymptotically stable if $d_2 < d_1$. $E_{11}(\check{I}^{(1)}, \check{U}^{(1)})$ is always unstable if it exists.

Theorem 6. If $\frac{d_2}{d_1} > 2-\xi$, then, equilibrium infected with wStri $E_{10}$ of system (5.2) is globally asymptotically stable in region $\{(I^{(1)}, U^{(1)}): I^{(1)} > 0, U^{(1)}\geq 0\}$.

Proof. Let $(I^{(1)}(t), U^{(1)}(t))$ be the trajectory of system (5.2) with the initial value $(I_0^{(1)}, U_0^{(1)})$. It follows from $\frac{d_2}{d_1} > 2-\xi$ that $\frac{d_2}{d_1} > 1$ and $\frac{d_2}{d_1} > 1-\xi$; therefore, the unstability of $E_{01}(0, \frac{b}{d_2})$, the nonexistence of the positive equilibrium $E_{11}$ and the local asymptotic stability of $E_{10}(\frac{b}{d_1}, 0)$ are obtained, where $E_{10}$ is a saddle point and its stable manifold is $U^{(1)}$-axis.

If $I_0^{(1)} > 0$ and $U_0^{(1)} = 0$, then for any $t > 0$, $U^{(1)}(t) = 0$ and $\frac{{\rm d}I^{(1)}(t)}{{\rm d}t} = bI^{(1)}-d_1(I^{(1)})^2$. Therefore, the system's trajectories will tend to $E_{10}$.

If $I_0^{(1)} > 0$ and $U_0^{(1)} > 0$, then we have $I^{(1)}(t) > 0$ and $U^{(1)}(t) > 0$ for all $t > 0$. To apply the Poincaré's method of tangential curves ^[30], we introduce the following function:

where $\alpha_1 > 0$, $\alpha_2 > 0$ and $\frac{\alpha_1}{\alpha_2} = \frac{d_1}{d_2}$. Taking the derivative of $H(t)$ with respect to the solution to system (5.2), we have the following:

Hence, $H(t)$ strictly increases with respect to $t$, and it follows that system (5.2) does not have any nontrivial periodic solution in region $\{(I^{(1)}, U^{(1)}): I^{(1)} > 0, U^{(1)}\geq 0\}$. Thus, the system's trajectories from the initial value $I_0^{(1)} > 0$ and $U_0^{(1)} > 0$ will be attracted to $E_{10}$.

Therefore, $E_{10}$ is globally asymptotically stable in region $\{(I^{(1)}, U^{(1)}): I^{(1)} > 0, U^{(1)}\geq 0\}$. □

For example, the parameters are taken as $b = 0.45$, $d_1 = 0.1$, $d_2 = 0.2$ and $\xi = 0.675$ because of $1-\frac{d_1}{d_2}(2-\xi) = 0.3375 > 0$, and the condition of Theorem 6 is satisfied. Thus, the equilibrium infected with wStri $E_{10} = (4.5, 0)$ is globally asymptotically stable. In such case, there exists no positive equilibrium, and the uninfected equilibrium $E_{01} = (0, 2.25)$ is a saddle point of system (5.2), whose stable manifold is $U^{(1)}$-axis. The numerical verification of Theorem 6 is shown in Figure 1.

At the end of this subsection, we provide the existence and stability of equilibria for system (5.1). There are always three equilibria $E^w_{0}(0, 0, 0, 0)$, $E^w_1(b/d_1, 0, b/d_1, 0)$, $E^w_2(0, b/d_2, 0, b/d_2)$, and there exists a positive equilibrium $E^w_3(\check{I}^{(1)}, \check{U}^{(1)}, \check{I}^{(2)}, \check{U}^{(2)})$ if and only if the condition $1-\xi < \frac{d_2}{d_1} < 1$ holds. For the stability of system (5.1), we have the following conclusion.

Theorem 7. $E^w_{0}(0, 0, 0, 0)$ is always unstable. $E^w_1$ is locally asymptotically stable if $1-\xi < \frac{d_2}{d_1}$. $E^w_2$ is locally asymptotically stable if $d_2 < d_1$. $E^w_3$ is always unstable if it exists. In addition, $E^w_1$ and $E^w_2$ are both locally asymptotically stable if $1-\xi < \frac{d_2}{d_1} < 1$.

5.2.   The wild N.lugens migration model

In this subsection, we study the migration dynamics of wild N.lugens with the long-winged morph (i.e., subsystem without N.lugens infected with wStri in system (2.3)) as follows:

System (5.4) reflects the southward and northward migration dynamic of the wild N.lugens. Denote the right hand function of system (5.4) by $(g_1, g_2)^T$, and $g_i (i = 1, 2)$ can be continuously extended to the origin $(0, 0)$ by the supplementary definition $g_1(U^{(1)}, 0) = bU^{(1)}-d_2(U^{(1)})^2$, $g_1(0, U^{(2)}) = 0$, $g_2(U^{(1)}, 0) = 0$, $g_2(0, U^{(2)}) = bU^{(2)}-d_2(U^{(2)})^2$, $g_1(0, 0) = g_1(0, 0) = 0$.

The interior $\{(U^{(1)}, U^{(2)}):U^{(1)} > 0, U^{(2)} > 0\}$ of the first quadrant for $U^{(1)}U^{(2)}$-plane is the positive invariant set of system (5.4), and its trajectory starting from this region is eventually bounded. It is easy to see that system (5.4) has a trivial equilibrium $E_{20}(0, 0)$, which is unstable by a direct calculation. Furthermore, if $\bar{a}_{21}+bc_2^2-\bar{b}_{12}c_2^3 > 0$, then there exists only a positive equilibrium $E_{21}(\hat{U}^{(1)}, \hat{U}^{(2)})$, where $c_2$ is the positive real root of equation (4.8), and $\hat{U}^{(1)}$ and $\hat{U}^{(2)}$ are the same as (4.9).

By a direct calculation, it is not hard to know that the trivial equilibrium $E_{20}(0, 0)$ is always unstable. We provide the global stability results for the positive equilibrium $E_{21}$ as follows. The Jacobi matrix corresponding to system (5.4) at $E_{21}$ is

and its characteristic equation is

where $B_3$ and $B_4$ are consistent with that of Eq (4.12). Therefore, we have the following result of the local stability for the positive equilibrium,

Lemma 2. If $B_4 > 0$, then the positive equilibrium $E_{21}(\hat{U}^{(1)}, \hat{U}^{(2)})$ of system (5.4) is locally asymptotically stable.

In fact, the positive equilibrium $E_{21}$ is not only locally asymptotically stable, but also globally asymptotically stable.

Theorem 8. If $B_4 > 0$, then the positive equilibrium $E_{21}(\hat{U}^{(1)}, \hat{U}^{(2)})$ of system (5.4) is globally asymptotically stable in the region $\{(U^{(1)}, U^{(2)}):U^{(1)} > 0, U^{(2)} > 0\}$.

Proof. We only need to show that system (5.4) does not have any nontrivial periodic solution in $\{(U^{(1)}, U^{(2)}):U^{(1)} > 0, U^{(2)} > 0\}$. According to Dulac's criteria (see, e.g., Theorem 2, Section 3.9 of Ref.^[31]), we choose the Dulac function $B(U^{(1)}, U^{(2)}) = \frac{1}{U^{(1)}U^{(2)}}$, where the divergence of the vector field $(g_1B, g_2B)$ satisfies the following:

for any $(U^{(1)}, U^{(2)})\in \{(U^{(1)}, U^{(2)}):U^{(1)} > 0, U^{(2)} > 0\}$. Thus, system (5.4) does not exist as a closed orbit in the interior of the first quadrant of the $U^{(1)}U^{(2)}$-plane. In addition, the condition $B_4 > 0$ implies that $E_{21}$ is locally asymptotically stable, and there exists no other equilibrium in the interior of the first quadrant for system (5.4). Therefore, the positive equilibrium $E_{21}(\hat{U}^{(1)}, \hat{U}^{(2)})$ is globally asymptotically stable in the region $\{(U^{(1)}, U^{(2)}):U^{(1)} > 0, U^{(2)} > 0\}$. □

We choose the parameters $b = 0.5$, $d_2 = 0.2$, $\xi = 0.675$, $\bar{b}_{12} = 0.75$, $\bar{b}_{21} = 0.85$, $\bar{a}_{12} = 0.7$ and $\bar{a}_{21} = 0.8$; then, we have $c_2 = 1.0405$, $\hat{U}^{(1)} = 2.2933$ and $\hat{U}^{(2)} = 2.2041$. Furthermore, we get $B_4 = 2.3761 > 0$. It follows from Theorem 8 that the positive equilibrium $E_{21} = (2.2933, 2.2041)$ is globally asymptotically stable, and the numerical simulation results as shown in Figure 2.

5.3.   The migration model of N.lugens infected with wStri

In this subsection, we will analyze the migration dynamics of N.lugens infected with wStri with long-winged morph (i.e., subsystem without the uninfected N.lugens in system (2.3)) as follows:

System (5.5) reveals the wStri-infected N.lugens migration dynamic after Wolbachia wStri completely invades the wild N.lugens. The study of system (5.5) is mainly aimed at proposing the artificial control measures for the wStri-infected N.lugens.

Similar to system (5.4), the first quadrant of the $I^{(1)}I^{(2)}$-plane is the positive invariant set of system (5.5). If $a_{21}+bc_1^2-b_{12}c_1^3 > 0$ holds, then system (5.5) only has a positive equilibrium $E_{31}(\bar{I}^{(1)}, \bar{I}^{(2)})$, where $c_1$ is the positive real root of (4.4), and $\bar{I}^{(1)}$ and $\bar{I}^{(2)}$ are consistent with that of formula (4.5) or (4.6).

The Jacobi matrix corresponding to system (5.5) at $E_{31}$ is as follows:

Its characteristic equation is as follows:

where $A_3$ and $A_4$ are completely consistent with that of Eq (4.11).

Lemma 3. If $A_4 > 0$, then the positive equilibrium $E_{31}(\bar{I}^{(1)}, \bar{I}^{(2)})$ of system (5.5) is locally asymptotically stable.

Theorem 9. If $A_4 > 0$, then the positive equilibrium $E_{31}(\bar{I}^{(1)}, \bar{I}^{(2)})$ of system (5.5) is globally asymptotically stable in the region $\{(I^{(1)}, I^{(2)}):I^{(1)} > 0, I^{(2)} > 0\}$.

Proof. We only need to show that there exists no closed orbit in the interior of the first quadrant of the $I^{(1)}I^{(2)}$-plane for system (5.5). Let $(h_1, h_2)^T$ be the right hand function of system (5.5). Similar to the proof of Theorem 8, to apply Dulac's criteria (see, e.g., Theorem 2, Section 3.9 of Ref.^[31]), we take the Dulac function $D(I^{(1)}, I^{(2)}) = \frac{1}{I^{(1)}I^{(2)}}$, for any $I^{(1)} > 0$ and $I^{(2)} > 0$, the divergence of the vector field $(h_1D, h_2D)$ is equal to the following:

According to the Dulac criterion, system (5.5) has no closed orbit in the interior of the first quadrant of the $I^{(1)}I^{(2)}$-plane. Furthermore, it follows from $A_4 > 0$ that the positive equilibrium $E_{31}$ is locally asymptotically stable. Additionally, there is no other equilibrium in the interior of the first quadrant. Therefore, $E_{31}(\bar{I}^{(1)}, \bar{I}^{(2)})$ is globally asymptotically stable in the region $\{(I^{(1)}, I^{(2)}):I^{(1)} > 0, I^{(2)} > 0\}$. □

For example, taking the parameters $b = 0.6$, $\xi = 0.675$, $d_1 = 0.1$, $b_{12} = 0.6$, $b_{21} = 0.85$, $a_{12} = 0.5$ and $a_{21} = 0.7$, we get $c_1 = 1.1106$, $\bar{I}^{(1)} = 5.0121$ and $\bar{I}^{(2)} = 4.5131$; then, we have $A_4 = 2.0410 > 0$. Hence, the positive equilibrium $E_{31} = (5.0121, 4.5131)$ is globally asymptotically stable from Theorem 9. The numerical simulation results are shown in Figure 3.

The result of Theorem 2, 4 and 9 indicate that although N.lugens infected with wStri are only released in area 1, they can coexist in both area 1 and area 2 after wStri completely invades the wild population.

6.   Numerical simulation

We need to estimate some biological parameters of N.lugens for a numerical simulation, such as $b, d_1, d_2$ and $\xi$, which can systematically reflect the biological characteristics of the N.lugens population dynamics and development, and are often used for population dynamics monitoring, prediction and bioassay ^[32]. N.lugens has different natural birth rates in the different rice varieties, ranging from $(0.1 \pm 0.016, 0.6)$ ^[32,33,34]. However, if N.lugens is not affected by any adverse factors, its birth rates should be high, and even close to $1$. The natural mortality of N.lugens is not only associated with artificial factors such as insecticide, but also with the meteorological conditions such as rainfall and temperature. Its daily mortality rates are approximately $3%$ to $5%$, and can reach $10%$ to $30%$ in case of severe weather ^[35,36]. Therefore, the natural mortality rates of wild N.lugens are within the ranges of $(0.03, 0.3)$.

The experimental data in Ref.^[2] indicated that the lifespan of female N.lugens infected with wStri was not shorter than that of the uninfected female N.lugens, and there was no significant difference in male survivorship among those two N.lugens types. Thus, we consider that the mortality rate of N.lugens infected with wStri has roughly the same range of values as that of the uninfected N.lugens. In addition, the CI intensity is $\xi = 0.675$ by ^[2].

When the N.lugens population migrates northward and returns southward, atmospheric circulation and the weather exert a direct influence on its takeoff and landing. In September, the peak of N.lugens infestation appears in the southern rice area under the influence of a low pressure circulation. It is easy to obtain the daily cumulative numbers of N.lugens infestation by using the data of the insect trap catches under lights. However, there is no literature reports on the successful landing rate and takeoff rate of N.lugens after its northward migration.

The above biological and migration parameters of N.lugens are summarized in Table 1.

6.1.   Verifying the existences and stabilities of two boundary equilibria of system (2.3)

We take parameters $b = 0.54$, $d_1 = 0.0015$, $d_2 = 0.0032$, $\xi = 0.675$, $a_{12} = 0.4222$, $a_{21} = 0.5453$, $\bar{a}_{12} = 0.543$, $\bar{a}_{21} = 0.4317$, $b_{12} = 0.568$, $b_{21} = 0.754$, $\bar{b}_{12} = 0.68$ and $\bar{b}_{21} = 0.64$. By solving equation (4.4), we obtain $c_1 = 1.0896$. Furthermore, we have $E_1 = (253.6439, 0,232.7944, 0)$, $A_1 = 2.5339 > 0$, $A_2 = 1.3586 > 0$ and $A_4 = 1.3144 > 0$. Therefore, from Theorem 4, equilibrium infected wStri $E_1$ is locally asymptotically stable, and the numerical simulation is shown in Figure 4(a). When parameters are taken as $b = 0.6438$, $d_1 = 0.0025$, $d_2 = 0.0011$, $\xi = 0.675$, $a_{12} = 0.4$, $a_{21} = 0.5$, $\bar{a}_{12} = 0.54$, $\bar{a}_{21} = 0.42$, $b_{12} = 0.56$, $b_{21} = 0.75$, $\bar{b}_{12} = 0.68$ and $\bar{b}_{21} = 0.65$, we obtain that $c_2 = 0.9602$, and $E_2 = (0,405.7542, 0,422.5527)$. In such a case we have $B_1 = 2.1020 > 0$, $B_2 = 0.8842 > 0$ and $B_4 = 1.681 > 0$. Thus, the condition of Theorem 5 is met, the uninfected equilibrium $E_2$ is locally asymptotically stable, and the numerical simulation result is shown in Figure 4(b).

6.2.   Existence of the positive equilibrium

Now, we numerically analyze the existence of the positive equilibrium $E_3$ of system (2.3). For simplicity of presentation, we still use the notation $J_{E_3} = (f_{ij})_{4\times 4}$ to represent the Jacobian matrix of system (2.3) at $E_3 = (\tilde{I}^{(1)}, \tilde{U}^{(1)}, \tilde{I}^{(2)}, \tilde{U}^{(2)})$, whose characteristic equation is $\lambda^4-C_1\lambda^3+C_2\lambda^2-C_3\lambda+C_4 = 0$. When the parameters are taken as $b = 0.6438$, $d_1 = 0.0025$, $d_2 = 0.0031$, $\xi = 0.675$, $a_{12} = 0.1$, $a_{21} = 0.34$, $\bar{a}_{12} = 0.1593$, $\bar{a}_{21} = 0.4317$, $b_{12} = 0.32$, $b_{21} = 0.64$, $\bar{b}_{12} = 0.2316$ and $\bar{b}_{21} = 0.5639$. System (2.3) has a positive equilibrium $E_3 = (29.4676,139.75, 21.9651,108.2649)$ and

The characteristic equation is $\lambda^4+2.9308\lambda^3+1.7108\lambda^2+0.2311\lambda-0.0172 = 0$, which has characteristic roots $\lambda_1 = -2.2036$, $\lambda_3 = -0.4414$, $\lambda_4 = -0.3382$ and $\lambda_2 = 0.0524$. Therefore, the positive equilibrium $E_3$ is unstable.

6.3.   Bistable phenomena

Further numerical simulation results indicate that system (2.3) can appear as a bistable phenomenon. For example, we take the parameters $b = 0.55$, $d_1 = 0.2$, $d_2 = 0.1$, $\xi = 0.675$, $a_{12} = 0.75$, $a_{21} = 0.8$, $\bar{a}_{12} = 0.6$, $\bar{a}_{21} = 0.8$, $b_{12} = 0.88$, $b_{21} = 0.93$, $\bar{b}_{12} = 0.85$ and $\bar{b}_{21} = 0.98$. By numerically solving equations (4.4) and (4.8), we get $c_1 = 1.0189$ and $c_2 = 1.0691$. Therefore, system (2.3) exists as an equilibrium infected wStri $E_1 = (2.1193, 0, 2.0799, 0)$. In addition, it is easy to obtain that $A_1 = 1.8903 > 0$, $A_2 = 0.4112 > 0$ and $A_4 = 2.2373 > 0$. From Theorem 4, it follows that $E_1$ is locally asmyptotically stable. On the other hand, system (2.3) also has an uninfected equilibrium $E_2 = (0, 3.4119, 0, 3.1913)$. Furthermore, we can obtain $B_1 = 2.0313 > 0$, $B_2 = 0.4283 > 0$ and $B_4 = 1.6263 > 0$. Thus, from Theorem 5, we know that $E_2$ is also locally asmyptotically stable. Therefore, system (2.3) appears as a bistable phenomenon. By randomly generating the initial values for numerical verification of 10 groups, we obtain the $I^{(1)}-U^{(1)}-I^{(2)}$-diagram for system (2.3), as shown in Figure 5. It is apparent from Figure 5 that system (2.3) has two stable equilibria.

7.   Discussion

7.1.   Impact of migration on dynamical properties of system

When N.lugens migration is not considered, system (5.1) reflects the population evolution of the wStri-infected N.lugens and the wild N.lugens. From Theorem 5.2, we know that equilibrium infected with wStri $E^w_1$ or the uninfected equilibrium $E^w_2$ are asymptotically stable under certain conditions, and $E^w_1$ and $E^w_2$ are simultaneously locally asymptotically stable as $1-\xi < \frac{d_2}{d_1} < 1$. This indicates that system (5.1) can exhibit a bistable phenomenon. In addition, the positive equilibrium of system (5.1) is always unstable, which shows that without considering migration, the two types of N.lugens must eventually not coexist in each patch. On the other hand, when migration is considered, system (2.3) reflects population evolution with N.lugens migrating in both patches. According to numerical simulation results from Subsection 6.2, we know that there exists a positive equilibrium for system (2.3) under certain conditions. Moreover, from subsection 6.3, we know that system (2.3) can exhibit bistable phenomenon.

We further qualitatively summarize the impact of migration parameters on equilibria of system (2.3). According to Theorems 4.1 and 4.3, we know that either decreasing the successful landing rate $\bar{a}_{ij}$ or increasing the take-off rate $\bar{b}_{ij}$ of the uninfected N.lugens is beneficial to the existence and stability of equilibrium infected with wStri $E_{1}$. Additionally, from Theorem 3 and Theorem 5, we know that either decreasing the successful landing rate $a_{ij}$ or increasing the take-off rate $b_{ij}$ of N.lugens infected wStri is conducive to the existence and stability of the uninfected equilibrium $E_{2}$. We will further analyze the controllable migration parameters threshold value for the existence and stability of each equilibrium in Subsection 7.3.

To numerically verifying the above theoretical results, we investigate the effect of migration coefficient $\bar{a}_{12}$ on the stability of the equilibrium infected with wStri $E_1$. The parameters are taken as $b = 0.6438$, $d_1 = 0.0025$, $d_2 = 0.0081$, $\xi = 0.675$, $a_{12} = 0.3$, $a_{21} = 0.2$, $\bar{a}_{21} = 0.82$, $b_{12} = 0.76$, $b_{21} = 0.55$, $\bar{b}_{12} = 0.28$, and $\bar{b}_{21} = 0.15$. Without considering migration, we have $1-\xi = 0.325 < \frac{d_2}{d_1} = 3.24$ and $d_2 > d_1$. It follows from Theorem 7 that $E^w_1$ is locally asymptotically stable, and $E^w_2$ and $E^w_3$ are not stable for system (5.1). When migration is considered, we find that increasing the migration parameter $\bar{a}_{12}$ to a threshold value $0.75$, equilibrium $E_1$ of system (2.3) is not stable, and $E_2$ is locally asymptotically stable, as shown in Figure 6. The numerical simulations show that when migration is not considered, equilibrium infected with wStri $E^w_1$ is stable; with the increase of migration coefficient $\bar{a}_{ij}$, the stability of $E^w_1$ can be changed from stable to unstable, which may lead $E^w_2$ to be stable.

7.2.   Impact of migration on control effects of N.lugens using Wolbachia

First, the influence of migration on the control effect of the wild N.lugens will be numerically analyzed by comparing the control degree of the wild N.lugens:

Notice that the control degree of the wild N.lugens at migration refers to the total control degree of both patches, where the concept of the control degree is consistent with that of the wild N.lugens introduced in Ref.^[6]. The control degree curves are plotted in Figure 7 after computing from system (2.3) (with migration) and system (5.1) (without migration). Here, Figure 7(b) illustrates that when migration is considered, all of the N.lugens infected with wStri are only released in area 1, namely the initial value $(I_0^{(1)}, U_0^{(1)}, I_0^{(2)}, U_0^{(2)}) = (230,530, 0, 20)$. When migration is not considered, the wStri-infected N.lugens is simultaneously released in two areas in proportion to the numbers of the uninfected N.lugens in both areas, so the initial value $(I_0^{(1)}, U_0^{(1)}, I_0^{(2)}, U_0^{(2)}) = (222,530, 8, 20)$. Numerical simulation results show that the control degree of the wild N.lugens with migration may be smaller than that without migration, which indicates that migration may lead to a decrease of the control effect of the wild N.lugens, which is consistent with the general understanding.

It is generally believed that effective control is achieved when the control degree reaches 95% ^[6]. The dotted lines in Figure 7 show that the goal of effective control is not achieved under migration. Therefore, in the presence of migration, the original strategy (i.e., quantity and timing of release) for releasing N.lugens infected with wStri may not achieve the goal of either complete or even effective control of the wild N.lugens.

Notice that the optimal state for controlling N.lugens using Wolbachia is stability of $E_1$. Next, we numerically analyze the effect of migration on this state. Taking $\bar{a}_{12}$ as an example, we fix $b = 0.6438$, $d_1 = 0.0025$, $d_2 = 0.0031$, $\xi = 0.675$, $a_{12} = 0.1222$, $a_{21} = 0.5453$, $\bar{a}_{21} = 0.4317$, $b_{12} = 0.4716$, $b_{21} = 0.7554$, $\bar{b}_{12} = 0.6867$ and $\bar{b}_{21} = 0.6439$. Here, we consider the birth rate of N.lugens in natural conditions, which is greater than that of N.lugens in resistant rice. Selecting $\bar{a}_{12} = 0.15, 0.25, 0.35, 0.45$ and $0.55$, respectively, the control degree curves under the different migration parameters $\bar{a}_{12}$ are plotted in Figure 8(a). Similarly, we analyze the influences of the successful landing rate $\bar{a}_{21}$ of a southward migration on the control degree of the wild N.lugens in both patches. Figure 8(b) shows the control degree curves of the wild N.lugens within $60$ days under the different values of $\bar{a}_{21}$. Although the stability of $E_1$ makes the control degree of the wild N.lugens close to 1 in a long time, it is not hard to find that increasing the successful landing rate $\bar{a}_{ij}$ of the wild N.lugens can decrease its control degree within a finite time.

7.3.   Measures to reduce the impact of migration on the control of N.lugens

A worthy question is how to take the effective measures to reduce the impact of migration on the control of wild N.lugens. The take-off rate of N.lugens is not easily controlled in practice, while the successful landing rate for N.lugens migration can be controlled through the relevant interference measures. Therefore, the successful landing rates of the two types of N.lugens are the controllable migration parameters. Next, we shall obtain the control measures by comparing the control degree of the wild N.lugens under the different successful landing rate values. First, the impact of the successful landing rate $\bar{a}_{ij}$ on the control degree have been calculated in subsection 7.2. From the numerical simulation results, we find that decreasing the successful landing rate of the wild N.lugens can improve its control effects. In future field trials, we can take measures to reduce the landing of the wild N.lugens in destination area and source area to improve the control effects of the wild N.lugens, such as spraying pesticides during their flight by unmanned aerial vehicle.

Second, by the similar numerical simulation method, the total control degree of the wild N.lugens in those two areas can be calculated under the different landing rates $a_{12}$ and $a_{21}$ of N.lugens infected with wStri. The control degree curves of the wild N.lugens are given in Figure 9(a) and (b). The numerical simulation results indicate that the control degree of the wild N.lugens increases when the successful landing rate of the northward or southward migration of N.lugens infected with wStri appropriately increases. Therefore, increasing the landing rate of N.lugens infected with wStri can improve the control degree of the wild N.lugens. Although it is not possible to directly increase the landing rate of N.lugens infected with wStri, we can achieve the same effect by releasing them in area 1 in autumn or in area 2 in summer. The numerical simulations demonstrate that this measure is feasible.

Moreover, we study the impacts of the initial release quantity and the CI intensity of the wStri-infected N.lugens on the control degree of wild N.lugens. Under the parameters $b = 0.6438$, $d_2 = 0.0031$, $d_1 = 0.0025$, $\xi = 0.675$, $a_{12} = 0.1222$, $a_{21} = 0.5453$, $\bar{a}_{12} = 0.54$, $\bar{a}_{21} = 0.4317$, $b_{12} = 0.4716$, $b_{21} = 0.7554$, $\bar{b}_{12} = 0.6867$ and $\bar{b}_{21} = 0.6439$, the control degree curves of the wild N.lugens with respect to the different initial release quantity of the wStri-infected N.lugens are shown in Figure 10(a). The numerical simulation shows that increasing the initial release quantity of N.lugens infected with wStri can effectively improve the control degree of the wild N.lugens. Similarly, when the CI intensity is selected as $\xi = 0.375$, $0.475$, $0.575$, $0.675$ and $\xi = 0.775$, respectively, and the remaining parameters are the same as that of Figure 10(a). The control degree curves of the wild N.lugens with respect to the different CI intensity are given in Figure 10(b). The numerical simulation results show that an increase of the CI intensity is beneficial to improve the control degree of the wild N.lugens. Therefore, to effectively control the wild N.lugens, the release measures can be taken in area 2 to improve the initial release quantity of the wStri-infected N.lugens with the higher CI intensity.

Previously, we found relevant control measures for migration by using numerical calculation method, namely reducing the landing rate $\bar{a}_{ij}$ or increasing the landing rate $a_{ij}$. However, the lower limit of $\bar{a}_{ij}$ reduction and the upper limit of $a_{ij}$ increase need to be determined theoretically. According to Theorem 2, we can take measures such as increasing the release quantity of wStri in area 1 such that the successful landing rate $a_{21}$ of the southward migration of N.lugens infected with wStri is greater than the threshold value $b_{12}c_1^3-bc_1^2$, or the successful landing rate $a_{12}$ of the northward migration of N.lugens infected with wStri is greater than the threshold value $\frac{1}{c_1^3}(b_{21}-bc_1)$ in area 2. Further, according to stability conditions of $E_1$, in order to ensure local asymptotic stability of $E_1$, the control measures we take should meet the following three conditions $d_2\bar{I}^{(1)}+d_2\bar{I}^{(2)}+\bar{b}_{12}c_1+\frac{\bar{b}_{21}}{c_1} > 2b(1-\xi)$, $\big(-b(1-\xi)+\bar{b}_{12}c_1+d_2\bar{I}^{(1)}\big)\big(-b(1-\xi)+\frac{\bar{b}_{21}}{c_1}+d_2\bar{I}^{(2)}\big) > \bar{a}_{12}\bar{a}_{21}$ and $\big(b+\frac{3\bar{a}_{21}}{c_1^2}\big)(b+3a_{12}c_1^2) > \big(b_{12}c_1^2+\frac{2a_{21}}{c_1}\big)\big(2a_{12}c_1+\frac{b_{21}}{c_1^2}\big)$. Therefore, the wild N.lugens can be completely eliminated by initially releasing a certain quantity of N.lugens infected with wStri, thereby achieving the ideal control state for the wild N.lugens in two areas.

From the existence conditions of the uninfected equilibrium $E_2$ in system (2.3), as long as the artificial control measures taken in area 2 ensure that the successful landing rate $\bar{a}_{12}$ of the wild N.lugens migrating northward does not exceed the threshold value $\frac{1}{c_2^3}(\bar{b}_{21}-bc_2)$, or those taken in area 1 ensure that the successful landing rate $\bar{a}_{12}$ of the wild N.lugens migrating southward does not exceed the threshold value $\bar{b}_{12}c_2^3-bc_2^2$, we can make system (2.3) not have equilibrium $E_2$. Even if the equilibrium $E_2$ exists, we can still take the control measures, such as increasing the landing rate $a_{ij}$ of N.lugens infected with wStri in destination area to make the inequality $(-b+d_1\hat{U}^{(1)}+b_{12}c_2)\big(-b+d_1\hat{U}^{(2)}+\frac{b_{21}}{c_2}\big) < a_{12}a_{21}$ holds, making equilibrium $E_2$ to be unstable.

7.4.   Migration strengthens the N.lugens population

The wild N.lugens migration model (5.4) has two important applications. First, the global asymptotic stability of the positive equilibrium for model (5.4) indicates that the wild N.lugens can survive in both areas long term in the process of south-north migration. Therefore, a successful migration can lead to the rapid development and smooth propagation of N.lugens after adapting to long-distance migration into a new habitat ^[37]. Second, it can be used to simulate the phenomenon of migration strengthening the N.lugens population. In order to know the relationship between migration and population reproduction, Shen and Cheng conducted the experiment of pair breeding on Shanyou 63 hybrid rice by using N.lugens to be emigrated in Guilin and N.lugens moving in An'qing ^[38]. The results showed that the amount of eggs laid and the reproductivity of N.lugens have significantly improved after a long distance migration.

The take-off rate data and the landing rate data of N.lugens are not found, based on the total labeled data of rice planthopper in Xiangxiang, Qidong, Linli and Guiyang observation sites in Hunan Province of China from September 1, 2021 to September 30, 2021 (data from Hunan Provincial plant protection and inspection station), we first estimate the parameters of model (5.4). Notice that the statistical data have a strong randomness. By using the Markov Chain Monte Carlo (MCMC) method to estimate the parameters and the initial value $(b, d_2, \bar{b}_{12}, \bar{a}_{21}, \bar{a}_{12}, \bar{b}_{21}, U_0^{(1)})$ in system (5.4), we obtain the estimation values $b = 0.6438, d_2 = 0.0031, \bar{b}_{12} = 0.2316, \bar{a}_{21} = 0.4317, \bar{a}_{12} = 0.4593, \bar{b}_{21} = 0.3639$, and $U_0^{(1)} = 130.0119$. We provide additional details about the estimation process in Appendix. In order to numerically verify that migration strengthens the N.lugens population, by substituting the estimated parameters into system (5.4), we get that the positive equilibrium $(\hat{U}^{(1)}, \hat{U}^{(2)}) = (260.8292,252.3654)$ and $B_4 = 2.4616 > 0$. Thus, from Theorem 8, the positive equilibrium is globally asymptotically stable. Therefore, the total number of N.lugens eventually approaches $\hat{U}^{(1)}+\hat{U}^{(2)} = 513.1945$. In addition, when there is no migration of N.lugens, the wild N.lugens population only interacts in their respective area; then, system (5.4) is further reduced to $\frac{{\rm d}U^{(1)}}{{\rm d}t} = bU^{(1)}-d_2(U^{(1)})^2$ and $\frac{{\rm d}U^{(2)}}{{\rm d}t} = bU^{(2)}-d_2(U^{(2)})^2$. It is easy to show that the number of N.lugens in south-north areas eventually approaches $\frac{2b}{d_2} = 415.3548$. Therefore, we have $\hat{U}^{(1)}+\hat{U}^{(2)} > \frac{2b}{d_2}$. This numerical simulation result shows that the number of wild N.lugens after migration is more than that of N.lugens without migration, as shown in Figure 11. The main experimental result of Ref.^[38] is verified.

8.   Conclusions

According to the release methods of N.lugens infected with wStri proposed in our previous work ^[6], the wild N.lugens in a paddy field or an area can be effectively controlled. However, N.lugens is an agricultural pest with a long distance migration habits. The Wolbachia spreading dynamic model, which takes the migration of N.lugens into account, can better reflect the reality in future field experiments; additionally, it is beneficial to understand how migration affects the control effect of Wolbachia to provide theoretical guidance for a accurate control of N.lugens and reduction of the risk of N.lugens occurrence.

The northern and southern rice areas were abstracted into two discrete patches. Based on the proposed population density proportional dependent migration coefficient, the discrete migration model reflecting the spread of Wolbachia in N.lugens is established.

The analysis results of equilibria show that the equilibrium infected with wStri and the uninfected equilibrium may be locally asymptotically stable under certain conditions. Based on numerical simulation, there may exist a positive equilibrium for system (2.3). From analyzing the subsystem without migration, we know that the positive equilibrium is always unstable and the subsystem can exhibit a bistable phenomenon.

The long-distance migration can increase the difficulty of controlling wild N.lugens. Numerical simulations results indicate that as the successful landing rate of the wild N.lugens increases, the equilibrium infected with wStri can change from stable state to unstable state, which may cause $E_2$ to be stable. Moreover, as the successful landing rate of the wild N.lugens decreases, the control degree of N.lugens increases, while the control degree of N.lugens increases as the successful landing rate of N.lugens infected with wStri increases. Additionally, some specific measures affecting the landing rate of N.lugens have also been suggested.

The global asymptotic stability of the positive equilibrium for the subsystem of the wild N.lugens migration shows that successful migration can lead to rapid development and smooth propagation of N.lugens after long distance migration. Further numerical simulation results suggested that the quantity of N.lugens population is increased after long distance migration, which is consistent with the experimental results found within literature ^[38].

Use of AI tools declaration

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

Acknowledgments

This work was supported by Hunan Provincial Natural Science Foundation of China [grant number 2023JJ50090, 2021JJ30309].

Conflict of interest

The authors declare there is no conflict of interest.

Appendix: Parameter estimation process based on MCMC

Since the observed data are the total labeled data of rice planthopper in four observation sites, which can be viewed as the number of wild N.lugens in area 1 to a certain extent. In order to use the MCMC method in Bayes statistical computation to estimate the six parameters and the initial value for model (5.4), let $X_i$ and $Y = (y_1, y_2, \cdots, y_{30})^T$ be the real data for wild N.lugens and the observed data, respectively. Because these statistics approximately reflect the Hunan Province labeling data for wild N.lugens in September, rather than the real number of wild N.lugens for area 1, there is strong randomness for these statistics, so we assume that the observed data statistic is a Poisson process. Therefore, given day $i$, under the condition that the real number of wild N.lugens is $x_i$, the probability of the statistical data $y_i$ is

and the observed data each day is independent.

We use Metropolis-Hastings algorithm in MCMC to estimate parameters for model (5.4). Since N.lugens returns to the south from September, we assume that the initial value of wild N.lugens in area 2 is $U_0^{(2)} = 200$. Let the likelihood function be $L(Y|\theta)$, where $\theta = (b, d_2, \bar{b}_{12}, \bar{a}_{21}, \bar{a}_{12}, \bar{b}_{21}, U_0^{(1)})\in \Theta$. In order to determine the posterior distribution, the seven parameters are nonnegative constant due to a practical problem is considered, we select non-informative prior distribution, and assuming the prior distribution is $p(\theta)\propto$ constant. Thus, the posterior distribution is

From the initial value $(b, d_2, \bar{b}_{12}, \bar{a}_{21}, \bar{a}_{12}, \bar{b}_{21}, U_0^{(1)}) = (0.5, 0.2, 0.06, 0.4, 0.5, 0.05, y_1)$, we choose a zero-centered uniform distribution as the proposed distribution, and a random walk Metropolis sampler is taken to produce the posterior distribution $p(\theta|Y)$ as stationary distribution of Markov chains ^[39,40], and implementing the following algorithm to generate random numbers (Markov chain):

1 Given the initial vector value $\theta^{(0)} = (b, d_2, \bar{b}_{12}, \bar{a}_{21}, \bar{a}_{12}, \bar{b}_{21}, U_0^{(1)})$;

2 Given the initial time $i = 0$;

3 Circulation

– Sampling $\theta'$ from the proposed distribution $g(\cdot|\theta_i)$;

– Extracting a random number $u$ from the uniform distribution $U(0, 1)$;

– If $u\leq {\rm min}\Big\{1, \frac{p(\theta'|Y)}{p(\theta|Y)}\Big\}$, then $\theta^{(i+1)} = \theta'$, else $\theta^{(i+1)} = \theta^{(i)}$.

– $i = i+1$;

– After a burn-in period, store $\theta^{(i+1)}$ after a cycle of each $s$ times.

4 Breaking circulation when $i$ is sufficiently large.

Here, we cycle MCMC $10000$ times, and take the burn-in period as $2000$, and calculate the average value of the parameters of the last $8000$ times, and get the estimated values of the seven parameters as follows: $b = 0.6438, d_2 = 0.0031, \bar{b}_{12} = 0.2316, \bar{a}_{21} = 0.4317, \bar{a}_{12} = 0.4593, \bar{b}_{21} = 0.3639$, and $U_0^{(1)} = 130.0119$.