Selection-recombination-mutation dynamics: Gradient, limit cycle, and closed invariant curve

How a living being appears and behaves is inherently determined by the variety of genes located on the loci of chromosomes. The entire genotype to phenotype map is very vast and complex, and is further complicated by the environmental feedback effects. While the evolutionary dynamics—e.g., replicator dynamics Taylor and Jonker (1978); Schuster and Sigmund (1983); Cressman and Tao (2014); Mukhopadhyay and Chakraborty (2021)—of simplified one-locus many-allele systems exhibits quite rich dynamical behaviour, the phenomenon of recombination during meiosis presents one with richer extension of the evolutionary dynamics in a population of sexually mating organisms. In the latter, the simplest non-trivial mathematical setup is that of a two-locus two-allele system. While two different loci can determine two different phenotypic traits (e.g., seed colour and seed colour in the revolutionary dihybrid cross experiments of Mendel O’Brien (1993); Mendel and and (1866)), there are plethora of examples where two loci control the same trait and their interaction has considerable effect on the phenotype. Some such examples that quickly comes to mind are comb-types on the head of chicken Bateson et al. (1909), flower colour in peas Dooner et al. (1991), wheat kernel color Carver (2009), blood group in humans Dean and Dean (2005), age-related hearing loss resistance in the Japanese wild-derived inbred MSM/Ms mice Yasuda et al. (2022), and eye colour of humans White and Rabago-Smith (2010).

Haldane Haldane (1931) and Wright Wright (1935) were among the firsts to mathematically explore the evolutionary outcome when selection acts on more than one locus assumed statistically independent, i.e., system is in linkage equilibrium. Later Lewontin and Kojima Lewontin and Kojima (1960) worked on general 2L2A model with both selection and recombination included, and showed that strong epistasis together with linkage disequilibrium can lead to significantly different outcomes. Subsequently, many extensions Bodmer and Felsenstein (1967); Karlin et al. (1970); Feldman and Libermann (1979); Hastings (1985); Ewens (1968) of the model were studied—starting from prediction and calculation of all possible fixed points and determination of stability condition for some special cases in simplified fitness model to the occurrence of cyclic motion. For continuous time model, Akin Akin (1982) and for discrete time model, Hastings Hastings (1981), Hofbauer and Iooss Hofbauer and Iooss (1984) showed that periodic orbits can be a possible outcome when system is in linkage disequilibrium; when the 2L2A model is in linkage equilibrium such complex behaviour can not occur and the system become gradient-like as shown by Nagylaki Nagylaki (1989). Nagylaki also showed that for multi-locus system under weak selection—selection strength sufficiently smaller compared to recombination—linkage disequilibrium decays close to zero within a few generations and the dynamics of the entire multilocus system is governed by the dynamical outcomes of time-continuous multi-locus system Nagylaki (1993); Nagylaki et al. (1999a); Nagylaki (2013a). Very recently Pontz et al. (2018), Pontz and coworkers studied the 2L2A model in detail under weak selection limit and classified all possible equilibrium structure and phase portraits for different payoff structures.

Mutation is the most important and indispensable ingredient of the evolutionary processes. It is important to note that the relationship between mutations at different loci is often complex and can be influenced by a variety of factors. It is possible for a mutation at one locus to affect the mutation rate at another locus on the same chromosome Nowak et al. (2004); Li et al. (2013) directly or indirectly. For example, a mutation in a gene that codes for a repair enzyme involved in maintaining the integrity of DNA could alter the mutation rate of other genes by affecting the efficiency of DNA repair Steele et al. (1998). Similarly, a mutation in a gene that regulates the transcription or translation of other genes can affect the rate of mutations in those genes. Mutation may not be always completely random: It can depend on environmental changes Galhardo et al. (2007) and fitness Baer (2008), and organisms can have evolved mechanisms which can influence the timing or genomic location of mutation Rando and Verstrepen (2007). In view of the above, in this paper, we distinguish between the cases of mutations that occur randomly (independent of what mutations occur at the other loci) from the cases where mutations can be treated to have occurred independently and randomly at a locus—we call the latter independent mutations.

