Modelling competition between hybridising subspecies

Journal of Theoretical Biology 486 (2020) 110072

Contents lists available at ScienceDirect

Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/jtb

Modelling competition between hybridising subspecies Nicholas J. Beeton a,∗, Geoffrey R. Hosack a, Andrew Wilkins b, Lawrence K. Forbes c, Adrien Ickowicz a, Keith R. Hayes a a

CSIRO, 3 Castray Esplanade, Battery Point, TAS 7004, Australia CSIRO, Queensland Centre for Advanced Technologies, 1 Technology Court, Pullenvale, QLD 4069, Australia c School of Mathematics and Physics, University of Tasmania, Australia b

a r t i c l e

i n f o

Article history: Received 19 July 2019 Revised 15 October 2019 Accepted 2 November 2019 Available online 7 November 2019 Keywords: Dynamical systems analysis Lotka-Volterra Climate change Invasive species Mate choice

a b s t r a c t The geographic niches of many species are dramatically changing as a result of environmental and anthropogenic impacts such as global climate change and the introduction of invasive species. In particular, genetically compatible subspecies that were once geographically separated are being reintroduced to one another. This is of concern for conservation, where rare or threatened subspecies could be bred out by hybridising with their more common relatives, and for commercial interests, where the stock or quality of desirable harvested species could be compromised. It is also relevant to disease ecology, where disease transmission is heterogeneous among subspecies and hybridisation may affect the rate and spatial spread of disease. We develop and investigate a mathematical model to combine competitive effects via the Lotka-Volterra model with hybridisation effects via mate choice. The species complex is structured into two classes: a subspecies of interest (named x), and other subspecies including any hybrids produced (named y). We show that in the absence of limit cycles the model has four possible equilibrium outcomes, representing every combination: total extinction, x-dominance (y extinct), y-dominance (x extinct), and at most a single coexistence equilibrium. We give conditions for which limit cycles cannot exist, then further show that the “total extinction” equilibrium is always unstable, that y-dominance is always stable, and that the other equilibria have stability depending on the model parameters. We demonstrate that both x-dominance and coexistence are achievable under a wide range of parameter values and initial conditions, which corresponds with empirical evidence of known competing-hybridising systems. We then briefly examine bifurcation behaviour. In particular, we note that a subcritical bifurcation is possible in which a “catastrophic” transition from x-dominance to y-dominance can occur, representing an invasion event. Finally, we briefly examine the common complication of time-varying carrying capacity, showing that such a case can make coexistence more likely. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction The geographic niches of many species are currently undergoing rapid contraction, expansion or wholesale movement through a mixture of global climate change (Walther, 2010) and the introduction of species into new environments, either accidentally (Simberloff et al., 2013) or deliberately for purposes such as aquaculture (Bostock et al., 2010). Sometimes these processes bring together genetically compatible but previously spatially separated subspecies, and thus raises the prospect of potential hybridisation. Hybridisation, as a natural or anthropogenically induced process, is an important ecological process with implications for managing threatened and endangered species (Allendorf et al., 2001).

∗

Corresponding author. E-mail address: [email protected] (N.J. Beeton).

https://doi.org/10.1016/j.jtbi.2019.110072 0022-5193/© 2019 Elsevier Ltd. All rights reserved.

The effects of hybridisation, however, are also important for managers seeking to use genetic control methods to suppress or replace populations of harmful species (Menge et al., 2005). The potential consequences of hybridisation between genetically compatible subspecies depend on the biology of the subspecies, the nature of their interaction and the effects of other forces that act to suppress or enhance the population of either species. In general, however, two outcomes are possible: 1. one or more subspecies are out-bred or out-competed, and go extinct in the shared niche; 2. all subspecies persist, despite competition or hybridisation pressures. Rhymer and Simberloff (1996) note that Outcome 1 may arise when hybrids outcompete one of the subspecies, which is a process of particular concern for a rare locally adapted subspecies that

2

N.J. Beeton, G.R. Hosack and A. Wilkins et al. / Journal of Theoretical Biology 486 (2020) 110072

confronts a more abundant invading subspecies. However, Outcome 1 can even result in the extinction of subspecies through outbreeding depression, where the combination of high hybridisation rates and reduced fitness of hybrids leads to loss of the locally adapted subspecies (Rhymer and Simberloff, 1996). Recognising which of the above two outcomes are possible and which are desirable or undesirable from a management standpoint is vital. In conservation settings, Outcome 1 may be undesirable for a rare or threatened subspecies, particularly if the shared niche represents a substantial portion of the subspecies’ overall niche. For example, the ruddy duck Oxyura jamaicencis is threatening the globally endangered white-headed duck Oxyura leucocephala in such a way, resulting in wide-scale eradication efforts across Europe (Muoz-Fuentes et al., 2007; Lovette, 2016). Conversely, Outcome 1 may be desirable if the subspecies facing extinction is an invasive pest. In commercial settings, the existence of a desirable harvested species could be threatened in Outcome 1, but Outcome 2 may also be problematic if it decreases the abundance or threatens the actual or perceived quality of a harvested species (Ford and Myers, 2008). In epidemiological settings, sympatric vector subspecies may have different levels of capacity to transmit vector– borne disease, such as human malaria by biting mosquitoes (Niang et al., 2019). If the vectorial capacity of hybrids is greater than that of endemic subspecies then Outcome 2 is most undesirable. Outcome 1, however, may also be problematic if hybridisation facilitates the spread of (for example) insecticide resistant alleles and contributes to the replacement of a susceptible subspecies with a resistant one (Anderson et al., 2018). Given the diversity of species for which this process is occurring or could potentially occur, developing a theoretical framework to quantify the conditions likely to produce each of the potential outcomes is a logical first step. To achieve this, we look to the well-established modelling literature for both competition and hybridisation. The competitive Lotka-Volterra model, as used to model two competing distinct species, has been well studied and tested against field data in both deterministic and stochastic settings (May and Mac Arthur, 1972). Less attention, however, has been given to its use in modelling competing subspecies, or hybridising species in a species complex. These systems are often modelled from a population genetics perspective, concentrating on the relevant alleles and treating the species as separate only in the heritability and fitness of each allele. This approach, for example, has been applied to different strains of Anopheles gambiae mosquitoes (Menge et al., 2005), and to the edible frog (Pelophylax esculentus) complex (Graf, 1986; Plótner and Grunwald, 1991). These types of studies may also be complemented by laboratory or field studies on competitive behaviour in the same or similar species complexes (Semlitsch and Ryer, 1992; Schneider et al., 20 0 0). The effects of hybridisation and competition within the Pelophylax esculentus species complex, and between the distantly related Atlantic salmon (Salmo salar) and brown trout (Salmo trutta), have been previously modelled (Quilodrán et al., 2014; 2015). In these studies, stage-structured models of density dependent competition were based on a version of the discrete Ricker model:





Ni,t+1 = Ni,t 1 + ri exp −(Ki −



 αi j N j,t )/Ki

,

j

where Ni,t is the number of animals of species i at timestep t, ri is the intrinsic growth rate, Ki the carrying capacity of species i and α ij the inter–species competition term (with αii = 1). The Pelophylax esculentus species complex has also been examined with density dependence modelled by a version (Currat et al., 2008; Quilo-

drán et al., 2018) of the closely related logistic map given by



Ni,t+1 = Ni,t 1 + ri (Ki −





αi j N j,t )/Ki .

j

In all of these models hybridisation occurs by a proportion of fecundity being allocated to interbreeding, depending on the relative abundances of each species and some weighting factor representing the increased difficulty of interbreeding. These approaches have allowed researchers to gain an understanding of introgression and invasion of hybridising species. The above discrete models usually (implicitly) assume some degree of synchronicity in the life histories of the subspecies: cases where there is a direct correspondence between discrete and continuous time mechanistic models are limited (Geritz and Kisdi, 2004). In this paper, we assess the ramifications of including both competition and hybridisation simultaneously. We develop a continuous-time ordinary differential equation (ODE) model that is appropriate for populations with overlapping generations that both compete and interbreed. Similar to the above studies that invoked discrete–time models, we construct the ODE model for each subspecies without age or stage structure. Fecundity and mortality levels therefore depend only on the population size and not its age or stage structure, and we can model the population of each subspecies and sex using a single variable only. The model incorporates mate choice, which is shown to have a critical role in determining the long term behaviour of the dynamic system. 2. Competitive Lotka-Volterra model We begin by considering the case of competition but no interbreeding. Denote the populations of the two species by x and y, changing over time t. All rate parameters are then defined in terms of the unit of time (e.g. per year if time t is defined in years). We shall use the competitive Lotka-Volterra model, with carrying capacities Kx and Ky , and interspecific competition set by α xy and α yx . The dynamics are:





