Bifurcation Analysis in Population Genetics Model with Partial Selfing

Abstract

A new model which allows both the effect of partial selfing selection and an exponential function of the expected payoff is considered. This combines ideas from genetics and evolutionary game theory. The aim of this work is to study the effects of partial selfing selection on the discrete dynamics of population evolution. It is shown that the system undergoes period doubling bifurcation, saddle-node bifurcation, and Neimark-Sacker bifurcation by using center manifold theorem and bifurcation theory. Numerical simulations are presented not only to illustrate our results with the theoretical analysis, but also to exhibit the complex dynamical behaviors, such as the period-3, 6 orbits, cascade of period-doubling bifurcation in period-2, 4, 8, and the chaotic sets. These results reveal richer dynamics of the discrete model compared with the model in Tao et al., 1999. The analysis and results in this paper are interesting in mathematics and biology.

1. Introduction

Evolutionary game theory was extended to the genetic model. The notion of an evolutionary stable strategy (ESS) which was proposed by Smith [1] (1982) is important. And the basic static solution concept of an ESS has been quite successful in predicting the equilibrium behavior of individuals in a population. A series of papers were devoted to the genetic matrix game models. But previous studies on the matrix game models all assumed that genetic mating was random. For example, Lessard [2] (1984) analyzed a frequency-dependent two-phenotype selection model of single-locus with multiallele. Cressman [3] (1996) studied a two-species evolutionary dynamics. Cressman et al. [4] (2003) discussed evolutionarily stable sets in the genetic model of natural selection. Tao et al. [5] (1999) investigated the discrete frequency dynamics of two phenotype diploid models. In fact, in genetic mating, there exists partial selfing selection (Rocheleau and Lessard [6]). From [6], each individual can reproduce by selfing or random outcrossing with constant probabilities, respectively, in a population. After considering partial selfing selection, we establish a nonlinear frequency dynamics, which becomes more realistic but more complicated.

We consider a one-locus two-allele model where genotypic fitness is an exponential function of the expected payoff and the frequency. Suppose that each individual of the population can reproduce by selfing with constant probability 𝛽. If 𝛽=0, then there is no partial selfing selection, which was discussed in [5]. If 𝛽≠0, then one-dimensional discrete-time dynamical systems are extended to two-dimensional systems, which lead to great effect on the genetic model. It was shown in [7] that a stable polymorphic equilibrium is not an ESS under some condition. In this paper, we focus on the bifurcation analysis of the genetic model. We find that the model admits the unstable period doubling bifurcation, the saddle-node bifurcation, and the Neimark-Sacker bifurcation, which are new interesting phenomena when considering partial selfing selection.

2. Model

Consider a diploid population with nonoverlapping generations and establish a genetic model with one locus and two alleles 𝐴₁ and 𝐴₂. The fitness of population individual is frequency-dependent (i.e., it depends on the genetic makeup of the population concerned). To show the results, the following assumptions will be made: (i)Mendelian segregation; (ii)sex ratio is independent of genotype; (iii)no gametic selection; (iv)the fecundity of offspring is equal; (v)no mutation or migration.

According to [6], each individual can reproduce by selfing or random outcrossing with probability 𝛽 or 1−𝛽  (0<𝛽<1), respectively, in the population. 𝑃_𝑖𝑗 and 𝑃^′_𝑖𝑗 are the frequencies of genotypes 𝐴_𝑖𝐴_𝑗 in the current generation and in the next generation, respectively. The genotypic frequency of genotype 𝐴_𝑖𝐴_𝑗 among adults in the current generation is 𝑃_𝑖𝑗. The corresponding selection values of genotypes 𝐴_𝑖𝐴_𝑗 are 𝐹_{𝐴𝑖𝐴𝑗}. 𝑝 and 𝑝^′ are the frequencies of allele 𝐴₁ in the current generation and in the next generation, respectively. 𝑞 and 𝑞^′ are the frequencies of allele 𝐴₂ in the current generation and in the next generation, respectively. 𝑝 and 𝑞 are genotypic frequencies of allele 𝐴₁ and allele 𝐴₂ among adults in the current generation, respectively. 𝑢_𝑖𝑗  (or 1−𝑢_𝑖𝑗) is the fraction of individuals of genotype 𝐴_𝑖𝐴_𝑗 which play pure strategy 𝑅₁  (or 𝑅₂). After selection but before mating, we have 𝑃₁₁=𝐹_𝐴1𝐴1𝑃₁₁𝐹_𝐴1𝐴1𝑃₁₁+𝐹_𝐴1𝐴2𝑃₁₂+𝐹_𝐴2𝐴2𝑃₂₂,𝑃₁₂=𝐹_𝐴1𝐴2𝑃₁₂𝐹_𝐴1𝐴1𝑃₁₁+𝐹_𝐴1𝐴2𝑃₁₂+𝐹_𝐴2𝐴2𝑃₂₂,𝑃₂₂=𝐹_𝐴2𝐴2𝑃₂₂𝐹_𝐴1𝐴1𝑃₁₁+𝐹_𝐴1𝐴2𝑃₁₂+𝐹_𝐴2𝐴2𝑃₂₂. (1) After mating and reproduction, there are 𝑃^′₁₁=𝛽𝑃₁₁+14𝑃₁₂+1−𝛽𝑃₁₁+12𝑃₁₂²,𝑃^′₁₂1=𝛽2𝑃₁₂+21−𝛽𝑃₁₁+12𝑃₁₂𝑃₂₂+12𝑃₁₂,𝑃^′₂₂=𝛽𝑃₂₂+14𝑃₁₂+1−𝛽𝑃₂₂+12𝑃₁₂². (2)