The effects of mutation was studied in one-locus two-allele model for both overlapping and non-overlapping generations, and the results are really interesting: Mutation can stop the extinction of cooperators by producing cyclic or chaotic outcomes and even can produce coexistence regions where both the cooperators and defectors can survive simultaneously Mobilia (2010); Toupo and Strogatz (2015); Mittal et al. (2020); Mukhopadhyay et al. (2021). Evolutionary dynamics were studied for the case of the 2L2A with mutation as well Park and Krug (2011); Qu et al. (2020); Bengtsson and Christiansen (1983); Jain and Seetharaman (2021); Bürger (1989); Christiansen and Frydenberg (1977); Karlin and McGregor (1971). Effect of weak selection in 2L2A model, in the absence of mutation, was studied in detail by Pontz et al. Pontz et al. (2018). If the selection is sufficiently weak then the linkage disequilibrium disappears after a few generations and dynamics become quite straightforward to find its equilibrium structure for the general fitness matrix. The system behaves like a gradient system and if all the fixed point are hyperbolic then for any initial condition dynamics ultimately approaches fixed points: No periodic or chaotic outcome can arise. On the other hand, the effect of both weak mutation and weak selection was studied for one-locus many-allele model by Hofbauer and Sigmund Hofbauer and Sigmund (1998). They have shown that the corresponding system becomes gradient system when the mutation matrix satisfies certain condition: When mutation probability dependents on the target gene only Hofbauer (1985), the system behaves like a gradient system, and no periodic or chaotic orbit is possible. In general, however, the system need not be a gradient system, and hence periodic and chaotic orbits can be possible outcome.

In such contexts, the study of the 2L2A model in presence of weak selection and weak mutation is still, to the best of our knowledge, remains to be reported. Specifically, we are interested in whether the 2L2A model is always gradient-like or there is a specific condition on mutation for that to happen. Does the answer depend on whether the mutations at the two locus are independent? Furthermore, if and when the epistasis induced dynamics in 2L2A system be altered due to presence of mutation is also a question of interest because epistasis or interaction between different loci has important effects on dynamics, like occurrence of cyclic motion Hastings (1981) and sustaining polymorphism Pontz et al. (2018).

Without further ado, we introduce the model in the next section (Sec. II) and then we present the gradient 2L2A system in Sec. III. Subsequently, we discuss the non-gradient dynamics in Sec. III before concluding in Sec. IV.

We consider a generation-wise non-overlapping, unstructured population of randomly mating diploid individuals. The population is assumed to evolve under the action of viability selection that acts on two diallelic, recombining loci. We, thus, have the standard two-locus two-allele (2L2A) model with viability selection Pontz et al. (2018); Karlin (1975); Karlin and Carmelli (1975); Sacker and von Bremen (2003); Hallgrímsdóttir and Yuster (2008); Lewontin and Feldman (1988); Franklin and Feldman (1977); Bodmer and Felsenstein (1967). Let $A_{1}$ and $A_{2}$ be the alleles at locus $A$, and $B_{1}$ and $B_{2}$ be the alleles at locus $B$. Let $G_{1}$, $G_{2}$, $G_{3}$, and $G_{4}$ represent the four possible gametes, viz., $A_{1}B_{1}$, $A_{1}B_{2}$, $A_{2}B_{1}$, and $A_{2}B_{2}$ respectively. Let the frequency of the $G_{i}$ gametes be $x_{i}$ (obviously, $\sum_{i=1}^{4}x_{i}=1$). In these notations, the frequencies of allele $A_{1}$ and $B_{1}$ are respectively $x_{1}+x_{2}$ and $x_{1}+x_{3}$ that we henceforth denote as $p$ and $q$ respectively. Defining $D\equiv x_{1}x_{4}-x_{2}x_{3}$, the measure of linkage disequilibrium, following identities follow:

Furthermore, let the fitness of the $G_{i}G_{j}$ genotype be $w_{ij}$. For simplicity, we assume that the fitness is independent of the fact which gamete is from mother and which one is from father, i.e., $w_{ij}=w_{ji}$. We also assume that $w_{14}=w_{23}$, the last assumption means that the corresponding fitnesses are the same regardless of whether $A_{1}$ and $B_{1}$ are on the same chromosome or opposite ones. Consequently, the following matrix suffices for fully representing the fitnesses of genotypes:

