On Large and Small Data Blow-Up Solutions in the Trojan Y Chromosome Model

Abstract

The Trojan Y Chromosome Strategy (TYC) is the only genetic biological control method in practice in North America for controlling invasive populations with an XX–XY sex determinism. Herein a modified organism, that is a supermale or feminised supermale, is introduced into an invasive population to skew the sex ratio over time, causing local extinction. We consider the three species TYC reaction diffusion model, and show that introduction of supermales above certain thresholds, and for certain initial data, solutions can blow-up in finite time. Thus, in order to have biologically meaningful solutions, one needs to restrict parameter and initial data regimes, in TYC type models.

1. Introduction

The detrimental effects of invasive species is well-documented [1,2,3,4,5,6,7]. Current control methods rely, primarily, on chemical treatment [8] and are environmentally detrimental. Biological control is an alternative means of control, where a species is released to control the invasive population by predation, competition, disease, or manipulating the mating system [9,10]. The recent advancement with gene drive technology makes this approach even more promising, from both modeling and practical standpoints [11,12,13,14]. The Trojan Y chromosome strategy (TYC) is a promising genetic biocontrol strategy which circumvents the detriment due to chemical control [15,16,17,18,19]. It involves introducing a YY male or YY male and feminised YY male into an invasive population with an XX–XY sex determinism system. The off-spring of the YY male or YY feminised male are only wild type males or YY males. This skews the sex ratio of subsequent generations towards all males, and extinction of the population may occur (see Figure 1). Note the TYC strategy has also been studied for ZZ–ZW systems [20].

A mathematical model for the TYC strategy was first pioneered in [16], and has been well investigated [11,19,21,22,23,24]. Prior work on TYC [25,26,27,28], shows that depending on the introduction rate of the supermale/feminized supermale, and initial data, eradication is always possible. These results rely on the assumption that solutions remain positive for all non-negative initial data. The TYC strategy is now in current practice and field studies that investigate the survivability and reproductivity of introduced supermales are reported in [14,29]. In current field trials, only the YY supermale population has been introduced into the wild [30]. Subsequently, what we have in practice then is a situation where only a YY supermale is introduced into a population of XX females and XY males. This leads to the following reaction diffusion system:

$$\begin{array}{ccc}\hfill \frac{\partial f}{\partial t}& =& D\Delta f+\frac{1}{2}\beta Lfm-\delta f,\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill \frac{\partial m}{\partial t}& =& D\Delta m+\frac{1}{2}\beta Lfm+\beta Lfs-\delta m,\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill \frac{\partial s}{\partial t}& =& D\Delta s+\mu -\delta s,\hfill \end{array}$$

(3)

specified over a bounded spatial domain $\mathsf{\Omega }\subset {\mathbb{R}}^n,n=1,2,3$ and subject to homogeneous Neumann boundary conditions, that is,

$$\nabla f·n=\nabla m·n=\nabla s·n=0.$$

The state variables $f,m,s$ are the population densities/numbers of individuals in each associated class: XX females, XY males, and YY males, respectively. The logistic term, $L=(1-(f+m+s)/K)$, where K is the carrying capacity of the ecosystem, attempts to penalize or encourage growth of populations when above or below K, respectively. The positive constants $\beta$ and $\delta$ represent the per capita birth and death rates, respectively; non-negative constant $\mu$ denotes the rate at which the YY males s are introduced. The kinetics describe that the female population increases due to mating between the XX females and XY males, whereas the male population increases due to mating between the XX females and XY males, as well as mating between the XX females and YY supermales.