In our paper, the fitness 𝐹_{𝐴𝑖𝐴𝑗} is taken to be an exponential function of the expected payoff as [6]. Assume that an individual in the population can use one of two strategies 𝑅₁ and 𝑅₂; each individual plays the same pure strategy throughout its lifetime. And the payoff matrix is 𝑎₁₁𝑎₁₂𝑎₂₁𝑎₂₂. (3) Let 𝑥 be the proportion of individuals using strategy 𝑅₁; then, we can get that the expected payoffs to 𝑅₁ and 𝑅₂ are 𝜑₁=𝑥𝑎₁₁+(1−𝑥)𝑎₁₂,𝜑₂=𝑥𝑎₂₁+(1−𝑥)𝑎₂₂. (4) Furthermore, we take the fitness function 𝐹_{𝐴𝑖𝐴𝑗}=𝑢_𝑖𝑗𝑒^𝜑1+1−𝑢_𝑖𝑗𝑒^𝜑2. (5) Since system (2) is on the simplex 𝑆₃={(𝑃₁₁,𝑃₁₂,𝑃₂₂)∣𝑃₁₁+𝑃₁₂+𝑃₂₂=1;𝑃₁₁,𝑃₁₂,𝑃₂₂≥0}, we can get the following system: 𝑃^′₁₁𝐹=𝛽^𝐴1𝐴1𝑃₁₁+(1/4)𝐹_𝐴1𝐴2𝑃₁₂𝐹+𝐹1−𝛽^𝐴1𝐴1𝑃₁₁+(1/2)𝐹_𝐴1𝐴2𝑃₁₂²𝐹²,𝑃^′₁₂()𝐹=𝛽1/2_𝐴1𝐴2𝑃₁₂𝐹×𝐹+21−𝛽_𝐴1𝐴1𝑃₁₁+12𝐹_𝐴1𝐴2𝑃₁₂×𝐹_𝐴2𝐴21−𝑃₁₁−𝑃₁₂+12𝐹_𝐴1𝐴2𝑃₁₂×𝐹²⁻¹, (6) where 𝑥=𝑃₁₁𝑢₁₁+𝑃₁₂𝑢₁₂+𝑃₂₂𝑢₂₂,𝜑_𝑖=𝑥𝑎_𝑖1+(1−𝑥)𝑎_𝑖2,𝐹_{𝐴𝑖𝐴𝑗}=𝑢_𝑖𝑗𝑒^𝜑1+1−𝑢_𝑖𝑗𝑒^𝜑2,𝐹=𝐹_𝐴1𝐴1𝑃₁₁+𝐹_𝐴1𝐴2𝑃₁₂+𝐹_𝐴2𝐴2𝑃₂₂. (7) From (1) and (2), we obtain 𝑝=𝑃₁₁+12𝑃₁₂=𝑃^′₁₁+12𝑃^′₁₂=𝑝^′. (8) Obviously, it is difficult to use 𝑝 and 𝑞 to denote the genotypic frequency 𝑃_𝑖𝑗 directly. If the population is polymorphic (i.e., 0<𝑝<1), we can use the fixation index 𝐹 (wright 1949 [8]) according to [6]. The genotypic frequencies can be written in the form 𝑃₁₁=𝑝²𝑃+𝑝𝑞𝐹,₁₂=2𝑝𝑞(1−𝐹),𝑃₂₂=𝑞²+𝑝𝑞𝐹, (9) where −1≤𝐹≤1. The value of 𝐹 varies from one generation to the next.

Equations (1) to (9) can deduce that𝑝^′𝑢=𝑝𝑝+𝑞𝐹₁₁𝑒^𝜑1+1−𝑢₁₁𝑒^𝜑2+(1−𝐹)𝑞𝑢₁₂𝑒^𝜑1+1−𝑢₁₂𝑒^𝜑2𝐹,𝐹^′=𝛽1−1−𝐹2𝑝1−1−𝐹_𝐴2𝐴2/𝐹_𝐴1𝐴2+𝑞𝑞+𝑝𝐹1−1−𝐹_𝐴1𝐴1/𝐹_𝐴1𝐴2,𝑝+𝑞𝐹 (10)where 𝑥_𝐴1=𝑝+𝑞𝐹𝑢₁₁+𝑞(1−𝐹)𝑢₁₂,𝑥_𝐴2=𝑞+𝑝𝐹𝑢₂₂+𝑝(1−𝐹)𝑢₁₂,𝑥𝑝,𝐹=𝑝²+𝑝𝑞𝐹𝑢₁₁+2𝑝𝑞(1−𝐹)𝑢₁₂+𝑞²+𝑝𝑞𝐹𝑢₂₂=𝑝𝑥_𝐴1+𝑞𝑥_𝐴2,𝜑_𝑖=𝑥𝑎_𝑖1+(1−𝑥)𝑎_𝑖2,𝐹_{𝐴𝑖𝐴𝑗}=𝑢_𝑖𝑗𝑒^𝜑1+1−𝑢_𝑖𝑗𝑒^𝜑2,𝐹=𝑥𝑒^𝜑1+(1−𝑥)𝑒^𝜑2. (11)

Throughout the paper, we discard the degenerate situations where all 𝑢_𝑖𝑗 are identical or where 𝛾=0 and suppose that 𝐴_𝑖𝐴_𝑗 is not the genotype in the parental generation.

3. Model Analysis and Basic Definitions

Recall Definition  1 in [7], which is a similar definition of phenotypic and genotypic equilibria under partial selfing selection according to [2]. The two types of the equilibria include all situations in which the population is in equilibrium.

Definition 1. A phenotypic equilibrium is an equilibrium where all pure strategies in current use have equal expected payoff. A genotypic equilibrium is a nonphenotypic equilibrium where the effective fitnesses of all alleles present in the current population are equal. The effective fitness of allele 𝐴_𝑖 is defined to be [𝑃_𝑖𝑖𝐹_{𝐴𝑖𝐴𝑖}+(1/2)𝑃₁₂𝐹_𝐴1𝐴2]/[𝑃_𝑖𝑖+(1/2)𝑃₁₂]. (In addition, in the case 𝑝=0, the effective fitness of allele 𝐴₁ is 𝐹_𝐴1𝐴2. In the case 𝑝=1, the effective fitness of allele 𝐴₂ is 𝐹_𝐴1𝐴2.)