x + αxy y dx = −μx + 1 − dt Kx



dy = −μy + dt

1−

  γx x,

αyx x + y Ky

 γy y.

(1)

(2)

This system of equations differs from the competitive LotkaVolterra model in that we have density-dependent and competitive effects acting only on fecundity rates for each species (γ x and γ y ) with mortality rates independent of density (μx and μy ). In the classical model, competition acts on both, though the two models are equivalent after scaling. Explicitly, let Kˆx = Kx (1 − μx /γx ), Kˆy = Ky (1 − μy /γy ), and define intrinsic growth rates rx = γx − μx and ry = γy − μy . Then we recover the original competitive LotkaVolterra model

dx = dt dy = dt



1−

1−



x + αxy y rx x, Kˆx

αyx x + y Kˆy



ry y.

Despite this equivalence, we keep our model in the form of Eqs. (1) and (2): the later addition of hybridisation to the model will require that fecundity and mortality be modelled separately. Biologically, the parameters μx , μy , γ x , γ y , Kx , and Ky are all positive, while α xy and α yx are non-negative. Also, the constraints μx < γ x and μy < γ y ensure the populations do not simply die out.

N.J. Beeton, G.R. Hosack and A. Wilkins et al. / Journal of Theoretical Biology 486 (2020) 110072

Of particular interest is the single coexistence equilibrium (x0 , y0 ) where both x0 and y0 are positive. This occurs when:



 μx x0 + αxy y0 = 1 − K and γx x

 μ αyx x0 + y0 = 1 − y Ky γy with solutions



−αxy x0 γx /Kx , −y0 γy /Ky

dyF = dt



λ +

γx x0 Kx

+

γy y0 Ky

 γγxy λ + (1 − αxy αyx ) x y 0 0 = 0. Kx Ky

The equilibrium can only be reached when x0 and y0 are positive, which constrains the model parameters. The characteristic γx x0 γy y0 polynomial coefficient + is therefore positive. The conKx Ky dition for stability, where the real parts of all eigenvalues are negative, then reduces to the classic condition α xy α yx < 1. 3. Sex and subspecies structured model To incorporate interbreeding between subspecies, we introduce sex structure to Eqs. (1) and (2) and model the breeding behaviour between males and females of each subspecies. We assume that each female will mate with a constant number of males; that is, females are always able to find a mate even when males are at low abundance. This assumption is used to model frequency-dependent transmission of sexually transmitted diseases in a pathogen transmission framework (McCallum et al., 2001) where reproduction is limited by the number of matings per female, as opposed to the traditional “mass-action” or density-dependent model where reproduction is limited by abundance of both males and females. We incorporate hybridisation by assuming that the probability per encounter of mating with the other subspecies is w < 1 times that of mating with a conspecific (e.g. Graf, 1986; Plótner and Grunwald, 1991). Each mating is presumed equally likely to produce offspring. As a result, every female xF of subspecies x will produce a constant number (γ x ) of offspring, so the total fecundity rate due to xF is γ x xF . The proportion of those offspring fathered by a male xM of xM subspecies x is . Following this logic provides a continxM + wyM gency table (Table 1) of total fecundity rates based on each pairing. We further assume a constant probability pF of any offspring being female, probability pxy of a female x and male y pairing being of subspecies x, and probability pyx of a female y and male x

Table 1 Fecundity rates based on subspecies of the father (xM vs yM ) and mother (xF vs yF ).

xF yF

(3)



xM

yM

xM γx x F xM + wyM wxM γy y F wxM + yM

wyM γx x F xM + wyM yM γy y F wxM + yM



dyM = dt

1−

1−

αyx x + y

(4)

 p F B y − μy y F ,

Ky

αyx x + y

(5)

 ( 1 − p F ) B y − μy y M ,

Ky

(6)

where x = xF + xM , y = yF + yM and

Bx =

wxM γy yF xM γx xF wyM γx xF + pxy + pyx xM + wyM xM + wyM wxM + yM

By =

yM γy yF wxM γy yF wyM γx xF + (1 − pxy ) + (1 − pyx ) . (8) wxM + yM xM + wyM wxM + yM

and the corresponding characteristic polynomial is 2





The Jacobian taken at this equilibrium is

−x0 γx /Kx J= −αyx y0 γy /Ky



x + αxy y dxM = 1− ( 1 − p F ) B x − μx x M , dt Kx





pairing being of subspecies x. We then obtain the following series of equations:

x + αxy y dxF = 1− p F B x − μx x F , dt Kx



 μy μx x0 = 1− K − αxy 1 − K , 1 − αxy αyx γx x γy y



 μy 1 μx y0 = 1− K − αyx 1 − K . 1 − αxy αyx γy y γx x 1

3

(7)

Summing Eqs. (3) and (4), and Eqs. (5) and (6) produces:





x + αxy y dx = 1− Bx − μx x, dt Kx dy = dt



1−

αyx x + y



Ky

By − μy y.

Note that the female-male relationship in both subspecies with interspecific density dependence and hybridisation is of the form

dF = p F f ( F , M ) − μF , dt dM = (1 − pF ) f (F , M ) − μM, dt where F is the number of females and M the number of males, with some function f depending on these numbers. Substituting F = pF Z and M = (1 − pF )Z into these equations results in them both being described by

dZ = f ( pF Z, (1 − pF )Z ) − μZ. dt As a result, substituting xF = pF x, xM = (1 − pF ) x, yF = pF y and yM = (1 − pF ) y consistently reduces the four equations to two equations. In other words, if the female and male populations for subspecies x are given initial conditions xF (t0 ) = pF X and xM (t0 ) = (1 − pF )X for some constant X, and similarly for subspecies y with some constant Y, their sex ratio will remain constant for all subsequent t. For alternative initial conditions, it can be shown (see Appendix A) that the sex ratio within a simple version of Eqs. (3)– (6) converges to the same constant ratio exponentially fast. We assume that such convergence has occurred and the sex ratio is stable. The two equations are then:





x + αxy y dx x = −μx x + pF 1 − γx x dt Kx x + wy  wy wx + pxy γx x + pyx γy y , x + wy wx + y



dy = −μy y + pF 1 − dt + (1 − pxy )

 αyx x + y  Ky

y γy y wx + y

(9)



wy wx γx x + (1 − pyx ) γy y . x + wy wx + y

(10)

4

N.J. Beeton, G.R. Hosack and A. Wilkins et al. / Journal of Theoretical Biology 486 (2020) 110072

4. Competition and one-way hybridisation model

4.1. Total extinction

An important first test of the biological relevance of a model that incorporates both competition and hybridisation between subspecies is whether it can reasonably represent interactions that are known to happen in the wild. In particular, there are multiple documented examples of closely related species that have strongly overlapping niches (suggesting substantial competition) and significant levels of hybridisation. Cases include the edible frog mentioned above, as well as the An. gambiae s.l. (sensu lato) species complex in mosquitoes, which includes the malaria-carrying subspecies An. gambiae s.s. (sensu stricto) and An. coluzzii. As the subspecies in these complexes are known to coexist in the wild over the long-term, a realistic model should include a stable coexistence equilibrium at relevant regions of model parameter space. In this paper, we will concentrate on the case pxy = pyx = 0 in Eqs. (9) and (10). Biologically, this means that all successful crossspecies matings will result in offspring of subspecies y (the choice is arbitrary - we could equally have done this for x). For example, subspecies x may be endangered with high priority for protection and so hybridisation is perceived as a threat to population recovery, or y may denote a suite of hybridising subspecies that interact with the focal subspecies x. Such a model gives the maximum advantage to one subspecies (y), and will be more likely to lead to exclusion of a subspecies (x) as opposed to coexistence. As a result, it provides a conservative test for the realism of our model: if we can show that coexistence is a reasonable outcome where pxy = pyx = 0, this will also be the case for more nuanced hybridisation models. Note that if we instead set pxy = 1 and pyx = 0, so that female subspecies determines the offspring subspecies, then we obtain the original competitive Lotka-Volterra model (Eqs. (1) and (2)). This is because each female in a subspecies produces a constant number of offspring, each of which shares its mother’s subspecies, so hybridisation will not affect the population numbers. We accordingly set pxy = pyx = 0 and redefine γ x and γ y to be the number of females produced per female (pF γ x and pF γ y ) instead of total offspring, resulting in our final model:

This is represented by the solution x = 0 and y = 0. At this point, the Jacobian is





x + αxy y dx = −μx + 1 − dt Kx



x x + wy

  γx x

= f (x, y ) x,



dy = −μy + dt = g(x, y ) y.

(11)

1−

 αyx x + y  Ky

γy +

wx γx x + wy





J=

f ( 0, 0 ) 0

0 g( 0 , 0 )





=

zγx − μx 0

0



γy + wzγx − μy

(13) x where z = lim . x→0+ , y→0+ x + wy The Jacobian is undefined at the critical point, and the limit as x → 0+ , y → 0+ can be seen to vary depending on the direction of approach. However, as we are only interested in the non-negative quadrant, this turns out not to matter: stability requires eigenx value λ1 = zγx − μx < 0 ⇒ z < μ γx and λ2 = γy + wzγx − μy < 0 ⇒ μ −γ

z < wy γx y . The second inequality requires z to be negative (since γ y > μy ) but this never occurs for all relevant solutions which approach (0,0) from the positive quadrant. As a result, this critical point is always either a source (all eigenvalues with positive real part) or a saddle (some eigenvalues with negative real part). Biologically, at least one of any small population of x and y will initially grow exponentially. 4.2. y-dominance (x extinction) In this case, x = 0 and



μ yc =Ky 1− γ y



y

. Assuming w > 0 (mating be-

tween x and y does occur, and produces y offspring) the Jacobian at this point is



J=

f ( 0, yc ) gx ( 0 , yc ) yc

0 gy ( 0 , yc ) yc





=

−μx gx ( 0 , yc ) yc

0



μy − γy

where gx and gy denote partial differentiation of g by x and y respectively. Each diagonal term is an eigenvalue of the matrix, so when w > 0 this point is always stable (all eigenvalues with negative real part) as μx is positive and γ y > μy . Biologically, for very small x, the chance of an x mating to produce another x is very small (when w > 0) so x disappears exponentially, while y assumes its equilibrium value of yc . Species y can then be considered as the established resident and x the (unsuccessful in the stable case) invader; in this case y benefits from priority effects (Weidlich et al., 2016). For w = 0, this point is a saddle as any small x will initially grow exponentially away from x = 0 with rate γx − μx . 4.3. x-dominance

y

Here,

(12)

When studying the dynamics of this system, we can without loss of generality set (for example) Kx = 1 and γx = 1 to nondimensionalise the time and population scales. In this case, we interpret any population numbers as a proportion of the carrying capacity of species x, and any time numbers in relation to the average generation time of species x. In other words, the behaviour of the system depends only on the generation times of the species involved (e.g. weeks for mosquitoes to years or decades for large mammals), and not on any intrinsic measure of time. Only the non-negative quadrant of the plane, x ≥ 0 and y ≥ 0, is physically meaningful. Only this quadrant needs to be considered: by dx dy Eq. (11), = 0 on the line x = 0, and by Eq. (12), = 0 at y = 0 dt dt (L’Hôpital’s rule can be used to confirm this at x = y = 0). Thus, solution trajectories cannot cross the x- or y-axes. Solving for f (x, y ) x = 0 = g(x, y ) y gives four types of equilibria for Eqs. (11) and (12), as discussed below.

J=

⎡



μ xc =Kx 1− γ x x

f x ( xc , 0 ) xc 0

=⎣



and y = 0. At this point, the Jacobian is

f y ( xc , 0 ) xc g( x c , 0 )

μx − γx

−μy +

0



1 − αyx

Kx Ky

⎤  fy (xc , 0 ) xc ⎦ μx 1− (γy + wγx ) γx



where fx and fy denote partial differentiation of f by x and y respectively. As μx < γ x , one eigenvalue is guaranteed to be negative, so stability requires



−μy +

1 − αyx

Kx Ky



1−

μx γx



(γy + wγx ) < 0

If αyx = 0, i.e. that x has no competitive effect on y, then the condition implies γ y < μy (false by assumption) and the point is a saddle. Biologically in this case, any small y will initially grow exponentially away from y = 0.

N.J. Beeton, G.R. Hosack and A. Wilkins et al. / Journal of Theoretical Biology 486 (2020) 110072

5

Otherwise the point is sometimes stable and sometimes a saddle depending on parameter values. As α yx is increased while keeping all other parameters fixed, there will come a point where this equilibrium becomes stable. Biologically, the competitive impact of x on y becomes large enough to drive the y population to extinction. When x is resident, it benefits from priority effects and will exclude an invasion of y. To illustrate: where αyx = 1 and Kx ≥ Ky , the inequality implies that

 γy γx γx > 1+ w μx μy γy

so γ x /μx > γ y /μy is a necessary but not sufficient condition for stability in this case. Biologically, for x to out-compete and drive a small population of y to extinction, two things need to happen: the growth rate for x over the time scale of a single lifespan (γ x /μx ) must be sufficiently higher than the corresponding rate for y, and the carrying capacity for x must be sufficiently large relative to that of y. These strict conditions are not surprising, considering that any x − y breeding is producing more y. 4.4. Coexistence Where x0 and y0 are both positive, we can determine the relationship between x0 and y0 using Eq. (11), and then substitute this relationship into Eq. (12) to determine the value of x0 and thus y0 . This process gives us the equations:

p( x 0 ) = 0

(14)

Kx (γx − μx )x0 − γx x20 Kx μx w + αxy γx x0

Fig. 1. Illustrative diagram of potential stable manifold configurations. In both figures the shaded region is T, the dotted curve is , and a and b are stable points within the manifold.

J21 =

−αyx Ky



γy +

wx0 γx x0 + wy0



w2 y0 γx + 1− (x0 + wy0 )2 J22 = −



αyx x0 + y0



Ky

y0

y0 γy y0 γx wx0 − (Ky w + x − αyx wx ) Ky Ky (x0 + wy0 )2

Eq. (18) will be used in the following section to qualitatively describe model behaviour in the phase plane. 5. Possible model behaviours

a3 = wγy γx2 − Dγx B + αyxCB2

To summarise our findings from Section 4 above, we have shown that in all cases the critical point at the origin (total extinction) is unstable (some or all eigenvalues with positive real part), and that there is one stable attractor on the y axis (where x goes extinct and y is dominant). The critical point on the x axis (where x is dominant) has stability which changes depending on parameter values. Stability of the remaining coexistence points, determined by Eqs. (14) and (15), may be determined using the eigenvalues of J in Eq. (18). As calculation of the coordinates and corresponding eigenvalues is difficult, further headway can more easily be made by a consideration of the possible phase-plane behaviour. To facilitate this, a trapping region is introduced.

a2 = −2 wγx γy P + D(BP − γx A ) + 2ABC αyx + wγx QB − RB2

Theorem 1. Any solution trajectories entering the region

y0 =

(15)

where p(x0 ) is a cubic, obtained after substantial algebra. For brevity, we define intermediate constants

A = Kx μx w,

B = αxy γx ,

P = Kx (γx − μx ),

Q = Ky (γy − μy )

C = wγx + γy , (16)

and then the further two constants

D = αyx wγy + C,

R = Q + Ky wγx .

(17)

We use these constants (16), (17) to define the coefficients

2

T = (x, y ) ∈ R2 : 0 ≤ x ≤ xc = Kx 1 −

2

a0 = RA + wQAP, from which it follows that the equilibrium value x0 satisfies Eq. (14) with

p( x ) = a 3 x 3 + a 2 x 2 + a 1 x − a 0 . The function p(x) has real coefficients, so there are either 1 or 3 distinct real solutions, though some of these may lie in the nonphysical x < 0 or y < 0 regions. Setting w = 0 gives the LotkaVolterra model result, as expected. When f (x0 , y0 ) = g(x0 , y0 ) = 0 the Jacobian takes the form



J=

J11 J21

where

J11

J12 J22



(18)



