Origin and persistence of polymorphism in loci targeted by disassortative preference: a general model

Abstract

The emergence and persistence of polymorphism within populations generally requires specific regimes of natural or sexual selection. Here, we develop a unified theoretical framework to explore how polymorphism at targeted loci can be generated and maintained by either disassortative mating choice or balancing selection due to, for example, heterozygote advantage. To this aim, we model the dynamics of alleles at a single locus A in a population of haploid individuals, where reproductive success depends on the combination of alleles carried by the parents at locus A. Our theoretical study of the model confirms that the conditions for the persistence of a given level of allelic polymorphism depend on the relative reproductive advantages among pairs of individuals. Interestingly, equilibria with unbalanced allelic frequencies were shown to emerge from successive introduction of mutants. We then investigate the role of the function linking allelic divergence to reproductive advantage on the evolutionary fate of alleles within the population. Our results highlight the significance of the shape of this function for both the number of alleles maintained and their level of genetic divergence. Large number of alleles are maintained with substantial replacement of alleles, when disassortative advantage slowly increases with allelic differentiation . In contrast, few highly differentiated alleles are predicted to be maintained when genetic differentiation has a strong effect on disassortative advantage. These opposite effects predicted by our model explain how disassortative mate choice may lead to various levels of allelic differentiation and polymorphism, and shed light on the effect of mate preferences on the persistence of balanced and unbalanced polymorphism in natural population.

Appendices

Birth-and-death model and convergence to dynamical system

The dynamical system (1) governing the population dynamics can be obtained as a large population limit of a stochastic individual based model (similarly as in Coron et al. (2018), Smadi et al. (2018) for instance).

More precisely, let K be a large parameter that gives the order size of the population. The microscopic population is represented by the process $(Z^K_{i}(t))_{i\in \{1,\ldots ,k\}}$ where $Z^K_{i}(t)$ gives the size, divided by K, of the population of type i individuals at time t. Each individual of type i reproduces with rate $\beta _i$, that is, after a time distributed as an exponential random variable of parameter $\beta _i$. As explained in the main text, the second parent is chosen uniformly at random among all individuals alive in the population. If the second parent is of type j, with probability $p_{ij}$ this mating event leads to the birth of a new individual, being of type i with probability one half, and of type j with probability one half when there is no mutation (otherwise the mutation is applied to the genotype i or j with the same probability). All the indivivuals have the same death rate, which is given by

$$\begin{aligned} d+c \sum _{i=1}^k Z^K_{i}(t) \end{aligned}$$

at time t. Under the assumption that the sequence of initial conditions $\big ((Z^K_{i}(0))_{i\in \{1,\ldots ,k\}}\big )_{K\in \mathbb {N}}$ converges (in probability) when K goes to infinity, the sequence of stochastic functions $\big ((Z^K_{i}(t))_{,i\in \{1,\ldots ,k\}}, t\in [0,T]\big )_{K\in \mathbb {N}}$ also converges (in probability for the uniform convergence) to the trajectory of the deterministic system (1).

Proof of the main results

This section gathers the theoretical results on the behavior of the solutions to system (2):

$$\begin{aligned} \dot{z}_i(t) = z_i(t) \left( b\sum _{j=1}^k (1+s_{ij}) \frac{z_j(t)}{z(t)} -d-cz(t) \right) , \quad \forall 1\le i\le k. \end{aligned}$$

Recall the matrix $M=(s_{ij})_{1\le i, j, \le k}$

1.1 Existence of positive equilibria

Our first result states conditions for the existence of a positive equilibrium, meaning that each coordinate of the vector is strictly positive.

Proposition B.1

If $\det (M)\ne 0$, the system admits a positive equilibrium if and only if $M^{-1}\textbf{1}>0$. In this case, the equilibrium is unique and is given by

$$\begin{aligned} Z^*=\frac{1}{c \textbf{1}^T M^{-1}\textbf{1}}\left( b+\frac{b}{\textbf{1}^T M^{-1}\textbf{1}}-d\right) M^{-1}\textbf{1}. \end{aligned}$$

(11)

If $\det (M)=0$, the system may have no equilibrium or an infinite number of equilibria.

Proof of Proposition B.1

$\underline{\textbf{Case}\,\det (M)\ne 0:}$ Firstly, we assume that there exists a positive equilibrium $Z^*=(z^*_i)_{i=1,..,k}$. From (2), the following equality holds for every $i\in \{1,..,k\}$,

$$\begin{aligned} \sum _{j=1}^ks_{ij}z^*_{j}= \frac{{z^*}(cz^*-b+d)}{b}>0, \end{aligned}$$

(12)

with $z^*=\sum _{i=1}^k z^*_{i}$. Note that this equation implies that $cz^*>b-d$, since $s_{ij}>0$. Using a matrix formulation, we have the following linear system:

$$\begin{aligned} M Z^*= \frac{{z^*}(cz^*-b+d)}{b} \textbf{1}. \end{aligned}$$

(13)

As $\det (M)\ne 0$, the matrix M has an inverse and

$$\begin{aligned} Z^* = \frac{{z^*}(cz^*-b+d)}{b} M^{-1} \textbf{1}>0. \end{aligned}$$

Since $\frac{{z^*}(cz^*-b+d)}{b} >0$, we thus deduce that $M^{-1} \textbf{1}>0$.

Moreover by summing all coordinates, we find

$$\begin{aligned} z^*=\frac{{z^*}(cz^*-b+d)}{b} \textbf{1}^T M^{-1}\textbf{1}\quad \Leftrightarrow \quad z^*=\frac{1}{c}\left( b+\frac{b}{\textbf{1}^T M^{-1}\textbf{1}}-d\right) , \end{aligned}$$

and finally, the equilibrium is unique and defined by (6).

