Evolutionary Games and Periodic Fitness

One thing that nearly all stability concepts in evolutionary game theory have in common is that they use a time-independent fitness matrix. Although this is a reasonable assumption for mathematical purposes, in many situations in real life it seems to be too restrictive. We present a model of an evolutionary game, driven by replicator dynamics, where the fitness matrix is a variable rather than a constant, i.e., the fitness matrix is time-dependent. In particular, by considering periodically changing fitness matrices, we model seasonal effects in evolutionary games. We discuss a model with a continuously changing fitness matrix as well as a related model in which the changes occur periodically at discrete points in time. A numerical analysis shows stability of the periodic orbits that are observed. Moreover, trajectories leading to these orbits from arbitrary starting points synchronize their motion in time. Several examples are discussed.


Introduction
Evolutionary game theory, started by the work of Maynard Smith and Price [7], studies the dynamic development of populations. Here, a population consists of interacting individuals of finitely many different types. Interactions between different types lead to different fitnesses for these types (e.g., number of offspring). Consequentially, the population distribution, i.e., the relative frequency of each of the types is subject to change over time. In evolutionary models, the rate and the direction of this change are determined by the dynamics. There are several dynamics in use (cf. Hofbauer & Sigmund [6]), the most common one being the replicator dynamics (Taylor & Jonker [10]) and the best response dynamics (Gilboa & Matsui [5]).
One thing that all dynamics in evolutionary games have in common is that they make use of time-independent fitness matrices. Although this is a reasonable assumption for mathe-matical purposes, in many situations in real life it seems to be too restrictive. For instance, if different types require different resources, then a possible environmental effect of having one type abundantly present could be that its resources run low, which in turn would lead to a lower fitness. Also, there may be external effects that influence the fitnesses of the different types in different ways. Little research has been done in this area; one example being the work by Broom [2], who discusses evolutionary games, where the fitness matrix converges to a fixed limit matrix.
In this paper, we present a model of an evolutionary game, driven by replicator dynamics, where the fitness matrix is a variable rather than a constant, i.e., the fitness matrix is timedependent, in the following way: We introduce periodically changing fitness matrices to model seasonal effects in evolutionary games. We present, by means of an example, the model of an evolutionary game with periodic fitnesses, i.e., the fitness matrix continuously changes in time in a periodic fashion. In the example, the population distribution is shown to converge to a periodic orbit. Moreover, trajectories leading to this orbit from arbitrary starting points synchronize their motion in time. Furthermore, we discuss games in which the fitness matrix still changes periodically, but only at discrete points in time. A similar result as for the continuously changing fitness matrix is obtained for an example where only two fitness matrices are used alternatingly over fixed time intervals.

The Model
An evolutionary game is determined by a fitness matrix, based on which a population distribution over different types will change. Here, the population distribution at time t can be described by the vector x(t) = (x 1 (t), x 2 (t), . . . , x n (t)), where x i (t) > 0 for all i (all types are present) and n i=1 x i (t) = 1. The fitness matrix is an n × n matrix A, that is to be interpreted as follows: The entry a ij is the fitness (the payoff, or the number of offspring) for an individual of type i when interacting with an individual of type j . So, given the population distribution x, the average fitness of an individual of type i is equal to e i Ax T and the average fitness of an individual in the population is xAx T . Here, the vector e i is a unit-vector with 1 in position i and 0 elsewhere. The emphasis of the research in these games is on stability. Loosely speaking, a population distribution y is stable if for the process {x(t) : t ≥ 0} we have that, if it ever gets close to y, then it will always stay close to y, or even converge to it. The most commonly used stability concept in evolutionary games is the so-called Evolutionarily Stable Strategy (or ESS) (cf. Maynard Smith and Price [7]). A population distribution (or a strategy) p is an ESS if for all strategies x = p we have Evolutionary stability is a refinement of the well-known Nash-equilibrium (cf. Nash [9]) for symmetric games. Condition (1) says that p should be a best reply against itself, while condition (2) addresses the stability of p. Namely, if x is also a best reply against p, then in order for the population distribution not to drift away in the direction of x, we need that p performs better against x than x against itself. The dynamics that are used most frequently are the replicator dynamics, introduced by Taylor & Jonker [10]. According to the replicator dynamics, the proportion of population members of type i, changes in time according to the following system of differential equations: So the replicator dynamics dictates, in a Darwinian way that the population distribution of those types that perform better than average (or have a higher amount of offspring than average) will grow, while the distribution of types that perform below average, will fall. It can be shown that any ESS p is an asymptotically stable point for the corresponding replicator dynamics, i.e., if x(0) is close enough to p, then the population distribution converges to p.
Both the original ESS-definition and the replicator dynamics assume that the fitness matrix is constant. In this paper, we investigate the replicator dynamics and a type of stability of the population distribution in evolutionary games with a time-dependent fitness matrix A(t). We will call this periodic stability. We do so by means of the following example. This example is based on the idea of having three types, two of which have a fitness that periodically depends on time, sometimes doing very good, sometimes very bad, while for another type the fitness is not directly affected by time at all: The parameters α, σ and ρ each have their own specific impact on the process: (i) σ determines the size or amplitude of the variation. Notice that if σ = 0, then the fitness matrix is time-independent. (ii) ρ determines the time it takes to complete one period. The smaller the value of ρ, the more time the population has to adapt to the changing environment. (iii) α is introduced to control the (time-independent) fitness of the type 2 individuals.
In the next section, we will analyze in detail the process for this example, which is illustrated by Fig. 1. This will eventually form the proof of the following proposition.

