Electronic Journal of Qualitative Theory of Differential Equations Proc. 8th Coll. QTDE, 2008, No. 6 1-10; http://www.math.u-szeged.hu/ejqtde/

NEIMARK-SACKER BIFURCATION IN A DISCRETE DYNAMICAL MODEL OF POPULATION GENETICS

ATTILA DÉNES

Abstract. In this paper we investigate a genetic model initiated by Gábor Tusnády. This model is a four-dimensional system of difference equations which describes the change of distribution of genotypes from generation to generation in the case of one locus and four alleles considering selection and mutation. During computer experiments Tusnády discovered cyclic behaviour in the evolution of genotypes. Using bifurcation theory we prove that the system indeed can have cycles and this occurrence is caused by a Neimark-Sacker bifurcation.

1. Introduction

G. Tusnády established a discrete population dynamical model which describes the change of distribution of genotypes from generation to generation considering selection, mutation and recombination.

The genetical program is stored in the cells of living beings in the form of chromosomes. One half of the chromosomes comes from the mother, the other half comes from the father, i. e., chromosomes appear in pairs. When such a diploid cell divides, each chromosome doubles and the two arising cells get the whole chromosome set. Our organism also contains haploid cells, in which only half as much chromosomes can be found: these are the gametes. These cells originate from diploid cells during meiosis which splits chromosome pairs. During fertilization gametes unite and the original chromosome number is re-established. The segments of the chromosomes that determine the different properties are called genes, their different variants are called alleles and their place in the chromosomes are called loci. The genotype is determined by the two alleles which are really present in the cell. The distribution of the genoytpes is affected mainly by selection, mutation and recombination. Selection means that different genotypes have different chances to create offspring. The succesfulness of a genotype is shown by the fitness. When a diploid cell divides and chromosomes double, their copying is not always absolutely precise: certain segments of the chromosome can change, an allele can change to another one. This phenomenon is called mutation. (For details see e. g. [1].)

In Section 2 we introduce Tusnády’s model in the case of one locus and four alleles. Using monograph [2], in Section 3 we delineate the Neimark-Sacker bifurcation and its conditions (Theorem A). In the proof of our main result (Theorem 4.1) we have to use a special procedure to check the genericity conditions of Theorem A. To make our paper self-contained we include this procedure from [2]. In Section 4 we apply the procedure presented in Section 3 to obtain our main result.

EJQTDE, Proc. 8th Coll. QTDE, 2008 No. 6, p. 1

This paper is in final form and no version of it is submitted for publication elsewhere.

2. The model

Let us model the above-mentioned notions mathematically in the case of one locus. Let n denote the number of alleles. Let x_{i} denote the density of the i-th allele in the given generation, x_{1} + ... + x_{n} = 1. Let Γ denote the set of genotypes. Then a γ ∈ Γ genotype is a vector with two independent allele components. Fitness is a mapping from Γ to the [0, 1] closed interval. Meiosis is a random mapping from Γ into the set of alleles.

Tusnády assumed that mutation depends on both parental genes and created the following model:

Let w(i, j) denote the fitness of the genotype ij, i. e., the probability of genotype ij creating offspring and {y_{ij}}_{i,j=1}^{n} denotes the distribution of genotypes in the next generation. If we consider selection only, we get the following selection model:

(2.1)  y_{ij} = \frac{w(i, j)x_{i}x_{j}}{∑_{p,q=1}^{n} w(p, q)x_{p}x_{q}} .

Let us now consider mutation as well. Let M_{ij}(k) denote the probability of the event that a gamete of type k issues from a cell of genotype ij. M_{ij}(k) contains mutation as well. Then, after meiosis, the new distribution of gametes is given by the following formula:

x′_{k} = ∑_{i,j=1}^{n} y_{ij}M_{ij}(k).

If we substitute the expression from the selection model (2.1) in the place of y_{ij} introducing the notation a(i, j, k) = w(i, j)M_{ij}(k), we get the following model:

(2.2)  x′_{k} = \frac{∑_{i,j=1}^{n} a(i, j, k)x_{i}x_{j}}{∑_{i,j,k=1}^{n} a(i, j, k)x_{i}x_{j}} .