On the other hand, if we assume that $M^{-1}\textbf{1}>0$, then we can define $Z^*$ by (6) and, using the above computation, it is straightforward to verify that it is a positive equilibrium of (2). This ends the proof for the case $\det (M)\ne 0$.

$\underline{{\textbf {Case}}\,\det (M)=0:}$

Let us verify whether there exists a vector X with positive coordinates such that $MX=\textbf{1}$.

Indeed, if such a vector X does not exist, then (12) can never be satisfied and there is no equilibrium of (2).

On the contrary, if we can find such a vector X, we set

$$\begin{aligned} Z^*=\frac{X}{c\ \textbf{1}^T X}\left( b+\frac{b}{\textbf{1}^T X}-d\right) . \end{aligned}$$

Then it is straightforward to verify that $MZ^*=z^*\frac{cz^*+d-b}{b}\textbf{1}$, with $z^*=\textbf{1}^TZ^*$. Hence, $Z^*$ satisfies (12) and is a positive equilibrium of (2). However, this equilibrium is not unique. Indeed, $\textbf{1}$ belongs to the image of M, and since M is a symmetric matrix, any vector Y in the kernel of M is orthogonal to $\textbf{1}$, i.e. $\textbf{1}^TY=0$. Then, for all vectors $Y\in Ker(M)$, $M(Z^*+Y)=z^*\frac{cz^*+d-b}{b}\textbf{1}$, with $z^*= \textbf{1}^T Z^*=\textbf{1}^T(Z^*+Y)$, and $Z^*+Y$ is also an equilibrium of (2).

We now construct examples for the two possibilities. Assume that M can be written as

$$\begin{aligned} M=\begin{pmatrix} A &{}\quad C \\ C^T &{}\quad 0 \end{pmatrix}, \text { with } C=\begin{pmatrix} s_{1k}\\ ..\\ s_{k-1 k} \end{pmatrix}, \end{aligned}$$

$\det (A)\ne 0$ and such that $A^{-1}\textbf{1}>0$. Since A is a symmetric matrix, there exists an orthonormal basis $(V_i)_{i=1,..,k-1}$ of $\mathbb {R}^{k-1}$ where $V_i$ is an eigenvector of A associated to its non-null eigenvalue $\lambda _i$ for all $i\in \{1,..,k-1\}$. Let us decompose C in this basis $C=\sum _{i=1}^{k-1}\alpha _i V_i$.

Our aim is to construct C, and thus to select $(\alpha _i)_{i=1,..,k-1}$, to find the examples. First, we have to select $(\alpha _i)_{i=1,..,k-1}$ such that $C=\sum _{i=1}^{k-1}\alpha _iV_i>0$. Then using the Schur complement, we have that $\det (M)=-(C^TA^{-1}C)\times \det (A)$ and thus $\det (M)=0$ if and only if

$$\begin{aligned} C^TA^{-1}C=\sum _{i=1}^{k-1}\frac{\alpha _i^2}{\lambda _i}=0. \end{aligned}$$

This gives us a way to construct a matrix M with $\det (M)=0$.

Finally, to ensure that there exists $Y>0$ with $MY=\textbf{1}$, let us consider vectors depending on a parameter $\eta \in \mathbb {R}$ such that

$$\begin{aligned} Y(\eta )=\begin{pmatrix} A^{-1} \textbf{1}\\ 0 \end{pmatrix} + \eta \begin{pmatrix} -A^{-1}C \\ 1 \end{pmatrix}. \end{aligned}$$

Let us study the image by M of this vector

$$\begin{aligned} MY(\eta )=\begin{pmatrix} \textbf{1}\\ C^T A^{-1} \textbf{1}\end{pmatrix}. \end{aligned}$$

Thus if $C^T A^{-1} \textbf{1}=\sum _{i=1}^{k-1}\frac{\alpha _i\omega _i}{\lambda _i}=1$ where $\textbf{1}=\sum _{i=1}^{k-1}\omega _i V_i$, then we can find $\eta >0$ sufficiently small, such that $MY(\eta )=\textbf{1}$ with $Y(\eta )>0$.

With this in mind, we can explicit an example for each possibility. Let

$$\begin{aligned} M= \begin{pmatrix} 0 &{}\quad 1 &{}\quad 1 &{}\quad \frac{2 \sqrt{2}+\sqrt{3}+1}{3 \sqrt{2}} \\ 1 &{}\quad 0 &{}\quad 1 &{}\quad \frac{2 \sqrt{2}-\sqrt{3}+1}{3 \sqrt{2}} \\ 1 &{}\quad 1 &{}\quad 0 &{}\quad \frac{2 \sqrt{2}-2}{3 \sqrt{2}} \\ \frac{2 \sqrt{2}+\sqrt{3}+1}{3 \sqrt{2}} &{}\quad \frac{2 \sqrt{2}-\sqrt{3}+1}{3 \sqrt{2}} &{}\quad \frac{2 \sqrt{2}-2}{3 \sqrt{2}} &{}\quad 0 \\ \end{pmatrix} \text { and} \quad X=\begin{pmatrix} \frac{1}{12} \left( \sqrt{2}+\sqrt{6}+4\right) \\ \frac{1}{12} \left( \sqrt{2}-\sqrt{6}+4\right) \\ \frac{1}{6} \left( 2-\sqrt{2}\right) \\ \frac{1}{2} \end{pmatrix}. \end{aligned}$$

M and X have non-negative coordinates, $\det (M)=0$ and $MX=\textbf{1}$. Thus, in this example, the associated system of equations (2) has an infinite number of positive stationary states.

Now, let $a>0$ and

$$\begin{aligned} M=\begin{pmatrix} 0 &{}\quad a &{}\quad a &{}\quad a \\ a &{}\quad 0 &{}\quad a &{}\quad a \\ a &{}\quad a &{}\quad 0 &{}\quad 4a \\ a &{}\quad a &{}\quad 4a &{}\quad 0 \end{pmatrix} \text { and} \quad X_0=\begin{pmatrix} -2 \\ -2 \\ 1 \\ 1 \end{pmatrix}. \end{aligned}$$