Positive initial data, $\left(f\right(x,0),m(x,0),s(x,0\left)\right)$, is assumed herein. Movement of species modeled via random diffusion, incorporated as reaction diffusion equations, has a long standing history in population dynamics [7,31]. Thus, we assume the species $f,m,s$ move via random diffusion, with equal speeds modeled via equal diffusion coefficients D. These could be thought of as purely movement of the species (D) or pure movement ($D_1$) plus some “eddy” diffusion ($D_1$), which the species could be subject to, due to the flow of the water body they are present in [31] (herein $D=D_1+D_2$). The bounded domain, $\mathsf{\Omega }$, represents the water body the species is present in, such as a lake or stream [32,33,34]. The Neumann or no flux boundary condition implies that no species enters or leaves through the boundary $\mathsf{\Omega }$ [31]. In the setting of a lake or pond, this implies that there are no inlets/outlets in or out of the water body in question. Thus, no species can enter the lake from surrounding water bodies or leave the lake, and go into those other water bodies [35]. Note Dirichlet or lethal boundary conditions could also be imposed, and this would imply the species die on hitting the boundary, due to rotenone poisoning sprayed around a pond for example [8,36].

The current manuscript shows that under certain initial data and parametric restrictions, the PDE model can produce negative (biologically spurious) solutions. Thus, we purport that in order for the TYC model to have practical applicability, such as for simulation purposes in the current field trials [30], one needs to restrict initial conditions and parameters. The current work focuses on showing:

• The three species reaction diffusion TYC model, (1)–(3), results in negative solutions in the male population m. These conditions include the choice of large initial conditions ($f_0+m_0+s_0>>K$) and/or large introduction rate of supermales s. This result is given via Lemma 2;
• The negativity of m may lead to finite time blow-up in $L^p$ of f and m. This happens for large initial data as seen via Theorem 1, but is also possible for small initial data, which is seen via Theorem 2, in the PDE model. Similar results are derived via Corollary A1, Corollary A2 for the ODE model;
• We numerically derive regions of invariance in the phase space, for the three species TYC PDE model, via Figure 2. Additionally, chemotaxis type mechanisms can damp blow-up solutions;
• We discuss the practical relevance of these results to biological control and possible restrictions.

Figure 2. In these simulations we consider the three species TYC PDE model (6)–(8) with Neumann boundary conditions. The phase space is partitioned into three regions showing various dynamics. The blue region in (a,b) indicates where positive solutions are guaranteed for f and m where $s\left(0\right)<s^{*}$ and $\gamma <{\gamma }^{*}$. The magenta region shows where negative solutions for m exist mathematically for $s^{**}>s\left(0\right)\ge s^{*}$ and ${\gamma }^{**}>\gamma \ge {\gamma }^{*}$. In this region, f converges asymptotically to zero in (a) and f goes to a non-negative steady state in (b). The red region in (a,b) indicates where negative solutions for m exist and f blows-up in finite time for $s\left(0\right)\ge s^{**}$ and $\gamma \ge {\gamma }^{**}$. In each simulation in (a), $\gamma =0$ and in (b), $s\left(0\right)=0$. Simulation uses $D=1$ and $\rho =18$.

Figure 2. In these simulations we consider the three species TYC PDE model (6)–(8) with Neumann boundary conditions. The phase space is partitioned into three regions showing various dynamics. The blue region in (a,b) indicates where positive solutions are guaranteed for f and m where $s\left(0\right)<s^{*}$ and $\gamma <{\gamma }^{*}$. The magenta region shows where negative solutions for m exist mathematically for $s^{**}>s\left(0\right)\ge s^{*}$ and ${\gamma }^{**}>\gamma \ge {\gamma }^{*}$. In this region, f converges asymptotically to zero in (a) and f goes to a non-negative steady state in (b). The red region in (a,b) indicates where negative solutions for m exist and f blows-up in finite time for $s\left(0\right)\ge s^{**}$ and $\gamma \ge {\gamma }^{**}$. In each simulation in (a), $\gamma =0$ and in (b), $s\left(0\right)=0$. Simulation uses $D=1$ and $\rho =18$.

2. Finite Time Blow up

2.1. Preliminaries

We aim to prove the possibility of finite time blow-up in (1)–(3). To this end we use standard techniques [37,38,39]. We first recall classical results guaranteeing non-negativity of solutions, local and global existence [38,40]:

Lemma 1.

Let us consider the following $m\times m$ - reaction diffusion system: for all $i=1,...,m,$