Tusnády investigated what can be said about the limit sets of the solutions of (2.2) obtained by iterating the mapping and whether there is a system in which there are solutions having more than one limit points. After a long numeric searching he found the following four-dimensional system:

a(2,4,1)=1042 a(2,4,2)=8
a(3,4,2)=113 a(1,2,3)=19
a(2,3,3)=9  a(1,3,4)=1078,

a(2,2,4)=414

where a(i, j, k) = a(j, i, k) and the coefficients not mentioned here are equal to zero. That is, our system is as follows:

(2.3)  F (x) = ( \frac{2084x_{2}x_{4}}{38x_{1}x_{2}+414x_{2}^{2}+2156x_{1}x_{3}+18x_{2}x_{3}+2100x_{2}x_{4}+226x_{3}x_{4}} , ; \frac{16x_{2}x_{4}+226x_{3}x_{4}}{38x_{1}x_{2}+414x_{2}^{2}+2156x_{1}x_{3}+18x_{2}x_{3}+2100x_{2}x_{4}+226x_{3}x_{4}} , ; \frac{38x_{1}x_{2}+18x_{2}x_{3}}{38x_{1}x_{2}+414x_{2}^{2}+2156x_{1}x_{3}+18x_{2}x_{3}+2100x_{2}x_{4}+226x_{3}x_{4}} , ; \frac{414x_{2}^{2}+2156x_{1}x_{3}}{38x_{1}x_{2}+414x_{2}^{2}+2156x_{1}x_{3}+18x_{2}x_{3}+2100x_{2}x_{4}+226x_{3}x_{4}} ) .

The limit set of this system is not one-dimensional. Probably its dimension is two and the behaviour of this system seems to be chaotic even for an ordinary observer.

EJQTDE, Proc. 8th Coll. QTDE, 2008 No. 6, p. 2

Figure 1. p = 8

Using bifurcation theory we prove that the appearance of the one-dimensional limit is not caused by errors of numeric calculating.

Let us examine the attractor of our system at different parameter values. We chose the coefficient p = a(2, 4, 2) = a(4, 2, 2) as a parameter. In Figure 1 the attractor of the system can be seen with the coefficients determined by Tusnády.

Figure 2 represents the attractor at p = 135, where the attractor is a stable closed curve, while in Figure 3, at p = 145 the attractor is a stable fixed point. The figures suggest that between the last two parameter values a Neimark-Sacker bifurcation takes place, as we can see a stable closed curve arising from the stable fixed point.

In Figure 4 we can see the whole dynamics of the system in the phase space. As the sum of the densities of the alleles is equal to one, we can reduce the system to three dimensions. The phase space is the four-dimensional simplex, i. e., a tetrahedron. When we change the value of the parameter p, a pair of complex eigenvalues of the Jacobian passes through the unit circle. To this complex pair corresponds a two-dimensional unstable manifold of the fixed point. The invariant closed curve appears on this manifold. The fixed point has a stable manifold as well; the unstable manifold attracts the solutions. At parameter value p = 139.455 the two eigenvalues have absolute value 1. The system has to satisfy some genericity conditions to assure the occurrence of the bifurcation.

EJQTDE, Proc. 8th Coll. QTDE, 2008 No. 6, p. 3 Figure 2. p = 135 Figure 3. p = 145 EJQTDE, Proc. 8th Coll. QTDE, 2008 No. 6, p. 4

Figure 4. Dynamics of the system

3. Neimark-Sacker bifurcation (see [2, Sections 4.7, 5.4])

If we change parameter p = a(2, 4, 2) = a(4, 2, 2), a complex pair of eigenvalues of the Jacobian of the system passes through the unit circle. However, the system has to satisfy some non-degenericity conditions to ensure a Neimark-Sacker bifurcation. First we formulate the Neimark-Sacker bifurcation theorem for two-dimensional systems.

Let us consider the following two-dimensional system:

x ↦ f(x, α), x ∈ ℝ^{2}, α ∈ ℝ^{1},