x0 γx wy0 =− 1− (K − (αxy − w )y0 ) Kx (x0 + wy0 )2 x

J12 = −

γx x20

Kx (x0 + wy0 )2

(w(Kx − x0 ) + αxy x )







a1 = wγy P + DAP + αyxCA + wQ (γx A − BP ) − 2ABR 2

y≥0

and Ky 1 −

μy γy



 μx , γx

− αyx x ≤ y ≤ Ky − αyx x



never leave that region, making T by definition a trapping region. Note that T is defined to contain its own boundary ∂ T. Proof. Trajectories never leave the non-negative quadrant. For x > xc and y ≥ 0, applying the assumption γ x > μx in Eq. (11) gives dx/dt < 0. For y > Ky − αyx x (and x ≥ 0 and y > 0) Eq. (12) implies that dy/dt < 0. For y < Ky (1 − μy /γy ) − αyx x (and x ≥ 0 and y > 0), Eq. (12) implies that dy/dt > 0. Therefore, all solution trajectories that enter T never leave it.  The trapping region and solution trajectories are shown in Fig. 1. The trapping region also contains all coexistence equilibria and the y-dominance equilibrium: Theorem 2. The trapping region contains all the physical equilibrium points corresponding to coexistence, and the only critical points

6

N.J. Beeton, G.R. Hosack and A. Wilkins et al. / Journal of Theoretical Biology 486 (2020) 110072

that lie outside the trapping region are the origin and possibly the x-dominance point (x, y ) = (xc , 0 ). Proof. By the proof of Theorem 1, all points in the strictly positive quadrant (x > 0 and y > 0) and outside T have dx/dt = 0 or dy/dt = 0 (or both), so can never be at equilibrium. The edges of the nonnegative quadrant may contain critical points. Along the line x = 0, there is the point at the origin, and the y-dominant stable equilibrium (x, y ) = (0, yc ) that lies within T. By Eq. (15) there are no other solutions of f = 0 = g with x0 = 0. Similarly, along y = 0, Eq. (15) implies that the only solutions of f = 0 = g are at the origin and the point (xc , 0).  Note that the point (xc , 0) may or may not lie within T depending on parameter values, and as such the condition x ≤ xc may not be required to help define T. Armed with these results, we can make inferences about the stability structure of the system. Theorem 3. If there are n distinct stable points in T, and no stable or unstable limit cycles in its interior int(T), there must be at least n − 1 saddles in int(T). Proof. By the stable manifold theorem, each stable point a ∈ T has an associated two-dimensional stable manifold A. As A must contain at least an open disc around a, at least some of this stable manifold must be in int(T). By the Poincaré-Bendixson theorem, every point in int(T) is either within the stable manifold of one of these stable points, or that of a stable limit cycle – but we have assumed that no such limit cycles exist. Then if more than one stable point exists (say, a and b), every pair of stable manifolds (corresponding manifolds A and B) must be separated by a one-dimensional boundary that does not self-intersect (i.e. a simple curve), which itself is an invariant set : any point that starts on  cannot leave towards either stable point, so must proceed along it. There are then two possibilities: 1. The curve is open, and terminated at both ends by ∂ T. As solutions along the curve must move away from the boundaries of the trapping region,  contains at least one critical point towards which all solutions in  equilibriate (see Fig. 1a). This is a saddle in int(T), as it is stable along  but unstable in any other direction: solutions that leave  will enter either A or B and head towards either a or b. 2. The curve is closed. If all solutions head in one direction along the curve,  is an unstable limit cycle (excluded by assumption). Any change of direction along the curve requires at least one stable equilibrium in  (that is, a saddle in int(T)) as in the previous case, with a corresponding unstable equilibrium (a source point in int(T)) (see Fig. 1b). If there are n distinct stable points in T, then int(T) will consist of n corresponding manifolds separated by at least n − 1 boundaries. As each will behave like  above, there must be at least n − 1 corresponding saddles and potentially additional source points in int(T).  The stable y-dominance point is known to lie in T. Theorem 3 states that if there were 2 or more distinct stable coexistence points (so 3 or more stable points in T including the y-dominance point) there must be 2 or more saddles in int(T), which is impossible, given there are at most 3 roots to the cubic p(x ) = 0. Hence, assuming no limit cycles in int(T), there can only be either none or one stable coexistence solutions. Along with the none or one stable x-dominance equilibrium, there are four possible stability configurations as shown in Fig. 2. For stable coexistence to be possible, p(x ) = 0 must therefore have 3 real roots, which requires its discriminant to be non-

Fig. 2. Phase portrait sketches. A red diamond indicates a stable attractor, while a black dot indicates an unstable attractor (a source or saddle). The blue lines show flow directions. Fig. 2a to 2 b show the only possible configurations of stable attractors in the model in the absence of limit cycles. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

negative. Hence

disc( p) < 0 ⇒ stable coexistence impossible. The discriminant is a function of the model parameters, so checking its sign allows easy assessment of the possibility of coexistence. As mentioned, these results require us to prove that no limit cycles exist within int(T). The following theorem demonstrates that this is true for at least some regions of parameter space. Elsewhere in parameter space, modifications to the proof may produce the desired result, or numerical study may demonstrate a lack of obvious limit cycles. Theorem 4. There can be no limit cycles anywhere in R2 if the three conditions 1. μx + μy < γy 2. wα yx γ y < μx 3.

αyx  Ky



wγx + γy <

γx Kx

+

2 wγy αyx Ky

.

on the parameters all hold simultaneously. Proof. Dulac’s Theorem (Edelstein-Keshet, 1988; Strogatz, 1994) states that, for two-dimensional systems of the general form

dx = F (x, y ) dt dy = G(x, y ) dt

N.J. Beeton, G.R. Hosack and A. Wilkins et al. / Journal of Theoretical Biology 486 (2020) 110072

and for any differentiable function D(x, y), there can be no limit cycles (no closed orbits in the phase plane) if the function

S(x, y ) =

∂   ∂   DF + DG ∂x ∂y

(19)

is not zero and does not change sign anywhere in the region of interest. It can be shown with some experimentation and algebra that the choice

D(x, y ) =

yμx /μy −1 (x + wy ) ≥ 0 x2

gives the inequality

S(x, y ) < D(x, y ) + +





γy +

wγx x x + wy

x x + wy

 x2

x + wy

 μ

x

γy



−μx + wαyx γy

 wγ α x yx Ky

μy  γy 

−1+

−

γx Kx

+

αyx γy Ky

−

2  wγy αyx . Ky

7

intraspecific competition coefficients to unity. Interspecific competition was considered symmetric and ranged from zero, which signified that some subspecies occupied different habitats, to a magnitude of 0.5 (in other case studies, even values greater than 1 are possible; see Begon et al., 2009, Chapter 4). Equating the interbreeding success rate of Quilodrán et al. (2015) to our probability per encounter of mating with the other subspecies w, we note that Quilodrán et al. (2015) investigated rates of interbreeding success that ranged from zero (reproductive isolation between certain combinations of subspecies) to one (panmictic mating). We also note that estimating mate choice parameters in non-random mating models is an active field of research and that obtaining such information from case studies is a difficult task (Carvajal-Rodrguez, 2018), and as such one that we have not attempted here. To graphically demonstrate qualitative long term behaviour of the dynamic system, we numerically simulate solutions across a range of parameter solutions using the programming language R

(20)

