Optimal Control and Analysis of a Modified Trojan Y-Chromosome Strategy

TL;DR

This paper develops and analyzes a new mathematical model for the Trojan Y Chromosome strategy, incorporating intraspecies competition and Allee effects, to optimize invasive species eradication methods.

Contribution

It introduces a novel model of the TYC strategy that accounts for previously neglected biological factors and applies optimal control analysis to improve eradication strategies.

Findings

Optimal control strategies can significantly enhance eradication success.

Intraspecies competition and Allee effects influence the effectiveness of the TYC strategy.

Mathematical analysis provides insights into biological control of invasive species.

Abstract

The Trojan Y Chromosome (TYC) Strategy is a promising eradication method that attempts to manipulate the female to male ratio to promote the reduction of the population of an invasive species. The manipulation stems from an introduction of sex-reversed males, called supermales, into an ecosystem. The offspring of the supermales is guaranteed to be male. Mathematical models have shown that the population can be driven to extinction with a continuous supply of supermales. In this paper, a new model of the TYC strategy is introduced and analyzed that includes two important modeling characteristics, that are neglected in all previous models. First, the new model includes intraspecies competition for mates. Second, a strong Allee effect is included. Several conclusions about the strategy via optimal control are established. These results have large scale implications for the biological control of invasive species.

Figures

Tables

Table 1. Table 1 : The lower (LCL) and upper (UCL) 95 % percent 95 95\% confidence interval for the simulations at increasing noise level. The LCL and UCL are used to create the regression lines in Fig.  7 .

Percent Noise	LCL	Mean	UCL
10	0.999914857	0.999959708	1.000004558
20	0.999896955	0.999979216	1.000061477
30	0.999835802	0.999986888	1.000137975
40	0.999776235	0.999966538	1.00015684
50	0.999770905	1.000003822	1.00023674
60	0.999666614	1.000169234	1.000471855
70	0.999749978	1.000115747	1.000481516
80	0.999739497	1.000118186	1.000496874
90	0.999570222	1.000001015	1.000431809
100	0.999672431	1.000142234	1.000612037

Equations

$$\dot{f}$$

$$\dot{m}$$

$$\dot{s}$$

$$\dot{f}$$

$$\dot{m}$$

$$\dot{s}$$

$${\rho }_1(m,s)=\frac{m}{m+s},{\rho }_2(m,s)=\frac{s}{m+s}.$$

$$\dot{f}$$

$$\dot{m}$$

$$\dot{s}$$

$$f\left(rL\left(\frac{f}{a}-1\right)m-1\right)=0m\left(rL\left(\frac{f}{a}-1\right)f-1\right)=0$$

$$m=f=\frac{1}{rL(f/a-1)}$$

$$g(f)=-2f^3+(2a+1)f^2-af-\frac{a}{r}.$$

$$\displaystyle J=\left(\begin{array}[]{ccc}\kappa+\mu-1&\kappa&c\\ \kappa+\mu&\kappa-1&c\\ 0&0&-1\end{array}\right)$$

$$\kappa$$

$$c$$

$$\mu$$

$${\lambda }_{1,2}$$

$${\lambda }_2(f_i)=3+2f_{i}r-\left(2+\frac{1}{a}\right)f_i^{2}r$$

$$q(r,a)=r\left(2+16a^3+ar+4a^{3}r+\sqrt{ar(3+6a+ar)}+4a^2\sqrt{ar(3+6a+ar)}\right)$$

$$-2(3+6a+2ar)\left(ar+\sqrt{ar(3+6a+ar)}\right)>0$$

$$J_0(\gamma )={\int }_0^T-(f+m)-\frac{1}{2}{\gamma }^{2}dt$$

$$U_0=\{\gamma |\gamma \text{mbox}measurable,0\le \gamma <\infty ,t\in [0,T],\forall T\}.$$

$$J_0({\gamma }^{\ast })=\underset{\gamma }{\max }{\int }_0^T-(f+m)-\frac{1}{2}{\gamma }^{2}dt$$

$$H_0=-(f+m)-\frac{1}{2}{\gamma }^2+{\lambda }_1{f}^{\prime }+{\lambda }_2{m}^{\prime }+{\lambda }_3{s}^{\prime }$$