$${\partial }_t{u}_i-d_i\Delta u_i=f_i(u_1,...,u_m)in{\mathbb{R}}_{+}\times \mathsf{\Omega },{\partial }_v{u}_i=0\mathrm{on}\partial \mathsf{\Omega },u_i\left(0\right)=u_{i0},$$

(4)

where $d_i\in (0,+\infty )$, $f=(f_1,...,f_m):{\mathbb{R}}^m\to {\mathbb{R}}^m$ is $C^1$ and $u_{i0}\in L^{\infty }\left(\mathsf{\Omega }\right)$. Then, there exists a $T>0$ and a unique classical solution of (4) on $[0,T).$ If $T^{*}$ denotes the greatest of these $T^{\prime }s$, then

$$\left[\underset{t\in [0,T^{*}),1\le i\le m}{sup}\left|\right|u_i\left(t\right){\left|\right|}_{L^{\infty }\left(\mathsf{\Omega }\right)}<+\infty \right]\Rightarrow [T^{*}=+\infty ].$$

If the nonlinearity ${\left(f_i\right)}_{1\le i\le m}$ is moreover quasi-positive, which means

$$\forall i=1,...,m,\forall u_1,...,u_m\ge 0,f_i(u_1,...,u_{i-1},0,u_{i+1},...,u_m)\ge 0,$$

then

$$[\forall i=1,...,m,u_{i0}\ge 0]\Rightarrow [\forall i=1,...,m,\forall t\in [0,T^{*}),u_i\left(t\right)\ge 0].$$

2.2. Negative Solutions

Remark 1.

Notice via (1)–(3), that $f_1(f,m,s)=\frac{1}{2}\beta Lfm-\delta f$, $f_2(f,m,s)=\frac{1}{2}\beta Lfm+\beta Lfs-\delta m$ and $f_3(f,m,s)=\mu -\delta s$. Thus $f_2(f,0,s)=\beta fs\left(1-(f+s)/K\right)$.

Even if $f,s>0$, $\left(1-(f+s)/K\right)$ does not have a definite sign, and so $f_2(f,0,s)=\beta fs\left(1-(f+s)/K\right)$ does not have a definite sign. This suggests negative solutions are a possibility in (2), and so m could possibly be negative.

We next state the following lemma.

Lemma 2.

Consider the TYC system given by (1)–(3). Then, there exists positive initial data $(f_0\left(x\right),m_0\left(x\right),s_0\left(x\right))\in L^{\infty }\left(\mathsf{\Omega }\right)$, and a time interval $[T_1,T_2]\in (0,\infty )$, s.t for solutions emanating from these data, $m(x,t)<0$ on $[T_1,T_2]$.

Proof. 

Consider (1)–(3), then we have,

$$f_1(0,m,s)=0,f_3(f,m,0)=\mu \ge 0.$$

Thus, $f,s\ge 0$ for all time. However,

$$f_2(f,0,s)=\beta fs\left(1-\frac{f+s}{K}\right),$$

and so if one chooses initial data $f_0\left(x\right),s_0\left(x\right)\gg K$, and $0<m_0\left(x\right)<<1$, we have by continuity of solutions $f,s>>K$ for $t\in [0,T^{*}]$, for some time $T^{*}$. Thus, for $t\in [0,T^{*}]$, $f_2(f,0,s)<0$, violating the necessary requirement for positive solutions from Lemma 1, and will yield negative solutions in m, for $t\in [T_1,T_2]\subset [0,T^{*}]$. □

2.3. Finite Time Blow-Up in the PDE Model

We see that for positive solutions to (1)–(3), classical methods [38], can be applied to yield global in time existence of solutions. In our case this is fairly easy as the equations could be added up, to derive an equation for a grouped variable $V=f+m+s$, as the diffusion coefficients are all the same, on which uniform estimates can then be made. However, these methods will not apply if $m<0$. We show next that the negativity of m, can, in turn, lead to finite time blow-up. We state the following theorem,

Theorem 1.

Consider the TYC system given by (1)–(3), with $\mu =0$. Then, there exists positive initial data $(f(x,0),m(x,0),s(x,0))\in L^{\infty }\left(\mathsf{\Omega }\right)$, such that solutions emanating from this data, can blow-up in finite time, that is