everywhere in int(T). The theorem then follows by making the coefficients of the three linearly independent functions on the right-hand side of the inequality in Eq. (20) all independently negative. This then guarantees that S(x, y) < 0, from which it follows by Dulac’s Theorem that no closed orbits can occur within int(T). As a result, limit cycles cannot occur in int(T). Further, every closed trajectory (as its own trapping region) contains a stationary point (by the Poincaré-Bendixson theorem). However, the only stationary points not in int(T) are on the x- or y-axes, so any closed trajectories in R2 outside int(T) would travel outside the non-negative quadrant, which is contradicted by previous results. As a result, limit cycles cannot occur anywhere in R2 .  6. Examples In Section 4, we showed that: (1) the y-dominance equilibrium (with x = 0) is always stable, regardless of the parameters (assuming w > 0); (2) the x-dominance equilibrium (with y = 0) is only stable if the competitive impact of x on y (the parameter α yx ) is large enough, the carrying capacity Kx is large enough, and/or the growth rate of x is suitably higher than y; and (3) in the absence of limit cycles, there is only one coexistence equilibrium (with x = 0 and y = 0) that is possibly stable. The four possible combinations of stable equilibria that satisfy all three results are depicted in Fig. 2. Estimation of density dependence in Lotka-Volterra systems from experimental or observational studies is a challenging task (Bender et al., 1984), and we are unaware of studies that have jointly estimated life history parameters of density independence, density dependence and the probability of hybridisation within species complexes. To investigate possible parameter values, we examine the parameterisation developed by Quilodrán et al. (2015) for a stage-structured discrete-time process model of a species complex of waterfrogs. In that study, a Ricker form of intraspecific and interspecific density dependence was applied to survival of larvae to adult stages. The Ricker model of density dependence is one of at least three forms of discretisation of the continuous-time logistic model (Turchin, 2003). Extending this Ricker from of discretisation to a multispecies Lotka–Volterra system, Hofbauer et al. (1987) argued that the Ricker-type discretisation method is the most defensible option because, in terms of defining the long term average abundances for the dynamic system, the parameters for both the Lotka-Volterra system and the Ricker form of discretisation share the same meaning. As a guide for parameterisation then, we note that Quilodrán et al. (2015) set

Fig. 3. Graphical illustration of domains of attraction for three sets of parameter combinations chosen to demonstrate qualitative differences in the phase plane. Yellow represents cases where either coexistence or y-dominance is possible; blue where x-dominance or y-dominance (but not coexistence) is possible; purple where only y-dominance is possible; white where all three are possible. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

8

N.J. Beeton, G.R. Hosack and A. Wilkins et al. / Journal of Theoretical Biology 486 (2020) 110072

Fig. 4. The stable (bold) and unstable (dashed) equilibria for y with varying w, with α = 0.95 and other parameters as in Fig. 3b.

(R Core Team, 2018) and package deSolve (Soetaert et al., 2010) The different combinations of parameters are chosen to illustrate insights about the model. For instance, Fig. 3 shows potential behaviours across a range of values for a competition parameter α = αxy = αyx (symmetric interspecific competition) and hybridisation parameter w, where the other parameters are fixed at γx = 1 (as defined in Section 4) and Kx = Ky = 1. Fig. 3a shows a case where mortality for x and y is equal: y is always present for α < 1 in this case, but the presence of x depends on both α and w. Reducing the mortality for subspecies x increases the region where x can exclude y as well as the region of coexistence, as seen in Fig. 3b. In a situation where x has a strong advantage over y in both mortality and fecundity (Fig. 3c), we lose the y-dominance only region but gain a region where all three types of result are possible. Boundaries between regions in Fig. 3 represent bifurcations, as the system’s topology changes immediately as the parameters’ values travel over the boundary. Subcritical bifurcations are of particular biological interest, as they represent a “catastrophic” change (see Fig. 4): a stable critical point (here, x-dominance) becomes unstable as the bifurcation parameter passes a threshold (here w, at about 0.12), and any trajectories starting near there will necessarily move to another critical point (here y-dominance). Biologically, this means any small-amplitude incursions of subspecies y will fail unless w exceeds the threshold, at which point any incursion will succeed. In the absence of limit cycles, x will then ultimately become extinct or extirpated in the long term. This represents a rateindependent hysteresis system: reducing w back below the threshold will not return x to long-term dominance, and this in fact cannot happen unless w is near enough to zero for a small-amplitude incursion to succeed. These subcritical bifurcations occur where a “y-dominance stable” region becomes an “x- and y-dominance stable” region or vice-versa. This can be demonstrated using the analysis presented in Boldin (2006, Theorem 2.2): the variable M in her paper is shown to be positive for a subcritical bifurcation. The value of M for our model is given below, and can be solved numerically for any set of parameter values.

M=



+ Ky γx (γx ((−1 + αxy αyx )Kx γx − 2αyx Kx γy + 2Ky γy )

(Ky γx + αyx Kx (−γx + μx ))μy + (αyx Kx − Ky )Ky2 γx3 μ2y )   γx2 (Ky γx + αyx Kx (−γx + μx ))2

In cases where more than one equilibrium is stable, the initial conditions will dictate which the system assumes as t → ∞. Fig. 5 shows several representative phase-plane portraits using parameters selected from Fig. 3. Note particularly that in Fig. 5d we have the unusual case of interspecific competition outstripping intraspecific competition (α > 1; see Begon et al., 2009, Chapter 4). 7. Time-varying k

γx (γy2 (−Ky γx + αyx Kx (γx − μx ))3

+Kx (γx − αxy αyx γx + αyx γy )μx )

Fig. 5. Phase-plane portraits showing directions of trajectories for parameter combinations chosen to demonstrate qualitatively different topologies in the domains of attraction. Here Kx = 1 = Ky , αxy = α = αyx and γx = 1. The other parameters are different for each figure, as noted. In each figure, a red dot indicates a stable attractor, while a black dot indicates an unstable attractor (a source or saddle). Filled-in areas represent domains of attraction for y-dominance (purple), x-dominance (blue) and coexistence (yellow). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

 (21)

Local carrying capacities are seldom constant over time in biology, and often vary wildly, particularly in the case of mosquitoes. Adding a time-varying parameter to our model makes the analyses described above impractical, but we can perform numerical simulations to gauge the potential behaviours. Here we use a square wave model for the variation of carrying capacity for each population, simulating a “high” and “low” period of time, each with

N.J. Beeton, G.R. Hosack and A. Wilkins et al. / Journal of Theoretical Biology 486 (2020) 110072

9

period P (e.g. annual seasonality, Stone et al., 2007):



β (t ) =

Kh t ∈ {[0, P ), [2P, 3P ), [4P, 5P ), . . .} Kl t ∈ {[P, 2P ), [3P, 4P ), [5P, 6P ), . . .}

(22)

where we replace Kx and Ky in Eqs. (11) and (12) with β (t)Kx and β (t)Ky , respectively. We can determine the behaviour of the system analytically as P approaches 0 by examining Eqs. (11) and (12). Theorem 5. A competing-hybridising model as defined in Eqs. (11)and (12) with β (t) defined as in Eq. (22)behaves like a constant-K 2Kh Kl model within T with β (t ) = Km as P → 0, where Km = ≤ 2Kl . Kh + Kl Proof. Though the introduction of discontinuous β (t) causes

dx dt

dy to also become discontinuous, they are continuous almost dt everywhere and bounded in the region of interest, and thus Riemann integrable. We can then use the Second Fundamental Theorem of Calculus to state that and





x (n + 2 )P − x(nP ) =

 (n+2)P dx(t ) dt dt nP

(23)

for positive integer n, and similarly for y. It follows that





x (n + 2 )P − x(nP ) =

 (n+1)P  (n+2)P dx(t ) dx(t ) dt + dt d t dt nP (n+1 )P

and we can then introduce the new definition of K(t):



Fig. 6. The phase-plane portrait given in Fig. 5c (see corresponding caption for description) with the addition of trajectories for time-varying K with Kh = 1, Kl = 0.2 for each of P = {20, 10, 1, 0} given in green, blue, red and white respectively. Each trajectory starts at (x, y ) = (0.015, 0.1 ), shown as a white dot. Constant-K model equilibria for both K = 1 and K = 0.2 are shown as in Fig. 5, with white dotted lines to illustrate the direct scaling relationship between the two models. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)



x (n + 2 )P − x(nP )

  (n+1)P x(t ) + αxy y(t ) = − μx + 1 − Kh Kx nP





x(t ) γx x(t ) dt x(t ) + wy(t )

  (n+2)P x(t ) + αxy y(t ) + − μx + 1 − Kl Kx (n+1 )P



x(t ) x(t ) + wy(t )

 γx x(t ) dt

(24)

As x, y and their derivatives are bounded within the trapping region T, x(t) and y(t) can be treated as constants there as P becomes small. Then Eq. (24) becomes

 



lim x (n + 2 )P − x(nP )

P→0

= lim



P→0





−μx + 1 −





+ −μx + 1 −

x + αxy y Kh Kx

x + αxy y Kl Kx





x x + wy

x x + wy

  γx xP

   γx xP

and so, dividing both sides by 2P and using the definition of a central derivative,



dx  dt 

=

t=(n+1 )P





− μx + 1 −