$$\begin{array}{rl}{\lambda }_1^{\prime }(t)=& {\lambda }_2\frac{r(m^2+2s^2)}{m+s}\left(f\left(\frac{f}{a}-1\right)+\left(\frac{2f}{a}-1\right)(f+m+s-1)\right)+\\ & {\lambda }_1\frac{m^{2}r}{m+s}\left(f\left(\frac{f}{a}-1\right)+\left(\frac{2f}{a}-1\right)(f+m+s-1)+1\right)+1\\ {\lambda }_2^{\prime }(t)=& \left(\frac{f}{a}-1\right)[{\lambda }_1\frac{fmr}{m+s}\left(m+2(f+m+s-1)-\frac{m(f+m+s-1)}{m+s}\right)+\\ & {\lambda }_2\frac{fr(m^2+2s^2)}{m+s}+{\lambda }_2\frac{fr(f+m+s-1)}{m+s}\left(2m-\frac{m^2+2s^2}{m+s}\right)]+1\\ {\lambda }_3^{\prime }(t)=& \left(\frac{f}{a}-1\right)\frac{fr}{m+s}[{\lambda }_2\{(f+m+s-1)\left(4s-\frac{(m^2+2s^2)}{m+s}\right)+\\ & (m^2+2s^2)\frac{}{}\}+{\lambda }_1\left(m^2-\frac{m^2(f+m+s-1)}{(m+s)}\right)]+{\lambda }_3\end{array}$$

$${\lambda }_1(T)={\lambda }_2(T)={\lambda }_3(T)=0$$

$$\frac{\partial H}{\partial \gamma }={\lambda }_3-\gamma$$

$${\gamma }^{\ast }(t)=\max (0,{\lambda }_3)$$

$$r\approx 17.8125,a\approx .06,0<t<\mathrm{12.96.}$$

$$df$$

$$dm$$

$$ds$$

Full text

1 Introduction

The detrimental effects of aquatic invasive species is well-documented [3, 4, 5, 6, 17, 19, 21, 30, 33]. Subsequently, a tremendous amount of effort by habitat controllers is devoted to designing effective eradication strategies such as chemical treatment, local harvesting, dewatering, ichthyocides, or a suitable combination [29]. However, these methods are known to negatively impact ecosystems, which may be already stressed by the presence of an aquatic invasive species [29].

The Trojan Y chromosome strategy (TYC) is a new eradication strategy which circumvents many of the known negative ecological impacts due to current practice [12, 32, 29]. TYC strategy involves an introduction of a sex-reversed male. The off-spring of the sex-reversed male, called a supermale, with a wild-type female is guaranteed to be male. Therefore, subsequent generations become male-dominant and this skews the sex ratio towards more males. The goal, is that through the gradual reduction in the female population, extinction of the population may occur (see Fig. 1). The supermale is not a genetically modified organism (GMO) and the TYC process is reversible, that is, if the introduction of supermales is stopped then the supermale population will die out [27]. The TYC strategy has seen tremendous experimental and theoretical interest [14, 28, 32, 7, 8, 9, 13, 22, 24, 23, 25, 26, 34, 35, 37].

The classical population model of the TYC strategy relates the populations of the wild-type XX females ($f$), wild-type XY males ($m$), and the YY supermale ($s$) populations over time. A mathematical model was first proposed by Gutierrez, Teem, and Parshad in [12, 32]:

[TABLE]

where $L=1-\frac{f+m+s}{K}$, $K$ is the carrying capacity, $\beta$ is the birth rate, $\delta$ is the death rate, and $\mu$ is the constant introduction rate. The parameters and populations are assumed to be nonnegative. Due to the nonlinearities, it is not necessary to assume that $\beta>\delta>0$, to obtain a persistent invasive population. It was shown in [32] that there exists a $\mu^{*}$ such that for all $\mu>\mu^{*}>0$ that $f,m\to 0$ in infinite time.

In the current manuscript, we seek to include two important and relevant modeling features. First, the influx of supermales leads to competition between wild-males and supermales for female mates. Second, if the female population is below a given threshold then the population loses fitness and extinction should occur [10, 15]. This latter phenomena is called the Allee effect. These two modeling features are introduced in Section 2. A stability analysis of the equilibrium solutions is also included in Section 2. In Section 3, we then investigate the influence that intraspecies competition and the Allee effect have on the optimal introduction rate that minimizes an objective function based on the total wild population and introduced super males. A stochastic model is introduced to examine the influence of noise on the birth and death rates have on the objective function’s value at the optimal introduction rate.