$$\underset{t\to T^{*}<\infty }{lim\; sup}{\left|\right|f\left|\right|}_p\to +\infty$$

and

$$\underset{t\to T^{*}<\infty }{lim\; sup}{\left|\right|m\left|\right|}_p\to \infty$$

for all $p\ge 1$.

Proof. 

Consider the equation for f expanded:

$$\frac{\partial f}{\partial t}=\Delta f+\frac{\beta }{2}fm-\frac{\beta }{2K}{f}^{2}m-\frac{\beta }{2K}f{m}^2-\frac{\beta }{2K}fms-\delta f.$$

Via Lemma 2 we know that for large $s(x,0)$ and $f(x,0)$, that is $s(x,0),f(x,0)\gg K$, that $m(x,t)<0$ for $t\in [T_1,T_2]$ and $x\in \mathsf{\Omega }$. Let $m=-\tilde{m}$, where $\tilde{m}>0$ for $t\in [T_1,T_2]$ and $x\in \mathsf{\Omega }$. By direct substitution,

$$\frac{\partial f}{\partial t}=\Delta f-\frac{\beta }{2}f\tilde{m}+\frac{\beta }{2K}{f}^2\tilde{m}-\frac{\beta }{2K}f{\tilde{m}}^2+\frac{\beta }{2K}f\tilde{m}s-\delta f.$$

Integrating over the spatial domain $\mathsf{\Omega }$ and rearranging yields,

$$\begin{array}{ccc}\hfill \frac{d}{dt}{\int }_{\mathsf{\Omega }}fdx& =& \frac{\beta }{2K}{\int }_{\mathsf{\Omega }}{f}^2\tilde{m}dx+\frac{\beta }{2K}{\int }_{\mathsf{\Omega }}f\tilde{m}sdx-\frac{\beta }{2}{\int }_{\mathsf{\Omega }}f\tilde{m}dx\hfill \\ & & -\frac{\beta }{2K}{\int }_{\mathsf{\Omega }}f{\tilde{m}}^{2}dx-\delta {\int }_{\mathsf{\Omega }}fdx.\hfill \end{array}$$

Since $m(x,t)<0$ for $(x,t)\in \mathsf{\Omega }\times [T_1,T_2]$ then $-{\delta }_2<m(x,t)<-{\delta }_1$, for positive constants ${\delta }_1$ and ${\delta }_2$. Thus via standard comparison as earlier and Hölder’s inequality we have,

$$\frac{d}{dt}{\int }_{\mathsf{\Omega }}fdx\ge \frac{\beta {\delta }_1}{2K}{\left({\int }_{\mathsf{\Omega }}fdx\right)}^2-\left(\frac{\beta {\delta }_2}{2}+\frac{\beta {\left({\delta }_2\right)}^2+\delta }{2K}\right){\int }_{\mathsf{\Omega }}fdx.$$

Define ${\displaystyle F\left(t\right)={\int }_{\mathsf{\Omega }}fdx}$, then,

$$\frac{d}{dt}F\left(t\right)\ge \frac{\beta {\delta }_1}{2K}{\left(F\left(t\right)\right)}^2-\left(\frac{\beta {\delta }_2}{2}+\frac{\beta {\left({\delta }_2\right)}^2+\delta }{2K}\right)F\left(t\right)$$

which yields the finite time blow-up of $F\left(t\right)$, for large enough initial data. That is for

$$F\left(0\right)={\int }_{\mathsf{\Omega }}f(x,0)dx\ge \frac{\left(\frac{\beta {\delta }_2}{2}+\frac{\beta {\left({\delta }_2\right)}^2+\delta }{2K}\right)}{\left(\frac{\beta {\delta }_2}{2K\delta }\right)}.$$