We include two further phenomena in the model: recombination and mutation. While we let $r$ denote the probability of recombination, the mutation from one gamete to another gamete is specified by the row stochastic matrix ${\sf Q}$ whose ($i,j$)th element $Q_{ij}$ is the probability that an offspring with gamete $G_{j}$ is born to a parent with gamete $G_{i}$. With this multiplicative Toupo and Strogatz (2015); Hofbauer and Sigmund (1998) mutation that takes place during DNA replication process, the evolution of gamete frequencies is given by the following equation Bengtsson and Christiansen (1983); Pritchett-Ewing (1981a); Christiansen and Frydenberg (1977); Karlin and McGregor (1971):

where $w_{i}\equiv\sum_{j}w_{ij}x_{j}$, $\bar{w}\equiv\sum_{j}w_{j}x_{j}$, $\sum_{j}Q_{ij}=1$, and $\theta_{1}=-\theta_{2}=-\theta_{3}=\theta_{4}=1$. Here, prime is the tag for immediately succeeding generation. We suppose that the mutation happens at the gametic stage Bengtsson and Christiansen (1983); Bergero et al. (2021).

As elaborated in the introduction to this paper, since our goal is to observe effect of both weak selection and weak mutation on the dynamics of the 2L2A selection-recombination model, we introduce a small parameter $s$ such that the fitnesses and mutations can be recast as

Thus, the limit of small $s$ renders the dynamics of Eq. (2) to be that under weak selection and weak mutation.

It is well-known Hofbauer and Sigmund (1998) that recombination drives 2L2A model to Wright manifold in the absence of selection and mutation: The Wright manifold is a linkage equilibrium manifold ($\Lambda_{0}$), where $D=0$, and a population in linkage equilibrium always remains in equilibrium, i.e. it is an invariant manifold. In the presence of weak selection $(s\ll r)$ and weak mutation, on using Eq. (3) in Eq. (2), we obtain

Interestingly, we find that the effect of mutation at this order is completely absent. Hence, the known results Nagylaki (2013b); Nagylaki et al. (1999b) for the case of the models without mutation (but with selection) holds in our case as well: Close to $\Lambda_{0}$, there exists a smooth invariant manifold, $\Lambda_{s}$, that is globally attracting for Eq. (2); for any initial values the linkage disequilibrium $D(t)$ becomes $\mathcal{O}(s)$ asymptotically in time.

Thus, in the weak selection and weak mutation limit, the linkage disequilibrium $D\to 0$, and consequently, the dynamical equation is further simplified because of the fact that it can now be specified using only two variables $p$ and $q$, for $t\geq t_{1}$. Specifically, in the set of Eqs. (1), $D$ is replaced by $\mathcal{O}(s)$ terms. Rescaling time $t(=0,1,2...)$ tagging generations as $\tau=st$, it is easy to see that as $s\to 0$, Eq. (2) approaches the following differential equations:

where, $\nu_{1}\equiv{(\epsilon_{31}+\epsilon_{32})},~{}\nu_{2}\equiv{(\epsilon_{41}+\epsilon_{42})},~{}\nu_{3}\equiv{(\epsilon_{13}+\epsilon_{14})},~{}\nu_{4}\equiv{(\epsilon_{23}+\epsilon_{24})},~{}\nu_{5}\equiv{(\epsilon_{21}+\epsilon_{23})},~{}\nu_{6}\equiv{(\epsilon_{41}+\epsilon_{43})},~{}\nu_{7}\equiv{(\epsilon_{12}+\epsilon_{14})},~{}\nu_{8}\equiv{(\epsilon_{32}+\epsilon_{34})}$, and

is the average fitness of the population wherein the marginal fitness of the gamete $i$,

Since, at $p=0$ and $p=1$, $\dot{p}\geq 0$ and $\dot{p}\leq 0$ respectively; and at $q=0$ and $q=1$, $\dot{q}\geq 0$ and $\dot{q}\leq 0$ respectively, Eq. (5) is forward invariant in the $p$-$q$ phase space—a closed unit square. It may also be noted that the dynamical equation remain invariant under addition (but not multiplication) of a constant matrix with the matrix ${\sf M}$ whose elements are $m_{ij}$’s. It may be remarked that had the mutation been taken as additive rather than multiplicative Mobilia (2010); Toupo and Strogatz (2015), the form of the final dynamical equation, in the limit of weak selection and weak mutation, remains same as Eqs. (5).