2 Modified TYC Model with Strong Allee Effect

In this paper, we investigate and propose a new model of the TYC strategy. Namely,

[TABLE]

where $L,\mu,\beta,$ and $\delta$ are as before. Again, the parameters and populations are assumed to be nonnegative and that $\beta>\delta>0$, that is, the birth supercedes the death rate.

Intraspecies competition between wild male and supermale populations for female mates is modeled through the nonnegative saturation term:

[TABLE]

The saturation terms provide the percentage of the total male population that is either wild-type or supermale. Clearly, the range of the saturations term is $[0,1].$ Notice, for a fixed wild-type population, as $s\rightarrow\infty$ then $\rho_{1}\rightarrow 0^{+}$ and $\rho_{2}\rightarrow 1^{-}$. In this situation, the birth of females approaches zero while male offspring only occurs from female and supermale progeny. Likewise, for a fixed supermale population, as $m\rightarrow\infty$ then $\rho_{1}\rightarrow 1^{-}$ and $\rho_{2}\rightarrow 0^{+}$; subsequently, male progeny only occurs from wild-type male and female mating. Therefore, $\rho_{1}$ and $\rho_{2}$ attempt to model the difficulty of wild-type female finding suitable mates from either wild-type or supermale populations.

The term $f/\alpha-1$ models a strong Allee effect and represents a loss of fitness in the female population when below the Allee threshold, $\alpha$. This effect is problem of undercrowding of a species and was first motivated by observations made by Allee in 1927 [1, 2, 20] . Since then, numerous evidences of this effect have been established [15, 31]. Notice, that when $f/\alpha-1<0$ then $\dot{f}<0$, provided $L>0$. Therefore, the female population will decrease toward extinction. Subsequently, a goal of the TYC strategy is to push, via the introduction of supermales, the female population below the Allee threshold.

In the forthcoming analysis, the equations are rescaled. The populations are scaled by the carrying capacity while the $t$ is scaled by the deathrate, that is $t\rightarrow t/\delta$. The rescaled equations are:

[TABLE]

where $L\rightarrow 1-(f+m+s)$, $r=\dfrac{K\beta}{2\delta}>1$, $a=\dfrac{\alpha}{K}\ll 1$, and $\gamma=\dfrac{\mu}{K\delta}.$

2.1 Equilibria and Stability Analysis.

A clear requirement of a valid intervention strategy is that in the absence of the strategy the invasive population would persist, while if the intervention strategy was employed control, ideally extinction, of the invasive species would be established. Therefore, an effective TYC strategy is one that provides an introduction rate $\gamma$ and initial supermale population such that the wild-population is driven to extinction, while in the supermale free case the population would persist.

It is clear that in the modified model exhibits the equilibrium solution $(0,0,\gamma)$, which is deemed the extinction state. Once $f<a$, the introduction rate of supermales, $\gamma$, should be set zero, causing $s\to 0$. Subsequently, we investigate the equilibrium solutions in the case where $\gamma=0$. A necessary requirement for a successful TYC strategy that there exists some time such that In the forthcoming equilibrium analysis we are interested in equilibrium solutions in the situation where $\gamma=0$. Obviously, when $\gamma=0$ then in equilibrium $s=0$, therefore the intraspecies terms $\rho_{1}(m,0)=1$ and $\rho_{2}(m,0)=0.$ In such case:

[TABLE]

We seek to analyze the stability and presence of nontrivial equilibrium solutions. Noting that in equilibrium $f/a-1\neq 0$ then

[TABLE]

Hence, all equilibrium solutions fall on the line $f=m$. The nontrivial equilibrium solutions are roots to the third degree polynomial:

[TABLE]

By Descartes rule of signs there always exists a one negative real root, which is neglected since this is not realistic. In addition, there are either two or zero positive real roots. If there are zero positive real roots then the only equilibrium solution is the extinction state and is globally attracting. In such case, the TYC strategy is not necessary. Therefore, we assume there are two real roots, $f_{1}$ and $f_{2}$, where $0<f_{1}<f_{2}.$

The Jacobian of our system of equations is

[TABLE]

where

[TABLE]