In this example, $\det (M)=0$, $Ker(M)=Vect(X_0)=Im(M)^T$, however $\textbf{1}^TX_0\ne 0$ hence $\textbf{1}$ does not belong to the image of M. In other words, the associated system of equations (2) has no positive stationary state. $\square$

1.2 Convergence toward the positive equilibria

We now aim at describing the dynamics of the population described by system (2). The following proposition gives conditions for a positive equilibrium to be stable. Moreover, it states that conditions for local stability and global stability are identical, i.e. a locally stable equilibrium attracts all trajectories starting with positive initial conditions.

Proposition B.2

Assume that $\det (M)\ne 0$ and $M^{-1}\textbf{1}>0$, and denote the unique positive equilibrium of (2) by $Z^*$, whose expression is given in (6).

• i) $Z^*$ is locally stable if and only if M has 1 positive eigenvalue and $k-1$ negative eigenvalues.

• ii) $Z^*$ is locally stable if and only if it is globally stable on $(\mathbb {R}^*_+)^k$.

Our proof relies on a study of the Jacobian matrix of the system computed at the positive equilibria and on the construction of a Lyapunov function. Actually, we prove that the criterion on the eigenvalues of M is equivalent to the fact that the Jacobian matrix of the system admits only negative eigenvalues.

The second point of the proposition will be proved using the following function

$$\begin{aligned} V({\textbf {z}}):=\sum _{\ell =1}^k \left( \frac{z_\ell }{z}-\frac{z_{\ell }^*}{z^*}\ln \left( \frac{z_{\ell }}{z}\right) \right) =1-\sum _{\ell =1}^k \frac{z_{\ell }^*}{z^*}\ln \left( \frac{z_{\ell }}{z}\right) , \end{aligned}$$

(14)

which is a Lyapunov function for the dynamical system (2).

Proof

$\underline{\mathrm{Proof\,of}\,i):}$ Let us first give the Jacobian matrix of the system (2) around equilibrium $Z^*$. For all $i\in \{1,..,k\}$, let

$$\begin{aligned} F_i(z)=z_i\left( b\left[ 1+\frac{\sum _{j=1}^ks_{ij}z_{j}}{z}\right] -d-cz\right) , \text { with } z=\sum _{i=1}^k z_{i}. \end{aligned}$$

Recall that $F_i(z^*)/z^*_i=0$. Then by differentiating the function $F_i$ w.r.t $z_i$ and $z_j$ for $j\ne i$, we find

$$\begin{aligned}&\frac{\partial F_i}{\partial z_i}(z^*)=z^*_i\left( b\frac{s_{ii}}{z^*}-b\frac{\sum _{\ell =1}^k s_{i\ell }z^*_\ell }{(z^*)^2}-c\right) =\frac{z^*_i}{z^*}(b-d-2cz^*+bs_{ii})\\&\frac{\partial F_i}{\partial z_j}(z^*)=z^*_i\left( b\frac{s_{ij}}{z^*}-b\frac{\sum _{\ell =1}^k s_{i\ell }z^*_\ell }{(z^*)^2}-c\right) =\frac{z^*_i}{z^*}(b-d-2cz^*+bs_{ij}). \end{aligned}$$

Hence the jacobian Matrix of the system can be written

$$\begin{aligned}{} & {} J^*:=D^2((b-d-2cz^*)I+bM), \text { with } D:=\begin{pmatrix} \sqrt{\frac{z^*_1}{z^*}} &{} 0 &{}\ldots &{} 0 \\ 0 &{} \sqrt{\frac{z^*_2}{z^*}} &{}\ldots &{} 0 \\ \ldots &{}\ldots &{}\ldots &{}\ldots \\ 0 &{}\ldots &{} 0 &{} \sqrt{\frac{z^*_k}{z^*}} \end{pmatrix}\text { and }\nonumber \\{} & {} I:= \begin{pmatrix} 1 &{}\quad 1 &{}\quad .. &{}\quad 1 \\ 1 &{}\quad 1 &{}\quad .. &{}\quad 1 \\ .. &{}\quad .. &{}\quad .. &{}\quad .. \\ 1 &{}\quad 1 &{}\quad .. &{}\quad 1 \end{pmatrix}. \end{aligned}$$

(15)

Then since for any $(A,B) \in \mathcal {M}_k(\mathbb {R})$, $\det (AB)=\det (BA)$, we have

$$\begin{aligned} \det (J^*-\lambda Id)&=\det (D^2((b-d-2cz^*)I+bM)-\lambda Id)\\&=\det (D((b-d-2cz^*)I+bM)D-\lambda Id). \end{aligned}$$

In other words, the eigenvalues of $J^*$ are the same as the ones of

$$\begin{aligned} B:=D((b-d-2cz^*)I+bM)D \end{aligned}$$

(16)

which is symmetric and has real eigenvalues. According to Proposition 2.1,

$$\begin{aligned} b-d-2cz^*=-(b-d)-2b/(\textbf{1}^TM^{-1}\textbf{1})<0. \end{aligned}$$

Hence, we can write the matrix B as

$$\begin{aligned} B=bDMD-b\left( \frac{2}{\textbf{1}^TM^{-1}\textbf{1}}+\frac{b-d}{b}\right) DID. \end{aligned}$$

Let us prove that B has negative eigenvalues if M has only 1 positive eigenvalue. To this aim, notice that as M and DMD are congruent, according to Sylvester’s law of inertia (Sylvester 1852), they have the same numbers of positive/negative eigenvalues.