Thus, the $L^1\left(\mathsf{\Omega }\right)$ norm of f blows-up in finite time. Since $L^p\left(\mathsf{\Omega }\right)\hookrightarrow L^1\left(\mathsf{\Omega }\right)$, for $p\ge 1$, we have that the $L^p$ norm of f blows-up for any p, for large enough initial conditions. This completes the proof of the blow-up of f. In order to establish the proof of blow-up in m, we can proceed by contradiction. That assumes ${\displaystyle F\left(t\right)={\int }_{\mathsf{\Omega }}fdx}$ blows up at time $T^{*}$, but ${\displaystyle G\left(t\right)={\int }_{\mathsf{\Omega }}mdx}$, does not and is bounded. So $F\left(T^{*}\right)=+\infty$, but $-\infty <G\left(T^{*}\right)<+\infty$. Let us integrate (2) in the time interval $[0,T^{*}]$ and then over the spatial domain $\mathsf{\Omega }$, we obtain

$$-\infty <G\left(T^{*}\right)={\int }_0^{T^{*}}{\int }_{\mathsf{\Omega }}\left(\frac{1}{2}\beta Lfm+\beta Lfs-\delta m\right)dxdt=-\infty .$$

This follows due to the negativity of m, and the blow-up of f, via Theorem 1. Thus, we have a contradiction, which implies the blow-up of m at $T^{*}$ as well. □

Remark 2.

Note, Theorem 1 shows that for (1)–(3), we do not have control of mass or $L^1\left(\mathsf{\Omega }\right)$ control.

The previous theorem proves that the finite time blow-up is a possibility even if the only introduction of s is through the initial condition. In the following theorem we prove that regardless of the initial condition size that there exists a threshold to the introduction rate, $\mu$, such that rates beyond this value will lead to finite time blow-up in f.

Theorem 2.

Consider the TYC system given by (1)–(3), with $\mu >0$. Then, for any positive initial data $(f(x,0),m(x,0),s(x,0))\in L^{\infty }\left(\mathsf{\Omega }\right)$, there exists a ${\mu }^{*}$, such that if $\mu >{\mu }^{*}$, then solutions emanating from these data, can blow-up in finite time, that is

$$\underset{t\to T^{*}<\infty }{lim\; sup}{\left|\right|f\left|\right|}_p\to +\infty$$

and

$$\underset{t\to T^{*}<\infty }{lim\; sup}{\left|\right|m\left|\right|}_p\to \infty$$

for all $p\ge 1$.

Proof. 

Let $\mu >\delta K$. Following similar estimates as in the previous theorem yields,

$$\frac{d}{dt}{\int }_{\mathsf{\Omega }}fdx=\frac{\beta }{2K}{\int }_{\mathsf{\Omega }}{f}^2\tilde{m}dx+\frac{\beta }{2K}{\int }_{\mathsf{\Omega }}f\tilde{m}sdx-\frac{\beta }{2}{\int }_{\mathsf{\Omega }}f\tilde{m}dx-\frac{\beta }{2K}{\int }_{\mathsf{\Omega }}f{\tilde{m}}^{2}dx-\delta {\int }_{\mathsf{\Omega }}fdx.$$

Via Holder’s inequality we have,

$$\begin{array}{ccc}\hfill \frac{d}{dt}{\int }_{\mathsf{\Omega }}fdx& \ge & C_1\frac{\beta {\delta }_1}{2K}{\left({\int }_{\mathsf{\Omega }}fdx\right)}^2+\frac{\beta {\delta }_2\mu }{2K\delta }{\int }_{\mathsf{\Omega }}fdx-\left(\frac{\beta {\delta }_2}{2}+\frac{\beta {\left({\delta }_2\right)}^2+\delta }{2K}\right){\int }_{\mathsf{\Omega }}fdx.\hfill \\ & =& C_1\frac{\beta {\delta }_1}{2K}{\left({\int }_{\mathsf{\Omega }}fdx\right)}^2+\left(\frac{\beta {\delta }_2\mu }{2K\delta }-\left(\frac{\beta {\delta }_2}{2}+\frac{\beta {\left({\delta }_2\right)}^2+\delta }{2K}\right)\right){\int }_{\mathsf{\Omega }}fdx\hfill \end{array}$$