Proposition 1
In the evolutionary game corresponding to the time-dependent fitness matrix A(t) from Example 1, for parameter values α = 0.88, ρ = 0.1 and σ = 1, the process {x(t) : t ≥ 0} converges to a periodically stable orbit.
Moreover, trajectories leading to this orbit from arbitrary starting points, synchronize their motion in time. 1

Analysis
In the analysis of Example 1, we first consider populations that do not contain any individuals of type 2. Also, we set the value of σ equal to 1, which means that the top-right and bottom-left entries of the matrix range from 1 to 3. At times t = k · 2π ρ , the fitness matrix only considering types 1 and 3 is and for times t = π ρ + k · 2π ρ we have This implies that at different points in time we are dealing with different fitness matrices, each having their own ESS. Thus, at different points in time, the population will be driven toward different attractors. The three matrices mentioned above correspond to the extremes in terms of the location and/or the average fitness of the corresponding ESS solution. At times corresponding with fitness matrix (3), the unique ESS would be p = ( 3 4 , 1 4 ) with average fitness 3 4 . At times corresponding to fitness matrix (4), the ESS would be p = ( 1 2 , 1 2 ) with average fitness 1 and at times for fitness matrix (5) the ESS is p = ( 1 4 , 3 4 ) with average fitness 3 4 . Consequently, if the initial population distribution of Example 1 does not contain any individuals of type 2, the population distribution will fluctuate (just like the time dependent ESS) between ( 1 4 , 0, 3 4 ) and ( 3 4 , 0, 1 4 ), whereas the average fitness of the ESS solution fluctuates between 3 4 and 1. In fact, at time t it is equal to 1 − 1 4 cos 2 (ρt). If the initial population distribution does contain individuals of type 2, we distinguish between 3 different situations for the values of α.
Case 1: α < 0.75. Then, for any value of t , the unique ESS of the game with fitness matrix A(t) combines types 1 and 3, while type 2 goes extinct. The ESS's oscillate in periods of 2π ρ between p = ( 3 4 , 0, 1 4 ) when t = k · 2π ρ and p = ( 1 4 , 0, 3 4 ) when t = π ρ + k · 2π ρ . Furthermore, at any time t the replicator dynamics drives the population distribution in the direction of the current ESS at time t . Case 2: α > 1. Now, at any time there is a unique ESS, namely e 2 = (0, 1, 0), and indeed the dynamics converge to e 2 . This has to do with the fact that, when α > 1, the fitness of type 2 individuals against any population distributionx that does not contain any individuals of type 2 is e 2 Ax T = α > 1 ≥xAx T . Case 3: 0.75 < α < 1. Now, depending on t , the process will be attracted to a point on the line segment between ( 1 4 , 0, 3 4 ) and ( 3 4 , 0, 1 4 ) or to e 2 . However, some care is needed, for in this case e 2 is no ESS as will be explained below. These time-dependent attractors allow for the possibility for all types to survive and some kind limit cycle to occur. One can observe that at certain time intervals, namely when α > 1 − 1 4 cos 2 (ρt), the type 2 individuals have the upper hand, whereas at the other time intervals types 1 and 3 flourish.
The value of ρ also plays an important role in the overall process. One can see ρ as the rate for the population to adapt to changing fitnesses. If ρ is high, then the population hardly has any time to adapt. In this case, the entire process behaves as if there was just one fitness matrix, the time average one. At the opposite, if ρ is small, the population has lots of time to adapt. In the extreme case where ρ approaches zero, some types may go extinct before their fitness values recover. The value ρ = 0.1 turns out to be good to observe a cyclic pattern.
Let us now examine the case with α = 0.88, ρ = 0.1 and σ = 1 in more detail. Note that the fitness matrix A(t) changes periodically with a period of 20π . During each period the process will be attracted to the line Nevertheless, a classical evolutionary game played with this fitness matrix would show a process converging to e 2 , because for this matrix a trajectory starting anywhere in the interior of the simplex, will eventually reach the line between ( 3 4 , 0, 1 4 ) and e 2 , and next it will move along this line toward e 2 . The fact that populations that consist almost exclusively of type 2 individuals can be invaded by mixtures of types 1 and 3, can be illustrated by examining the so called force field of the replicator dynamics for this matrix A(0). In Fig. 2a, we see that in the interior of the population space the directional derivatives seem to be pointing at e 2 . In Fig. 2c, we take a closer look at the corner near e 2 , where we have stretched up this corner point to make the lines [e 3 , e 2 ] and [e 1 , e 2 ] parallel on the left and right, respectively. We have done so by switching from coordinates (x, y, z) to ( x−z 1−y , y) for y < 1. Next, in Fig. 2d, we have normalized the vertical components of all the arrows and neglected the horizontal parts. This shows that in any neighborhood of e 2 there are points for which the process would drift away from e 2 . Yet, the point e 2 is asymptotically stable because when the process would drift away from e 2 then it will slowly move in the direction of the line [e 3 , e 2 ] until it reaches the line segment [( 3 4 , 0, 1 4 ), e 2 ], as is indicated in Fig. 2b; from that moment onward the process will continue to e 2 directly.