According to Definition 1, a phenotypic equilibrium is an equilibrium where all pure strategies in current use have equal expected payoff. So, a polymorphic phenotypic equilibrium exists if and only if 𝜑₁=𝜑₂. Then, we have 𝜑₂−𝜑₁=𝑥𝑎₂₁−𝑎₁₁+(1−𝑥)𝑎₂₂−𝑎₁₂=𝛾(𝑥−𝜉), (12) where 𝛾=𝑎₁₂−𝑎₂₂+𝑎₂₁−𝑎₁₁, and 𝜉=(𝑎₁₂−𝑎₂₂)/(𝑎₁₂−𝑎₂₂+𝑎₂₁−𝑎₁₁).

Let (𝑝^∗,𝐹^∗) denote polymorphic phenotypic equilibrium; then, we have the Jacobian matrix at (𝑝^∗,𝐹^∗) of system (10) 𝐴^∗=1−𝑘₁₁𝛾𝑘₁₂𝛾𝑘₂₁𝛾𝛽2−𝑘₂₂𝛾, (13) where 𝑢=𝑢₁₁+𝑢₂₂−2𝑢₁₂,𝑎=𝑢₁₁−𝑢₁₂,𝑏=𝑢₂₂−𝑢₁₂,𝐹^∗=𝛽,𝑘2−𝛽₂=𝑝^∗2𝑢₂₂−𝑢₁₂+𝑞^∗2𝑢₁₁−𝑢₁₂,𝑘=𝑢₁₂−𝑢₁₁+𝑝^∗𝑢,𝑧=𝑢₂₂−𝑢₁₂−𝑝^∗𝑘𝑢,₁₁=𝑝^∗1−𝑝^∗𝑧+𝐹^∗𝑘2−𝐹^∗𝑧+𝐹^∗𝑘,𝑘₁₂=𝑝^∗2𝑞^∗2𝑢𝑧+𝐹^∗𝑘,𝑘₂₁=𝛽21−𝐹^∗𝑝^∗𝑞^∗𝑢+𝐹^∗𝑘₂2−𝐹^∗𝑧+𝐹^∗𝑘,𝑘₂₂=𝛽21−𝐹^∗𝑝^∗𝑞^∗𝑢𝑝^∗𝑞^∗𝑢+𝐹^∗𝑘₂. (14)

We studied the system (10) and obtained Theorem 3.1 (i) and (ii) in [7]. In [7], the authors showed that a stable polymorphic equilibrium is not an ESS under some condition. In this paper, we focus on the bifurcation of system (10). For convenience, we present some lemmas and theorems of the stability of polymorphic phenotypic equilibrium in [7] at first.

Lemma 2. 𝑥^∗=𝑥(𝑝^∗,𝐹^∗)=𝜉(0<𝜉<1) is an interior ESS if and only if 𝛾>0.

Lemma 3. If 𝑘₁₁>0, then (𝛽/2)𝑘₁₁+𝑘₂₂>0. If 𝑘₁₁<0, then ((2+𝛽)/4)𝑘₁₁+𝑘₂₂>0.

Theorem 4. Suppose that (𝑝^∗,𝐹^∗) is a polymorphic phenotypic equilibrium.(i)Suppose that 𝑘₁₁>0. If 𝑥^∗ is not an ESS, then (𝑝^∗,𝐹^∗) is unstable. If 𝑥^∗ is an ESS, then (𝑝^∗,𝐹^∗) is stable for 0<𝛾<𝛾_𝑐𝛾2. (ii)Suppose that 𝑘₁₁<0. If 𝑥^∗ is an ESS, then (𝑝^∗,𝐹^∗) is unstable. If 𝑥^∗ is not an ESS, then (𝑝^∗,𝐹^∗) is stable for 𝛾_𝑐𝛾1<𝛾<0.(iii)Suppose that 𝑘₁₁=0. If 𝑎≠𝑏, then (𝑝^∗,𝐹^∗) is unstable for 𝛾≠(𝛽±2)/2𝑘₂₂. If 𝑎=𝑏, then (𝑝^∗,𝐹^∗) is stable for 𝛾∈(0,4(𝛽+2)/𝛽𝑎²(1−𝐹^∗2)), unstable for 𝛾∈(−∞,−4/𝑎²(1−𝐹^∗2)𝐹^∗)∪(−4/𝑎²(1−𝐹^∗2)𝐹^∗,0)∪(4(𝛽+2)/𝛽𝑎²(1−𝐹^∗2),+∞).