Let us first assume that M as well as DMD have only 1 positive eigenvalue. As DMD is symmetric, there exists an orthonormal basis $(V_i)_{i=1,..,k}$ of $\mathbb {R}^k$ formed with eigenvectors of DMD, associated to its eigenvalues $(\lambda _i)_{i=1,..,k}$. In the following we denote by $\lambda _1>0$ the positive eigenvalue. Using the definition of D in (15), we obtain that

$$\begin{aligned} D^2\textbf{1}= \frac{Z^*}{z^*} = \left( \frac{1}{\textbf{1}^TM^{-1}\textbf{1}} M^{-1}\textbf{1}\right) \end{aligned}$$

combining (13) and (6). Then

$$\begin{aligned} DMD(D\textbf{1})=DM(D^2\textbf{1})=DM\left( \frac{1}{\textbf{1}^TM^{-1}\textbf{1}} M^{-1}\textbf{1}\right) =\frac{1}{\textbf{1}^TM^{-1}\textbf{1}} D\textbf{1}. \end{aligned}$$

(17)

Since by assumption $\frac{1}{\textbf{1}^TM^{-1}\textbf{1}} >0$, we deduce that it is the unique positive eigenvalue of DMD and that $D\textbf{1}$ is a positive eigenvector associated to the positive eigenvalue of DMD.

From (17), we have

$$\begin{aligned} V_1=D\textbf{1}\text { and } \lambda _1=\frac{1}{\textbf{1}^TM^{-1}\textbf{1}}. \end{aligned}$$

Now, let us consider any $Y\in \mathbb {R}^k$ and its decomposition in the eigenvectors’ basis: there exist $(\gamma _1,..,\gamma _k)\in \mathbb {R}^k$ such that $Y=\sum _{i=1}^k\gamma _i V_i,$ and

$$\begin{aligned} Y^TBY&=bY^T(DMD)Y-b\left( \frac{2}{\textbf{1}^TM^{-1}\textbf{1}}+\frac{b-d}{b}\right) Y^T(DID)Y \nonumber \\&=b\sum _{i=1}^k\lambda _i\gamma _i^2-b\left( \frac{2}{\textbf{1}^TM^{-1}\textbf{1}}+\frac{b-d}{b}\right) Y^T(D\textbf{1})(D\textbf{1})^TY\nonumber \\&=\gamma _1^2\left( b\lambda _1-b\left( \frac{2}{\textbf{1}^TM^{-1}\textbf{1}}+\frac{b-d}{b}\right) \right) +b\sum _{i=2}^k\gamma _i^2\lambda _i\nonumber \\&=-b\gamma _1^2\left( \frac{1}{\textbf{1}^TM^{-1}\textbf{1}}+\frac{b-d}{b}\right) +b\sum _{i=2}^k\gamma _i^2\lambda _i< 0, \end{aligned}$$

(18)

since the eigenvalues $\lambda _2,.., \lambda _k$ of DMD are negative by assumption. In other words, B is a negative definite matrix. Recalling that B and the Jacobian matrix $J^*$ have the same eigenvalues, we deduce that the equilibrium is locally stable.

Conversely, if M has more than 1 positive eigenvalue, it is also the case of DMD. Then, assuming for example that $\lambda _2>0$, we have, by an application of (18) with $Y=V_2$,

$$\begin{aligned} V_2^TBV_2=b\lambda _2>0. \end{aligned}$$

Hence B has at least 1 positive eigenvalue. The stationary state is unstable.

$\underline{\mathrm{Proof\,of}\,ii): }$ We are now ready to prove the second point of the proposition. We only have to prove the direct implication, as the other one is obvious. To this aim, let us assume that the equilibrium is locally stable and thus all eigenvalues of $J^*$ are negative. Recall the definition of the function V in (14). By differentiating, we find

$$\begin{aligned} \frac{d}{dt} V(z(t))&= \sum _{\ell } \frac{\dot{z}_\ell }{z}-\frac{z^*_\ell }{z^*}\frac{\dot{z}_\ell }{z_\ell }-\frac{z_\ell }{z^2}\dot{z} + \frac{z^*_\ell }{z^*} \frac{\dot{z}}{z}\\&=\sum _{\ell } \left( \frac{z_\ell }{z}-\frac{z_{\ell }^*}{z^*}\right) \left( b\left[ 1+\frac{\sum _{j=1}^k s_{\ell j}z_{j}}{z}\right] -d-cz\right) \\&\quad -\frac{\dot{z}}{z} \sum _{\ell } \left( \frac{z_\ell }{z}-\frac{z_{\ell }^*}{z^*}\right) \\&= b\sum _{\ell } \left( \frac{z_\ell }{z}-\frac{z_{\ell }^*}{z^*}\right) \sum _{j=1}^k s_{\ell j}\frac{z_{j}}{z}, \end{aligned}$$

since $\sum _\ell z_\ell /z=\sum _{\ell }z^*_\ell /z^*=1$. In addition with the fact that $\sum _{j=1}^k s_{\ell j}\frac{z_{j}^*}{z^*}$ is a constant, we find

$$\begin{aligned} \frac{d}{dt} V(z(t))&= b\sum _{\ell } \left( \frac{z_\ell }{z}-\frac{z_{\ell }^*}{z^*}\right) \sum _{j=1}^k s_{\ell j}\left( \frac{z_{j}}{z}-\frac{z_{j}^*}{z^*}\right) \\&=b Y^T M Y, \end{aligned}$$

with

$$\begin{aligned} Y=\begin{pmatrix} \frac{z_{1}}{z}-\frac{z_{1}^*}{z^*}\\ \frac{z_{2}}{z}-\frac{z_{2}^*}{z^*} \\ .. \\ \frac{z_{k}}{z}-\frac{z_{k}^*}{z^*} \end{pmatrix}. \end{aligned}$$