Clearly, $\lambda_{3}=-1$ is an eigenvalue, which indicates exponentially decay in the supermale population. The remaining two eigenvalues are determined by the characteristic equation of the submatrix $J_{33}$. Namely,

[TABLE]

To investigate the sign of $\lambda_{2}$ we recall that $f_{i}$ is a root of Eq. $(\ref{EquilEq})$. This fact, is used to determine an expression for $\lambda_{2}$ as a function of $f_{i}$:

[TABLE]

Notice, that this function is a concave down quadratic with a maximum location located at $f=1/(2+1/a)>0$. Define $f_{+}$ is the positive root of the quadratic function $\lambda_{2}(f)$. Then $\lambda_{2}>0$ if $f_{i}<f_{+}$ and $\lambda_{2}<0$ if $f_{i}>f_{+}.$ Notice that $\lambda_{2}(1)=3-r/a$. Since $a\ll 1$ implies that $r/a\gg 1$ then $\lambda_{2}(1)<0$. Subsequently, we have $0<f_{+}<1$.

Notice that $g(0),g(1)<0$. Now, $g(f_{+})>0$ provided that

[TABLE]

Fig. 2 shows the region where $q(r,a)>0$ and $q(r,a)<0$. While the Allee threshold is difficult to know precisely for any biological system [10], it is reasonable to expect values less than $5\%$ of the carrying capacity [10, 20]. In such case, $g(f_{+})>0$ then by the Intermediate Value Theorem we have $f_{1}<f_{+}<f_{2}$. Subsequently, $(f_{1},f_{1},0)$ and $(f_{2},f_{2},0)$ are a saddle and a sink, respectively. We call $(f_{2},f_{2},0)$ the sustained state. An effective TYC strategy will push populations away from the sustained state and toward the basin of attraction, ideally into, of the extinction state.

3 Optimal Control Analysis

3.1 Optimal control analysis

The goal of this section is to investigate the mechanisms in our TYC system of equations, that, if controlled, could lead to optimal levels of both densities. We assume that the scaled introduction rate $\gamma$ is not known a priori and enter the system as a time-dependent control, such that $0\leq\gamma(t)<\infty$. Consider the objective function

[TABLE]

subject to the governing equations $(\ref{feq1})$-$(\ref{seq1})$ and initial conditions. Optimal strategies are derived for the objective function, where we minimize both female and male populations while also minimizing the introduction rate $\gamma$. Optimal controls are searched for within the set $U_{0}$, namely,

[TABLE]

The goal is to seek an optimal $\gamma^{*}(t)$ such that,

[TABLE]

We use the Pontryagin’s maximum principle to derive the necessary conditions on the optimal control [16]. The Hamiltonian for $J_{0}$ is given by

[TABLE]

We use the Hamiltonian to find a differential equation of the adjoint $\lambda_{i},i=1,2,3.$ Namely,

[TABLE]

with the transversality condition given by

[TABLE]

In consideration of the optimality conditions, the Hamiltonian function is differentiated with respect to control variable $\gamma$ resulting in:

[TABLE]

A compact way of writing the optimal control $\gamma^{*}$ is

[TABLE]

3.2 Numerical simulations

In this section, we will numerically simulate the optimal control for the modified TYC model with strong Allee effect. The following unscaled parameters used for simulation are provided from a least squares approximation of population experiments of guppy fish in [18], namely, $\beta=0.0057,\delta=0.0648,K=405,\alpha=24$, and time interval of $(0,200)$. In scaled variables:

[TABLE]

In Fig. 3, the numerically determined optimal $\gamma(t)$ is shown in conjunction with the scaled populations using the optimal control $\gamma(t)$. This indicates, that a graduate introduction of supermales is ideal. Notice that the maximum introduction rate is located near the inflection point of the female population at $t\approx 4.25$. Therefore, the acceleration of the declining population can be used as an indicator as to when to begin rapidly reducing the introduction rate.