Proof. Consider only the case 𝑘₁₁>0 since the proof of the case 𝑘₁₁<0 is analogous. From the Jacobian matrix 𝐴^∗ at (𝑝^∗,𝐹^∗), we have det𝐴^∗𝛽=−2𝑘₁₁+𝑘₂₂𝛽𝛾+2,tr𝐴^∗=2+𝛽2−𝛾𝑘₁₁+𝑘₂₂. (15) If 𝑥^∗ is not an ESS, then 𝛾<0 from Lemma 2. The inequality 𝑘₁₁𝛾<0 implies that tr𝐴^∗>det𝐴^∗+1. So, (𝑝^∗,𝐹^∗) is unstable. If 𝑥^∗ is an ESS, then 𝛾>0. By the discussion on the monotone function, we can obtain the following inequalities: det𝐴^∗<1,tr𝐴^∗<det𝐴^∗+1. (16) So, (𝑝^∗,𝐹^∗) is stable.
Now, we discuss the stability of (𝑝^∗,𝐹^∗) when 𝑘₁₁=0. If 𝑢=0, then 𝑘₁₁=0. If 𝑘₁₁=0, we have the following.
Case  1. (𝑝^∗,𝐹^∗)=(𝑝₂,𝛽/(2−𝛽)), where 𝑝₂=(2𝑏−𝑢𝐹^∗)/2𝑢(1−𝐹^∗) and 𝐹^∗=𝛽/(2−𝛽).
Case  2. (𝑝^∗,𝐹^∗)=(𝑝₁,𝛽/(2−𝛽)), where 𝑝₁=(𝑏−𝑎𝐹^∗)/𝑢(1−𝐹^∗) and 𝐹^∗=𝛽/(2−𝛽).
We only prove Case 1 since Case 2 is analogous.
Write 𝑎₁=𝑘₁₂𝛾,𝑏₁=𝛽2−𝑘₂₂𝑎𝛾,₁₁=−2𝑢1−𝐹^∗2𝑝^∗𝑞^∗𝑎𝑧𝛾,₁₂=𝑝^∗2𝑞^∗2𝑢2𝑘+𝑧+𝐹^∗𝑘𝑝^∗𝑞^∗𝑎𝑢𝛾𝛾,₂₂=𝑝^∗𝑞^∗𝑢21−2𝑝^∗𝑧+𝐹^∗𝑘−𝑝^∗𝑞^∗𝑢1−𝐹^∗𝑏𝛾,₁₁=−𝛽,𝑢1−𝐹^∗2𝑝^∗𝑞^∗𝑢+𝑏𝑝^∗2+𝑎𝑞^∗2𝐹^∗𝛾, (17) and let 𝑃=^{1−𝑎1/(1−𝑏1)}₀₁. Using the translation 𝑝𝐹𝜇𝜈+𝑝=𝑃²𝐹^∗, (18) system (10) becomes 𝜇^′1=𝜇+2×𝑎₁₁𝑎𝜇−¹1−𝑏₁𝜈²+𝑎₁₂𝜈𝑎𝜇−¹1−𝑏₁𝜈+𝑎₂₂𝜈²+𝑎₁21−𝑏₁⋅𝑏₁₁𝑎𝜇−¹1−𝑏₁𝜈²+𝑏₁₂𝜈𝑎𝜇−¹1−𝑏₁𝜈+𝑏₂₂𝜈²𝜈+h.o.t,_′=𝑏₁𝜈+12𝑏₁₁𝑎𝜇−¹1−𝑏₁𝜈²+𝑏₁₂𝜈𝑎𝜇−¹1−𝑏₁𝜈+𝑏₂₂𝜈²+h.o.t. (19) By center manifold theory, we can obtain the following reduced system which is locally homeomorphic with system (10): 𝜇^′1=𝜇+2𝑎₁₁+𝑎₁1−𝑏₁𝑏₁₁𝜇²+𝑂𝜇³, (20) where (1/2)(𝑎₁₁+(𝑎₁/(1−𝑏₁))𝑏₁₁)   =  𝑝^∗𝑞^∗𝑢𝑧(1−𝐹^∗)²𝛾{(1/(1/𝑝^∗𝑞^∗𝑢[𝑝^∗𝑞^∗𝑢+𝐹^∗(𝑎𝑞^∗2+𝑏𝑝^∗2)𝛾]𝐹^∗(1−𝐹^∗))+1)−1}. If 𝑎≠𝑏, then 𝑧≠0. So, we have (1/2)(𝑎₁₁+(𝑎₁/(1−𝑏₁))𝑏₁₁)≠0. If 𝑎=𝑏, we can calculate the reduced system by similar analysis, 𝑥^′=𝑥+𝑘₃𝑥³+𝑂𝑥⁴, (21) where 𝑘₃=−(1−𝐹^∗)²𝑎²𝛾/(1+(1/4)𝑎²(1−𝐹^∗2)𝐹^∗𝛾). It is easy to obtain the results. This completes the proof.

4. Bifurcation Analysis

Based on the analysis in Section 3, we discuss the period doubling bifurcations, the saddle-node bifurcation, and the Neimark-Sacker bifurcation of the positive fixed point (𝑃^∗,𝐹^∗) in this section. We choose parameter 𝛾 as a bifurcation parameter to study the period doubling bifurcations and the Neimark-Sacker bifurcation and parameter 𝛿 as a bifurcation parameter to study the saddle-node bifurcation by using center manifold theorem and bifurcation theory in [9, 10]. Suppose that 𝑘_𝑖𝑗 is the same as that in (14).

4.1. Period Doubling Bifurcation

In the analysis of period doubling bifurcations, we take 𝛾 as the bifurcation parameter and prove that there is period doubling bifurcation at the fixed point (𝑝^∗,𝐹^∗) for 𝛾=𝛾_𝑐𝛾2. When 𝛾=𝛾_𝑐𝛾2, the characteristic polynomial of Jacobian matrix at (𝑝^∗,𝐹^∗) is 𝜆²+4−𝛽²𝑘₁₁22+𝛽𝑘₁₁+4𝑘₂₂𝛽𝑘𝜆−𝛽+2₁₁+8𝑘₂₂22+𝛽𝑘₁₁+4𝑘₂₂=0. (22) We have the following characteristic values of (22): 𝜆₁=−1,𝜆₂𝛽𝑘=𝑠=𝛽+2₁₁+8𝑘₂₂22+𝛽𝑘₁₁+4𝑘₂₂. (23) Let 𝑐𝐶=¹¹𝑐₁₂𝑐₂₁𝑐₂₂=𝛽+2𝑘₁₁4𝑘₁₂−4𝑘₂₁𝛽+2𝑘₁₁, (24) And let 𝐷 denote 𝐶⁻¹. Then, 𝑑𝐷=¹¹𝑑₁₂𝑑₂₁𝑑₂₂=1𝛽+2²𝑘²₁₁+16𝑘₁₂𝑘₁₂𝛽+2𝑘₁₁−4𝑘₁₂4𝑘₂₁𝛽+2𝑘₁₁. (25)