Now, we recall that the stationary state is assumed to be locally stable, i.e. all eigenvalues of $J^*$ (or B) are negative. Since B and $(b-d-2cz^*)I+bM$ are congruent, these matrices are both negative definite. Furthermore, for all $X\in \mathbb {R}^k\setminus \{0\}$, $X^T((b-d-2cz^*)I+bM)X<0$. Since $\textbf{1}^TY=0$, we find

$$\begin{aligned} \frac{d}{dt} V(z(t))=b Y^T M Y =Y^T((b-d-2cz^*)I+bM)Y<0, \end{aligned}$$

as soon as $Y\ne 0$, and $\frac{d}{dt} V(z(t))=0$ if and only if $Y=0$, i.e. ${z_i}/{z}={z_i^*}/{z^*}$ for all $1 \le i \le k$. Moreover, we can prove that the set $\{z\in \mathbb {R}^k, Y=0\}$ is a positive invariant set for the dynamical system (2). Indeed, for all $i\in \{1,..,k\}$,

$$\begin{aligned} \frac{d}{dt}\left( \frac{z_i}{z}\right) \big |_{Y=0}= & {} \frac{F_i(z)}{z}-\sum _{j=1}^k \frac{F_j(z)}{z^2}\big |_{Y=0}\\= & {} b\frac{z_i}{z}\left( \sum _{j=1}^k s_{ij} \frac{z_j}{z} -\sum _{(j,l)} s_{j\ell }\frac{z_j}{z}\frac{z_\ell }{z}\right) \big |_{Y=0}=0, \end{aligned}$$

where we used (12). Theorem 1 of LaSalle (1960) then ensures that any bounded positive solution to (2) converges toward this invariant set. Finally, we prove that the solutions of (2) remain bounded by studying the dynamics of the total population size z, which satisfies

$$\begin{aligned} \frac{d}{dt}z=z\left( b+b\sum _{(i,j)}s_{ij}\frac{z_i}{z}\frac{z_j}{z}-d-cz\right) . \end{aligned}$$

We easily obtain that

$$\begin{aligned} \frac{d}{dt}z\le z\left( b+b\sum _{(i,j)}s_{ij}-d-cz\right) . \end{aligned}$$

and thus classical results on the logistic equation ensure that the total population size remains bounded through time.

Lasalle’s Theorem then ensures that $\frac{z_i(t)}{z(t)}$ converges to $\frac{z_i^*}{z^*}$. To conclude the proof it remains to prove that z(t) converges to $z^*$. Let us write $\frac{z_i(t)}{z(t)}=\frac{z_i^*}{z^*}+\epsilon _i(t)$ and according to (12)

$$\begin{aligned} \frac{d}{dt}z=z(cz^*+b\mathcal {E}(t)-cz), \end{aligned}$$

with

$$\begin{aligned} \mathcal {E}(t)=\sum _i(cz^*-b+d)\epsilon _i(t)+\sum _{i,j}s_{ij}\epsilon _i(t)\epsilon _j(t)\underset{t\rightarrow \infty }{\rightarrow } 0. \end{aligned}$$

For any $\varepsilon >0$, as soon as t is large enough, we deduce that

$$\begin{aligned} z(cz^*-\varepsilon -cz)\le \frac{d}{dt}z\le z(cz^*+\varepsilon -cz). \end{aligned}$$

Using classical results on logistic equations, we deduce that z(t) converges to $z^*$. This ends the proof of Proposition 2.1. $\square$

1.3 Complete study of the three allele case

Consider three alleles $A_1$, $A_2$, and $A_3$ whose selection matrix is given by

$$\begin{aligned} M=\begin{pmatrix} 0 &{}\quad s_{12} &{}\quad s_{13} \\ s_{12} &{}\quad 0 &{}\quad s_{23} \\ s_{13} &{}\quad s_{23} &{}\quad 0 \end{pmatrix}\quad \text {with } s_{12}>0, s_{13}>0 \text { and } s_{23}>0. \end{aligned}$$

(19)

Since the selective advantage parameters are positive, $\det (M)>0$ and one obtains with a simple computation that the condition for existence of a three-allele positive equilibrium $M^{-1}\textbf{1}>0$ can be written as (8). Moreover, the three-alleles equilibrium equals

$$\begin{aligned} \begin{pmatrix} z^*_1 \\ z^*_2 \\ z^*_3 \end{pmatrix}=\frac{z^*}{\textbf{1}^T M^{-1} \textbf{1}} M^{-1}\textbf{1}= \frac{z^*}{\textbf{1}^T M^{-1} \textbf{1}} \frac{1}{2s_{12}s_{13}s_{23}} \begin{pmatrix} s_{23}(s_{12} + s_{13} - s_{23}) \\ s_{13}(s_{12} + s_{23} - s_{13}) \\ s_{12}(s_{13} + s_{23} - s_{12}) \end{pmatrix}, \end{aligned}$$

where

$$\begin{aligned}{} & {} z^*=\frac{1}{c}\left( b-d+\frac{b}{\textbf{1}^T M^{-1} \textbf{1}}\right) \quad \text {and} \textbf{1}^T M^{-1} \textbf{1}=\frac{1}{2s_{12}s_{13}s_{23}} \sum _{(i,j,k) \in \{1,2,3\}^3}s_{ij}(s_{ik}+s_{jk}-s_{ij}). \end{aligned}$$

In this case, the behavior of the dynamical system (2) can be fully characterized.

Proposition B.3

Let us consider a system of three alleles, whose interactions are characterized by matrix (19).

• Assume (8), then the co-existence equilibrium given above is globally stable on $(\mathbb {R}^*_+)^3$.

• otherwise, assume for example that $s_{23}\ge s_{12}+s_{13}$, then the two-alleles equilibrium

  $$\begin{aligned} \left( 0,\frac{b-d+b{s_{23}}/{2}}{2c},\frac{b-d+b{s_{23}}/{2}}{2c}\right) , \end{aligned}$$

  is globally stable on $(\mathbb {R}^*_+)^3$.