=

−μx +

1−

x + αxy y x + αxy y − 2Kh Kx 2Kl Kx

(x + αxy y ) Km Kx

 

2Kh Kl dy  where Km = , and similarly for Kh + Kl dt 



x x + wy

t=(n+1 )P

x x + wy

  γx x

  γx x (25)

. Letting Kh = cKl

for c ≥ 1, it is easily seen that Km = 2Kl c/(1 + c ) ≤ 2Kl . For any given time t and positive constant  , for P <  there exists a positive integer m such that (m + 1 )P ≤ t < (m + 2 )P . We can then complete the proof by choosing n = m in Eq. (25) and taking the limit as  → 0. 

In the other extreme, where P → ∞, the system will originally equilibriate at the relevant stable point for carrying capacity Kh , and then from this starting point it will re-equilibriate at the relevant stable point for carrying capacity Kl . The system will then oscillate between these two solutions. Fig. 6 demonstrates these cases: where P is large in green, small in red, and the theoretical limit where P = 0 in white. An intermediate value where P = 10 is also shown in blue: despite starting in the basin of attraction of (i.e. the region in which all trajectories converge to) y-dominance for the constant-Kh model, the resulting trajectory ends up oscillating between the two constant-K coexistence steady states. This behaviour does not occur for larger P (P = 20 in green), but does occur for smaller P (P = 1 in red and P = 0 in white), demonstrating that the period can determine which equilibrium is reached by the model. By substituting x → Kx, y → Ky, Kx → KKx and Ky → KKy into Eqs. (11) and (12), it can be shown that the system is scaleindependent as long as both subspecies are scaled proportionately. Using this result, it can then be proved (see Appendix B) that if the basins of attraction of the x-dominance (where it exists) and y-dominance equilibria are concave and monotonically increasing for the constant-Kh model, then the corresponding basins of attraction in the varying-K model will be subsets of these. Note that the result in Appendix Appendix B is not limited to Eq. (23), and in fact describes any time-varying β (t) such that β (t) ≤ Kh for all t. In brief, any decrease in the carrying capacity of both species will result in coexistence in some cases where it otherwise wouldn’t (given the conditions hold true). The concavity condition is satisfied for the parameters used in Fig. 6, but a more general proof for parameter space is much more difficult and is not presented here. There is a substantial literature on periodically time-varying parameters in related models (Cushing, 1980; 1986; Rinaldi et al., 1993) that suggests that such a conclusion would be reasonable.

10

N.J. Beeton, G.R. Hosack and A. Wilkins et al. / Journal of Theoretical Biology 486 (2020) 110072

8. Discussion In this paper, we have explored the dynamics of competitive and hybridising interactions between two coexisting subspecies. We have used a simple two-compartment ODE model following well-established modelling assumptions for competition and mate choice. The model predicts various outcomes. There is a ydominance state, in which the subspecies x goes extinct but y persists. This equilibrium is always stable, but as with all stable equilibria in this model, may only be reached under particular initial conditions. There is also an x-dominance state in which conversely x persists but y goes extinct, despite including all hybrids. This equilibrium is sometimes stable depending on a restricted range of parameter values. Finally, there is a single coexistence equilibrium in which both subspecies persist, which only sometimes exists and is also sometimes stable. As a result, assuming the absence of limit cycles, any given modelled competing-hybridising system with known or estimated parameters will have one of four combinations of possible equilibrium outcomes: • • • •

y-dominance only y-dominance and x-dominance y-dominance and coexistence y-dominance, x-dominance and coexistence

We are able to show that each combination of possible outcomes occurs for plausible parameter values, and further that the domain of attraction for each possible outcome is sufficiently large to make it plausible depending on the given set of parameter values. The assumption of a constant carrying capacity as modelled here is seldom backed up by empirical evidence (del Monte-Luna et al., 2004), especially in the case of mosquitoes as mentioned, which are known to have seasonal population cycles (Minakawa et al., 2001) with consequences for disease transmission (Wyse et al., 2007; Reiner et al., 2015). We have therefore extended the model to incorporate a square-wave periodic oscillation in the carrying capacity of both subspecies, and have shown that such a modification can make long-term coexistence more likely. Many of the case studies mentioned in this paper support coexistence as the stable state — this is supported both by the lines of evidence above, and simply that in many cases of subspecies extinction, there may not be evidence in the fossil or written record that the extinct subspecies ever existed. The parametrisation of dynamic process models from experimental or observational data is a difficult task (Bender et al., 1984). Here, we have investigated parametrisations for our system based on the study by Quilodrán et al. (2015) for a species complex of waterfrogs, which used a discrete–time model with interspecific density dependence described by a Ricker functional form. The comparison between studies is supported by the fact that the parameters of a Ricker model derived as a discrete approximation to a continuous-time ODE system of Lotka-Volterra equations have an identical role as the continuous-time parameters in defining the long term time averages of the species abundances, which are equivalent for both continuous or discrete time formulations (Hofbauer et al., 1987). Hosack and Trenkel (2019) used an analogous argument to relate the parameters of an ODE system with Gompertz–type density dependence to a discrete time approximation, where the deterministic equilibrium was equivalent whether the model was formulated in continuous or discrete time. The discrete–time process model was then embedded in a Bayesian state space model that enabled estimation of both density dependent and density independent parameters while allowing for both process and observation uncertainty. Developing such statistical approaches for the estimation of parameters within dynamic models of species complexes will become increasingly possible as pub-

lic data sources with long term time series of species abundances are developed. For example, the development of Findable, Accessible, Interoperable and Reusable (FAIR; Wilkinson et al., 2016; 2018) databases for the abundance of arthropods (Rund et al., 2019a) including mosquito disease vectors (Rund et al., 2019b), such as the Culex pipiens species complex (Harbach, 2012) that transmit West Nile virus (Fonseca et al., 2004), provides an empirical base to support the statistical estimation of a species complex dynamic model from spatio-temporal abundance data. An advantage of developing robust mechanistic models is that their study can uncover management strategies or experiments that are more likely to succeed in relevant systems than naïve experimentation. For example, a two-subspecies system that is expected to have multiple stable equilibria could be manipulated by artificially changing the population sizes of one or both subspecies: this would be equivalent to changing the initial condition, ideally to the basin of attraction of a desirable outcome. Failing this, the parameters of the system could be artificially changed: for example, increasing or decreasing the amount of suitable habitat in order to change the carrying capacity (Kx or Ky ). Performing a sensitivity analysis on the model could inform which parameters would be most effective to change. Adopting an adaptive management approach of monitoring both before and after interventions would also help improve parameter estimates, making the model more robust for future management. Our model framework was designed to be as general as possible, using assumptions that are likely to be reasonable or parsimonious in a wide range of cases. To recap, these are: 1. each female mates with a constant number of males (i.e. females can reliably find males even at low density) 2. hybridisation is controlled by mate choice, with all hybrids treated identically 3. competition is controlled by a Lotka-Volterra model 4. the sex ratio of births and of the population are fixed 5. only fecundity, not mortality, is affected by population density Naturally, many specific case studies will be substantially different from what we have described. Perhaps most importantly, the fixing of sex ratios and binary treatment of hybrids removes our ability to model finer-scale dynamics such as changes in F1 and multiple generational hybrids or differential introgression across semipermeable species boundaries (Harrison and Larson, 2014). However, where any such new mechanisms can be mathematically described and potentially parameterised, our model can readily be modified to accommodate such additional detail. For example, asymmetric mating permissiveness has been documented for species groups such as the An. gambiae species complex (Marsden et al., 2011; Niang et al., 2015). In this case, w could be elaborated to wi . Compared to models that only account for competitive interactions among members of this species complex (Tene Fossog et al., 2015), the addition of asymmetric mating leads to rich and complex model behaviour. Far more complex models and datasets can also be considered, however, allowing the use of stochastic modelling, spatiotemporal modelling and statistical inference where supported by available data on population genetics, abundance and associated life history parameters. This will allow future models to represent behaviour such as the spatial emergence of hybridisation zones.

Declaration of Competing Interest None.

N.J. Beeton, G.R. Hosack and A. Wilkins et al. / Journal of Theoretical Biology 486 (2020) 110072

11