As a basis of comparison, the optimal $\mu$ in unscaled variables for the classical model was determined in [18] (see Fig. 5), which is shown in Fig. 4 in scaled variables for convenience. Without modeling intraspecies competition and the Allee effect yields a larger objective function value at the optimal introduction rate. In particular, the objective values for the classical and modified models are $-501.9203$ and $-315.1675$, respectively. The reduced value of the objective function in the modified model is clear by inspection of the plots of the integrands of the corresponding objective functions given in Fig. 5. Therefore, the inclusion of intraspecies competition and the Allee effect greatly influence the overall cost of the TYC strategy and indicate the strategy is less costly than previously considered in [18]. In the case, of the TYC model that includes intraspecies competition for mates and does not consider the Allee effect obtains an objective value of $-381.3074$, which is still considerably lower than the classical model. However, this suggests that intraspecies competition for mates has a greater influence on the choice of $\gamma$ than the Allee effect.

In Fig. 6 we examine the influence of increasing the value of the dimensionless variable $r$. As $r$ increases we see that optimal control is maintains the same overall shape. In particular, the initial introduction remains the same, however, the initial drop in the introduction rate is less severe in the case of higher reproductive rates.

3.3 Stochastic Model $\&$ Sensitivity Analysis

There are numerous environmental influences that may cause perturbations to the birth and death rates of the invasive species. Here, we assume that $\beta$ and $\delta$ fluctuate around average values [11, 36]. Therefore, the birth and death rates may be treated as random variables such that $\beta\mapsto\beta+\sigma_{\beta}\dot{W}_{\beta}$ and $\delta\mapsto\delta+\sigma_{\delta}\dot{W}_{\delta}$ [36], namely,

[TABLE]

where $W_{\beta}$ and $W_{\delta}$ are independent Brownian motions with intensities $\sigma_{\beta}$ and $\sigma_{\delta}$, respectively.

Here, we examine the influence of noise on the objective value in simulations with the optimal control given Fig. 5. The noise on the birth, $d\beta$, and death, $d\delta$, rates are taken from a normal distribution. An Euler–Maruyama method is employed to determine a numerical solution to the stochastic model. In each simulation we determine the objective function value (Eq. 12) and compare the percent difference between the objective function value determined in case of no noise. The results for increasing percent noise are shown in Fig. 7. A linear regression through the $95\%$ confidence interval is shown as a reference and the intervals are given in the Table 1. Notice, that for large amount of noise in the death and birth rates that the objective function is no more than $8\%$ away from the no noise situation. The slopes of the upper and lower $95\%$ confidence intervals are relatively small and, subsequently, means that the value of objective function is changed only slightly for perturbations in the birth and death rates. This provides experimental evidence that the determined optimal introduction rate is fairly robust to noise the death and birth rates.

4 Conclusions and Future Work

The mathematical analysis and improvement of models of the TYC eradication strategy are essential to understanding the efficacy of the strategy as a control. This is especially important in light of recent field studies of reproductivity and survivability of introduced Trojan supermale populations of brook trout (Salvelinus fontinalis), in the Big Lost River basin in south-central Idaho [14]. In this paper, a modification to the classical model of the TYC strategy is proposed that includes intraspecies competition between males and supermales for female mates and a strong Allee effect in the female population. It is shown that the dynamical system exhibits an extinction and recovery equilibrium solution for realistic parameter values of the carrying capacity and birth and death rates. In such case, the equilibrium solutions are shown to be asymptotically stable and hence the goal of an effective TYC strategy is to push the wild-type populations toward the basin of attraction of the extinction state.

An optimal introduction rate of supermales is determined through optimal control theory and was chosen to minimize an objective function that measures the total amount of wild population and introduced supermales. The optimal introduction rate for the classical model was given in [18]. Here, we compare the influence of including the intraspecies competition for mates and the Allee effect on the optimal introduction rate. It is determined that the optimal introduction rate yields a significantly smaller objective function value as compared to the classical model. This suggests that it is important to include intraspecies competition and the Allee effect to appropriately determine the overall cost of the TYC eradication strategy. In addition, this indicates that the overall cost of the strategy is smaller than previously predicted [18]. Lastly, a stochastic model is proposed to investigate the influence of noise in the birth and death rates on the objective function’s value. These results show the sensitivity in the objective function’s value in light of perturbations to the birth and death rates. Numerical results indicate that the optimal introduction rate is indeed robust to noise in the birth and death rate. This indicates that the optimal introduction rate is not greatly influenced by noise in death and birth rates.

Acknowledgements

JL and RP would like to acknowledge valuable support from the NSF via DMS-1715377 and DMS-1839993. MB, SB, HC, and TG would like to acknowledge valuable support from the NSF via DMS-1715044.