Skip to main content
Log in

Pulled, pushed or failed: the demographic impact of a gene drive can change the nature of its spatial spread

  • Published:
Journal of Mathematical Biology Aims and scope Submit manuscript

Abstract

Understanding the temporal spread of gene drive alleles—alleles that bias their own transmission—through modeling is essential before any field experiments. In this paper, we present a deterministic reaction-diffusion model describing the interplay between demographic and allelic dynamics, in a one-dimensional spatial context. We focused on the traveling wave solutions, and more specifically, on the speed of gene drive invasion (if successful). We considered various timings of gene conversion (in the zygote or in the germline) and different probabilities of gene conversion (instead of assuming 100\(\%\) conversion as done in a previous work). We compared the types of propagation when the intrinsic growth rate of the population takes extreme values, either very large or very low. When it is infinitely large, the wave can be either successful or not, and, if successful, it can be either pulled or pushed, in agreement with previous studies (extended here to the case of partial conversion). In contrast, it cannot be pushed when the intrinsic growth rate is vanishing. In this case, analytical results are obtained through an insightful connection with an epidemiological SI model. We conducted extensive numerical simulations to bridge the gap between the two regimes of large and low growth rate. We conjecture that, if it is pulled in the two extreme regimes, then the wave is always pulled, and the wave speed is independent of the growth rate. This occurs for instance when the fitness cost is small enough, or when there is stable coexistence of the drive and the wild-type in the population after successful drive invasion. Our model helps delineate the conditions under which demographic dynamics can affect the spread of a gene drive.

This is a preview of subscription content, log in via an institution to check access.

Access this article

Price excludes VAT (USA)
Tax calculation will be finalised during checkout.

Instant access to the full article PDF.

Fig. 1
Fig. 2
Fig. 3
Fig. 4
Fig. 5
Fig. 6

Similar content being viewed by others