Acknowledgments The authors would like to thank Sean McElwain, Dan Pagenham and two anonymous reviewers for their thoughtful comments. This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors. Appendix A. Sex ratio convergence To demonstrate the rate of sex ratio convergence, we construct a simple model of the interactions between two sexes in a single subspecies. We first remove density dependence as this term only affects xF + xM and does not affect the ratio between sexes. Once this is removed, the sex ratios for each subspecies according to Eqs. (3)–(6) can be shown to satisfy the equations

 dqx Bx  = pF − qx ( 1 − pF ) dt xM

(A.1)

 dqy By  = pF − qy ( 1 − pF ) dt yM

(A.2)

x y where qx = F and qy = F . The equilibria of these equations can xM yM be easily seen to be

qx = qy =

pF 1 − pF

which corresponds to females comprising a proportion of pF , as expected. These equilibria are stable, with the rate of convergence By Bx proportionally related to and respectively. If these are kept xM yM constant, assuming that the effect of hybridisation on birth rates remains the same for the time period of interest, the solutions to Eqs. (A.1) and (A.2) are









pF Bx pF qx (t ) = + Ix − exp − (1 − pF )t 1 − pF 1 − pF xM    B  pF pF y qy (t ) = + Iy − exp − (1 − pF )t , 1 − pF 1 − pF yM

where qx (0 ) = Ix and qy (0 ) = Iy . The expected time scale over pF which the transition from Ix or Iy to occurs is 1 − pF

xM t≈ Bx ( 1 − pF )

or

yM . By ( 1 − pF )

1

γx Ix (1 − pF )

or

1

γy Iy (1 − pF )

where Kx and Ky are replaced with β (t)Kx and β (t)Ky (see example Eq. (23)) where 0 < β (t) ≤ 1 for all t, then the corresponding basins of attraction for this modified model will be subsets of those of the original model. Proof. Consider the constant-K model described in Eqs. (11) and (12) with a heteroclinic orbit described as in Fig. B.1. Let the heteroclinic orbit through P be defined by the function y = f (x ). Then for the constant-β K model described in the Figure with β ≤ 1, the corresponding heteroclinic orbit through P will satisfy the formula

1

β

y= f

1  β

x

.

Where the initial conditions Ix and Iy are already near the equi1 libria, this can be further simplified to t ≈ where γ reprepF γ sents the relevant birth rate. In biological terms, this is approximately the number of generations needed to replace all of the females in the population. In terms of the given examples, this may be a few months for anopheline mosquitoes to potentially a decade or more for edible frogs. Appendix B. Proof of larger coexistence basin of attraction with periodically decreasing carrying capacity Theorem 6. Assume the basins of attraction of the x-dominance (where it exists) and y-dominance equilibria are concave and increasing monotonically for the model described by Eqs. (11)and (12). Then

(B.1)

and can thus be defined by the function

y = g( x ) = β f

1  β

x

(B.2)

but then taking derivatives of both sides in Eq. (B.2) results in



1 g (x ) = f

x

β



(B.3)

and as f(x) is concave, its derivative is monotonically decreasing, so f (x ) ≥ f β1 x and

g (x ) ≤ f (x )

If we assume that w is very small, Eqs. (7) and (8) give us By γy yF Bx γx xF ≈ and ≈ (this is exact where w = 0 and no hyxM xM yM yM bridisation occurs). Then the time scale for transition is

t≈

Fig. B.1. An exaggerated example of the basin of attraction for the x-dominance equilibrium (Q), such as seen in Fig. 5d, is defined by the x-axis and the heteroclinic orbit between the origin and the nearest unstable coexistence equilibrium P. This construction is replicated on the phase plane for a similar model where the original carrying capacities Kx and Ky are both multiplied by some factor β < 1, generating new points P and Q .

(B.4)

Then as f (0 ) = g(0 ) = 0 and g (x) ≤ f (x) for all x > 0, g(x) ≤ f(x) over the same domain, and the basin of attraction for Q is a subset of the basin for Q. We can repeat the same argument with the axes reversed for the basin of attraction of the y-dominance equilibrium. By the definition of a basin of attraction, if y > f(x) ≥ g(x) for any point in time t in the constant-β K model, the trajectory of x and y will never reach the equilibrium point Q . But this also holds true at any point in time in the varying-β (t)K model as β (t) ≤ 1, so the (likely far more complicated) basin of attraction of the varyingβ (t)K model is also a subset of the constant-K model, completing our proof.  References Allendorf, F.W., Leary, R.F., Spruell, P., Wenburg, J.K., 2001. The problems with hybrids: setting conservation guidelines. Trends Ecol. Evolut. 16 (11), 613–622. doi:10.1016/S0169-5347(01)02290-X. Anderson, C.J., Oakeshott, J.G., Tay, W.T., Gordon, K.H.J., Zwick, A., Walsh, T.K., 2018. Hybridization and gene flow in the mega-pest lineage of moth, Helicoverpa. Proc. Natl. Acad. Sci. 115 (19), 5034–5039. https://www.pnas.org/content/115/19/ 5034.full.pdf. Begon, M., Mortimer, M., Thompson, D.J., 2009. Population Ecology: A Unified Study of Animals and Plants, 3rd ed. Wiley-Blackwell.

12

N.J. Beeton, G.R. Hosack and A. Wilkins et al. / Journal of Theoretical Biology 486 (2020) 110072

Bender, E.A., Case, T.J., Gilpin, M.E., 1984. Perturbation experiments in community ecology: theory and practice. Ecology 65 (1), 1–13. Boldin, B., 2006. Introducing a population into a steady community: the critical case, the center manifold, and the direction of bifurcation. SIAM J. Appl. Math. 66 (4), 1424–1453. Bostock, J., McAndrew, B., Richards, R., Jauncey, K., Telfer, T., Lorenzen, K., Little, D., Ross, L., Handisyde, N., Gatward, I., Corner, R., 2010. Aquaculture: global status and trends. Philosoph. Trans. R. Soc. B Biol. Sci. 365, 2897–2912. doi:10.1098/ rstb.2010.0170. Carvajal-Rodr-guez, A., 2018. Non-random mating and information theory. Theor. Popul. Biol. 120, 103–113. doi:10.1016/j.tpb.2018.01.003. Currat, M., Ruedi, M., Petit, R.J., Excoffier, L., 2008. The hidden side of invasion: massive introgression by local genes. Evolution 62 (8), 1908–1920. doi:10.1111/ j.1558-5646.20 08.0 0413.x. Cushing, J.M., 1980. Two species competition in a periodic environment. J. Math. Biol. 10 (4), 385–400. doi:10.10 07/bf0 0276097. Cushing, J.M., 1986. Periodic Lotka-Volterra competition equations. J. Math. Biol. 24 (4), 381–403. doi:10.1007/bf01236888. Edelstein-Keshet, L., 1988. Mathematical Models in Biology. Reprint. Soc. Ind. Appl. Math.. Fonseca, D.M., Keyghobadi, N., Malcolm, C.A., Mehmet, C., Schaffner, F., Mogi, M., Fleischer, R.C., Wilkerson, R.C., 2004. Emerging vectors in the Culex pipiens complex. Science 303 (5663), 1535–1538. Ford, J.S., Myers, R.A., 2008. A global assessment of salmon aquaculture impacts on wild salmonids. PLoS Biol. 6 (2), e33. doi:10.1371/journal.pbio.0 060 033. Geritz, S.A.H., Kisdi, E., 2004. On the mechanistic underpinning of discrete-time population models with complex dynamics. J. Theor. Biol. 228 (2), 261–269. doi:10.1016/j.jtbi.2004.01.003. Graf, J., 1986. Population genetics of the Rana esculenta complex: a model.. In: Roc˘ ek, Z. (Ed.), Studies in Herpetology: Proceedings of the European Herpetological Meeting (3rd Ordinary General Meeting of the Societas Europaea Herpetologica), Prague, 1985. Charles University, Prague. Harbach, R.E., 2012. Culex pipiens: species versus species complex – taxonomic history and perspective. J. Am. Mosq. Control Assoc. 28 (4s), 10–23. Harrison, R.G., Larson, E.L., 2014. Hybridization, introgression, and the nature of species boundaries. J. Hered. 105, 795–809. Hofbauer, J., Hutson, V., Jansen, W., 1987. Coexistence for systems governed by difference equations of Lotka-Volterra type. J. Math. Biol. 25, 553–570. Hosack, G.R., Trenkel, V.M., 2019. Functional group based marine ecosystem assessment for the Bay of Biscay via elasticity analysis. PeerJ 7, e7422. doi:10.7717/ peerj.7422. Lovette, I., 2016. Cornell lab of Ornithology’s Handbook of Bird Biology. John Wiley & Sons, Inc.. Marsden, C.D., Lee, Y., Nieman, C.C., Sanford, M.R., Dinis, J., Martins, C., Rodrigues, A., Cornel, A.J., Lanzaro, G.C., 2011. Asymmetric introgression between the M and S forms of the malaria vector, Anopheles gambiae, maintains divergence despite extensive hybridization. Mol. Ecol. 20 (23), 4983–4994. May, R.M., Mac Arthur, R.H., 1972. Niche overlap as a function of environmental variability. Proc. Natl. Acad. Sci. 69 (5), 1109–1113. https://www.pnas.org/ content/69/5/1109.full.pdf. McCallum, H., Barlow, N., Hone, J., 2001. How should pathogen transmission be modelled? Trends Ecol. Evolut. 16 (6), 295–300. doi:10.1016/S0169-5347(01) 02144-9. Menge, D.M., Guda, T., Zhong, D., Pai, A., Zhou, G., Beier, J.C., Gouagna, L., Yan, G., 2005. Fitness consequences of Anopheles gambiae population hybridization. Malar. J. 4 (1), 44. doi:10.1186/1475- 2875- 4- 44. Minakawa, N., Githure, J.I., Beier, J.C., Yan, G., 2001. Anopheline mosquito survival strategies during the dry period in western Kenya. J. Med. Entomol. 38 (3), 388–392. del Monte-Luna, P., Brook, B.W., Zetina-Rejón, M.J., Cruz-Escalona, V.H., 2004. The carrying capacity of ecosystems. Global Ecol. Biogeogr. 13 (6), 485–495. Muoz-Fuentes, V., Vil, C., Green, A.J., Negro, J.J., Sorenson, M.D., 2007. Hybridization between white-headed ducks and introduced ruddy ducks in Spain. Mol. Ecol. 16 (3), 629–638. https://onlinelibrary.wiley.com/doi/pdf/10.1111/j. 1365-294X.2006.03170.x. Niang, A., Epopa, P.S., Sawadogo, S.P., Maïga, H., Konaté, L., Faye, O., Dabiré, R.K., Tripet, F., Diabaté, A., 2015. Does extreme asymmetric dominance promote hybridization between Anopheles coluzzii and Anopheles gambiae s.s. in seasonal malaria mosquito communities of West Africa? Paras. Vect. 8 (1), 586.