Denote the right-hand side of (10) by 𝑝𝐹, and use the translation 𝑝𝐹𝑥𝑦+𝑝=𝐶^∗𝐹^∗,𝛾=𝛾_𝑐𝛾2+𝜆; (26) system (10) becomes 𝑥^′𝑦^′=𝑥𝑦+𝐻𝐺−100𝑠𝑥,𝑦,𝜆𝑥,𝑦,𝜆, (27) where 𝐻𝑥,𝑦,𝜆=𝑥+𝑑₁₁𝑓𝑥,𝑦,𝜆−𝑝^∗+𝑑₁₂𝑔𝑥,𝑦,𝜆−𝐹^∗,𝐺𝑥,𝑦,𝜆=−𝑠𝑦+𝑑₂₁𝑓𝑥,𝑦,𝜆−𝑝^∗+𝑑₂₂𝑔𝑥,𝑦,𝜆−𝐹^∗. (28) It is easy to obtain that 𝐻()𝐺()0,0,𝜆0,0,𝜆=0,𝜕𝐻𝑥,𝑦,𝜆𝜕𝑥𝜕𝐻𝑥,𝑦,𝜆𝜕𝑦𝜕𝐺𝑥,𝑦,𝜆𝜕𝑥𝜕𝐺𝑥,𝑦,𝜆𝜕𝑦_{(𝑥,𝑦,𝜆)=(0,0,0)}𝜕=0,₂𝐻𝑥,𝑦,𝜆𝜕𝑥𝜕𝜆_{(𝑥,𝑦,𝜆)=(0,0,0)}=𝑑₁₁−𝑐₁₁𝑘₁₁+𝑐₂₁𝑘₂₁+𝑑₁₂𝑐₁₁𝑘₂₁−𝑐₂₁𝑘₂₂𝑘=−𝛽+2₁₁+4𝑘₂₂²𝛽+2²𝑘₁₁+16𝑘₂₂. (29) According to Lemma 3, there are 𝛽+2𝑘₁₁+4𝑘₂₂≠0,𝛽+2²𝑘₁₁+16𝑘₂₂>0. (30) So, we have 𝜕𝑐=−2²𝐻𝑥,𝑦,𝜆𝜕𝑥𝜕𝜆>0. (31) In order for system (27) to undergo period doubling bifurcation, we require that the following 𝐿 is not zero [10]: 𝐿=−2𝐻_𝑥𝑥𝑥−3𝐻²_𝑥𝑥−6𝐺1−𝑠^𝑥𝑥𝐻_𝑥𝑦. (32) If 𝑘₁₁≠0, then |𝑠|<1 for 𝑘₁₁>0, |𝑠|>1 for 𝑘₁₁<0.

From the previous analysis and the theorem in [10], we have the following result.

Proposition 5. Suppose that (𝑝^∗,𝐹^∗) is a polymorphic phenotypic equilibrium. If 𝑘₁₁≠0 and 𝐿≠0, then system (10) undergoes a period doubling bifurcation at (𝑝^∗,𝐹^∗) for 𝛾=𝛾_𝑐𝛾2.

Moreover, we have the following.(i)If 𝑘₁₁>0 and 𝐿>0  (𝐿<0), then period doubling bifurcation happens for 𝛾<𝛾_𝑐𝛾2(𝛾>𝛾_𝑐𝛾2). And 2-periodic points are saddle points (asymptotically stable nodes).(ii)If 𝑘₁₁<0 and 𝐿>0  (𝐿<0), then period doubling bifurcation happens for 𝛾<𝛾_𝑐𝛾2(𝛾>𝛾_𝑐𝛾2). And 2-periodic points are unstable nodes (saddle points).

Example 6. Let 𝑢₁₁=0.6, 𝑢₁₂=0.9, 𝑢₂₂=0.3, 𝛽=0.4, and 𝜉=0.6515625. There are two polymorphic phenotypic equilibria (𝑝^∗₁,𝐹^∗)=(0.6944444,0.25) and  (𝑝^∗₂,𝐹^∗)=(0.75,0.25). For (𝑝^∗₁,𝐹^∗), we have 𝑘₁₁(𝑝^∗₁,𝐹^∗)=0.4475911459×10⁻³>0 and 𝐿(𝑝^∗₁,𝐹^∗)=−4.839394634<0. Then, period doubling bifurcation happens for 𝛾>𝛾_𝑐𝛾2=149.7775520, and 2-periodic points are asymptotically stable nodes. For (𝑝^∗₂,𝐹^∗), we have 𝑘₁₁(𝑝^∗₂,𝐹^∗)=−0.1318359375×10⁻³<0 and 𝐿(𝑝^∗₂,𝐹^∗)=9.126558784>0. Then, period doubling bifurcation happens for 𝛾<𝛾_𝑐𝛾2=186.1395138, and 2-periodic points are unstable nodes.
The effects of partial selfing selection on the dynamics of population evolution are shown further in Figures 1, 2, and 3. More complex dynamical behaviors of the genetic system are exhibited by numerical simulations. In Figure 1, the partial selfing selection leads period doubling bifurcation to emerge earlier, and leads chaos to emerge earlier. In Figure 3, the model exhibits the complex dynamical behaviors, such as the period-3, 6 orbits, cascade of period-doubling bifurcation in period-2, 4, 8, and the chaotic sets. By choosing 𝛽 as a bifurcation parameter, we show that the complex dynamical behaviors such as the period-3, 4, 6 orbits, cascade of period-doubling bifurcation in period-2, 4, 8, and the chaotic sets can occur as 𝛽 crosses some critical values in Figure 2.

(a)

(b)

Figure 1 