Continuously Changing Fitness Matrix
When we examine the process obtained by the replicator dynamics on the continuously changing fitness matrix A(t), then we observe that no matter where we start in the interior of the simplex, the process always converges to the same cyclic trajectory. This is illustrated in Fig. 1. Again the parameters used are: α = 0.88, ρ = 0.1, and σ = 1.
In order to show that the observed periodic orbit is indeed periodically stable we used the tool ARIADNE (cf. Collins et al. [3]) that was developed for the analysis of dynamic systems using rigorous numerical methods. (See the Appendix for an overview of rigorous numerics and the methods used for the analysis.) We considered the replicator dynamics for the specific parameter values given, which has a forcing period of T = 2π ρ . An approximation to the time-T return map r over an initial domain D was computed, along with an error bound in the uniform norm. The interval Newton operator was then used to prove the existence of a fixed-point of r, yielding a period-T orbit of the replicator dynamics.
One of the eigenvalues of the derivative matrix of the return map was computed to be approximately 0.95, whereas the other two were calculated to be of order 10 −15 , which is comparable to the machine epsilon. This validates the observation made during the simulations, that in the region of interest, there is very fast convergence to an invariant manifold, followed by slow convergence to the fixed-point of the return map. The eigenvalue 0.95 of the derivative matrix indicates that a perturbation of the system changes the position of the

Changes in Fitness at Periodic Discrete Points in Time
After having observed these periodically stable population distributions based on the continuously changing fitness matrix, we checked whether cyclic population distributions can also be obtained by discretely changing fitness matrices. This question was answered affirmatively for the model where A(0) is used as fitness matrix for t ∈ [2k π ρ , (2k + 1) π ρ ], while A( π ρ ) is used for t ∈ [(2k + 1) π ρ , (2k + 2) π ρ ], k = 0, 1, 2, 3, . . . , which means that the fitness matrix changes only twice per period at fixed moments in time. As we can see in Fig. 3, this still leads to a cyclic pattern. This is quite remarkable, because when looking at the replicator dynamics for each of the fitness matrices separately, we find that each of those processes has a unique attractor, namely e 2 . However, when these processes are combined as one, then this alternating process never reaches e 2 , because the trajectories leading to e 2 are along different sides of the simplex for the separate processes. Also, in this case, no matter what population distribution we start with, the process will always lead to the unique periodically stable orbit. Obviously, this observation depends very much on the two fitness matrices and the specific values of the parameters.

