that both $r_{f}$ and $r_{m}$ approach $0$ as $t\rightarrow+\infty$. We denote by $$\label{eqintroduce3}
u(t)=I_{f}+I_{m} ~\hbox{ and }~v(t)=U_{f}+U_{m}.$$ Assuming equal determination case, which means that $\delta={1}/{2}$, $I_{f}=I_{m}$ and $U_{f}=U_{m},$ then system can be obtained by setting $b_{1}={ b_{I}}/{2}$ and $b_{2}={ b_{U}}/{2}$. In order to obtain the spatiotemporal dynamics of , Huang et al. [@huangmugen1; @huangmugen2] studied the following reaction-diffusion system: $$\label{Eq_Huang}
\left\{\begin{array}{ll}
\displaystyle{\frac{\partial u}{\partial t}= d_{1}\Delta u+u(b_{1}-\delta_{1}(u+v)),}&t>0,~~ x\in \Omega, \vspace{7 pt}\\
\displaystyle{\frac{\partial v}{\partial t}=d_{2}\Delta v+v\left(\frac{b_{2}v}{u+v}-\delta_{2}(u+v)\right),}&t>0, ~~ x\in \Omega,\vspace{7 pt}\\
\displaystyle{\frac{\partial u}{\partial\nu}=\frac{\partial v}{\partial\nu}=0,} &t>0,~~  x\in\partial\Omega,\vspace{7 pt}\\
u(0,x)=u_{0}(x), ~~ v(0,x)=v_{0}(x), &x\in \Omega.
\end{array}\right.$$ In , $d_{1}$ and $d_{2}$ are the diffusion rates, $\Delta$ denotes the Laplace operator in the spatial variable $x$, and $\nu$ denotes the unit outward normal vector to the boundary of $\Omega$. We mention that is obtained from a delay differential equation model in [@BOZHENG] after ignoring the delay factor and incorporating the spatial inhomogeneity. Similarly, there has been several mathematical models formulated to describe the [*Wolbachia*]{} spreading dynamics [@huangmugen2018JTB; @yujianshe2018siam; @yujianshe-zhengbo2019; @BOZHENG-JBD]. These models focused on studying the subtle relation between the threshold releasing level for [*Wolbachia*]{}-infected mosquitoes and several important parameters including the CI intensity and the fecundity cost of [*Wolbachia*]{} infection.

We also note that female Aedes 