𝑢11=1, 𝑢12=0.5, 𝑢22=0, and 𝜉=0.85. If 𝛾<0, then all interior trajectories evolve to the stable monomorphic equilibria (𝑝∗=0 and 𝑝∗=1). (a) For 𝛽=0, there is unique polymorphic phenotypic equilibrium (𝑝∗,𝐹∗)=(0.85,0.25), which is stable for 0<𝛾<𝛾𝑐𝛾2=31.37254902. (b) For 𝛽=0.4, there is unique polymorphic phenotypic equilibrium (𝑝∗,𝐹∗)=(0.85,0.25), which is stable for 0<𝛾<𝛾′𝑐𝛾2=25.09803922.

(a)

(b)

(c)

(d)

Figure 2 

𝑢11=1, 𝑢12=0.5, 𝑢22=0, and 𝜉=0.85. When 𝛾=27, (𝑝∗,𝐹∗) is stable for 𝛽∈(0,𝛽∗)(𝛽∗=0.27875). Stable 2-periodic points emerge at 𝛽=𝛽∗. When 𝛾=38, stable 2-periodic points become stable 4-periodic points with the increase of 𝛽. When 𝛾=41, there is chaos for 0.47<𝛽<0.49. When 𝛾=69, stable 3-periodic points become stable 6-periodic points and then become stable 3-periodic points again.

Figure 3 

𝑢11=1, 𝑢12=0.5, 𝑢22=0, 𝜉=0.85, and 𝛽=0.4. As 𝛾 is sufficiently large, there are stable 3-periodic points (at 𝛾≃56.5), and then stable 6-periodic points occur at 𝛾≃67.5. And stable 3-periodic points occur at 𝛾≃201 again.

4.2. Saddle-Node Bifurcation

We take 𝜉 as the bifurcation parameter for studying the saddle-node bifurcation by using center manifold theorem.

If 𝜉^∗=𝑢₂₂−((𝑢𝐹^∗−2𝑏)²/4𝑢(1−𝐹^∗)), then there is unique polymorphic phenotypic equilibrium, and (𝑝^∗,𝐹^∗)=((2𝑏−𝑢𝐹^∗)/2𝑢(1−𝐹^∗),𝐹^∗). The Jacobian matrix of (𝑝^∗,𝐹^∗) is 𝐴^∗=1𝑘^∗₁₂𝛾0𝛽2−𝑘^∗₂₂𝛾, (33) where 𝑘^∗₁₂=(2𝑏−𝑢𝐹^∗)²(2𝑎−𝑢𝐹^∗)²(𝑏−𝑎)𝐹^∗/32𝑢(1−𝐹^∗)² and 𝑘^∗₂₂=(𝛽/32)(((2𝑏−𝑢𝐹^∗)(2𝑎−𝑢𝐹^∗))/𝑢²(1−𝐹^∗)²)[4𝑎𝑏+2(𝑎−𝑏)²𝐹^∗−𝑢²𝐹^∗2].

The characteristic values of 𝐴^∗ are 𝜆₁=1,𝜆₂=𝛽2−𝑘^∗₂₂𝛾, (34) where 𝛾≠(𝛽±2)/2𝑘^∗₂₂ can make sure that |𝜆₂|≠1.

Using the translation 𝑥𝑦=𝑘1−^∗12𝛾1−𝛽/2+𝑘^∗₂₂𝛾01𝑝−𝑝^∗𝐹−𝐹^∗, (35) we have 𝑥^′𝑦^′=100𝜆^∗𝑥𝑦+𝑓₁𝜉̃𝑔𝑥,𝑦,₁𝑥,𝑦,𝜉̃, (36) where 𝑓₁=𝑎^∗̃1𝜉+2𝑎₁₁𝑥−𝐵₀𝑦²+𝑎₁₂𝑦𝑥−𝐵₀𝑦+𝑎₂₂𝑦²+𝑎₀₁𝑥−𝐵₀𝑦𝜉̃+𝑎₀₂𝑦𝜉̃+𝑎₀₃𝜉̃²+12𝐵₀𝑏₁₁𝑥−𝐵₀𝑦²+𝑏₁₂𝑥−𝐵₀𝑦𝑦+𝑏₂₂𝑦²+𝑏₀₁𝑥−𝐵₀𝑦𝜉̃+𝑏₀₂𝑦𝜉̃+𝑏₀₃𝜉̃²𝑔+h.o.t,₁=𝑏^∗̃1𝜉+2𝑏₁₁𝑥−𝐵₀𝑦²+𝑏₁₂𝑦𝑥−𝐵₀𝑦+𝑏₂₂𝑦²+𝑏₀₁𝑥−𝐵₀𝑦𝜉̃+𝑏₀₂𝑦𝜉̃+𝑏₀₃𝜉̃²𝜆+h.o.t,_∗=𝛽2−𝑘^∗₂₂𝛾,𝐵₀=𝑘^∗₁₂𝛾1−𝛽/2+𝑘^∗₂₂𝛾,𝑎^∗=𝑎−𝑏2𝑝^∗𝑞^∗𝐹^∗𝛾1−𝛽/21−𝛽/2+𝑘^∗₂₂𝛾,𝑏^∗=𝛽𝛾8𝑢4𝑎𝑏−𝑢²𝐹^∗2+2𝐹^∗(𝑎−𝑏)²,𝑎₁₁=−2𝑢1−𝐹^∗2𝑝^∗𝑞^∗𝑎𝑧𝛾,₁₂=𝑝^∗2𝑞^∗2𝑢2𝑘+𝑧+𝐹^∗𝑘𝑝^∗𝑞^∗𝑎𝑢𝛾𝛾,₂₂=𝑝^∗𝑞^∗𝑢21−2𝑝^∗𝑧+𝐹^∗𝑘−𝑝^∗𝑞^∗𝑢1−𝐹^∗𝑏𝛾,₁₁=−𝛽𝑢1−𝐹^∗2𝑝^∗𝑞^∗𝑢+𝑏𝑝^∗2+𝑎𝑞^∗2𝐹^∗𝛾. (37)