The result follows for ${\mu }^{*}$ chosen s.t,

$$\mu >{\mu }^{*}={\displaystyle \frac{\left({\displaystyle \frac{\beta {\delta }_2}{2}}+{\displaystyle \frac{\beta {\left({\delta }_2\right)}^2+\delta }{2K}}\right)}{\left({\displaystyle \frac{\beta {\delta }_2}{2K\delta }}\right)}}.$$

(5)

□

Remark 3.

The thresholds provided in (5) are not sharp. In fact, numerical experiments suggest the critical values of initial condition size or introduction may be much smaller than (5) for finite time blow-up. However, this result shows that small data blow-up is possible in the three species TYC system, if the introduction rate μ is sufficiently large.

3. Numerical Experiments for TYC PDE Model

We perform numerical experiments to show that the f variable has a potential to blow-up in finite time in a three species classical TYC PDE model. We also show via simulations that blow-up phenomenon in the f variable can be curtailed when spatial damping mechanisms are introduced.

3.1. Scaling of TYC PDE Model

We non-dimensionalize the PDE form of TYC models. We let $f\to {\displaystyle \frac{f}{K}}$, $m\to {\displaystyle \frac{m}{K}}$, $s\to {\displaystyle \frac{s}{K}}$, $\tau \to \delta t$, $\gamma ={\displaystyle \frac{\mu }{K\delta }}$ and $\rho ={\displaystyle \frac{\beta K}{2\delta }}$. The dimensionless variables $\rho$ is a ratio of the two time scales in the TYC model, that is, the birth and death rates. D is the diffusion co-efficient. The dimensionless spatio-temporal three species TYC model is given as,

$$\begin{array}{ccc}\hfill \frac{\partial f}{\partial t}& =& D\Delta f+\rho mfL-f,\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill \frac{\partial m}{\partial t}& =& D\Delta m+\rho mfL+2\rho sfL-m,\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill \frac{\partial s}{\partial t}& =& D\Delta s+\gamma -s,\hfill \end{array}$$

(8)

specified over the scaled spatial domain $\mathsf{\Omega }=(0,1)$ and $t\in (0,\infty )$. We, prescribe both homogeneous and Neumann boundary conditions in simulations seen in Figure 3.

3.2. Spatial Damping

We aim to investigate spatial pressures, via numerical simulations, that might inhibit or damp the blow-up. For example, what happens if f moves towards lower concentrations of s? We now consider the spatially explicit version of the TYC model with a chemotaxis term [41] given by,

$$\begin{array}{ccc}\hfill \frac{\partial f}{\partial t}& =& \Delta f+\chi \nabla ·(f\nabla s)+\frac{1}{2}\beta Lfm-\delta f,\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill \frac{\partial m}{\partial t}& =& \Delta m+\frac{1}{2}\beta Lfm+\beta Lfs-\delta m,\hfill \end{array}$$

(10)

$$\begin{array}{ccc}\hfill \frac{\partial s}{\partial t}& =& \Delta s+\mu -\delta s,\hfill \end{array}$$

(11)

specified over the domain $(x,t)\subset \mathsf{\Omega }\times (0,\infty )$ and subject to a Neumann and mixed boundary conditions on the boundary $\partial \mathsf{\Omega }$, that is,

$$\nabla f·n=\nabla m·n=0,s+\nabla s·n=0,$$

where L, $\beta$, $\delta$, and $\mu$ are as defined previously and $\chi$ is the chemotatic parameter. Here, we use mixed boundary conditions for s, as homogenous Neumann boundary conditions would lead to blow-up upon a simple integration of (9). Our goal is to investigate the effect of the chemotaxis pressure herein. See Figure 4 and Figure 5.

4. Discussion and Conclusions

This paper proves and provides numerical experiments indicating negativity of solutions and subsequent finite time blow-up in the three species reaction diffusion TYC model. The cause of the negativity of solutions is a result of the form of the reaction terms incorporating mating between females and males, as well as the mating between females and supermales. Note, although classical methods to prove global existence [37,38,40] require quasi-positivity, see Lemma 2, there are not many examples in the literature (to the best of our knowledge) where loss of positivity, leads to global non-existence such as via finite time blow-up.