If a dynamical system is gradient-like, then the non-convergent solutions like periodic or chaotic solutions can be ruled out the system. So it is very useful to find the condition on the mutation rates for which the 2L2A model becomes gradient-like.

It is of notational convenience to represent Eqs. (5) in the form:

here, $z_{1}\equiv p$, $z_{2}\equiv q$, $z_{3}\equiv(1-p)$ and $z_{4}\equiv(1-q)$ correspond to the four coordinates in the allele frequency space. The explicit expression of $f_{i}$’s is obvious when compared with Eqs. (5). The phase space $\mathbb{P}$, hence, is ${\Sigma}^{2}\times\Sigma^{2}$ embedded in $\mathbb{R}^{4}$; here, $\Sigma^{2}$ is one dimensional simplex. We need to consider the interior of $\mathbb{P}$, ${\rm int}\mathbb{P}$ (say); at any point ${\bm{z}}\in{\rm int}\mathbb{P}$, let the tangent space be denoted by $\mathbb{T}_{\bm{z}}\mathbb{P}$.

Now for ${\bm{z}}\in{\rm int}\mathbb{P}$ and $\bm{\eta}\in\mathbb{T}_{\bm{z}}\mathbb{P}$, the inner product

where (and henceforth) sum over repeated indices—Einstein’s summation convention—has been imposed. $g_{ij}$ is the metric in ${\rm int}\mathbb{P}$. Considering the metric $g_{ij}={\delta_{ij}}/[z_{i}(1-z_{i})]$—reminiscent of the Shahshahani gradient Hofbauer and Sigmund (1998); Shahshahani (1979); Sigmund (1985)—and recalling Eq. (8), it is obvious that

if there exists a continuous and differentiable scalar function $U({\bm{z}})$ such that $f_{i}(\bm{z})={\partial U({\bm{z}})}/{\partial z_{i}}$. Actually, since $\sum_{i=1}^{4}\eta_{i}=0$ because $\bm{\eta}\in\mathbb{T}_{\bm{z}}\mathbb{P}$, even if a more general condition, viz., $f_{i}(\bm{z})={\partial U({\bm{z}})}/{\partial z_{i}}+\phi({\bm{z}})$ ($\phi$ is a scalar function) is satisfied, Eq. (10) holds good (see Appendix A).

In summary, what we have found is that if $\dot{z}_{i}={f}_{i}({\bm{z}})$ is gradient system, then so is $\dot{z}_{i}=\tilde{f}_{i}({\bm{z}})$. Consequently, for Eqs. (5) to be gradient system, one requires

This condition, on putting the expression of $f_{i}$’s, ultimately takes the form

As $p$ and $q$ are independent variables, the system is gradient system for all possible allowed values of $p$ and $q$ if

Evidently, the condition on mutation for 2L2A model to be gradient-like is completely different from the one-locus four-allele model Hofbauer and Sigmund (1998).

In the light of condition (13), Eqs. (5) take the form

It is interesting to note that with $V(p,q)\equiv\exp(\bar{m})p^{2\nu_{1}}(1-p)^{2\nu_{3}}q^{2\nu_{5}}(1-q)^{2\nu_{7}}$ , Eq. (14) can be recast as:

the interior fixed points $(p^{*},q^{*})$ correspond to ${\partial V}/{\partial p}={\partial V}/{\partial q}=0$. In fact, it can be verified (see Appendix B) that $V(p,q)$ is local Lyapunov function for the gradient system (14) because

implying $\dot{V}=0$ only at any $(p^{*},q^{*})\in(0,1)^{2}$ and $\dot{V}>0$ always for all possible $(p,q)\in(0,1)^{2}$ other than the fixed points.