By the center manifold theory, we know that the stability of (0,0) near 𝜉̃=0 can be determined by a one-parameter family of equations on a center manifold, which can be represented as follows: 𝑊^𝑐(0)=𝑥,𝑦,𝜉̃∈𝑅³𝜉̃(),()∣𝑦=ℎ𝑥,,ℎ0,0=0,𝐷ℎ0,0=0∈𝑅³ (38) for 𝑥 and 𝜉̃ sufficiently small. Assume that a center manifold has the form ℎ𝑥,𝜉̃=𝑙₀𝜉̃+𝑙₁𝑥𝜉̃+𝑙₂𝑥²+𝑙₃𝜉̃²+𝑂|𝑥|+|𝜉̃|³. (39) The center manifold must satisfy 𝑁ℎ𝑥,𝜉̃=ℎ𝐴^∗𝑥+𝑓₁,𝜉̃−𝜆^∗ℎ𝑥,𝜉̃−𝑔₁𝑥,ℎ𝑥,𝜉̃,𝜉̃=0. (40) So, we obtain the map restricted to the center manifold 𝑥→𝑓𝑥,𝜉̃=𝑥+𝑐₁𝜉̃+𝑐₂𝑥²+𝑂|𝑥|+|𝜉̃|², (41) where 𝑐₁=((𝑎−𝑏)/2)𝑝^∗𝑞^∗𝐹^∗𝛾((1−𝛽/2)/(1−(𝛽/2)+𝑘^∗₂₂𝛾)), and 𝑐₂=−((𝑎−𝑏)/4)𝑝^∗𝑞^∗(1−𝐹^∗)𝑢𝛾(𝛽/(1−(𝛽/2)+𝑘^∗₂₂𝛾)).

By condition (41) and the saddle-node bifurcation theorem in [9], we can state the following result.