Proof

The first point easily derives from Proposition 2.1.

For the second point, we exhibit a Lyapunov function, which is similar to the one of the previous proof. Let us consider

$$\begin{aligned}{} & {} V({\textbf {z}})=\sum _{\ell =1}^3 \left( \frac{z_\ell }{z}-\frac{z_{\ell }^*}{z^*}\ln \left( \frac{z_{\ell }}{z}\right) \right) , \quad \text {with} {\textbf {z}}^*=\left( 0,\frac{b-d+b{s_{23}}/{2}}{2c},\frac{b-d+b{s_{23}}/{2}}{2c}\right) . \end{aligned}$$

Then, by a direct computation, we find

$$\begin{aligned} \begin{aligned} \dot{V}({\textbf {z}})&=-\frac{b}{2z^2}\left[ s_{12}z_1z+s_{23}z_3z+s_{13}z_1z+s_{23}z_2z-4s_{12}z_1z_2-4s_{13}z_1z_3-4s_{23}z_2z_3\right] \\&=-\frac{b}{2z^2}\left[ (s_{12}+s_{13})z_1^2+s_{23}(z_2-z_3)^2+z_1(z_2+z_3)(s_{23}-s_{12}-s_{13})\right. \\&\qquad \qquad \left. +2z_1(-s_{13}z_3+s_{12}z_3-s_{12}z_2+s_{13}z_2) \right] . \end{aligned} \end{aligned}$$

Recall that $s_{23}\ge s_{12}+s_{13}$, and without loss of generality, we may assume that $s_{12} \le s_{13}$ (otherwise, exchange the roles of $z_2$ and $z_3$ in the computations). Then

$$\begin{aligned} \begin{aligned} \dot{V}({\textbf {z}})&\le -\frac{b}{2z^2}\left[ (s_{13}-s_{12})(z_1+z_2-z_3)^2+2s_{12}z_1^2+(s_{23}-s_{13}+s_{12})(z_2-z_3)^2 \right] \\&\le 0. \end{aligned} \end{aligned}$$

Therefore V is a Lyapunov function and $\dot{V}({\textbf {z}})=0$ if and only if $z_1=0$ and $z_2=z_3$. We then conclude with Lasalle Theorem (Theorem 1 of LaSalle (1960)) and an argument similar to the one at the end of the previous proof. $\square$

1.4 Successive introductions of types

Assume that k types coexist. In view of previous sections, if $\det (M)\ne 0$, there is coexistence of the k types if and only if $M^{-1}\textbf{1}>0$ and the second eigenvalue of M is negative (or equivalently the eigenvalues of $J^*$ are all negative).

If a mutant characterized by parameters $S=(s_{k+1,i})_{i=1,..,k}$ and $\sigma =s_{k+1,k+1}$ appears, this mutant may invade if the Jacobian matrix of (2) around equilibrium $(z^*,0)$ is unstable. This Jacobian matrix can be written

$$\begin{aligned} J_{(z^*,0)}=\begin{pmatrix} J^* &{} \star \\ 0 \,.. \, 0 &{} \frac{\partial F_{k+1}}{\partial z_{k+1}}(z^*,0) \end{pmatrix}. \end{aligned}$$

As previously written, the matrix $J^*$ is the Jacobian matrix of the resident system of size k at equilibrium $z^*$ and its eigenvalues are all negative, thus the equilibrium $(z^*,0)$ is unstable if the last eigenvalue is positive, i.e.

$$\begin{aligned} \frac{\partial F_{k+1}}{\partial z_{k+1}}(z^*,0)=b-d-cz^*+\frac{b}{z^*}\sum _{i=1}^{k}s_{i,k+1}z^*_i=-\frac{b}{\textbf{1}^T M^{-1}\textbf{1}}(1-S^T M^{-1}\textbf{1})>0,\nonumber \\ \end{aligned}$$

(20)

which gives (9).

In what follows, we will be interested in proving the following lemma which states that if the mutant may invade, the new community will be composed of $k+1$ types if and only if the coexistence state with $k+1$ types exists.

Proposition B.4

Let us consider a stable resident population with k types given by a matrix M, i.e M satisfies $M^{-1}\textbf{1}>0$ and its second eigenvalue is negative.

Let us consider a mutant type arising in this resident population characterized by S and $\sigma$. Denote the new fitness matrix by

$$\begin{aligned} \bar{M}=\begin{pmatrix} M&{}S\\ S^T &{}\sigma \end{pmatrix}. \end{aligned}$$

If $\bar{M}^{-1}\textbf{1}>0$, i.e. the equilibrium with $k+1$ species exists, then it is globally asymptotically stable.

Proof

Let us first give conditions under which the second eigenvalue of $\bar{M}$ is positive. As $\bar{M}$ is a symmetric matrix and M a principal sub-matrix of $\bar{M}$, Proposition D.1 implies that

$$\begin{aligned} \lambda _1^{\bar{M}}\ge \lambda _1^{M}\ge \lambda _2^{\bar{M}}\ge \lambda _2^M\ge \cdots \ge \lambda _k^M\ge \lambda _{k+1}^{\bar{M}}, \end{aligned}$$

(21)

where $(\lambda _i^M)_{i=1,\ldots ,k}$ denote the eigenvalues of M and $(\lambda _i^{\bar{M}})_{i=1,\ldots ,k+1}$ the ones of $\bar{M}$. Since $\lambda _2^M<0$, $\bar{M}$ can have at most 2 positive eigenvalues. On the other hand, using Proposition D.2, we compute

$$\begin{aligned} \det ({\bar{M}})=(\sigma -S^TM^{-1}S)\det (M). \end{aligned}$$

Since $\det (M)$ has the sign of $(-1)^{k-1}$, we deduce the following equivalences:

$$\begin{aligned} \lambda _2^{\bar{M}}<0 \quad \Longleftrightarrow \quad \det (\bar{M}) \text { has the sign of } (-1)^k\quad \Longleftrightarrow \quad \sigma -S^TM^{-1}S<0.\qquad \end{aligned}$$

(22)

Let us now prove that this last condition holds as soon as $\bar{M}^{-1}\textbf{1}>0$. For all $i\in \{1,..,k+1\}$ we denote by $M_i$ the matrix resulting from the suppresion of the i$^{th}$ column and the i$^{th}$ line of $\bar{M}$, by $S_i$ the i$^{th}$ column of $\bar{M}$ without the i$^{th}$ coefficient, and by $\sigma _i$ the coefficient $\bar{M}_{ii}$. Notice that $M_{k+1}=M$, $S_{k+1}=S$ and $\sigma _{k+1}=\sigma$. Using these notations and Proposition D.2, we notice that for all i

$$\begin{aligned} \det (\bar{M})=(\sigma _i-S_i^TM_i^{-1}S_i)\det (M_i), \end{aligned}$$

(23)

Let us now denote $\bar{M}^{-1}\textbf{1}$ by $(u_i)_{i=1,..,k+1}$. Then using Cramer’s rule along with (23), we have for all $i\in \{1,..,k+1\}$

$$\begin{aligned} u_i=\frac{1-S_i^TM_i^{-1}\textbf{1}}{\sigma _i-S_i^TM_i^{-1}S_i}. \end{aligned}$$

(24)

Hence, we deduce that if $\bar{M}^{-1}\textbf{1}>0$, then $u_{k+1}>0$, and since $1-S^TM^{-1}\textbf{1}<0$, then $\sigma -S^TM^{-1}S<0$. Therefore according to (22) and Proposition 2.1, the equilibrium with $k+1$ species is globally stable. This ends the proof. $\square$

Now, let us study the case where the equilibrium with $k+1$ species does not exist, i.e. there exists $i\in \{1,..,k+1\}$ such that $u_i<0$. We are not capable of deciphering which alleles will go extinct due to mutant invasion. Numerical simulations reveal that alleles $1\le i\le k$ that verify $u_i u_{k+1}>0$ might go extinct in the equilibrium. Notice that in some cases, all these alleles disappear, while in other cases, only a subset does.

Proofs concerning the genetic model

This section gathers theoretical results associated with the specific model in Sect. 3.1.

1.1 Globally stable equilibrium in a particular case

In this section we focus on the case where $L=3$, such that there are 8 possible alleles (see Fig. 6 for a representation).

We first prove the convergence of the population state in a particular case for which there is no positive stable equilibrium. This example will be important for the proof of Proposition 3.1 below.

Assume that 4 alleles are present in the population:

$$\begin{aligned} (0,1,1) \quad (1,0,1) \quad (0,0,0)\quad (1,1,1) \end{aligned}$$

up to permutation of the sites. Then, the selection matrix M equals

$$\begin{aligned} M=\begin{pmatrix} 0 &{}\quad 2^{\alpha } &{}\quad 2^{\alpha } &{}\quad 1 \\ 2^{\alpha } &{}\quad 0 &{}\quad 2^{\alpha } &{}\quad 1 \\ 2^{\alpha } &{}\quad 2^{\alpha } &{}\quad 0 &{}\quad 3^{\alpha } \\ 1 &{}\quad 1 &{}\quad 3^{\alpha } &{}\quad 0\end{pmatrix}. \end{aligned}$$

Our aim is to prove that the equilibrium

$$\begin{aligned} \left( 0,0,\frac{b-d+3^{\alpha }b}{2c},\frac{b-d+3^{\alpha }b}{2c}\right) \end{aligned}$$

(25)

is globally positively stable in $(\mathbb {R}^+)^4$, i.e. all trajectories starting from $(\mathbb {R}^+)^4$ converge to this equilibrium.

To this aim, let us compute the function V defined in Eq. (14) with $z^*$ defined by (25), i.e. $V(z)=1-\frac{1}{2}\ln \left( \frac{z_3}{z}\right) -\frac{1}{2}\ln \left( \frac{z_4}{z}\right)$. We get

$$\begin{aligned} \dot{V}(z_1,z_2,z_3,z_4)&=-\frac{b}{2z^2}\left[ 2^{\alpha }(z_1-z_2)^2+(z_1+z_2)z_3(1+3^{\alpha }-2^{\alpha +1})\right. \nonumber \\&\quad +(z_1+z_2)z_4(3^{\alpha }-3)+\left. (z_1+z_2)^2+3^{\alpha }(z_3-z_4)^2+2^{\alpha }(z_1+z_2)(z_4-z_3)\right] \nonumber \\&\le - \frac{b}{2z^2}\left[ (z_1+z_2)^2+3^{\alpha }(z_3-z_4)^2+2^{\alpha }(z_1+z_2)(z_4-z_3)\right] . \end{aligned}$$

(26)

Finally, note that if $\alpha \le 2$ then $\alpha \le \frac{\ln (4)}{\ln (4/3)}$, and then the two dimensional function $(x,y)\mapsto x^2+3^{\alpha }y^2+2^{\alpha }xy$ is positive on $\mathbb {R}^2$. So V is a Lyapunov function and its derivative is null only if $z_1=z_2=0$ and $z_3=z_4$. We then conclude by using an argument similar to the one of the end of proof of Proposition B.2. This finally proves the convergence of all positive trajectories to Equilibrium (25).

1.2 Proof of Proposition 3.1

Proposition 3.1 arises from two facts: we are able to exhibit stable sets of alleles (see Fig. 6) and we can characterize the impact of mutations on the number of coexisting alleles in the population.