The existence of the Lyapunov function leads to the conclusion for Eq. (2): For suffciently small $s$ and if all the equilibria of Eqs. (15) are hyperbolic, every solution of Eq. (2) converges to a fixed point. This is essentially as extension of a theorem due to Nagyalaki Nagylaki et al. (1999a) and can be succinctly understood as follows. Substituting $w_{ij}=1+sm_{ij}$ and $Q_{ij}=s\epsilon_{ij}$ for all $i\neq j$ in equation (2), and assuming linkage equilibrium ($D=0$) and condition (13), we arrive at

which reduce to Eqs. (15) as $s\rightarrow 0$. Here, $\bar{w}=1+s\bar{m}$. It is straightforward to note that the fixed points of Eqs. (15) and Eqs. (17), and their stability properties for small enough $s$ are same. This means that if the fixed points of Eqs. (15) are hyperbolic, then due to the system being a gradient system, the initial conditions will be attracted to some stable fixed points; and so is going to be the fate of the initial conditions under the map given by Eq. (2) whose corresponding fixed points are in $\Lambda_{s}$ manifold.

The search for which kind of replicator-mutator 2L2A system can be a gradient system interestingly stops at the most simple case of random point mutation, happening independent of what is happening elsewhere (e.g., at another locus): We call such mutations independent. It can be expressed in our setup as follows: Let the probability of mutation from $A_{i}$ allele to $A_{j}$ allele be given by $\mu_{ij}^{A}$ and the probability of mutation from $B_{i}$ allele to $B_{j}$ allele is given by $\mu_{ij}^{B}$. Naturally $\mu_{ii}^{A}$ and $\mu_{ii}^{B}$ correspond to the probabilities of no mutation from $A_{i}$ and $B_{i}$ alleles respectively. So the corresponding mutation matrix ${\sf Q}$ representing the probability of mutation from one gamete to another gamete, when the mutation at the two locus are independent, can be represented by the following matrix:

Now note that condition (13) for the gradient system in more convenient notions is

This is trivially satisfied by the aforementioned mutation matrix for independent mutation because of the fact that $\sum_{j=1}^{2}\mu^{A}_{ij}=\sum_{j=1}^{2}\mu^{B}_{ij}=1$—the normalization condition of probability. In conclusion, in the presence of independent point mutations at the two different loci, the 2L2A system always behaves as a gradient system in the weak selection limit and no complicated dynamics like oscillatory or chaotic dynamics can appear. It may be noted, however, that converse of this result need not be true: Mathematically speaking, mutations that are not independent (mentioned in Sec. I) may or may not satisfy the gradient condition.

Now let us see some of the example where non-convergent dynamical outcomes are possible in 2L2A model in the presence of mutation in the weak selection limit. As is clear from the preceding discussions, we have to consider non-independent mutations as they can result in non-gradient dynamics, and hence oscillatory outcomes.

Let us consider a simple fitness matrix well studied in literature Lewontin and Kojima (1960); Hastings (1981):

It basically says that a genotype homozygous at both locus has the fitness $k$, a genotype homozygous at $A$ ($B$) locus and heterozygous at $B$ ($A$) locus has the fitness $L_{1}$ ($L_{2}$), and a genotype heterozygous at both locus has the fitness $K$. Our aim is to illustrate the non-fixed point type dynamics that can appear in such systems in the presence of mutation.

To this end, let us take a mutation matrix that does not satisfies condition (13), and hence the corresponding dynamics given by Eqs. (5) is not a gradient system. Specifically, let us take $\epsilon_{31}=\epsilon_{24}=\epsilon_{1}$, $\epsilon_{12}=\epsilon_{43}=\epsilon_{2}$, and other mutation probabilities as zero. Furthermore, in order to achieve analytical tractability, let us choose $\nu_{1}+\nu_{2}-\nu_{3}-\nu_{4}=0$ and $\nu_{5}+\nu_{6}-\nu_{7}-\nu_{8}=0$ (note that our choice of ${\sf Q}$ satisfies these) so that the internal fixed point $(p^{*},q^{*})$ is conveniently located at $(1/2,1/2)$.

While performing linear stability analysis of Eqs. (5) about fixed point $(1/2,1/2)$, the Jacobian comes out to be