Proposition 7. Suppose that (𝑝^∗,𝐹^∗) is a polymorphic phenotypic equilibrium. If 𝛽∈(0,min{4𝑏/(𝑎+3𝑏),4𝑎/(3𝑎+𝑏)}), (𝑎−𝑏)𝑢𝛾≠0, 𝑎𝑏>0, and 𝛾≠(𝛽±2)/2𝑘^∗₂₂, then system (10) undergoes a saddle-node bifurcation at (𝑝^∗,𝐹^∗) for 𝜉=𝜉^∗ where 𝑝^∗=(2𝑏−𝑢𝐹^∗)/2𝑢(1−𝐹^∗), 𝜉^∗=𝑢₂₂−((𝑢𝐹^∗−2𝑏)²/4𝑢(1−𝐹^∗), and 𝑘^∗₂₂=(𝛽/32)((2𝑏−𝑢𝐹^∗)(2𝑎−𝑢𝐹^∗)/𝑢²(1−𝐹^∗)²)[4𝑎𝑏+2(𝑎−𝑏)²𝐹^∗−𝑢²𝐹^∗2].

Moreover, if 𝑢>0, then two new polymorphic phenotypic equilibria are created for 𝜉>𝜉^∗ and disappear for 𝜉<𝜉^∗. If 𝑢<0, then two new polymorphic phenotypic equilibria are created for 𝜉<𝜉^∗ and disappear for 𝜉>𝜉^∗. (see Figure 4).

Figure 4 

𝑢11=0.5, 𝑢12=1, 𝑢22=0, 𝛽=0.6, and 𝛾=10. A saddle-node bifurcation happens for (𝑝∗,𝐹∗,𝜉∗), where 𝑝∗=0.7916666665,𝐹∗=0.4285714286,and𝜉∗=0.5372023808. Two new polymorphic phenotypic equilibria are created for 𝜉<𝜉∗ and disappear for 𝜉>𝜉∗.

In Proposition 7, the conditions 𝛽∈(0,min{4𝑏/(𝑎+3𝑏),4𝑎/(3𝑎+𝑏)}) and 𝑎𝑏>0 guarantee that 𝑝^∗=(2𝑏−𝑢𝐹^∗)/2𝑢(1−𝐹^∗)∈(0,1). If there is no partial selfing selection (i.e., 𝛽=0), then the saddle-node bifurcation does not occur. This means that the parameter 𝛽 of partial selfing selection leads to more complex dynamical behavior of the genetic system.

4.3. Neimark-Sacker Bifurcation

We take 𝛾 as the bifurcation parameter and prove the existence of Neimark-Sacker bifurcation. The characteristic polynomial of Jacobian matrix at (𝑝^∗,𝐹^∗) is 𝜆²+𝑀𝛾𝜆+𝑁𝛾=0, (42) where 𝑀𝛾=𝑘₁₁+𝑘₂₂𝛾−𝛽+22,𝑁𝛾𝛽=−2𝑘₁₁+𝑘₂₂𝛽𝛾+2. (43) If the following condition is satisfied 𝑀𝛾²−4𝑁𝛾<0, (44) then the eigenvalues of the characteristic equation are complex conjugate. When 𝛾=𝛾₀=(𝛽−2)/(𝑘₁₁𝛽+2𝑘₂₂), (44) is equivalent to 𝑘₁₁<0. (45) If condition (45) is satisfied, we can calculate that 𝑑𝜆𝛾𝑑=𝑑𝛾𝛾=𝛾₀1=−2𝛽2𝑘₁₁+𝑘₂₂𝜆≠0,_𝑗𝛾₀≠1,𝑗=1,2,3,4,𝜆𝛾₀=𝑁𝛾₀^1/2. (46)

Next, we study the normal form of (10) when 𝛾=𝛾₀.

Use the translation 𝑝𝐹=−1𝑘₂₁−𝑘₁₁𝛽+24²𝑘₁₁+𝑘₂₂−𝛽+24𝑘₁₁𝑘₂₁×𝑥𝑦+𝑝01_∗𝐹^∗; (47) system (10) becomes 𝑥^′𝑦^′=−𝑀𝛾₀2−4−𝑀𝛾₀²24−𝑀𝛾₀²2−𝑀𝛾₀2×𝑥𝑦+𝐻𝐺,𝑥,𝑦𝑥,𝑦 (48) where 𝑚=−(1/𝑘₂₁)−𝑘₁₁[((𝛽+2)/4)²𝑘₁₁+𝑘₂₂],𝑛=−((𝛽+2)/4)(𝑘₁₁/𝑘₂₁), 𝐻=1𝑥,𝑦𝑚𝑓𝑥,𝑦−𝑝^∗−𝑛𝑚𝑔𝑥,𝑦−𝐹^∗+𝑀𝛾₀2𝑥+4−𝑀𝛾₀²2𝐺=𝑔𝑦,𝑥,𝑦𝑥,𝑦−𝐹^∗−4−𝑀𝛾₀²2𝑀𝛾𝑥+⁰2𝑦. (49)

In order for system (48) to undergo Neimark-Sacker bifurcation, we require that the following discriminatory quantity 𝐶₀ is not zero [9]: 𝐶₀()=−Re1−2𝜆𝜆²𝜉1−𝜆¹¹𝜉₂₀−12𝜉₁₁²−𝜉₀₂²+Re𝜆𝜉₂₁, (50) where 𝜉₂₀=18𝐻_𝑥𝑥−𝐻_𝑦𝑦+2𝐺_𝑥𝑦+𝑖𝐺_𝑥𝑥−𝐺_𝑦𝑦−𝐻_𝑥𝑦,𝜉₁₁=14𝐻_𝑥𝑥+𝐻_𝑦𝑦+𝑖𝐺_𝑥𝑥+𝐺_𝑦𝑦,𝜉₀₂=18𝐻_𝑥𝑥−𝐻_𝑦𝑦−2𝐺_𝑥𝑦+𝑖𝐺_𝑥𝑥−𝐺_𝑦𝑦+2𝐻_𝑥𝑦,𝜉₂₁=1𝐻16^𝑥𝑥𝑥+𝐻_𝑥𝑦𝑦+𝐺_𝑥𝑥𝑦+𝐺_𝑦𝑦𝑦+𝑖𝐺_𝑥𝑥𝑥+𝐺_𝑥𝑦𝑦−𝐻_𝑥𝑥𝑦−𝐻_𝑦𝑦𝑦. (51) From previous analysis and the theorem in [9], we have the following theorem.

Proposition 8. Suppose that (𝑝^∗,𝐹^∗) is a polymorphic phenotypic equilibrium of system (10). If 𝑘₁₁<0 and 𝐶₀≠0, then system (10) undergoes a Neimark-Sacker bifurcation at (𝑝^∗,𝐹^∗) for 𝛾₀=(𝛽−2)/(𝑘₁₁𝛽+2𝑘₂₂).

Moreover, if 𝑑=−(1/2)((𝛽/2)𝑘₁₁+𝑘₂₂)<0 and 𝐶₀<0, then an attracting invariant closed curve bifurcates from (𝑝^∗,𝐹^∗) for 𝛾<𝛾₀. If 𝑑=−(1/2)((𝛽/2)𝑘₁₁+𝑘₂₂)<0 and 𝐶₀>0, then a repelling invariant closed curve bifurcates from (𝑝^∗,𝐹^∗) for 𝛾>𝛾₀.

Example 9. Let 𝑢₁₁=0.8, 𝑢₁₂=0.4, 𝑢₂₂=1, 𝛽=0.75, and 𝜉=0.776. We have that (𝑝^∗,𝐹^∗)=(0.8,0.6) is a polymorphic phenotypic equilibrium which satisfies 𝑘₁₁<0. According to 𝐶₀=−3.1284538<0, and 𝑑=−(1/2)((𝛽/2)𝑘₁₁+𝑘₂₂)<0, we have that the invariant closed curve bifurcates from (𝑝^∗,𝐹^∗) for 𝛾<𝛾₀=−65.76178451  being attracting. Let 𝑢₁₁=0.5, 𝑢₁₂=1, 𝑢₂₂=0, 𝜉=0.549287, and 𝛽=38/69. We have that (𝑝^∗,𝐹^∗)=(0.79,0.38) is a polymorphic phenotypic equilibrium which satisfies 𝑘₁₁<0. According to 𝐶₀=2.5988353>0 and 𝑑=−(1/2)((𝛽/2)𝑘₁₁+𝑘₂₂)<0, the invariant closed curve bifurcates from (𝑝^∗,𝐹^∗) for 𝛾>𝛾₀=−34.72413828  being repelling.

5. Conclusion

A discrete population genetics model with partial selfing selection is investigated. We assume that each individual can reproduce by selfing or random outcrossing with probability 𝛽 or 1−𝛽  (0<𝛽<1), respectively, in the population. Some population, such as human population, does not mate at random. So, the partial selfing selection model is reasonable.

In this paper, the conditions for the stability of polymorphic phenotypic equilibria are obtained by the Jury conditions and the center manifold theorem. ESS is not the necessary condition of the stability of the polymorphic phenotypic equilibria. The theoretical analysis and numerical simulations present the existence of stable and unstable period doubling bifurcations, saddle-node bifurcation, and Neimark-Sacker bifurcation. Numerical simulations exhibit more complex dynamical behavior of the genetic system under partial selfing selection.

Acknowledgments

The author is thankful to the referees for helpful suggestions for the improvement of this paper. The research of the first author was supported by the Scientific Research Fund of Sichuan Provincial Education Department under Grant 12ZB082.