In fact, we obtain a stochastic process which jumps across stable states of the population depending on the arrival of mutant alleles. This process is characterized by a transition matrix which explicits the possible transitions and is represented on Fig. 8.

Fig. 8

Transition matrix for the allelic composition of the population. The red color corresponds to stable population states. Blue dashed arrows show the initial distribution after the first mutation

Proof

The set of possible alleles is $\{0,1\}^3$. We first prove that the two following types of genetic compositions are stable: (Stable population 1) Only the two most-differentiated alleles, (Stable population 2) four alleles that are all at distance 2 from each other.

We already know that (stable population 1) is stable from the end of Sect. 3.1. For (stable population 2), we study the conditions for invasion of a mutant in a population composed of four alleles all at distance 2 from each other. The new mutant will necessarily be at distance 1 from 3 of the resident types, and at distance 3 from the forth one. Then following the notations of Eq. (9) in Sect. 2.3, we get that $S=(1, 1, 1, 3^{\alpha })^T$. Therefore the new mutant can invade only if

$$\begin{aligned} \frac{1-S^TM^{-1}\textbf{1}}{\textbf{1}^TM^{-1}\textbf{1}} =3\times 2^{\alpha -2} \left( 1-\frac{1+3^{\alpha -1}}{2^{\alpha }}\right) <0. \end{aligned}$$

In particular for $\alpha \in [1,2]$, this gives that the mutant cannot invade, hence the population with four equidistant alleles is stable.

Now, we consider the different possible orders of mutations arising in the population to prove the proposition. Let us assume without loss of generality that the population starts with only individuals with allele (0, 0, 0). The second arising allele in the population can be at distance 1 (with probability 3/7), 2 (with probability (3/7), or 3 (with probability 1/7), from (0, 0, 0). In any case, Proposition B.3 implies that this mutant will invade and the population will then exhibit two different alleles at equilibrium.

• If this first mutant is at distance 3 from the resident allele (0, 0, 0), then the population has reached the configuration (Stable population 1) and no new mutation will be able to invade.

• If the first mutant is at distance 2 from the resident allele (0, 0, 0), then there are 3 possibilities for the second mutant (each arising with probability 1/3):

  • (i) The second mutant can be at distance 1 from each of the two pre-existing alleles. Then Eq. (9) shows that it cannot invade, if $\alpha >1$. The population thus stays with the two resident alleles, at a distance 2 from each other.

  • (ii) it can be at distances 1 and 3 from the two first alleles, in which case only the mutant and its opposite type will remain in the population (according to Proposition B.3 and since $3^\alpha >2^\alpha +1$), which corresponds to a configuration (Stable population 1).

  • (iii) The second mutant can be at distance 2 from the two pre-existing alleles. From Eq. (8), the population will then converge to an equilibrium with the three first alleles.

    If the third mutant is in the fourth position at distance 2 from the three pre-existing alleles then it will also invade (see Sect. 2.3). This mutant arises with probability 1/5 and brings the population to (Stable population 2) (as illustrated in Fig. 6).

    If the third mutant is in a position such that two most-differentiated alleles co-exist in the population then only these two alleles will remain (see Sect. C.1 for details). This occurs with probability 3/5.

    Lastly, if (with probability 1/5) the third mutant is at distance 1 from the three first types then it cannot invade, according to Eq. (9).

• Finally, if the first mutant is at distance 1 from (0, 0, 0), then there are two possibilities for the second mutant:

  • (i) The second mutant is at distance 2 and 3 from the two resident alleles, respectively (with probability 1/3). Only the two most-differentiated alleles will then remain from Proposition B.3 and this population state is stable.

  • (ii) The second mutant is at distances 1 and 2 (with probability 2/3) from the pre-existing alleles, respectively. Thus, one of the two pre-existing alleles will disappear and the population will exhibit only two alleles, at distance 2 from each other, as already proved by Proposition B.3.

Finally, the transition graph of the genetic composition of the population is given in Fig. 8. This gives that the probability that the population ends up with 4 equidistant genotypes is equal to

$$\begin{aligned} \frac{3}{7}\times \frac{\frac{1}{3}}{\frac{2}{3}}\times \frac{\frac{1}{5}}{\frac{4}{5}}+\frac{3}{7}\times \frac{2}{3}\times \frac{\frac{1}{3}}{\frac{2}{3}}\times \frac{\frac{1}{5}}{\frac{4}{5}}=\frac{5}{56}. \end{aligned}$$

$\square$

Complements of linear algebra

Proposition D.1

(Eigenvalue Interlacing Theorem). Let $A\in \mathbb {R}^{n\times n}$ be a symmetric matrix and B be a principal sub-matrix of size $m<n$. If $\lambda ^A_1\ge \lambda ^A_2 \ge \cdots \ge \lambda ^A_n$ the eigenvalues of A and $\lambda ^B_1\ge \cdots \ge \lambda ^B_m$ the eigenvalues of B. Then for any $1\le k\le m$

$$\begin{aligned} \lambda ^A_k\ge \lambda ^B_k\ge \lambda ^A_{k+n-m}. \end{aligned}$$

In particular if $m=n-1$ then

$$\begin{aligned} \lambda ^A_1\ge \lambda ^B_1\ge \lambda ^A_2\ge \lambda ^B_2\ge \cdots \ge \lambda ^A_{n-1}\ge \lambda ^B_{n-1}\ge \lambda ^A_n. \end{aligned}$$

Proposition D.2

(Schur complement). Let M be a matrix defined by blocs as $M=\begin{pmatrix} A&{}B\\ C&{}D\end{pmatrix}$ If D is invertible then the Schur complement of D is defined by $A-BD^{-1}C$. Moreover

$$\begin{aligned} \det (M)=\det (A)det( A-BD^{-1}C) \end{aligned}$$

For proof of these two results, see Example 7.5.3 and Sect. 6.2 in Meyer (2000).