where the smooth function f has at α = 0 the fixed point x = 0 with simple eigenvalues μ_{1,2} = e^{±iθ_{0}}, 0 < θ_{0} < π. By the Implicit Function Theorem, the system has a unique fixed point x_{0}(α) for all sufficiently small |α| in some neighbourhood of the origin, since μ = 1 is not an eigenvalue of the Jacobian matrix. With a parameter-dependent coordinate shift we can place the fixed point at the origin, therefore we can assume that x = 0 is the fixed point for |α| sufficiently small. Thus, the system can be written as

(3.1)  x ↦ A(α)x + F (x, α),

where F is a smooth vector function and its components F_{1,2} have Taylor expansions in x starting with at least quadratic terms, F (0, α) = 0 for all sufficiently small |α|. The Jacobian matrix A(α) has two multipliers

μ_{1,2} = r(α)e^{±iφ(α)},

where r(0) = 1, φ(0) = θ_{0}. Thus r(α) = 1 + β(α) for some smooth β(α) function, where β(0) = 0. Suppose that β′(0) ≠ 0. Then we can use β as a new parameter

EJQTDE, Proc. 8th Coll. QTDE, 2008 No. 6, p. 5 and express the multipliers in terms of β: μ_{1}(β) = μ(β), μ_{2}(β) = \overline{μ}(β), where μ(β) = (1 + β)e^{θ(β)} and θ(β) is a smooth function such that θ(0) = θ_{0}.

To formulate the theorem we need another transformation. Omitting the technical details it is enough to say that introducing a complex variable and a new parameter system (3.1) can be transformed for all sufficiently small |α| into the following form:

z ↦ μ(β)z + g(z, \overline{z}, β),

where β ∈ ℝ^{1}, z ∈ ℂ^{1}, μ(β) = (1 + β)e^{iθ(β)} and g is a smooth complex-valued function of z, \overline{z} and β, whose Taylor expansion with respect to (z, \overline{z}) contains quadratic and higher-order terms:

g(z, \overline{z}, β) = ∑_{k+l≥2} \frac{1}{k!l!}g_{kl}(β)z^{k}\overline{z}^{l},

k, l = 0, 1, ….

Now we can state the following theorem for Neimark-Sacker bifurcation in two dimensions:

Theorem A (Generic Neimark-Sacker bifurcation ([2])). For any generic twodimensional one-parameter system

x ↦ f(x, α),

having at α = 0 the fixed point x_{0} = 0 with complex multipliers μ_{1,2} = e^{±iθ_{0}} there is a neighbourhood of x_{0} in which a unique closed invariant curve bifurcates from x_{0} as α passes through zero.

The system has to satisfy the following genericity conditions:

(1) r′(0) ≠ 0, where μ_{1,2}(α) = r(α)e^{±iφ(α)}, r(0) = 1, φ(0) = θ_{0},

(2) e^{±ikθ_{0}} ≠ 1, k = 1, 2, 3, 4,

(3) a(0) ≠ 0, where a(0) = Re (\frac{e^{−iθ_{0}} g_{21}}{2} )−Re (\frac{(1−2e^{iθ_{0}} )e^{−2iθ_{0}}}{2(1−e^{iθ_{0}} )} g_{20}g_{11})−\frac{1}{2}|g_{11}|^{2}−

\frac{1}{4}|g_{02}|^{2}, where g_{ij} = g_{ij}(0).

Coefficient a(0) in condition (3) determines the direction of the appearance of the invariant curve in a generic system exhibiting Neimark-Sacker bifurcation. If a(0) is negative, the bifurcation is supercritical, i. e., a stable closed invariant curve bifurcates from a stable fixed point while the fixed point becomes unstable. If a(0) is positive, the bifurcation is subcritical, i. e., an unstable closed curve disappears as we pass through the critical value.

As for systems with dimension higher than 2 essentially the same takes place: there exists a two-dimensional invariant manifold on which the system exhibits the bifurcation, while the behaviour off the manifold is “trivial”.

During the proof of Theorem 2.1 we followed the procedure described in [2]. Using this method we avoid the transformation of the system into its eigenbasis. We use the eigenvectors belonging to the critical eigenvalues of the Jacobian A of the system and of its transpose A^{T} to “project” the system into the critical eigenspace and its complement. We used Mathematica to perform the necessary calculations.