Niang, E.H.A., Konateé, L., Faye, O., Diallo, M., Dia, I., 2019. Vector bionomics and malaria transmission in an area of sympatry of An. arabiensis, An. coluzzii and An. gambiae. Acta Trop. 189, 129–136. Plótner, J., Grunwald, C., 1991. A mathematical model of the structure and the dynamics of Rana ridibunda/esculenta-♂♂-populations (Anura, Ranidae). Z. Zool. Syst. Evolut. 29, 201–207. Quilodrán, C.S., Austerlitz, F., Currat, M., Montoya-Burgos, J.I., 2018. Cryptic biological invasions: a general model of hybridization. Sci. Rep. 8 (1), 2414. doi:10.1038/ s41598- 018- 20543- 6. Quilodrán, C.S., Currat, M., Montoya-Burgos, J.I., 2014. A general model of distant hybridization reveals the conditions for extinction in Atlantic salmon and brown trout. PLoS ONE 9 (7), 1–10. doi:10.1371/journal.pone.0101736. Quilodrán, C.S., Montoya-Burgos, J.I., Currat, M., 2015. Modelling interspecific hybridization with genome exclusion to identify conservation actions: the case of native and invasive Pelophylax waterfrogs. Evol. Appl. 8 (2), 199–210. doi:10. 1111/eva.12245. R Core Team, 2018. R: A language and environment for statistical computing. R Foundation for Statistical Computing. Vienna, Austria. Reiner, R.C., Geary, M., Atkinson, P.M., Smith, D.L., Gething, P.W., 2015. Seasonality of Plasmodium falciparum transmission: a systematic review. Malar. J. 14 (1), 343. Rhymer, J.M., Simberloff, D., 1996. Extinction by hybridization and introgression. Annu. Rev. Ecol. Syst. 27 (1), 83–109. Rinaldi, S., Muratori, S., Kuznetsov, Y., 1993. Multiple attractors, catastrophes and chaos in seasonally perturbed predator-prey communities. Bull. Math. Biol. 55 (1), 15–35. doi:10.1016/s0 092-8240(05)80 060-6. Rund, S.S.C., Braak, K., Cator, L., Copas, K., Emrich, S.J., Giraldo-Calderón, G.I., Johansson, M.A., Heydari, N., Hobern, D., Kelly, S.A., et al., 2019. MIReAd, a minimum information standard for reporting arthropod abundance data. Sci. Data 6 (1), 40. Rund, S.S.C., Moise, I.K., Beier, J.C., Martinez, M.E., 2019. Rescuing troves of hidden ecological data to tackle emerging mosquito-borne diseases. J. Am. Mosq. Control Assoc. 35 (1), 75–83. Schneider, P., Takken, W., McCall, P., 20 0 0, . Interspecific competition between sibling species larvae of Anopheles arabiensis and An. gambiae. Med. Vet. Entomol. 14. 165–70. doi: 10.1046/j.1365-2915.20 0 0.0 0204.x. Semlitsch, R.D., Ryer, H.-U., 1992. Performance of tadpoles from the hybridogenetic Rana esculenta complex: interactions with pond drying and interspecific competition. Evolution 46 (3), 665–676. Simberloff, D., Martin, J.-L., Genovesi, P., Maris, V., Wardle, D.A., Aronson, J., Courchamp, F., Galil, B., Garca-Berthou, E., Pascal, M., Pyek, P., Sousa, R., Tabacchi, E., Vil, M., 2013. Impacts of biological invasions: what’s what and the way forward. Trends Ecol. Evolut. 28 (1), 58–66. doi:10.1016/j.tree.2012.07.013. Soetaert, K., Petzoldt, T., Setzer, R.W., 2010. Solving differential equations in R: package deSolve. J. Stat. Softw. 33 (9), 1–25. doi:10.18637/jss.v033.i09. Stone, L., Olinky, R., Huppert, A., 2007. Seasonal dynamics of recurrent epidemics. Nature 446, 533. doi:10.1038/nature05638. Strogatz, S.H., 1994. Nonlinear Dynamics and Chaos. Perseus Books Publishing, LLC. Tene Fossog, B., Ayala, D., Acevedo, P., Kengne, P., Ngomo Abeso Mebuy, I., Makanga, B., Magnus, J., Awono-Ambene, P., Njiokou, F., Pombi, M., et al., 2015. Habitat segregation and ecological character displacement in cryptic African malaria mosquitoes. Evol. Appl. 8 (4), 326–345. Turchin, P., 2003. Complex Population Dynamics. Princeton University Press, Princeton, New Jersey. Walther, G.-R., 2010. Community and ecosystem responses to recent climate change. Philosoph. Trans. R. Soc. B Biolog. Sci. 365, 2019–2024. doi:10.1098/rstb.2010. 0021. Weidlich, E.W.A., von Gillhaussen, P., Delory, B.M., Blossfeld, S., Poorter, H., Temperton, V.M., 2016. The importance of being first: exploring priority and diversity effects in a grassland field experiment.. Front Plant Sci. 7, 2008. Wilkinson, M.D., Dumontier, M., Aalbersberg, I.J., Appleton, G., Axton, M., Baak, A., Blomberg, N., Boiten, J.-W., da Silva Santos, L.B., Bourne, P.E., et al., 2016. The FAIR guiding principles for scientific data management and stewardship. Sci. Data 3. Wilkinson, M.D., Sansone, S.-A., Schultes, E., Doorn, P., da Silva Santos, L.O.B., Dumontier, M., 2018. A design framework and exemplar metrics for FAIRness. Sci Data 5. Wyse, A.P.P., Bevilacqua, L., Rafikov, M., 2007. Simulating malaria model for different treatment intensities in a variable environment. Ecol Modell 206 (3–4), 322–330.