Notes

  1. \(\partial _{xx}:^2 n_{\textrm{DD}} = \partial _{xx}^2 n p_{\textrm{D}} = \partial _{x} ( p_{\textrm{D}} \ \partial _{x} n + n \ \partial _{x} p_{\textrm{D}} ) = p_{\textrm{D}} \ \partial _{xx}^2 n + 2 \ \partial _{x} p_{\textrm{D}} \ \partial _{x} n + n \ \partial _{xx}^2 p_{\textrm{D}} \)

  2. \(\partial _{xx}:^2 n_{\textrm{D}} = \partial _{xx}^2 n p_{\textrm{D}} = \partial _{x} ( p_{\textrm{D}} \ \partial _{x} n + n \ \partial _{x} p_{\textrm{D}} ) = p_{\textrm{D}} \ \partial _{xx}^2 n + 2 \ \partial _{x} p_{\textrm{D}} \ \partial _{x} n + n \ \partial _{xx}^2 p_{\textrm{D}} \).

  3. \(\partial _{xx}:^2 n_{\textrm{D}} = \partial _{xx}^2 n p_{\textrm{D}} = \partial _{x} ( p_{\textrm{D}} \ \partial _{x} n + n \ \partial _{x} p_{\textrm{D}} ) = p_{\textrm{D}} \ \partial _{xx}^2 n + 2 \ \partial _{x} p_{\textrm{D}} \ \partial _{x} n + n \ \partial _{xx}^2 p_{\textrm{D}} \).

  4. \(M_z'(p_{\textrm{D}}) = \mathscr {A}_z ( 1 - 2 p_{\textrm{D}}) - s \le 0 \) therefore \(M_z(1) \le M_z(p_{\textrm{D}}) \le M_z(0)\).

  5. \(M_g'(p_{\textrm{D}}) = \mathscr {A}_g ( 1 - 2 p_{\textrm{D}}) - s \le 0 \) therefore \(M_g(1) \le M_g(p_{\textrm{D}}) \le M_g(0)\).

References

Download references

Acknowledgements

This work is funded by ANR-19-CE45-0009-01 “TheoGeneDrive”. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programm (grant agreement No 865711). We are grateful to the INRAE MIGALE bioinformatics facility (MIGALE, INRAE, 2020. Migale bioinformatics Facility, DOI: 10.15454/1.5572390655343293E12) for providing computing resources.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Léna Kläy.

Ethics declarations

Conflict of interest

The authors declare that they have no conflict of interest.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Appendices

Appendix

A Model with partial conversion: growth term details

To obtain the global growth term for each genotypes in models (4) and (5), we calculate type proportions among the offspring for each possible couple, and then we sum the corresponding terms. The calculations follow standard lines of population genetics, differing only by the timing of gene conversion. When conversion occurs in the zygote, the parameter c appears in between the gametes and the offspring production, whereas when conversion occurs in the germline, it appears before gametes production. These equations in densities are consistent with the one obtained in frequency in the literature (when conversion occurs in the zygote (Deredec et al. 2008; Unckless et al. 2015), or in the germline (Deredec et al. 2008; Rode et al. 2019)).

1.1 A.1 Conversion occurring in the zygote

Table 8 Growth term details when conversion occurs in the zygote

1.2 A.2 Conversion occurring in the germline

Table 9 Growth term details when conversion occurs in the germline

B System rewritten with variables \((n, p_{\textrm{D}})\)

Below, we present the details of the reformulation from models (6), (17) and (27) in terms of total population density n and drive allele frequency \(p_{\textrm{D}}\).

1.1 B.1 Model with perfect conversion

We rewrite model (6) with variables:

$$\begin{aligned} n = n_{\textrm{WW}} + n_{\textrm{DD}}, \quad \quad \quad \quad p_{\textrm{D}} = \dfrac{n_{\textrm{D}}}{n_{\textrm{W}} + n_{\textrm{D}}} = \dfrac{n_{\textrm{DD}}}{n_{\textrm{WW}} + n_{\textrm{DD}}} . \end{aligned}$$
(35)

where n is the total population density and \(p_{\textrm{D}}\) is the drive allele frequency, or equivalently in this model, the frequency of drive homogeneous individuals.Footnote 1

Equation on n:

$$\begin{aligned} \begin{aligned} \partial _t n - \partial _{xx}^2 n&= \Big (r \ (1-n) + 1 \Big ) \ \Big ( (1-s) \ \dfrac{n_{\textrm{DD}}^2 + 2 \ n_{\textrm{WW}} n_{\textrm{DD}}}{n_{\textrm{WW}}+n_{\textrm{DD}}} + \ \dfrac{n_{\textrm{WW}}^2}{n_{\textrm{WW}}+n_{\textrm{DD}}} \Big ) - n, \\&= \Big (r \ (1-n) + 1 \Big ) \ \Big ( (1-s) \ p_{\textrm{D}}^2 + (1-s) \ 2 \ p_{\textrm{D}} (1-p_{\textrm{D}}) \\&\quad + (1-s ) (1-p_{\textrm{D}})^2+ s (1-p_{\textrm{D}})^2 \Big ) n - n, \\&= \Big (r \ (1-n) + 1 \Big ) \ \Big ( (1-s) (p_{\textrm{D}} + 1 - p_{\textrm{D}})^2 + s (1-p_{\textrm{D}})^2 \Big ) n - n, \\&= \Big (r \ (1-n) + 1 \Big ) \ \Big ( (1-s) + s (1-p_{\textrm{D}})^2 \Big ) n - n, \\&= \Big (r \ (1-n) + 1 \Big ) \ \Big ( (1-s) \ p_{\textrm{D}} \ (2-p_{\textrm{D}}) + (1-p_{\textrm{D}})^2 \Big ) n - n. \end{aligned} \end{aligned}$$
(36)

Equation on \(\mathbf {n_{\textrm{DD}}}\)

(37)

Equation on \(\mathbf {p_{\textrm{D}}= \dfrac{ n_{DD}}{n}}\):

(38)

Combining equations on n and \(p_{\textrm{D}}\), we obtain model (12).

1.2 B.2 Model with partial conversion

1.2.1 B.2.1 Conversion in the zygote

We rewrite model (17) with variables:

$$\begin{aligned} n = n_{\textrm{W}} + n_{\textrm{D}}, \quad \quad \quad \quad p_{\textrm{D}} = \dfrac{n_{\textrm{D}}}{n_{\textrm{W}} + n_{\textrm{D}}}. \end{aligned}$$
(39)

where n is the total population density and \(p_{\textrm{D}}\) is the drive allele frequency.

Equation on \(\textbf{n}\):

$$\begin{aligned} \begin{aligned} \partial _t n - \partial _{xx}^2 n&= \dfrac{ r \ (1-n)+1 }{n} \ \Big ( (1-s) \ n_{\textrm{D}}^2 + [ 2 \ c \ (1-s) + 2 \ (1-sh) \ (1-c) ] \\&\qquad n_{\textrm{D}} n_{\textrm{W}} + n_{\textrm{W}}^2 \Big ) - n, \\&= \dfrac{ r \ (1-n)+1 }{n} \ \Big ( (1-s) \ (np_{\textrm{D}})^2 + [ 2 \ c \ (1-s) + 2 \ (1-sh) \ (1-c) ] \\&\qquad p_{\textrm{D}} (1-p_{\textrm{D}}) n^2 + (1-p_{\textrm{D}})^2 n^2 \Big ) - n, \\&= \big ( r \ (1-n)+1 \big ) \ \Big ( (1-s) \ p_{\textrm{D}}^2 + [ 2 \ c \ (1-s) + 2 \ (1-sh) \ (1-c) ] \\&\qquad p_{\textrm{D}} \ (1-p_{\textrm{D}}) + (1-p_{\textrm{D}})^2 \Big ) n - n,\\&= \big ( r \ (1-n)+1 \big ) \ \Big ( (1-s) \ p_{\textrm{D}}^2 + 2 \ p_{\textrm{D}} \ (1-p_{\textrm{D}}) \ [ c \ (1-s) + (1-c) \ (1-sh) ] \\&\qquad + (1-p_{\textrm{D}})^2 \Big ) n - n. \end{aligned} \end{aligned}$$
(40)

Equation on \(\mathbf {n_{\textrm{D}}}\)Footnote 2

(41)

Equation on \(\mathbf {p_{\textrm{D}}= \dfrac{n_{\textrm{D}}}{n}}\):

(42)

Combining equations on n and \(p_{\textrm{D}}\), we obtain model (21).

1.2.2 B.2.2 Conversion in the germline

We rewrite model (27) with variables:

$$\begin{aligned} n = n_{\textrm{W}} + n_{\textrm{D}}, \quad \quad \quad \quad p_{\textrm{D}} = \dfrac{n_{\textrm{D}}}{n_{\textrm{W}} + n_{\textrm{D}}}. \end{aligned}$$
(43)

where n is the total population density and \(p_{\textrm{D}}\) is the drive allele frequency.

Equation on n:

$$\begin{aligned} \begin{aligned} \partial _t n - \partial _{xx}^2 n&= \dfrac{ r \ (1-n)+1 }{n} \ \Big ( (1-s) \ n_{\textrm{D}}^2 + 2 \ (1-sh) \ n_{\textrm{D}} n_{\textrm{W}} + n_{\textrm{W}}^2 \Big ) - n, \\&= \dfrac{ r \ (1-n)+1 }{n} \ \Big ( (1-s) \ (np_{\textrm{D}})^2 + 2 \ (1-sh) \ p_{\textrm{D}} (1-p_{\textrm{D}}) n^2 \\&\quad + ((1-p_{\textrm{D}})n)^2 \Big ) - n, \\&= \big ( r \ (1-n)+1 \big ) \ \Big ( (1-s) \ p_{\textrm{D}}^2 + 2 \ (1-sh) \ p_{\textrm{D}} \ (1-p_{\textrm{D}}) + (1-p_{\textrm{D}})^2 \Big ) n - n. \\ \end{aligned} \end{aligned}$$
(44)

Equation on \(\mathbf {n_{\textrm{D}}}\)Footnote 3

(45)

Equation on \(\mathbf {p_{\textrm{D}}= \dfrac{n_{\textrm{D}}}{n}}\):

$$\begin{aligned} \begin{aligned} \partial _t p_{\textrm{D}}&= \dfrac{ n \ (\partial _t n_{\textrm{D}}) - n_{\textrm{D}} \ (\partial _t n)}{n^2} = \dfrac{ n \ (\partial _t n_{\textrm{D}}) - n \ p_{\textrm{D}} \ (\partial _t n)}{n^2} = \dfrac{1}{n} \big ( \ \partial _t n_{\textrm{D}} - p_{\textrm{D}} \ (\partial _t n) \big ), \\&= \dfrac{1}{n} \Big [ 2 \ \partial _x n \ \partial _x p_{\textrm{D}} + n \ \partial _{xx}^2 p_{\textrm{D}} \ + \big ( r \ (1-n)+1 \big ) \Big ( \ (1-s) \ p_{\textrm{D}}\\&\quad \quad + (1-sh) \ (1+c) \ (1-p_{\textrm{D}}) \\&\quad \quad - (1-s) \ p_{\textrm{D}}^2 - 2 \ (1-sh) \ p_{\textrm{D}} \ (1-p_{\textrm{D}}) - (1-p_{\textrm{D}})^2 \Big ) \ p_{\textrm{D}} \ n \ \Big ], \\&= 2 \ \partial _x \log (n) \ \partial _x p_{\textrm{D}} + \partial _{xx}^2 p_{\textrm{D}} \ + \big ( r \ (1-n)+1 \big ) \Big ( (1-s) \ p_{\textrm{D}} + (1-sh)(1+c) \\&\qquad - 2 \ (1-sh) \ p_{\textrm{D}} - (1-p_{\textrm{D}}) \Big ) \ p_{\textrm{D}} \ (1-p_{\textrm{D}}), \\&= 2 \ \partial _x \log (n) \ \partial _x p_{\textrm{D}} + \partial _{xx}^2 p_{\textrm{D}} \ + \big ( r \ (1-n)+1 \big ) \Big ( (2h-1) \ s \ p_{\textrm{D}} \\&\qquad + (1-sh) (1+c) - 1 \Big ) \ p_{\textrm{D}} \ (1-p_{\textrm{D}}) . \end{aligned} \end{aligned}$$
(46)

Combining equations on n and \(p_{\textrm{D}}\), we obtain model (29).

C Proofs for model (12) with perfect conversion in the zygote

1.1 C.1 Numerical evidence for the continuity when \(r \rightarrow 0\)

In Fig. 7, we plot the speed of the traveling wave solutions of model (12) for a range of r and s values, and for \(r=0\). A positive speed correspond to drive invasion.

We observe continuity in the speed value when \(r \rightarrow 0\) away from \(s = \frac{1}{2}\), meaning that the case \(r=0\) is relevant to approximate very small intrinsic growth rates.

Fig. 7
figure 7

Wave speed values in model with perfect conversion in the zygote (12), regarding parameters r the intrinsic growth rate (log scale in between 0.01 and 10, plus the exact value \(r=0\) in the bottom color line) and s the fitness disadvantage for drive (normal scale). Below the pure drive persistence line (light green), a well-mixed population containing only drive homozygous individuals will necessarily go extinct (color figure online)

1.2 C.2 Proof of the statements in Tables 3 and 4 when perfect conversion occurs in the zygote

In this section we prove the statements of Tables 3 and 4 on the two models of interest:

\(\underline{{\textbf {r}} = \infty }\)

$$\begin{aligned} \partial _t p - \partial _{xx}^2 p = \dfrac{s \ p \ (1-p) \ \big (p- \dfrac{2s -1}{s} \big )}{1-s+s(1-p)^2} = f^{\infty }(p). \end{aligned}$$
(47)

\(\underline{{\textbf {r}} = 0}\)

$$\begin{aligned} \left\{ \begin{array}{ll} \partial _t n_{\textrm{DD}} - \partial _{xx}^2 n_{\textrm{DD}} &{} = \ (1-s) \ \dfrac{n_{\textrm{WW}} n_{\textrm{DD}}}{n_{\textrm{WW}}+n_{\textrm{DD}}} - s \ n_{\textrm{DD}} = f^0(n_{\textrm{DD}}, n_{\textrm{WW}}) \\ \\ \partial _t n_{\textrm{WW}} - \partial _{xx}^2 n_{\textrm{WW}} &{} = \ \dfrac{- \ n_{\textrm{WW}} n_{\textrm{DD}}}{n_{\textrm{WW}}+n_{\textrm{DD}}}. \end{array} \right. \end{aligned}$$
(48)

Monostable / Bistable

\(\underline{{\textbf {r}} = \infty }\)

\(0<s<0.5\) :

The equation admits two admissible steady states 0 and 1.

As \((f^{\infty })'(0)>0\) and \((f^{\infty })'(1)<0\), the only stable state is \(p=1\).

\(0.5<s<1\) :

The equation admits three admissible steady states 0, \(\dfrac{2s-1}{s}\) and 1.

As \((f^{\infty })'(0)<0\), \((f^{\infty })'(\frac{2s-1}{s})>0\) and \((f^{\infty })'(1)<0\), both \(p=0\) and \(p=1\) are stable states.

\(\underline{{\textbf {r}} = 0}\)

\(0<s<0.5\) :

The system admits \((n_{\textrm{DD}}=0, n_{\textrm{WW}} \in [0,1])\) as admissible steady states. The Jacobian matrix, when switching to n and \(p_{\textrm{D}}\) variables, indicates that the only stable state is \((n=0, p_{\textrm{D}}=1)\), i.e. \((n_{\textrm{DD}}=0, n_{\textrm{WW}}=0)\).

\(0.5<s<1\) :

The system admits \((n_{\textrm{DD}}=0, n_{\textrm{WW}} \in [0,1])\) as admissible steady states. The Jacobian matrix, when switching to n and \(p_{\textrm{D}}\) variables, indicates that the stable states are \((n=0, p_{\textrm{D}}=1)\) and \((n \in [0,1], p_{\textrm{D}}=0)\), i.e. \((n_{\textrm{DD}}=0, n_{\textrm{WW}} \in [0,1])\).

Existence of critical traveling waves

\(\underline{{\textbf {r}} = \infty }\)

  • The existence of traveling waves for the scalar equation (13) in both monostable and bistable cases is a classical result in the theory of reaction-diffusion equations, see for instance the seminal works in Aronson and Weinberger (1975, 1978).

\(\underline{{\textbf {r}} = 0}\)

\(0<s<0.5\) :

We apply the results of Appendix D with \(\beta _1 = 1\) and \(\beta _2 = 1 - s\). Therefore system (48) admits a traveling wave when \(0<s<0.5\).

\(0.5<s<1\) :

There is no drive propagation due to the gene drive clearance: the drive allele density decreases uniformly in space (details in Sect. C.2.1). Regarding the wild-type alleles, their dynamic is given by the heat equation implying only diffusion and no growth. It cannot admit traveling wave solutions.

Pulled/pushed waves and speed values

For both models, the speed of the linearized problem around zero density of drive allele is given by \( 2 \sqrt{1-2\,s} = 2 \sqrt{(f^{\infty })'(0)} = 2 \sqrt{ \partial _p f^0(0,1)}\).

\(\underline{{\textbf {r}} = \infty }\)

\(0<s \lesssim 0.35\) :

Numerically, we observe that the speed of the wave is equal to the minimal speed of the linearized problem: the wave is pulled (detail in Sect. C.2.2)

\(0.35 \lesssim s<0.5\) :

Numerically, we observe that the speed of the wave is strictly above the minimal speed of the linearized problem: the wave is pushed (detail in Sect. C.2.2).

\(0.5<s<1\) :

As the system is bistable, the wave is necessarily pushed. The numerical approximation \(s \approx 0.70\) indicating whether the drive of the wild-type population will invade the environment was already determined in the work of Tanaka et al. (2017).

\(\underline{{\textbf {r}} = 0}\)

\(0<s<0.5\) :

We apply the results of Appendix D with \(\beta _1 = 1\) and \(\beta _2 = 1 - s\). Therefore system (48) admits a traveling wave with speed \(v = 2 \sqrt{1-2s}\) when \(0<s<0.5\). This value corresponds to the KPP speed, by definition the wave is pulled.

\(0.5< s < 1 \) :

No wave (see above, in Existence of critical traveling waves).

1.2.1 C.2.1 Gene drive clearance for \(s \in (0.5,1)\) when \(r = 0\)

Consider the model with perfect conversion in the zygote (6). The densities \(n_{\textrm{WW}}\) and \(n_{\textrm{DD}}\) dynamics are qualitatively given in Fig. 8 for \(r=0\).

Fig. 8
figure 8

Qualitative dynamics of the drive homozygotes density \(n_{\textrm{DD}}\) (red line) and the wild-type homozygotes density \(n_{\textrm{WW}}\) (blue line) in space (color figure online)

When \(s>0.5\) and \(r=0\) we observe gene drive clearance (in Fig. 8B). More precisely, we have the following estimate, deduced from (15):

$$\begin{aligned} \partial _t n_{\textrm{DD}} - \partial _{xx}^2 n_{\textrm{DD}} \le \ (1-2s) \ n_{\textrm{DD}}, \end{aligned}$$
(49)

Therefore, \(n_{\textrm{DD}}\) is exponentially decaying in time, uniformly in space. The dynamics of the wild type then boil down to the standard heat equation, there cannot exist a traveling wave.

1.2.2 C.2.2 Numerical approximation of s threshold value for the pulled/pushed wave when \(r=+ \infty \)

In order to determine an approximation of the threshold value at which the wave switches from a pulled wave to a pushed wave, we used the recent continuation procedure published in Avery et al. (2022). Figure 9 presents the value of the wave speed obtained via the latter continuation numerical scheme (Holzer 2022), for a wide range of s values. Notice the transition between pulled fronts (plain red) and pushed fronts (plain green). For the sake of clarity, the value of the minimal speed of the linearized problem \(v = 2 \sqrt{1-2s}\) is shown in red for \(s\in (0,\frac{1}{2})\). Notice that the speed of the pushed front changes sign approximately at \(s\approx 0.70\), in agreement with the theoretical criterion (9).

Fig. 9
figure 9

Value of the wave speed when \(r=+ \infty \) obtained via the numerical scheme in Holzer (2022), for a wide range of s values. The transition between pulled fronts (plain red) and pushed fronts (plain green) is approximately 0.35. For the sake of clarity, the value of the minimal speed of the linearized problem \(v = 2 \sqrt{1-2\,s}\) is shown in red for \(s\in (0,\frac{1}{2})\) (color figure online)

D Critical travaling wave for an SI similar model.

Consider the following epidemiological model,

$$\begin{aligned} \left\{ \begin{array}{ll} \partial _t S - \partial _{xx}^2 S &{} = \ - \beta \ \dfrac{\ S \ I}{S + I},\\ \\ \partial _t I - \partial _{xx}^2 I &{} = \ \beta \ \dfrac{ S \ I}{S + I} - \gamma I.\\ \end{array} \right. \end{aligned}$$
(50)

where S is the density of susceptible individuals, I is the density of infected individuals, \(\gamma \) is the mortality of infected individuals and \(\beta \) is the transmission coefficient. This model has already been studied in the literature, see (Zhou et al. 2019) and references therein. In particular, the existence of a minimal traveling wave has been established in the latter reference.

Models (15), (25) and (33) are very similar to the above SI model (50), except that the coefficient \(\beta \) is different in the first and the second equation of the system. We write this new system with two coefficients \(\beta _1, \beta _2\):

$$\begin{aligned} \left\{ \begin{array}{ll} \partial _t S - \partial _{xx}^2 S &{} = \ - \beta _1 \ \dfrac{ \ S \ I}{S + I}, \\ \\ \partial _t I - \partial _{xx}^2 I &{} = \ \beta _2 \ \dfrac{ S \ I}{S + I} - \gamma I. \end{array} \right. \end{aligned}$$
(51)

1.1 D.1 Existence of critical traveling wave solutions

We are able to establish the following Theorem by adapting the proof in Zhou et al. (2019).

Theorem

Suppose that \(\beta _1>0\), and \(\beta _2 > \gamma \), then system (51) admits a positive and bounded traveling wave solution with profile \((S^*, I^*)\), and speed \(v = 2 \sqrt{\beta _2 - \gamma }\). Furthermore, both \(S^*\) and \(I^*\) are positive, and bounded by 1 and \(\frac{\beta _2 - \gamma }{\gamma }\) respectively.

By adapting further the elements of (Zhou et al. 2019), it would be possible to prove that the profile \(S^*\) is increasing, whereas the profil \(I^*\) is unimodal, as shown in Fig. 10.

Fig. 10
figure 10

Qualitative shape of the solution (\(I^*\) in red, \(S^*\) in blue)

1.2 D.2 Proof of the theorem

We proceed as follows:

  1. 1.

    Although the system does not satisfy the comparison principle due to a lack of monotonicity, the construction of traveling waves is performed through a construction of sub-solutions (\({\underline{S}}\), \({\underline{I}}\)) and super-solutions (\({\bar{S}}\),\({\bar{I}}\)) for the system.

  2. 2.

    Using Schauder’s fixed point theorem, we prove the existence of a critical traveling wave solution \((S^*, I^*)\) with speed v such that \({\underline{S}}(z) \le S^*(z) \le {\bar{S}}(z)\) and \({\underline{I}}(z) \le I(z) \le {\bar{I}}(z)\) for all z in \(\mathbb {R}\).

  3. 3.

    Finally, we conclude with the positivity of the critical traveling wave solution thanks to the strong maximum principle.

1.2.1 D.2.1 Construction of sub- and super-solutions

We are seeking sub- and super-solutions, respectively (\({\underline{S}}\), \({\underline{I}}\)), (\({\bar{S}}\),\({\bar{I}}\)). Because of the non-monotonic coupling in the system, the following set of cross-relationships must be satisfied:

  1. 1.

    \( -v \ {{\bar{S}}}' - {\bar{S}}'' \ge - \beta _1 \ \dfrac{{\bar{S}} \ I}{{\bar{S}}+I}\quad \forall {\underline{I}} \le I \le {\bar{I}} \);

  2. 2.

    \( -v \ {\underline{S}}' - {\underline{S}}'' \le - \beta _1 \dfrac{{\underline{S}} \ I}{{\underline{S}}+I} \quad \forall {\underline{I}} \le I \le {\bar{I}} \);

  3. 3.

    \( -v \ {\bar{I}}' - {\bar{I}}'' \ge \beta _2 \ \dfrac{S \ {\bar{I}}}{S + {\bar{I}}} - \gamma \ {\bar{I}} \quad \forall {\underline{S}} \le S \le {\bar{S}} \);

  4. 4.

    \( -v \ {\underline{I}}' - {\underline{I}}'' \le \beta _2 \ \dfrac{S \ {\underline{I}}}{S + {\underline{I}}} - \gamma \ {\underline{I}}\quad \forall {\underline{S}} \le S \le {\bar{S}} \).

Inequalities 1,2,3 and 4 are valid in a weak sense. As \({\underline{S}}\) is a piece-wise differentiable function, the quantity \({\underline{S}}\) consists of functions on each sub-interval with a Dirac mass at the point of \(\mathcal {C}^1\) discontinuity. However, the Dirac mass has the good sign in this case (the transition has a convex shape), hence the second derivative \({\underline{S}}''\) is non-negative in the sense of a measure. All signs are correct: \({\underline{I}}\) has a convex transition also, and \({\bar{I}}\) has a concave transition. This key point is equivalent to the standard principle in the theory of parabolic equations (widely used for reaction-diffusion equations involving comparison techniques): "the maximum of sub-solutions is a sub-solution" (here, \({\underline{S}}\), \({\underline{I}}\)) and "the minimum of super-solutions is a super-solution" (here, \({\bar{S}}\), \({\bar{I}}\)).

To define our sub and super-solutions, it is useful to introduce the the following family of functions \(\mathcal I(z) = e^{- \lambda ^* z}\), where \(\lambda ^*\) is solution of the following dispersion equation:

$$\begin{aligned} (\lambda ^*)^2 - v \lambda ^* + ( \beta _2 - \gamma ) = 0. \end{aligned}$$
(52)

They are solutions of the linearized problem

$$\begin{aligned} v \mathcal I' + \mathcal I'' + ( \beta _2 - \gamma )\ \mathcal I = 0. \end{aligned}$$
(53)

For the critical speed \(v = 2 \sqrt{(\beta _2 - \gamma )} \), the corresponding double root is \( \lambda ^* = \dfrac{v}{2} = \sqrt{(\beta _2 - \gamma )} \).

Lemma

There exist two large enough constants \(L_1 > 0\) and \(L_2 > 0\), such that the functions \( {\bar{S}}\), \({\underline{S}}\) \({\bar{I}}\), \( {\underline{I}} \) defined below satisfy the conditions 1. 2. 3. and 4.:

$$\begin{aligned}{} & {} {\bar{S}} = 1. \end{aligned}$$
(54)
$$\begin{aligned}{} & {} {\underline{S}} = \left\{ \begin{array}{ll} 0 \quad \quad &{} \forall z \le L_1 \log (L_1), \\ \\ 1 - L_1 e^{- \dfrac{z}{L_1} } \quad \quad &{} \forall z > L_1 \log (L_1). \\ \end{array} \right. \end{aligned}$$
(55)
$$\begin{aligned}{} & {} {\bar{I}} = \left\{ \begin{array}{ll} M \quad \quad &{} \forall z \le \dfrac{1}{\lambda ^*}, \\ \\ e M \lambda ^* z e^{- \lambda ^* z} \quad \quad &{} \forall z > \dfrac{1}{\lambda ^*}. \\ \end{array} \right. \end{aligned}$$
(56)
$$\begin{aligned}{} & {} {\underline{I}} = \left\{ \begin{array}{ll} 0 \quad \quad &{} \forall z \le \big ( \dfrac{L_2}{e M \lambda ^*} \big ) ^2,\\ \\ \big ( e M \lambda ^* z - L_2 \sqrt{z} \big ) e^{-\lambda ^* z} \quad \quad &{} \forall z > \big ( \dfrac{L_2}{e M \lambda ^*} \big ) ^2. \end{array} \right. \end{aligned}$$
(57)

with \(M = \dfrac{\beta _2 - \gamma }{\gamma } = \dfrac{(\lambda ^*)^2}{\gamma } \).

Proof

Before we proceed with the proof, we introduce the following set of conditions, which are more restrictive than 1,2,3,4, but may appear more useful at some point in the calculations. When satisfied, they clearly imply 1,2,3,4.

  1. (i)

    \( -v \ {\bar{S}}' - {\bar{S}}'' \ge 0 \);

  2. (ii)

    \( - v \ {\underline{S}}' - {\underline{S}}'' \le - \beta _1 \ {\bar{I}} \);

  3. (iii)

    \( -v \ {\bar{I}}' - {\bar{I}}'' \ge (\beta _2 - \gamma ) \ {\bar{I}} \);

  4. (iv)

    \( -v \ {\underline{I}}' - {\underline{I}}'' \le \beta _2 \ \dfrac{{\underline{S}} \ {\underline{I}}}{{\underline{S}} + {\underline{I}}} - \gamma \ {\underline{I}} \).

We verify each of the four conditions on \( {\bar{S}}\), \({\underline{S}}\) \({\bar{I}}\), \( {\underline{I}} \):

Condition 1: \( -v \ {{\bar{S}}}' - {\bar{S}}'' \ge - \beta _1 \ \dfrac{{\bar{S}} \ I}{{\bar{S}}+I} \).

The constant function \( {\bar{S}} = 1 \) satisfies the more restrictive condition (i) \( -v \ {\bar{S}}' - {\bar{S}}'' \ge 0 \).

Condition 2: \({\underline{S}}' - {\underline{S}}'' \le - \beta _1 \dfrac{{\underline{S}} \ I}{{\underline{S}}+I} \).

Let us take \(L_1\) sufficiently large such that \( \dfrac{1}{\lambda ^*} < L_1 \log (L_1)\).

  • For \( z > L_1 \log (L_1) \):

Since \( {\bar{I}} = e \ M \ \lambda ^* \ z \ e^{\lambda ^*z}\) and \( \ \ {\underline{S}} = 1 - L_1 e^{- \dfrac{z}{L_1} }\), the condition (ii) \( - v \ {\underline{S}}' - {\underline{S}}'' \le - \beta _1 \ {\bar{I}} \) holds for a sufficiently large \(L_1 > 0 \):

$$\begin{aligned} - v \ {\underline{S}}' - {\underline{S}}'' = (\dfrac{1}{L_1} - v ) \ e^{- \dfrac{z}{L_1}}&\le - \beta _1 \ e \ M \ \lambda ^* \ z \ e^{- \lambda ^*z} = - \beta _1 \ {\bar{I}}, \end{aligned}$$
(58)
$$\begin{aligned} \iff \beta _1 \ e \ M \ \lambda ^* \ z \ e^{ ( \dfrac{1}{L_1} - \lambda ^* ) z}&\le ( v - \dfrac{1}{L_1} ) . \end{aligned}$$
(59)
  • For \( z \le L_1 \log (L_1) \):

The condition 2. is verified since \({\underline{S}} = 0\).

Condition 3: \( -v \ {\bar{I}}' - {\bar{I}}'' \ge \beta _2 \ \dfrac{S \ {\bar{I}}}{S + {\bar{I}}} - \gamma \ {\bar{I}} \).

  • For \( z \le \dfrac{1}{\lambda ^*} \):

With \({\bar{I}} = M = \dfrac{\beta _2 - \gamma }{\gamma } \):

$$\begin{aligned} - v \ {\bar{I}}' - {\bar{I}}'' {=} 0 {=} \beta _2 \ \dfrac{M}{1+M} - \gamma \ M = \beta _2 \ \dfrac{{\bar{S}} \ {\bar{I}}}{{\bar{S}} {+} {\bar{I}}} - \gamma \ {\bar{I}} \ge \beta _2 \ \dfrac{S \ {\bar{I}}}{S + {\bar{I}}} - \gamma \ {\bar{I}} \quad \quad \quad \quad \forall S \le {\bar{S}}. \end{aligned}$$
(60)
  • For \( z > \dfrac{1}{\lambda ^*} \):

Since \({\bar{I}}\) is proportional to \(z e^{- \lambda ^* z}\), and \(\lambda ^*\) is precisely the double root of the characteristic equation (52), we deduce that condition (iii) \( -v \ {\bar{I}}' - {\bar{I}}'' - (\beta _2 - \gamma ) \ {\bar{I}} \ge 0 \) is verified.

Condition 4: \( -v \ {\underline{I}}' - {\underline{I}}'' \le \beta _2 \ \dfrac{S \ {\underline{I}}}{S + {\underline{I}}} - \gamma \ {\underline{I}} \). Let us take \(L_2\) sufficiently large such that \( L_2 > e M \lambda ^* \sqrt{L_1 \log \left( L_1\right) }\).

  • For \( z \le \Big ( \dfrac{L_2}{e \ M \ \lambda ^*} \Big )^2 \):

\({\underline{I}} = 0\) so condition 4. is satisfied.

  • For \( z > \Big ( \dfrac{L_2}{e \ M \ \lambda ^*} \Big )^2 \):

The choice of \(L_2\) implies \(z > L_1 \log (L_1)\), which means \( {\underline{S}} = 1 - L_1 e^{- \dfrac{z}{L_1}}\).

We can reformulate condition (iv) as follows:

$$\begin{aligned} -v \ {\underline{I}}' - {\underline{I}}'' \le \beta _2 \ \dfrac{{\underline{S}} \ {\underline{I}}}{{\underline{S}} + {\underline{I}}} - \gamma \ {\underline{I}} \iff -v \ {\underline{I}}' - {\underline{I}}'' - (\beta _2 - \gamma ) \ {\underline{I}} \le - \beta _2 \ \dfrac{ {\underline{I}}^2}{{\underline{S}} + {\underline{I}}}. \end{aligned}$$
(61)

With \(L_3 = e M \lambda ^*\), and:

$$\begin{aligned}&{\underline{I}} = [ L_3 \ z - L_2 \ \sqrt{z} ] \ e^{- \lambda ^* z}, \end{aligned}$$
(62)
$$\begin{aligned}&{\underline{I}}' = [ L_3 - L_2 \ \dfrac{1}{2 \ \sqrt{z}} ] \ e^{-\lambda ^* z} - [ L_3 z - L_2 \ \sqrt{z} ] \ \lambda ^* \ e^{-\lambda ^* z},\end{aligned}$$
(63)
$$\begin{aligned}&{\underline{I}}'' = [ L_2 \ \dfrac{1}{4 \ z \ \sqrt{z}} ] \ e^{- \lambda ^* z} - 2 \ [ L_3 - L_2 \ \dfrac{1}{2 \ \sqrt{z}} ] \ \lambda ^* \ e^{ - \lambda ^* z} + [ L_3 z - L_2 \ \sqrt{z} ] \ (\lambda ^*)^2 \ e^{ - \lambda ^* z} \end{aligned}$$
(64)

On the one hand, we obtain the following identities:

$$\begin{aligned} -v \ {\underline{I}}' - {\underline{I}}'' - (\beta _2 - \gamma ) \ {\underline{I}}&= e^{ - \lambda ^* z} \Big [ L_3 \Big ( -v + v \ \lambda ^* \ z + 2 \ \lambda ^* - (\lambda ^*)^2 \ z - (\beta _2 - \gamma ) z \Big ) \end{aligned}$$
(65)
$$\begin{aligned}&\quad + \ L_2 \Big ( v \ \dfrac{1}{2 \ \sqrt{z}} - v \ \sqrt{z} \ \lambda ^* - \dfrac{1}{4 \ z \ \sqrt{z}} - \lambda ^* \ \dfrac{1}{ \sqrt{z}} + (\lambda ^*)^2 \ \sqrt{z} \nonumber \\&\quad + (\beta _2 - \gamma ) \ \sqrt{z} \Big ) \Big ] ,\end{aligned}$$
(66)
$$\begin{aligned}&= e^{- \lambda ^* z} \Big [ L_3 \Big ( ( 2 \lambda ^* - v ) - z \ \Big ( - v \lambda ^* + (\lambda ^*)^2 + (\beta _2 - \gamma ) \Big ) \Big )\end{aligned}$$
(67)
$$\begin{aligned}&\quad + \ L_2 \Big ( \sqrt{z} \ \Big ( - v \ \lambda ^* + (\lambda ^*)^2 + (\beta _2 - \gamma ) \Big ) + \dfrac{1}{2 \ \sqrt{z}} ( v - 2 \lambda ^*) \nonumber \\&\quad - \dfrac{1}{4 \ z \ \sqrt{z}} \Big ) \Big ] ,\end{aligned}$$
(68)
$$\begin{aligned}&= - L_2 \ e^{- \lambda ^* z} \ \dfrac{1}{4 \ z \ \sqrt{z}} \end{aligned}$$
(69)

On the other hand, we have:

$$\begin{aligned} - \beta _2 \ \dfrac{ {\underline{I}}^2}{{\underline{S}} + {\underline{I}}} = - \beta _2 \ \dfrac{[ L_3 \ z - L_2 \ \sqrt{z} ]^2 \ e^{ - 2 \lambda ^* z} }{ 1 - L_1 e^{- \dfrac{z}{L_1}} + [ L_3 \ z - L_2 \ \sqrt{z} ] \ e^{ - \lambda ^* z}}. \end{aligned}$$
(70)

We resume with the reformulation (61), which is now equivalent to the following:

$$\begin{aligned}&- L_2 \ e^{-\lambda ^* z} \ ( 1 - L_1 e^{- \dfrac{z}{L_1} } + [ L_3 \ z - L_2 \ \sqrt{z} ] \ e^{- \lambda ^* z} ) \le - 4 \beta _2 \ [ L_3 \ z - L_2 \ \sqrt{z} ]^2 \ e^{- 2 \lambda ^* z} \ z \ \sqrt{z}\end{aligned}$$
(71)
$$\begin{aligned}&\quad \iff 4 \beta _2 \ [ L_3 \ z - L_2 \ \sqrt{z} ]^2 \ e^{ - \lambda ^* z} \ z \ \sqrt{z} - L_2 \ [ L_3 \ z - L_2 \ \sqrt{z} ] \ e^{ - \lambda ^* z} \le \ L_2 \ ( 1 - L_1 e^{- \dfrac{z}{L_1}} )\nonumber \\ ,\end{aligned}$$
(72)
$$\begin{aligned}&\quad \iff 4 \beta _2 \ e^{ - \lambda ^* z} \ z^3 \sqrt{z} \ (L_3)^2 + e^{ - \lambda ^* z} \Big ( ( 1 - 8 \beta _2 \ z^2) \ L_3 \ L_2 \ z + (L_2)^2 \sqrt{z} (1 - 4 \beta _2 \ z^2 ) \Big ) \nonumber \\&\quad \le \ L_2 \ ( 1 - L_1 e^{- \dfrac{z}{L_1} } ). \end{aligned}$$
(73)

We may increase \(L_2\) such that \( 1 - 4 \beta _2 \ \Big ( \dfrac{L_2}{e M \lambda ^*} \Big )^4 \le 0\). Then, since \( z > (\dfrac{L_2}{e M \lambda ^*})^2\):

$$\begin{aligned} 1 - 8 \beta _2 \ z^2 \le 1 - 4 \beta _2 \ z^2 < 1 - 4 \beta _2 \ \Big ( \dfrac{L_2}{e M \lambda ^*} \Big ) ^4 \le 0. \end{aligned}$$
(74)

Since \(( 1 - 8 \beta _2 \ z^2)\) and \(( 1 - 4 \beta _2 \ z^2)\) are negative terms, we need to show \( L_2 \ ( 1 - L_1 e^{- \dfrac{z}{L_1} } ) \ge 4 \beta _2 \ e^{- \lambda ^* z} \ z^3 \sqrt{z} \ (L_3)^2 \).

Let \(g(z) = \beta _2 \ e^{ - \lambda ^* z} \ z^3 \sqrt{z} \ (L_3)^2\) be a \( \mathscr {C}^1([0;- \infty [)\) function. Since \( \lim \limits _{z \rightarrow 0} (g(z)) = 0 \) and \( \lim \limits _{z \rightarrow + \infty } (g(z)) = 0 \) there exists a constant C (which is independent from \(L_2\)) such that \(g(z) < C \ \ \forall z \ge 0 \). We finally increase \(L_2\) so that condition (iv) is verified. \(\square \)

1.2.2 D.2.2 Existence and positivity of a critical traveling wave solution

Now, exactly as in Zhou et al. (2019), we are in a position to define a set of functions

$$\begin{aligned} \Gamma = \{(S,I)\in B_\mu (\mathbb R,\mathbb {R}^2)\ |\ \underline{S}\le S\le \bar{S},\ \underline{I}\le I\le \bar{I}\}, \end{aligned}$$

where \(B_\mu (\mathbb R,\mathbb {R}^2)\) is the set of two-component continuous functions with each component growing at infinity slower than \(e^{\mu |z|}\), as well as an operator \(F:\Gamma \rightarrow C(\mathbb {R},\mathbb {R}^2)\) that will satisfy the assumptions of the Schauder fixed point theorem and whose fixed point in \(\Gamma \) will precisely be the solution (SI) we seek. Note that the inequalities 1., 2., 3., 4. (beginning of Sect. D.2.1) are precisely what we use to prove that \(F(\Gamma )\subset \Gamma \). Details can be found in Zhou et al. (2019).

The positivity of both S and I comes from the use of the strong maximum principle, again exactly as in Zhou et al. (2019).

E Study of the reaction term when \( r = + \infty \) in section 3.2

We are searching for conditions implying a pulled monostable wave, using criterion (10).

1.1 E.1 Conversion occurring in the zygote

We rewrite limit equation (22):

(75)

With \( \mathscr {A}_z:= s \ \big [ 2 (1-c) (1-h) - 1 \big ] \in [-s, s]\):

$$\begin{aligned} \partial _t p_{\textrm{D}} - \partial _{xx}^2 p_{\textrm{D}} \ = \dfrac{ \big ( - \mathscr {A}_z \ p_{\textrm{D}} + \frac{1}{2} ( \mathscr {A}_z - s) + c (1-s) \big ) \ (1-p_{\textrm{D}}) \ p_{\textrm{D}} }{-\mathscr {A}_z p_{\textrm{D}}^2 + (\mathscr {A}_z - s) \ p_{\textrm{D}} + 1}. \end{aligned}$$
(76)

Note that the mean fitness \( M_z(p_{\textrm{D}}) = -\mathscr {A}_z p_{\textrm{D}}^2 + (\mathscr {A}_z - s) \ p_{\textrm{D}} + 1 \in [1-s, 1]\)Footnote 4. When \( \mathscr {A}_z \ne 0 \), equation (76) can be rewritten:

$$\begin{aligned} \partial _t p_{\textrm{D}} - \partial _{xx}^2 p_{\textrm{D}} \ = \dfrac{ - \mathscr {A}_z \ ( p_{\textrm{D}} - {p_{\textrm{D}}^*}_z ) \ (1-p_{\textrm{D}}) \ p_{\textrm{D}} }{ -\mathscr {A}_z p_{\textrm{D}}^2 + (\mathscr {A}_z - s) \ p_{\textrm{D}} + 1} \quad \text {with} \quad {p_{\textrm{D}}^*}_z := \frac{1}{2} + \dfrac{ 2 c \ (1-s) - s }{ 2 \mathscr {A}_z}. \end{aligned}$$
(77)

Let us introduce \(s_1:= \dfrac{c}{1-h(1-c)}\) and \(s_{2,z}:= \dfrac{c}{2c + h(1-c)}\). Note that \( \mathscr {A}_z > 0 \iff s_1 < s_{2,z}\). We draw the reaction term regarding the sign of \(\mathscr {A}_z \) and the s values in Fig. 11.

When \(\mathscr {A}_z < 0\) and \( s \in (s_{2,z}, s_1)\), equation (22) admits two stable steady states (bistability). The final proportion will then strongly depend on the initial condition. On the other hand, when \(\mathscr {A}_z > 0\) and \( s \in (s_1, s_{2,z})\), the only possible equilibrium state is a coexistence state: the final proportion \(p_{\textrm{D}}\) will be strictly in between 0 and 1.

Independently of the sign of \(\mathscr {A}_z\), if \(s < \min (s_1,s_{2,z})\) the only stable steady state is \(p_{\textrm{D}}=1\) meaning that for an initial condition outside of the steady states, we expect that the drive always invades the population. If \(s > \max (s_1,s_{2,z})\), the only stable steady state is \(p_{\textrm{D}}=0\) meaning that for an initial condition outside of the steady states, we expect that the wild-type always invades the population.

Fig. 11
figure 11

Reaction term of equation (76) regarding the sign of \( \mathscr {A}_z = s [2 (1-c)(1-h) - 1] \) and the s values. The s threshold values are \(s_1 = \frac{c}{1-h(1-c)}\) and \(s_{2,z} = \frac{c}{2c + h(1-c)}\). The dots on the x axis correspond to the steady states: \(p_{\textrm{D}}=0\) and \(p_{\textrm{D}}=1\) in black, \( {p_{\textrm{D}}^*}_z = \frac{1}{2} + \frac{ 2 c \ (1-s) - s }{ 2 \mathscr {A}_z} \) in red when it exists (color figure online)

In case of bistability, the wave is always pushed (Hadeler and Rothe 1975); we can dismiss condition in the research of pulled monostable waves. We use criterion (10) on monostable cases, i.e. drive invasion , wild type invasion , and coexistence state (the numbers refer to the subgraphs in Fig. 11).

1.1.1 E.1.1 Monostable drive invasion

When \(\mathscr {A}_z=0\) and \( s< s_1 = s_{2,z} = \dfrac{2c}{2c + 1} \iff s < 2 c \ (1-s)\)

From equation (76), we have for all \(p_{\textrm{D}} \in [0,1]\):

$$\begin{aligned} \sigma (0) - (1-p_{\textrm{D}}) \ \sigma (p_{\textrm{D}})= & {} \Big ( c(1-s) - \dfrac{s}{2} \Big ) - \dfrac{ \Big ( c(1-s) - \dfrac{s}{2} \Big ) \ (1-p_{\textrm{D}}) }{1 - s p_{\textrm{D}}} \nonumber \\= & {} \Big ( c(1-s) - \dfrac{s}{2} \Big ) \ \dfrac{ (1- s) \ p_{\textrm{D}}}{ 1 - s p_{\textrm{D}}} \ge 0, \end{aligned}$$
(78)

where \(\sigma \) is the selection term defined by equation (8). Criterion (10) is verified.

When \(\mathscr {A}_z > 0 \) and \( {p_{\textrm{D}}^*}_z > 1 \)

From equation (77), we have for all \(p_{\textrm{D}} \in [0,1]\):

$$\begin{aligned} \sigma (0) - (1-p_{\textrm{D}}) \ \sigma (p_{\textrm{D}})= & {} \mathscr {A}_z {p_{\textrm{D}}^*}_z - \dfrac{ - \mathscr {A}_z \ ( p_{\textrm{D}} - {p_{\textrm{D}}^*}_z ) \ (1-p_{\textrm{D}}) }{ -\mathscr {A}_z p_{\textrm{D}}^2 + (\mathscr {A}_z - s ) \ p_{\textrm{D}} + 1} \nonumber \\= & {} \mathscr {A}_z p_{\textrm{D}} \dfrac{ - (\mathscr {A}_z {p_{\textrm{D}}^*}_z + 1) \ p_{\textrm{D}} + ( \mathscr {A}_z + 1 - s ) \ {p_{\textrm{D}}^*}_z + 1}{ -\mathscr {A}_z p_{\textrm{D}}^2 + (\mathscr {A}_z - s ) \ p_{\textrm{D}} + 1}.\nonumber \\ \end{aligned}$$
(79)

Note that \( -\mathscr {A}_z p_{\textrm{D}}^2 + (\mathscr {A}_z - s ) \ p_{\textrm{D}} + 1> (1-s) > 0\) and \( \mathscr {A}_z p_{\textrm{D}} > 0 \). The affine term \( - (\mathscr {A}_z {p_{\textrm{D}}^*}_z + 1) \ p_{\textrm{D}} + ( \mathscr {A}_z + 1 - s ) \ {p_{\textrm{D}}^*}_z + 1 \) decreases with \(p_{\textrm{D}}\). In order to show that it is positive for all \(p_{\textrm{D}} \in [0,1]\), we just need to verify that this it is true for \(p_{\textrm{D}}=1\):

$$\begin{aligned} - (\mathscr {A}_z {p_{\textrm{D}}^*}_z + 1) + ( \mathscr {A}_z + 1 - s ) \ {p_{\textrm{D}}^*}_z + 1= & {} ( 1-s ) \ {p_{\textrm{D}}^*}_z \ge 0 \quad \Rightarrow \quad \sigma (0) - (1-p_{\textrm{D}}) \ \sigma (p_{\textrm{D}}) \nonumber \\{} & {} \quad \ge 0 \forall p_{\textrm{D}} \in [0,1]. \end{aligned}$$
(80)

Criterion (10) is verified.

When \(\mathscr {A}_z <0\) and \({p_{\textrm{D}}^*}_z<0\) and \(s< s_{2,z} \iff 0 < c - 2 sc - sh + sch \)

We consider equation (79) with \(-\mathscr {A}_z p_{\textrm{D}}^2 + (\mathscr {A}_z - s ) \ p_{\textrm{D}} + 1> 1-s > 0 \) and \(\mathscr {A}_z p_{\textrm{D}} < 0\). The affine term \( - (\mathscr {A}_z {p_{\textrm{D}}^*}_z + 1) \ p_{\textrm{D}} + ( \mathscr {A}_z + 1 - s ) \ {p_{\textrm{D}}^*}_z + 1\) decreases with \(p_{\textrm{D}}\). In order to show that it is negative for all \(p_{\textrm{D}} \in [0,1]\), we introduce a condition implying the negativity for \(p_{\textrm{D}}=0\):

$$\begin{aligned}{} & {} ( \mathscr {A}_z + 1 - s ) \ {p_{\textrm{D}}^*}_z + 1 = \dfrac{ ( \mathscr {A}_z + 1 - s ) \ (\mathscr {A}_z + 2 c (1- s) -s + 2 \mathscr {A}_z )}{ 2 \mathscr {A}_z} < 0 \end{aligned}$$
(81)
$$\begin{aligned}{} & {} \iff \Big ( 1 - 2s [1 - (1-h)(1-c)] \Big ) \Big ( c - 2 sc - sh + sch \Big ) + s \ \big [ 2 (1-c) (1-h) - 1 \big ] > 0\nonumber \\ \end{aligned}$$
(82)

Criterion (10) is verified when condition (82) is true.

1.1.2 E.1.2 Monostable wild-type invasion

In case of a monostable wild-type invasion, we need to consider the wild-type proportion \(p_{\textrm{W}} = 1-p_{\textrm{D}} \in [0,1]\) and rewrite the equation (76):

$$\begin{aligned} - \partial _t p_{\textrm{W}} + \partial _{xx}^2 p_{\textrm{W}} \ {}&= \dfrac{ \big ( - \mathscr {A}_z \ (1-p_{\textrm{W}}) + \frac{1}{2} \ (\mathscr {A}_z - s) + c (1-s) \big ) \ (1-p_{\textrm{W}}) \ p_{\textrm{W}} }{-\mathscr {A}_z (1-p_{\textrm{W}})^2 + ( \mathscr {A}_z - s) (1-p_{\textrm{W}}) + 1} \nonumber \\ \iff \partial _t p_{\textrm{W}} - \partial _{xx}^2 p_{\textrm{W}} \ {}&= \dfrac{ \big ( - \mathscr {A}_z \ p_{\textrm{W}} + \frac{1}{2} \ (\mathscr {A}_z + s) - c (1-s) \big ) \ (1-p_{\textrm{W}}) \ p_{\textrm{W}} }{-\mathscr {A}_z p_{\textrm{W}}^2 + ( \mathscr {A}_z + s) \ p_{\textrm{W}} + (1-s)} \end{aligned}$$
(83)

When \( \mathscr {A}_z \ne 0\), equation (83) can be rewritten:

$$\begin{aligned} \partial _t p_{\textrm{W}} - \partial _{xx}^2 p_{\textrm{W}} \,= & {} \dfrac{ - \mathscr {A}_z \ (p_{\textrm{W}} - {p_{\textrm{W}}^*}_z) \ (1-p_{\textrm{W}}) \ p_{\textrm{W}} }{-\mathscr {A}_z p_{\textrm{W}}^2 + ( \mathscr {A}_z + s) \ p_{\textrm{W}} + (1-s)} \quad \text {with} \quad {p_{\textrm{W}}^*}_z \nonumber \\= & {} \frac{1}{2} - \dfrac{ 2 c \ (1-s) - s }{ 2 \mathscr {A}_z} = 1 - {p_{\textrm{D}}^*}_z \end{aligned}$$
(84)

When \( \mathscr {A}_z = 0\) and \( s< s_1 = s_{2,z} = \dfrac{2c}{2c + 1} \iff 2 c \ (1-s) < s \)

From equation (83) we have for all \(p_{\textrm{W}} \in [0,1]\):

$$\begin{aligned} \sigma (0) - (1-p_{\textrm{W}}) \ \sigma (p_{\textrm{W}}) = \big ( \dfrac{s}{2} - c(1-s) \big ) \Big ( \dfrac{1 }{ 1 - s } - \dfrac{1-p_{\textrm{W}}}{ s p_{\textrm{W}} + 1 - s } \Big ) \ge \dfrac{ p_{\textrm{W}} \ \big ( \dfrac{s}{2} - c(1-s) \big ) }{ 1 - s } \ge 0 \end{aligned}$$
(85)

Criterion (10) is verified.

When \( \mathscr {A}_z > 0\) and \( {p_{\textrm{W}}^*}_z = 1 - {p_{\textrm{D}}^*}_z > 1 \)

From equation (84), we have for all \(p_{\textrm{W}} \in [0,1]\):

$$\begin{aligned} \begin{aligned} \sigma (0) - (1-p_{\textrm{W}}) \ \sigma (p_{\textrm{W}})&= \dfrac{ \ \mathscr {A}_z \ {p_{\textrm{W}}^*}_z}{1-s} - \dfrac{ - \mathscr {A}_z \ ( p_{\textrm{W}} - {p_{\textrm{W}}^*}_z) \ (1-p_{\textrm{W}}) }{ -\mathscr {A}_z p_{\textrm{W}}^2 + ( \mathscr {A}_z + s) \ p_{\textrm{W}} + (1-s)} \\&= \ \mathscr {A}_z p_{\textrm{W}} \ \dfrac{ \ ( - {p_{\textrm{W}}^*}_z+ 1 - s ) \ p_{\textrm{W}} + {p_{\textrm{W}}^*}_z(\mathscr {A}_z + 1) + ( 1 - s ) }{( 1 - s )(-\mathscr {A}_z p_{\textrm{W}}^2 + ( \mathscr {A}_z + s) \ p_{\textrm{W}} + (1-s) )}. \end{aligned} \end{aligned}$$
(86)

Note that \( ( 1 - s )(-\mathscr {A}_z p_{\textrm{W}}^2 + ( \mathscr {A}_z + s) \ p_{\textrm{W}} + (1-s) )> (1-s)^2 > 0 \) and \( \mathscr {A}_z p_{\textrm{W}}> 0 \). As \( {p_{\textrm{W}}^*}_z> 1 \), the affine term \(( - {p_{\textrm{W}}^*}_z+ 1 - s ) \ p_{\textrm{W}} + {p_{\textrm{W}}^*}_z(\mathscr {A}_z + 1) + ( 1 - s ) \) decreases with \(p_{\textrm{W}}\). In order to show that it is positive for all \(p_{\textrm{W}} \in [0, 1]\), we just need to verify that this it is true for \(p_{\textrm{W}} = 1\):

$$\begin{aligned}{} & {} ( - {p_{\textrm{W}}^*}_z+ 1 - s ) + {p_{\textrm{W}}^*}_z(\mathscr {A}_z + 1) + ( 1 - s ) = 2 \ ( 1 - s ) + \mathscr {A}_z {p_{\textrm{W}}^*}_z \ge 0 \nonumber \\{} & {} \quad \quad \Rightarrow \quad \sigma (0) - (1-p_{\textrm{W}}) \ \sigma (p_{\textrm{W}}) \ge 0 \quad \forall p_{\textrm{W}} \in [0,1]. \end{aligned}$$
(87)

Criterion (10) is verified.

When \(\mathscr {A}_z <0\) and \({p_{\textrm{W}}^*}_z = 1 - {p_{\textrm{D}}^*}_z < 0 \)

We consider equation (86) with \( ( 1 - s )(-\mathscr {A}_z p_{\textrm{W}}^2 + ( \mathscr {A}_z + s) \ p_{\textrm{W}} + (1-s) )> (1-s)^2 > 0 \) and \( \mathscr {A}_z p_{\textrm{W}} < 0 \). The affine term \(( - {p_{\textrm{W}}^*}_z + 1 - s ) \ p_{\textrm{W}} + {p_{\textrm{W}}^*}_z (\mathscr {A}_z + 1) + ( 1 - s ) \) is strictly positive for \(p_{\textrm{W}}=1\), therefore criterion (10) is not verified.

1.1.3 E.1.3 Monostable coexistence state

When \(\mathscr {A}_z > 0 \) and \( 0< {p_{\textrm{D}}^*}_z = 1 - {p_{\textrm{W}}^*}_z < 1 \)

In the coexistence case, we have to verify that both waves, the drive invasion wave going to the right and the wild-type invasion wave going to the left, are pulled waves (see Fig. 3).

For the drive invasion wave we consider equation (79) with \( 0< {p_{\textrm{D}}^*}_z < 1 \) and \(p_{\textrm{D}} \in [0,{p_{\textrm{D}}^*}_z]\) (the term drive wave implies that the proportion of wild type increases after the wave passes; therefore the global stable steady state \({p_{\textrm{D}}^*}_z\) is also the maximum proportion). Once again, we need to prove that the affine term \( - (\mathscr {A}_z {p_{\textrm{D}}^*}_z + 1) \ p_{\textrm{D}} + ( \mathscr {A}_z + 1 - s ) \ {p_{\textrm{D}}^*}_z + 1 \) is positive. As it decreases with \(p_{\textrm{D}} \in [0,{p_{\textrm{D}}^*}_z]\), we determine its sign for \(p_{\textrm{D}} = {p_{\textrm{D}}^*}_z\):

$$\begin{aligned} - (\mathscr {A}_z {p_{\textrm{D}}^*}_z + 1) \ {p_{\textrm{D}}^*}_z + ( \mathscr {A}_z + 1 - s ) \ {p_{\textrm{D}}^*}_z + 1&= - \mathscr {A}_z \ ({p_{\textrm{D}}^*}_z)^2 + ( \mathscr {A}_z - s ) \ {p_{\textrm{D}}^*}_z + 1 \ge 1-s \ge 0 \nonumber \\&\Rightarrow \ \sigma (0) - (1-p_{\textrm{D}}) \sigma (p_{\textrm{D}}) \ge 0 \quad \forall p_{\textrm{D}} \in [0,{p_{\textrm{D}}^*}_z].\nonumber \\ \end{aligned}$$
(88)

Criterion (10) is verified for the drive wave.

For the wild-type invasion wave, we consider equation (86) with \( 0< {p_{\textrm{W}}^*}_z = 1 - {p_{\textrm{D}}^*}_z < 1 \) and \(p_{\textrm{W}} \in [0,{p_{\textrm{W}}^*}_z]\) (the term wild-type wave implies that the proportion of wild type increases after the wave passes; therefore the global stable steady state \({p_{\textrm{W}}^*}_z\) is also the maximum proportion). Once again, we need to prove that the affine term \(( - {p_{\textrm{W}}^*}_z + 1 - s ) \ p_{\textrm{W}} + {p_{\textrm{W}}^*}_z (\mathscr {A}_z + 1) + ( 1 - s )\) is positive. As it decreases with \(p_{\textrm{W}} \in [0,{p_{\textrm{W}}^*}_z]\), we determine its sign for \(p_{\textrm{W}} = {p_{\textrm{W}}^*}_z\):

$$\begin{aligned} \begin{aligned}&( - {p_{\textrm{W}}^*}_z + 1 - s ) \ {p_{\textrm{W}}^*}_z + (\mathscr {A}_z + 1) \ {p_{\textrm{W}}^*}_z + ( 1 - s ) = - ( {p_{\textrm{W}}^*}_z )^2 + (\mathscr {A}_z + 2 - s) \ {p_{\textrm{W}}^*}_z + 1 - s \\&\ge \min (1, \mathscr {A}_z + 2 ( 1- s)) \ge 0 \quad \Rightarrow \quad \sigma (0) - (1-p_{\textrm{W}}) \sigma (p_{\textrm{W}}) \ge 0 \quad \forall p_{\textrm{W}} \in [0,{p_{\textrm{W}}^*}_z]. \end{aligned} \end{aligned}$$
(89)

Criterion (10) is verified for the wild-type wave.

1.2 E.2 Conversion occurring in the germline

We rewrite limit equation (30):

$$\begin{aligned} \partial _t p_{\textrm{D}} - \partial _{xx}^2 p_{\textrm{D}} \ = \dfrac{ \Big ( - (1-2h) \ s \ p_{\textrm{D}} + [ (1-sh) (1+c) - 1 ] \Big ) \ p_{\textrm{D}} \ (1-p_{\textrm{D}}) }{- s (1-2h) p_{\textrm{D}}^2 - 2 s h p_{\textrm{D}} + 1} \end{aligned}$$
(90)

With \( \mathscr {A}_g:= s \ (1-2h) \in [ -s, s ]\):

$$\begin{aligned} \partial _t p_{\textrm{D}} - \partial _{xx}^2 p_{\textrm{D}} \ = \dfrac{ \big ( - \mathscr {A}_g \ p_{\textrm{D}} + \frac{1}{2} ( \mathscr {A}_g - s ) + c (1-sh) \big ) \ p_{\textrm{D}} \ (1-p_{\textrm{D}}) }{- \mathscr {A}_g p_{\textrm{D}}^2 + (\mathscr {A}_g - s) \ p_{\textrm{D}} + 1}. \end{aligned}$$
(91)

Note that the mean fitness \( M_g(p_{\textrm{D}}) = -\mathscr {A}_g p_{\textrm{D}}^2 + (\mathscr {A}_g - s) \ p_{\textrm{D}} + 1 \in [1-s, 1]\)Footnote 5.When \(\mathscr {A}_g \ne 0\), equation (91) can be rewritten:

$$\begin{aligned} \partial _t p_{\textrm{D}} - \partial _{xx}^2 p_{\textrm{D}} \,= & {} \dfrac{- \mathscr {A}_g \ ( p_{\textrm{D}} - {p_{\textrm{D}}^*}_g \big ) \ (1-p_{\textrm{D}}) \ p_{\textrm{D}}}{- \mathscr {A}_g p_{\textrm{D}}^2 + (\mathscr {A}_g - s) \ p_{\textrm{D}} + 1} = f_g(p_{\textrm{D}}) \quad \text {with} \nonumber \\{} & {} \quad {p_{\textrm{D}}^*}_g:= \dfrac{1}{2} + \dfrac{2 c \ (1- sh) - s }{ 2 \ \mathscr {A}_g}. \end{aligned}$$
(92)

Let us introduce \(s_1:= \dfrac{c}{1-h(1-c)}\) and \(s_{2,g}:= \dfrac{c}{2ch + h(1-c)} = \dfrac{ c }{h (1 + c) }\). We draw the reaction term regarding the sign of \( \mathscr {A}_g \) and the s values (in Fig. 12).

Fig. 12
figure 12

Reaction term of equation (91) regarding the sign of \( \mathscr {A}_g = s \ (1-2h) \) and the s values. The s threshold values are \(s_1 = \frac{c}{1-h(1-c)}\) and \(s_{2,g} = \frac{c}{2ch + h(1-c)} = \frac{ c }{ h (1 + c)}\). The dots on the x axis correspond to the steady states: \(p_{\textrm{D}}=0\) and \(p_{\textrm{D}}=1\) in black, \( {p_{\textrm{D}}^*}_g = \frac{1}{2} + \frac{2 c \ (1- sh) - s }{ 2 \ \mathscr {A}_g} \) in red when it exists (color figure online)

In case of bistability, the wave is always pushed (Hadeler and Rothe 1975); we can dismiss condition in the research of pulled monostable waves. We use criterion (10) on monostable cases, i.e. drive invasion , wild type invasion , and coexistence state (the numbers refer to the subgraphs in Fig. 12).

1.2.1 E.2.1 Monostable drive invasion

When \(\mathscr {A}_g=0 \iff h = \dfrac{1}{2}\) and \( s< s_1 = s_{2,g} = \dfrac{2 c }{c + 1} \iff \dfrac{s}{2} \ (c + 1) < c \)

From equation (91), we have for all \( p_{\textrm{D}} \in [0,1]\):

$$\begin{aligned} \sigma (0) - (1-p_{\textrm{D}}) \sigma (p)= & {} \Big ( c - \dfrac{s}{2} \ (c+1) \Big ) - \dfrac{ (c - \dfrac{s}{2} \ (c+1) ) \ (1-p_{\textrm{D}}) }{1 - s p_{\textrm{D}}} \nonumber \\= & {} \Big ( c - \dfrac{s}{2} \ (c+1) \Big ) \ \dfrac{ (1- s) \ p_{\textrm{D}} }{ 1 - s p_{\textrm{D}}} \ge 0, \end{aligned}$$
(93)

where \(\sigma \) is the selection term defined by equation (8). Criterion (10) is verified.

When \(\mathscr {A}_g > 0 \) and \( {p_{\textrm{D}}^*}_g > 1 \)

From equation (92), we have for all \(p_{\textrm{D}} \in [0,1]\):

$$\begin{aligned} \sigma (0) - (1-p_{\textrm{D}}) \sigma (p)= & {} \mathscr {A}_g {p_{\textrm{D}}^*}_g - \dfrac{ - \mathscr {A}_g \ ( p_{\textrm{D}} - {p_{\textrm{D}}^*}_g ) \ (1-p_{\textrm{D}}) }{ -\mathscr {A}_g p_{\textrm{D}}^2 + (\mathscr {A}_g - s ) \ p_{\textrm{D}} + 1} \nonumber \\= & {} \mathscr {A}_g p_{\textrm{D}} \dfrac{ - (\mathscr {A}_g {p_{\textrm{D}}^*}_g + 1) \ p_{\textrm{D}} + ( \mathscr {A}_g + 1 - s ) \ {p_{\textrm{D}}^*}_g + 1}{ -\mathscr {A}_g p_{\textrm{D}}^2 + (\mathscr {A}_g - s ) \ p_{\textrm{D}} + 1}.\nonumber \\ \end{aligned}$$
(94)

Note that \(-\mathscr {A}_g p_{\textrm{D}}^2 + (\mathscr {A}_g - s ) p_{\textrm{D}} + 1> (1-s) > 0\) and \( \mathscr {A}_g p_{\textrm{D}} > 0 \). The affine term \(- (\mathscr {A}_g {p_{\textrm{D}}^*}_g + 1) \ p_{\textrm{D}} + ( \mathscr {A}_g + 1 - s ) \ {p_{\textrm{D}}^*}_g + 1\) decreases with \(p_{\textrm{D}}\). In order to show that it is positive for all \(p_{\textrm{D}} \in [0,1]\), we just need to verify that this it is true for \(p_{\textrm{D}}=1\):

$$\begin{aligned}{} & {} - (\mathscr {A}_g {p_{\textrm{D}}^*}_g + 1) + ( \mathscr {A}_g + 1 - s ) \ {p_{\textrm{D}}^*}_g + 1 = ( 1-s ) \ {p_{\textrm{D}}^*}_g \ge 0 \quad \Rightarrow \quad \sigma (0) - (1-p_{\textrm{D}}) \sigma (p) \nonumber \\{} & {} \quad \ge 0 \quad \forall p_{\textrm{D}} \in [0,1]. \end{aligned}$$
(95)

Criterion (10) is verified.

When \(\mathscr {A}_g <0\) and \({p_{\textrm{D}}^*}_g <0\) and \(s< s_{2,g} \iff 0 < c - s h (1 + c )\)

We consider equation (94) with \(-\mathscr {A}_g p_{\textrm{D}}^2 + (\mathscr {A}_g - s ) \ p_{\textrm{D}} + 1> 1-s > 0 \) and \(\mathscr {A}_g p_{\textrm{D}} < 0\). The affine term \( - (\mathscr {A}_g {p_{\textrm{D}}^*}_g + 1) \ p_{\textrm{D}} + ( \mathscr {A}_g + 1 - s ) \ {p_{\textrm{D}}^*}_g + 1\) decreases with \(p_{\textrm{D}}\). In order to show that it is negative for all \(p_{\textrm{D}} \in [0,1]\), we introduce a condition implying the negativity for \(p=0\):

$$\begin{aligned}{} & {} ( \mathscr {A}_g + 1 - s ) \ {p_{\textrm{D}}^*}_g + 1 = \dfrac{(1-2sh) (\mathscr {A}_g + 2 c (1- sh)-s) + 2 \mathscr {A}_g }{ 2 \mathscr {A}_g} < 0 \end{aligned}$$
(96)
$$\begin{aligned}{} & {} \iff ( 1 - 2 sh) (c - sh ( c + 1 ) ) + s (1 - 2h) > 0 \end{aligned}$$
(97)

Criterion (10) is verified when condition (97) is true.

1.2.2 E.2.2 Monostable wild-type invasion

In case of a monostable wild-type invasion, we need to consider the wild-type proportion \(p_{\textrm{W}} = 1-p_{\textrm{D}} \in [0,1]\) and rewrite the equation (91):

$$\begin{aligned} - \partial _t p_{\textrm{W}} + \partial _{xx}^2 p_{\textrm{W}} \ {}&= \dfrac{ \big ( - \mathscr {A}_g \ (1-p_{\textrm{W}}) + \frac{1}{2} \ (\mathscr {A}_g - s) + c (1-sh) \big ) \ (1-p_{\textrm{W}}) \ p_{\textrm{W}} }{-\mathscr {A}_g (1-p_{\textrm{W}})^2 + ( \mathscr {A}_g - s) (1-p_{\textrm{W}}) + 1} \nonumber \\ \iff \partial _t p_{\textrm{W}} - \partial _{xx}^2 p_{\textrm{W}} \ {}&= \dfrac{ \big ( - \mathscr {A}_g \ p_{\textrm{W}} + \frac{1}{2} \ (\mathscr {A}_g + s) - c (1-sh) \big ) \ (1-p_{\textrm{W}}) \ p_{\textrm{W}} }{-\mathscr {A}_g p_{\textrm{W}}^2 + ( \mathscr {A}_g + s) \ p_{\textrm{W}} + (1-s)} \end{aligned}$$
(98)

When \( \mathscr {A}_g \ne 0\), equation (98) can be rewritten:

$$\begin{aligned} \partial _t p_{\textrm{W}} - \partial _{xx}^2 p_{\textrm{W}} \,= & {} \dfrac{ - \mathscr {A}_g \ (p_{\textrm{W}} - {p_{\textrm{W}}^*}_g) \ (1-p_{\textrm{W}}) \ p_{\textrm{W}} }{-\mathscr {A}_g p_{\textrm{W}}^2 + ( \mathscr {A}_g + s) \ p_{\textrm{W}} + (1-s)} \quad \text {with} \nonumber \\{} & {} \quad {p_{\textrm{W}}^*}_g = \frac{1}{2} - \dfrac{ 2 c \ (1-sh) - s }{ 2 \mathscr {A}_g} = 1 - {p_{\textrm{D}}^*}_g \end{aligned}$$
(99)

When \( \mathscr {A}_g = 0 \iff h=\dfrac{1}{2}\) and \( s_1 = s_{2,g} = \dfrac{2c}{c + 1}< s \iff c < \dfrac{s}{2} \ (c+1) \)

From equation (98) we have for all \(p_{\textrm{W}} \in [0,1]\):

$$\begin{aligned} \sigma (0) - (1-p_{\textrm{W}}) \sigma (p_{\textrm{W}}) = \big ( \dfrac{s}{2} \ (c+1) - c \big ) \Big ( \dfrac{1}{ 1 - s } - \dfrac{1-p_{\textrm{W}}}{ s p_{\textrm{W}} + 1 - s } \Big ) \ge \dfrac{ p_{\textrm{W}} \ \big ( \dfrac{s}{2} \ (c+1) - c \big ) }{ 1 - s } \ge 0 \end{aligned}$$
(100)

Criterion (10) is verified.

When \( \mathscr {A}_g > 0\) and \( {p_{\textrm{W}}^*}_g = 1 - {p_{\textrm{D}}^*}_g > 1 \)

From equation (99), we have for all \(p_{\textrm{W}} \in [0,1]\):

$$\begin{aligned} \begin{aligned} \sigma (0) - (1-p_{\textrm{W}}) \sigma (p_{\textrm{W}})&= \dfrac{ \ \mathscr {A}_g \ {p_{\textrm{W}}^*}_g}{1-s} - \dfrac{ - \mathscr {A}_g \ ( p_{\textrm{W}} - {p_{\textrm{W}}^*}_g ) \ (1-p_{\textrm{W}}) }{ -\mathscr {A}_g p_{\textrm{W}}^2 + ( \mathscr {A}_g + s) \ p_{\textrm{W}} + (1-s)} \\&= \ \mathscr {A}_g p_{\textrm{W}} \ \dfrac{ \ ( - {p_{\textrm{W}}^*}_g + 1 - s ) \ p_{\textrm{W}} + {p_{\textrm{W}}^*}_g (\mathscr {A}_g + 1) + ( 1 - s ) }{( 1 - s )(-\mathscr {A}_g p_{\textrm{W}}^2 + ( \mathscr {A}_g + s) \ p_{\textrm{W}} + (1-s) )}. \end{aligned} \end{aligned}$$
(101)

Note that \( ( 1 - s )(-\mathscr {A}_g p_{\textrm{W}}^2 + ( \mathscr {A}_g + s) \ p_{\textrm{W}} + (1-s) )> (1-s)^2 > 0 \) and \( \mathscr {A}_g p_{\textrm{W}}> 0 \). As \( {p_{\textrm{W}}^*}_g > 1 \), the affine term \(( - {p_{\textrm{W}}^*}_g + 1 - s ) \ p_{\textrm{W}} + {p_{\textrm{W}}^*}_g (\mathscr {A}_g + 1) + ( 1 - s ) \) decreases with \(p_{\textrm{W}}\). In order to show that it is positive for all \(p_{\textrm{W}} \in [0, 1]\), we just need to verify that this it is true for \(p_{\textrm{W}} = 1\):

$$\begin{aligned}{} & {} ( - {p_{\textrm{W}}^*}_g + 1 - s ) + {p_{\textrm{W}}^*}_g (\mathscr {A}_g + 1) + ( 1 - s ) = 2 \ ( 1 - s ) + \mathscr {A}_g {p_{\textrm{W}}^*}_g \ge 0 \nonumber \\{} & {} \quad \quad \Rightarrow \sigma (0) - (1-p_{\textrm{W}}) \sigma (p_{\textrm{W}}) \ge 0 \quad \forall p_{\textrm{W}} \in [0, 1]. \end{aligned}$$
(102)

Criterion (10) is verified.

When \(\mathscr {A}_g <0\) and \({p_{\textrm{W}}^*}_g = 1 - {p_{\textrm{D}}^*}_z < 0\)

We consider equation (101) with \( ( 1 - s )(-\mathscr {A}_g p_{\textrm{W}}^2 + ( \mathscr {A}_g + s) \ p_{\textrm{W}} + (1-s) )> (1-s)^2 > 0 \) and \( \mathscr {A}_g p_{\textrm{W}}< 0 \). The affine term \(( - {p_{\textrm{W}}^*}_g + 1 - s ) \ p_{\textrm{W}} + {p_{\textrm{W}}^*}_g (\mathscr {A}_g + 1) + ( 1 - s ) \) is strictly positive for \(p_{\textrm{W}}=1\), therefore criterion (10) is not verified.

1.2.3 E.2.3 Monostable coexistence state

When \(\mathscr {A}_g > 0 \) and \( 0< {p_{\textrm{D}}^*}_g = 1 - {p_{\textrm{W}}^*}_g < 1 \)

In the coexistence case, we have to verify that both sub-traveling waves, the drive invasion wave going to the right and the wild-type invasion wave going to the left, are pulled waves (see Fig. 3).

For the drive invasion wave we consider equation (94) with \( 0< {p_{\textrm{D}}^*}_g < 1 \) and \(p_{\textrm{D}} \in [0,{p_{\textrm{D}}^*}_g]\) (the term drive wave implies that the proportion of wild type increases after the wave passes; therefore the global stable steady state \({p_{\textrm{D}}^*}_g\) is also the maximum proportion). Once again, we need to prove that the affine term \( - (\mathscr {A}_g {p_{\textrm{D}}^*}_g + 1) \ p_{\textrm{D}} + ( \mathscr {A}_g + 1 - s ) \ {p_{\textrm{D}}^*}_g + 1 \) is positive. As it decreases with \(p_{\textrm{D}} \in [0,{p_{\textrm{D}}^*}_g]\), we determine its sign for \(p_{\textrm{D}} = {p_{\textrm{D}}^*}_g\):

$$\begin{aligned} \begin{aligned}&- (\mathscr {A}_g {p_{\textrm{D}}^*}_g + 1) \ {p_{\textrm{D}}^*}_g + ( \mathscr {A}_g + 1 - s ) \ {p_{\textrm{D}}^*}_g + 1 = - \mathscr {A}_g \ ({p_{\textrm{D}}^*}_g)^2\\&\quad \quad + ( \mathscr {A}_g - s ) \ {p_{\textrm{D}}^*}_g + 1 \ge 1-s \ge 0 \\&\quad \Rightarrow \quad \sigma (0) - (1-p_{\textrm{D}}) \sigma (p_{\textrm{D}}) \ge 0 \quad \forall p_{\textrm{D}} \in [0,{p_{\textrm{D}}^*}_g]. \end{aligned} \end{aligned}$$
(103)

Criterion (10) is verified for the drive wave.

For the wild-type invasion wave, we consider equation (101) with \( 0< {p_{\textrm{W}}^*}_g = 1 - {p_{\textrm{D}}^*}_g < 1 \) and \(p_{\textrm{W}} \in [0,{p_{\textrm{W}}^*}_g]\) (the term wild-type wave implies that the proportion of wild type increases after the wave passes; therefore the global stable steady state \({p_{\textrm{W}}^*}_g\) is also the maximum proportion). Once again, we need to prove that the affine term \(( - {p_{\textrm{W}}^*}_g + 1 - s ) \ p_{\textrm{W}} + {p_{\textrm{W}}^*}_g (\mathscr {A}_g + 1) + ( 1 - s )\) is positive. As it decreases with \(p_{\textrm{W}} \in [0,{p_{\textrm{W}}^*}_g]\), we determine its sign for \(p_{\textrm{W}} = {p_{\textrm{W}}^*}_g\):

$$\begin{aligned} \begin{aligned}&( - {p_{\textrm{W}}^*}_g + 1 - s ) \ {p_{\textrm{W}}^*}_g + (\mathscr {A}_g + 1) \ {p_{\textrm{W}}^*}_g + ( 1 - s ) = - ( {p_{\textrm{W}}^*}_g )^2 + (\mathscr {A}_g + 2 - s) \ {p_{\textrm{W}}^*}_g + 1 - s \\&\ge \min (1, \mathscr {A}_g + 2 ( 1- s)) \ge 0 \quad \Rightarrow \quad \sigma (0) - (1-p_{\textrm{W}}) \sigma (p_{\textrm{W}}) \ge 0 \quad \forall p_{\textrm{W}} \in [0,{p_{\textrm{W}}^*}_g]. \end{aligned} \end{aligned}$$
(104)

Criterion (10) is verified for the wild-type wave.

F Heatmap supplementary materials

1.1 F.1 Effect of fitness disadvantage (s) and dominance coefficient (h) on drive dynamics, for r \( = + \infty \)

Below, we compute heatmaps indicating the stability regime of systems (4) and (5) when \(r=+\infty \), depending on the values of (h, s) and for a fixed value of c.

Fig. 13
figure 13

Effect of fitness disadvantage (s) and dominance coefficient (h) on drive dynamics for system (4) (when conversion occurs in the zygote) when \( r = + \infty \). Parameters for Fig. 4 (\(c=0.25\) and \(h=0.1\)) and Fig. 5 (\(c=0.75\) and \(h=0.1\)) are materialized by dotted lines

Fig. 14
figure 14

Effect of fitness disadvantage (s) and dominance coefficient (h) on drive dynamics for system (5) (when conversion occurs in the germline) when \( r = + \infty \). Parameters for Fig. 6A (\(c=0.25\) and \(h=0.3\)) and Fig. 6B (\(c=0.25\) and \(h=0.75\)) are materialized by dotted lines

Such a figure has already been computed in Rode et al. (2019) for \(c=0.85\), and conversion occurring in the germline.

1.2 F.2 Heatmap lines

1.2.1 F.2.1 Pure drive line

Consider model (6) with \(n_{\textrm{WW}} = 0\). A well-mixed population containing only drive homozygous individuals will persist in the environment if its equilibrium state \(n_{\textrm{DD}}^*\) is strictly positive, i.e. if:

$$\begin{aligned} n_{\textrm{DD}}^* = \min \Big ( 0, 1 - \frac{s}{r \ (1-s)} \Big )> 0 \quad \iff \quad r > \dfrac{s}{1-s} \end{aligned}$$
(105)

In case of partial conversion, calculations give the same threshold (consider models (4) and (5) with \(n_{\textrm{DW}} = 0\) and \(n_{\textrm{WW}} = 0\)).

1.2.2 F.2.2 Composite persistence line

Similarly, in case of coexistence, a well-mixed population will persist in the environment only if its equilibrium state \(n^*\) is strictly positive. Using Mathematica, we compute this population density equilibrium when conversion occurs in the zygote (\(n_z^*\)) or in the germline (\(n_g^*\)) based on systems (21) and (29). We obtain the following:

$$\begin{aligned} n_z^* = \min \Big ( 0, 1 - \dfrac{1 - M_z({p_{\textrm{D}}^*}_z)}{ r M_z({p_{\textrm{D}}^*}_z)} \Big ) \quad \text {and} \quad n_g^* = \min \Big (0, 1 - \dfrac{1 - M_g({p_{\textrm{D}}^*}_g)}{ r M_g({p_{\textrm{D}}^*}_g)} \Big ), \end{aligned}$$
(106)

where the mean fitness \(M_z\) and \(M_g\) (already defined in Appendix E) are given by:

$$\begin{aligned} M_z(p_{\textrm{D}}) = -\mathscr {A}_z p_{\textrm{D}}^2 + (\mathscr {A}_z - s) \ p_{\textrm{D}} + 1 \quad \text {and} \quad M_g(p_{\textrm{D}}) = -\mathscr {A}_g p_{\textrm{D}}^2 + (\mathscr {A}_g - s) \ p_{\textrm{D}} + 1, \end{aligned}$$
(107)

and the proportions \({p_{\textrm{D}}^*}_z\) and \({p_{\textrm{D}}^*}_g\) (already defined in Appendix E) are given by:

$$\begin{aligned} {p_{\textrm{D}}^*}_z = \frac{1}{2} + \dfrac{ 2 c \ (1-s) - s }{ 2 \mathscr {A}_z} \quad \text {and} \quad {p_{\textrm{D}}^*}_g := \dfrac{1}{2} + \dfrac{2 c \ (1- sh) - s }{ 2 \ \mathscr {A}_g} \end{aligned}$$
(108)

with,

$$\begin{aligned} \mathscr {A}_z = s \ \big [ 2 (1-c) (1-h) - 1 \big ] \quad \text {and} \quad \mathscr {A}_g = s \ (1-2h). \end{aligned}$$
(109)

Finally, the threshold values for r are given by:

$$\begin{aligned} n_z^*>0 \quad \iff \quad r > \dfrac{1-M_z({p_{\textrm{D}}^*}_z)}{M_z({p_{\textrm{D}}^*}_z)} \end{aligned}$$
(110)

when conversion occurs in the zygote and,

$$\begin{aligned} n_g^*>0 \quad \iff \quad r > \dfrac{1-M_g({p_{\textrm{D}}^*}_g)}{M_g({p_{\textrm{D}}^*}_g)} \end{aligned}$$
(111)

when conversion occurs in the germline.

Rights and permissions

Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Kläy, L., Girardin, L., Calvez, V. et al. Pulled, pushed or failed: the demographic impact of a gene drive can change the nature of its spatial spread. J. Math. Biol. 87, 30 (2023). https://doi.org/10.1007/s00285-023-01926-4

Download citation

  • Received:

  • Revised:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1007/s00285-023-01926-4

Mathematics Subject Classification

Navigation