EJQTDE, Proc. 8th Coll. QTDE, 2008 No. 6, p. 6

We write the system in the form

(3.2)  G(x) = Ax + H(x), x ∈ ℝ,

where H(x) = O(∥x∥^{2}) is a smooth function with Taylor expansion

H(x) = \frac{1}{2}B(x, x) + \frac{1}{6}C(x, x, x) + O(∥x∥^{4}),

where B(x, x) and C(x, x, x) are multilinear functions. In coordinates:

(3.3)  B_{i}(x, y) = ∑_{j,k=1}^{n} \frac{∂^{2}H_{i}(ξ)}{∂ξ_{j}∂ξ_{k}} |_{ξ=0} x_{j}y_{k}

and

(3.4)  C_{i}(x, y, z) = ∑_{j,k,l=1}^{n} \frac{∂^{3}H_{i}(ξ)}{∂ξ_{j}∂ξ_{k}∂ξ_{l}} |_{ξ=0} x_{j}y_{k}z_{l},

where i = 1, 2, …, n.

A has a simple pair of complex eigenvalues on the unit circle: μ_{1,2} = e^{±iθ_{0}} Let q ∈ ℂ^{n} denote the complex eigenvector corresponding to μ_{1}:

(3.5)  Aq = e^{iθ_{0}}q, A\overline{q} = e^{−iθ_{0}} \overline{q}.

Introduce also the adjoint eigenvector p ∈ ℂ^{n} which has the property

(3.6)  A^{T} p = e^{−iθ_{0}}p, A^{T} \overline{p} = e^{iθ_{0}} \overline{q}.

and satisfies the normalization ⟨p, q⟩ = 1.

The critical real eigenspace T^{c} corresponding to μ_{1,2} is two-dimensional and is spanned by Re q, Im q. The real eigenspace T^{su} corresponding to all other eigenvalues of A is (n − 2)-dimensional. y ∈ T^{su} if and only if ⟨p, y⟩ = 0. As y ∈ ℝ^{n} is real and p ∈ ℂ^{n} is complex, the condition ⟨p, y⟩ = 0 implies two real constraints on y. We decompose x ∈ ℝ^{n} as

x = zq + \overline{z}\overline{q} + y,

where z ∈ ℂ^{1} and zq + \overline{z}\overline{q} ∈ T^{c}, y ∈ T^{s}u.

We have

\{ \begin{matrix}z & = ⟨p, x⟩ \\ y & = x − ⟨p, x⟩q − ⟨\overline{p}, x⟩\overline{q}. \end{matrix}

After a series of transformations we get the following form for the map:

(3.7)  z̃ = e^{iθ_{0}}z + \frac{1}{2}g_{20}z^{2} + g_{11}z\overline{z} + \frac{1}{2}g_{02}\overline{z}^{2} + \frac{1}{2}g_{21}z^{2}\overline{z} + ⋯,

where

(3.8)  g_{20} = ⟨p, B(q, q)⟩, g_{11} = ⟨p, B(q, \overline{q})⟩, g_{02} = ⟨p, B(\overline{q}, \overline{q})⟩

EJQTDE, Proc. 8th Coll. QTDE, 2008 No. 6, p. 7 and

g_{21} = ⟨p, C(q, q, \overline{q})
+ 2⟨p, B(q, (E − A)^{−1}B(q, \overline{q})))⟩ + ⟨p, B(\overline{q}, (e^{2iθ_{0}}E − A)^{−1}B(q, q))⟩
+ \frac{e^{−iθ_{0}}(1 − 2e^{iθ_{0}})}{1 − e^{iθ_{0}}}⟨p, B(q, q)⟩⟨p, B(q, \overline{q})⟩
− \frac{2}{1 − e^{−iθ_{0}}} |⟨p, B(q, \overline{q})⟩|^{2} − \frac{e^{iθ_{0}}}{e^{3iθ_{0}} − 1}|⟨p, B(\overline{q}, \overline{q})⟩|^{2}.

If e^{ikθ_{0}} ≠ 1 for k = 1, 2, 3, 4, we can transform (3.7) into the following form:

z̃ = e^{iθ_{0}z}(1 + d(0)|z|^{2}) + O(|z|^{4}).

Here the real number a(0) = Re d(0) determines the direction of bifurcation of a closed invariant curve, and can be obtained from the following formula:

a(0) = Re (\frac{e^{−iθ_{0}}g_{21}}{2} ) − Re (\frac{(1 − 2e^{iθ_{0}})e^{−2iθ_{0}}}{2(1 − e^{iθ_{0}})} g_{20}g_{11})
− \frac{1}{2}|g_{11}|^{2} − \frac{1}{4}|g_{02}|^{2}.

Using this formula with the above-definded coefficients, we obtain the following invariant expression:

\begin{aligned}a(0) & = \frac{1}{2}Re \{e^{−iθ_{0}}[⟨p, C(q, q, \overline{q})⟩ + 2⟨p, B(q, (E − A)^{−1}B(q, \overline{q}))⟩ \\ & + ⟨p, B(\overline{q}, (e^{2iθ_{0}}E − A)^{−1}B(q, q))⟩]\}. \end{aligned} ; (3.9)

This formula allows us to verify the nondegeneracy of the nonlinear terms at a nonresonant Neimark-Sacker bifurcation of n-dimensional maps with n ≥ 2.

4. Main result

Theorem 4.1. Tusnády’s system (2.3) undergoes a Neimark-Sacker bifurcation at parameter value 139.455, i. e., a stable invariant closed curve bifurcates from a stable fixed point while the fixed point becomes unstable.

Proof:

To check condition (1) of Theorem A we calculated the value of r′(0) numerically and obtained r′(0) = −0.00217713. (Actually, to check the transversality directly, we draw the graph of the interpolating function to r, see Figure 5.)

The critical eigenvalues are 0.561391 + 0.827552i and 0.561391 − 0.827552i, i. e., they are not complex roots of unity of order four or less and condition (2) of Theorem A is satisfied.

We calculate the value of a(0) in (3.9) for system (2.3). First we use the transformation G(x) = F (a − x) − a to place the fixed point at the origin.

EJQTDE, Proc. 8th Coll. QTDE, 2008 No. 6, p. 8

Figure 5. Interpolating function to r

We determine the eigenvector q and the adjoint eigenvector p introduced in (3.5) and (3.6):

q = (− 0.67126 − 0.0108908i, −0.0793651 − 0.049971i,
0.0181909 + 0.0608618i, 0.732434 + i)

p = ( 0.0134851 − 0.0321706i, −0.19597 − 0.00610824i,
0.97727, −0.0340428 + 0.0642307i)

Now we have to compute the functions B and C introduced in (3.3) and (3.4). Here we only show B_{1} and C_{1}; the rest of the functions look similar to these.

B_{1}(x, y) = ( \begin{matrix}0.281071x_{1}y_{1} + 1.3569x_{2}y_{1} − 4.99837x_{3}y_{1} + 0.337012x_{4}y_{1}+ \\ 1.3569x_{1}y_{2} − 30.0461x_{2}y_{2} + 42.9415x_{3}y_{2} − 1.51484x_{4}y_{2}− \\ 4.99837x_{1}y_{3} + 42.9415x_{2}y_{3} + 208.644x_{3}y_{3} + 7.85539x_{4}y_{3}+ \\ 0.337012x_{1}y_{4} − 1.51484x_{2}y_{4} + 7.85539x_{3}y_{4} − 2.8086x_{4}y_{4} \end{matrix} )

C_{1}(x, y, z) =