In the three species reaction diffusion model (in the most realistic case $\mu =0$) blow-up seems to occur only for sufficiently large initial conditions, about three times the magnitude of carrying capacity. A necessary condition on the data guaranteeing globally existing solutions for initial conditions between carrying capacity and three times carrying capacity, is yet to be proven. Therefore, caution and discretion should be taken in choosing the parameter regime and initial condition size when utilizing the TYC model for predictions of the efficacy of the TYC strategy. Additionally note, in [26], population experiments with guppy fish were performed and subsequently used to determine the best parameters, in the least squares sense, to the mating model, with no supermales. The best fit parameters suggested $\rho \approx 17.8125.$ Thus, we use this value throughout the numerical simulations.

Revisions to the TYC model may remove this inconsistency present in current models. A promising revision and modification to the classical TYC model are models that include the strong Allee effect [11,23,42] and intraspecies competition [11] for females among males and supermales, and finite time effects [43]. Sexual selection, might also play a role herein. For example, if f selects m over the mutant male s, and so moves away from s, blow-up can be attenuated, see Figure 4 and Figure 5.

Numerical experiments in these works indicate that negativity of solutions remains a mathematical possibility for certain parameter regime and initial data choice—but blow-up is not noticed in simulations. We conjecture that the intraspecies competition terms remove this unrealistic dynamic since large populations will attenuate the growth rates of the female and male populations. An interesting numerical observation is that a blow-up solution in the 3 species model can be attenuated in the 4 species one [16,36], if enough feminised supermales are introduced. Alternatively, one can phrase this as introducing supermales (in 3 species model) accentuates blow-up, but introducing feminized supermales (in 4 species model) attenuates blow-up. This remains unproven and is the focus of our current and future research efforts.

Clearly, additional refinement of the TYC model is still required to eliminate the possibility of negative solutions and/or blow-up solutions. An interesting question to study is if the choice of very different diffusion coefficients in the species could lead to attenuation of negativity and blow-up—at least in certain regimes of initial conditions. Additionally, alternatives to the classical logistic term are plentiful. For instance, one such alternative is to consider a logistic term of the form $\exp (1-(f+m+s\left)\right)$. Hence, populations are still penalized with a dampened growth rate when populations exceed the carrying capacity. However, preliminary analysis of this type of logistic term have proven difficult for mathematical analysis and, moreover, can generate stable nontrivial equilibrium solutions, such that the total population size asymptotically approaches values above the carrying capacity. Another alternative, is based on the work of mating models, see [44] and references therein. Another interesting direction, although challenging could be relaxing of the quasi-positivity condition, via Lemma 2, to the reaction terms being

$$\forall i=1,...,m,\forall u_1,...,u_m\ge 0,f_i(u_1,...,u_{i-1},0,u_{i+1},...,u_m)\ge -C,$$

$C>0$, and derive bounds on the C, perhaps in terms of initial conditions or other parameters in the system, to prove global existence or blow-up. In particular it might be interesting to see how the C in question scales with the $ϵ$, from recent results, where $|f_i{\left(u\right)|\le |u|}^{2+ϵ}$, is sufficient for global existence [45,46].

In all, we show that earlier TYC models produce unphysical (blow-up) solutions, if initial conditions are chosen much larger than carrying capacity. This is due to a modeling flaw in their formulation. It is important to note that the TYC strategy remains the most effective use of genetic bio-control in the United States. Preliminary results from its use in field trials are very promising [12]. It also yields as fast an extinction, if we restrict initial conditions to less than carrying capacity. Note, modeling forecasts using classical TYC have maintained these restrictions [47].

It, thus, remains an open and paramount problem in invasive species control to continue to determine a complete mathematical model for the TYC strategy. That is one which is valid in a full parameter regime, and applicable in various realistic ecological scenarios, that also yields realistic solutions. Additionally, parametrizing these new models to field data, such as [12], as it becomes available, should also be a primary focus of current and future work on invasive species control.