Other Patterns with Periodic Fitness Matrices
We have also explored another periodic fitness matrix. We used the following fitness matrix, based on the rock-paper-scissors game, and examined the processes that arise for several choices of values for the parameters: Figure 4 gives the periodically stable orbits for several settings of the parameters ρ 1 , ρ 2 , and ρ 3 . It should be stressed that in each of these cases the orbits will be reached from any initial population distribution.
We also observed that processes starting at different points of the simplex quickly start following synchronized patterns. This can be seen in Fig. 5, where this is illustrated for the periodically stable orbit of Fig. 4c.

Conclusions
In conclusion, we have discovered that it is possible to reach a periodically stable population pattern when following the replicator dynamics based on a periodically changing fitness matrix. Because the fitness matrix changes continuously in time, the force field of the dynamics that makes the population move into the direction of a particular population distribution, changes with it. In our stepwise periodic model, we have sharp changes of the fitness matrix at fixed times, and the number of attraction points is finite. In fact, for the particular case examined, there is just one attraction point all the time, but the combined process never reaches it.
We want to remark that we could only employ the Ariadne tool for showing stability of the orbits observed, because the fitness matrix was changing in a deterministic way. Such would not have been possible when it was changing in stochastic way, because then there would not necessarily be any closed loop.
A challenging question is to predict the existence of periodically stable orbits based on the periodic fitness matrix. A second question is to find closed form expressions for the periodic orbits observed.
Movies for the dynamic processes discussed can be observed at http://www.youtube.com/ watch?v=C65Z7fcLA4s.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.

Appendix: Rigorous Numerics
Rigorous numerical methods are a class of tools in which numerical approximations to exact mathematical objects are computed with bounds on the approximation errors, allowing for rigorous conclusions to be made about the calculated quantities. In order to control roundoff error, use is typically made of the round-up and round-down modes of floating-point arithmetic built in to modern microprocessors, and specialized algorithms are used to handle truncation errors.
A fundamental tool of rigorous numerics is interval analysis (cf. Moore [8]), in which upper and lower bounds on the value of a real quantity are specified. In other words, a computation of a real quantity x returns an intervalx = [x, x] containing x. Typically, the values x and x are machine-precision floating-point numbers. Interval arithmetic can be implemented using rounded floating-point arithmetic, for example, addition is defined as Similarly, we can describe a function f : R n → R using an approximation g over a bounded domain D with a uniform error bound . In ARIADNE [3], we represent a function using a polynomial function model, which is a tuple D, s, p, e where D = [a 1 , b 1 ] × · · · × [a n , b n ] is the domain of the approximation, s is an affine scaling function with s i : [−1, +1] → [a i , b i ], p is a polynomial p(x) = α∈N n c α x α 1 1 · · · x αn n , and e is a uniform error bound. The values a i , b i , c α and e are machine-precision floatingpoint numbers. A polynomial modelf = D, s, p, e represents a function f on D if sup x∈D |f (x) − p(s −1 (x))| ≤ e. The scaling is useful to improve the conditioning of the representing polynomial. This representation is essentially the same as the Taylor model representation of Berz and Makino [1].
The flow φ of a differential equationẋ = f (x) is defined by ∂φ ∂t (x, t) = f φ(x, t), t and φ(x, 0) = x In order to rigorously prove the existence of a periodic orbit by numerical methods, it is necessary to be able to compute the function φ(x, T ) in a rigorous way, and to rigorously validate the existence of a solution of φ(x, T ) − x = 0. A validated approximation to the function φ(x, T ) can be computed in ARIADNE by evaluating flow stepsφ(x, h) using Picard iteration on function models: Repeated iteration of the Picard operator must refine the approximationφ n to φ if the time step h is sufficiently small. An approximationr to the time-T return map r(x) = φ(x, T ) over an initial domain D is computed by composing flow steps.
In order to prove existence of a fixed-point of the time-T return map, a modified version of the interval Newton operator (cf. Moore [8]) is provided. Suppose r(x) ∈ p(x) + E for x ∈ X, where p is differentiable and E is a fixed box, and define N(X, x, p, E) = x − Dp(X) −1 where [Dp(X)] is an interval matrix containing {Dp(x) | x ∈ X} and x is a fixed element of x. It can be shown, using Theorem 6.1 from Collins et al. [4], that if N(X, x, p, E) ⊂ X, then r has a zero in X, and any zero of r in X must lie in N(X, x, p, E). Applying the interval Newton operator tor(x) − x allows a proof of the existence of a periodic orbit.