( \begin{matrix}0.413396x_{1}y_{1}z_{1} + 3.40342x_{2}y_{1}z_{1} − 13.5565x_{3}y_{1}z_{1} + 0.961551x_{4}y_{1}z_{1}+ \\ 3.40342x_{1}y_{2}z_{1} + 6.19157x_{2}y_{2}z_{1} + 28.8248x_{3}y_{2}z_{1} − 1.50999x_{4}y_{2}z_{1}− \\ 13.5565x_{1}y_{3}z_{1} + 28.8248x_{2}y_{3}z_{1} − 369.354x_{3}y_{3}z_{1} + 10.0174x_{4}y_{3}z_{1}+ \\ 0.961551x_{1}y_{4}z_{1} − 1.50999x_{2}y_{4}z_{1} + 10.0174x_{3}y_{4}z_{1} + 0.136465x_{4}y_{4}z_{1}+ \\ 3.40342x_{1}y_{1}z_{2} + 6.19157x_{2}y_{1}z_{2} + 28.8248x_{3}y_{1}z_{2} − 1.50999x_{4}y_{1}z_{2}+ \\ 6.19157x_{1}y_{2}z_{2} − 696.409x_{2}y_{2}z_{2} + 234.526x_{3}y_{2}z_{2} − 13.6853x_{4}y_{2}z_{2}+ \\ 28.8248x_{1}y_{3}z_{2} + 234.526x_{2}y_{3}z_{2} + 2845.51x_{3}y_{3}z_{2} − 41.4316x_{4}y_{3}z_{2}− \\ 1.50999x_{1}y_{4}z_{2} − 13.6853x_{2}y_{4}z_{2} − 41.4316x_{3}y_{4}z_{2} − 0.0275582x_{4}y_{4}z_{2}− \\ 13.5565x_{1}y_{1}z_{3} + 28.8248x_{2}y_{1}z_{3} − 369.354x_{3}y_{1}z_{3} + 10.0174x_{4}y_{1}z_{3}+ \\ 28.8248x_{1}y_{2}z_{3} + 234.526x_{2}y_{2}z_{3} + 2845.51x_{3}y_{2}z_{3} − 41.4316x_{4}y_{2}z_{3}− \\ 369.354x_{1}y_{3}z_{3} + 2845.51x_{2}y_{3}z_{3} + 8360.91x_{3}y_{3}z_{3} + 642.894x_{4}y_{3}z_{3}+ \\ 10.0174x_{1}y_{4}z_{3} − 41.4316x_{2}y_{4}z_{3} + 642.894x_{3}y_{4}z_{3} + 0.598556x_{4}y_{4}z_{3}+ \\ 0.961551x_{1}y_{1}z_{4} − 1.50999x_{2}y_{1}z_{4} + 10.0174x_{3}y_{1}z_{4} + 0.136465x_{4}y_{1}z_{4}− \\ 1.50999x_{1}y_{2}z_{4} − 13.6853x_{2}y_{2}z_{4} − 41.4316x_{3}y_{2}z_{4} − 0.0275582x_{4}y_{2}z_{4}+ \\ 10.0174x_{1}y_{3}z_{4} − 41.4316x_{2}y_{3}z_{4} + 642.894x_{3}y_{3}z_{4} + 0.598556x_{4}y_{3}z_{4}+ \\ 0.136465x_{1}y_{4}z_{4} − 0.0275582x_{2}y_{4}z_{4} + 0.598556x_{3}y_{4}z_{4} − 18.9188x_{4}y_{4}z_{4} \end{matrix} )

With the help of these functions we can calculate the values of the scalar products in the formula for a(0).

For the value of a(0) at parameter value 139.455 we get −13.9966. As a(0) ≠ 0, the system is non-degeneric, i. e., there is a Neimark-Sacker bifurcation taking place at parameter value 139.455. As a(0) is negative, the bifurcation is supercritical, i. e., a stable closed invariant curve bifurcates from the fixed point while the fixed point becomes unstable.

5. Acknowledgement

The author is grateful for the suggestions and comments by the referee that helped to improve the paper.

References

[1] Hofbauer, S., Sigmund, K., Evolutionary Games and Population Dynamics, Cambridge University Press, 1998.

[2] Kuznetsov, Y. A., Elements of applied bifurcation theory, Springer-Verlag, 1998.

[3] Tusnády, G., Mutáció és szelekció, Magyar Tudomány 7 (1997), 792–805.

[4] Wiggins, S., Global Bifurcations and Chaos, Springer-Verlag, New York, 1988.

(Received September 13, 2007)

Bolyai Institute, University of Szeged, Aradi vértanúk tere 1, H-6720 Szeged, Hungary