where $\alpha={(k-K)}/{4}-(\nu_{1}+\nu_{2}+\nu_{3}+\nu_{4})/{2},\beta=(\nu_{1}-\nu_{2}-\nu_{3}+\nu_{4})/2,\gamma=(\nu_{5}-\nu_{6}-\nu_{7}+\nu_{8})/2$, and $\delta={(k-K)}/{4}-(\nu_{5}+\nu_{6}+\nu_{7}+\nu_{8})/2$; and the corresponding eigenvalues are

For the Hopf bifurcation to occur one must have $4\beta\gamma<-(\alpha-\delta)^{2}$ and it happens at $\alpha+\delta=0$ which can be recast to read $(k-K)=\sum_{i=1}^{8}\nu_{i}$. For illustrative purpose, if we now fix $\epsilon_{2}=0.2$ and $(k-K)=1$, and work with $\epsilon_{1}$ as a variable parameter then at $\epsilon_{1}=0.3$, the Hopf bifurcation should occurs (see Fig. 1).

We observe that in the setting we are working, the condition of the Hopf bifurcation is independent of any constraints on $L_{1}$ and $L_{2}$. This has an interesting implication that we now point out: The case $K>L_{1},L_{2}$ and $k>L_{1},L_{2}$ corresponds to strong epistasis in literature Hastings (1981) where it was shown that for oscillatory behaviour to appear in 2L2A model (without mutation), strong epistasis is a necessary condition. However, when non-independent mutation is in action, the limit cycle can appear even in the absence of the strong epistasis because the corresponding Hopf bifurcation condition is dependent only on the values of $k$ and $K$.

In this context, we recall that the 2L2A gradient system closely approximates the corresponding discrete-time dynamics of 2L2A model in the weak selection limit: Each initial condition approaches a stable fixed point. Somewhat analogous effect is seen corresponding to the limit cycles in discrete-time 2L2A non-gradient system: For the same parameter values and sufficiently small $s$, an invariant closed curve is generated. However, in general, the orbits on the invariant curve need not necessarily be periodic with finite number of periodic points as another possibility is that the orbit on the curve can be everywhere dense (may be due to quasiperiodcity or chaos). As long as $s$ is sufficiently small (we can ignore $\mathcal{O}(s^{2})$ and higher order terms), it is easily seen that the time-discrete equation (which is effectively the Euler forward discretization of the continuous-time equation) under the same parameter values chosen above (while discussing the Hopf bifurcation), undergoes the Neimark–Sacker bifurcation. We do not present the calculations for the sake of avoiding trivial repetition. A numerical evidence of this is showcased in Fig. 1.

In conclusion, we have mathematically investigated the two-locus two-allele selection-recombination-mutation dynamics through the replicator equations in the limit of weak selection and weak mutation. We note that when the selection and mutation are sufficiently weak, one can ignore linkage disequilibrium completely and the system becomes two dimensional continuous time differential equation which faithfully approximates the evolutionary dynamics when the dynamics is gradient-like. Our search for the cases of mutation for which the dynamics is gradient-like has led to the result that whenever the point mutation rates at one locus is not affected by what happens at the other, the dynamics is always gradient-like. When the dynamics is not gradient-like, we found that stable limit cycle or stable closed invariant curve can appear. The dynamics on the invariant curve can, in principle, sustain non-convergent solutions like periodic, quasiperiodic, and chaotic orbits. We have also discussed how the dynamical manifestation of epistasis can completely changed in presence of mutation.

We are tempted to claim that our results imply that, within the paradigm of replicator dynamics of evolutionary process in the weak selection and mutation limit, an observation of oscillatory gametic or allelic frequencies is indirect indication of the existence of mutations that are not independent. An immediate future work that builds on the present paper is to find how the aforementioned results manifest themselves in a finite population where stochastic effects can no longer be ignored. Furthermore, many characteristics like our body size Reed et al. (2003), height Lettre et al. (2008), and the shape of different organs like the ear Adhikari et al. (2015) are controlled by more than two loci of a chromosome. So understanding the evolutionary dynamical outcome of multi-locus systems Zhivotovsky and Pylkov (1998); Wray and Goddard (2010), in the presence of mutation, is still another natural direction to investigate. While our main result—occurrence of gradient system in the presence of independent mutation—is valid for multi-locus two-allele systems (see Appendix A), the condition for general multi-locus multi-allele systems’ becoming gradient-like remains an open problem.