EvoDyn-3s : A Mathematica computable document to analyze evolutionary dynamics in 3-strategy games

EvoDyn-3s generates phase portraits of evolutionary dynamics, as well as data for the analysis of their equilibria. The considered evolutionary dynamics are ordinary differential equations based on adaptive processes taking place in a population of players who are randomly and repeatedly matched in couples to playa2-playersymmetricnormal-formgamewiththreestrategies. EvoDyn-3s calculatestherestpointsof thedynamicsusingexactarithmetic,andrepresentsthem.ItalsoprovidestheeigenvaluesoftheJacobian of the dynamics at the isolated rest points, which are useful to evaluate their local stability. The user only needs to specify the 3 × 3 payoff matrix of the game and choose the dynamics. © 2018 The Authors. Published by Elsevier B


Motivation and significance
Social and biological interactions among agents who may adopt different actions are usually modeled as games.Frequently, these games are studied from an adaptive or evolutionary perspective, leading to systems of ordinary differential equations known as evolutionary game dynamics.A paradigmatic example is the replicator dynamics [1][2][3][4], which have become a standard reference case when analyzing adaptive processes.
In general, many of the basic properties of evolutionary game dynamics can be illustrated using games with three strategies [5,6], and applications to the evolution of cooperation are often studied in this framework [7].To understand and analyze evolutionary game dynamics with three strategies, one of the most useful and intuitive tools are phase portraits, which are geometric representations of the trajectories of the dynamical system in the phase plane.Generating this type of graph typically requires expensive software or customized programming, and this can present a considerable barrier for many researchers studying evolutionary processes.EvoDyn-3s has been designed to overcome this barrier.
Specifically, EvoDyn-3s provides high-quality print-ready vector phase portraits for a diverse group of evolutionary game dynamics with three strategies.Using it does not require any programming and most of its functionality can be run with the free Wolfram CDF Player. 1 EvoDyn-3s also calculates and presents the equilibria of the selected dynamics and performs an eigenvalue analysis of the linearized dynamics using exact arithmetic, a feature that is not available in other programs and can be very useful for theoretical analysis.
The software most closely related to EvoDyn-3s is Dynamo [8,9].Dynamo is also open-source software that runs on Mathematica, and has also been designed to create phase diagrams and other images related to dynamical systems from evolutionary game theory.Dynamo can be used to generate graphs for single-population games with 3 or 4 strategies and for some multipopulation games.Dynamo is more flexible and general than EvoDyn-3s, but significantly less user-friendly, and it uses numerical approximations rather than exact arithmetic.Another software somewhat related to EvoDyn-3s is PDToolbox [10], which is a set of functions coded in Matlab for analyzing some evolutionary dynamics, as well as finite-population agent-based models related to those dynamics.The Python package egtplot [11] creates phase diagrams for the replicator dynamics.While narrower in scope and less userfriendly than EvoDyn-3s, egtplot has the commendable feature of running on an open-source platform.Lastly, ABED [12] (Agent-Based Evolutionary Dynamics) is also free and completely opensource software for simulating adaptive processes, but in finite populations.It provides a complementary approach to the analysis of evolutionary dynamics followed in EvoDyn-3s, in the sense that many of the adaptive processes considered and implemented in ABED (following an agent-based approach) can be approximated -for sufficiently large populations -by differential equations corresponding to the evolutionary dynamics implemented in EvoDyn-3s.

Software description
EvoDyn-3s is a computable document written in Mathematica language.The document contains the executable program, detailed instructions on how to use it, and the source code.There is no need to compile the code.The program can be used directly, by simply opening the computable document with the free Wolfram CDF Player or with Mathematica.
Fig. 1 shows the interface of EvoDyn-3s.The left part contains a series of input boxes and various controls that are used to set the values of all parameters.The right part shows the main output of the program: a phase portrait in the 2-dimensional simplex, a table showing all the isolated rest points and the eigenvalues of the Jacobian of the chosen dynamic at each of the isolated rest points (if the Jacobian is defined), and another table showing the components of rest points (if any exist).
The effect of changing the value of any parameter (except the payoff matrix) on the phase portrait and on the computation of rest points and eigenvalues is immediate, i.e. there is no need to compile or rerun the program.2Thus, for example, the user can gradually move any parameter slider and immediately appreciate how this change affects the output of the program.

Software architecture
EvoDyn-3s conducts the following high-level operations, which are sketched in Fig. 2: • Creates the system of differential equations using the following input provided by the user: payoff matrix, baseline dynamic, probability of random strategy µ, and -only for dynamics Logit and Single-match imitative logit -parameter η.The generated system of differential equations is of the form (ẋ 1 , ẋ2 , ẋ3 ) = f (x 1 , x 2 , x 3 ), where x i represents the fraction of the population using strategy i.
• Solves the system of differential equations numerically for various initial conditions and represents the solutions.The number of initial conditions and the length of the computed trajectories can be controlled by the user.
• Represents a series of orbits and arrows showing the direction of movement in the 2-dimensional simplex • Colors the background of the simplex according to the speed of the dynamic, using the color gradient selected by the user.
• Computes the rest points of the dynamic using exact arithmetic.
• Computes numerical approximations to the isolated rest points.Shows them in a table and represents them in the 2-dimensional simplex.
• For dynamics where the Jacobian is defined, computes the eigenvalues of the Jacobian of the dynamic at the exact isolated rest points using exact arithmetic.Computes numerical approximations to the eigenvalues and shows them in a table.
• Shows the components of rest points in a table and represents them in the 2-dimensional simplex.

Software functionalities
EvoDyn-3s generates phase portraits of evolutionary game dynamics, colors the background according to their speed (i.e. the modulus of the derivative vector), calculates the rest points, and provides the eigenvalues of the Jacobian at isolated rest points in order to analyze their local stability (in dynamics where the Jacobian is defined).
The user can analyze any 3 × 3 game by setting the values a ij of the 3 × 3 payoff matrix, which represent the payoff that a player using strategy i ∈ {1, 2, 3} obtains when interacting with a player using strategy j ∈ {1, 2, 3}.There is also a list of predefined games that the user can choose from.The next parameter in the interface is the dynamic.A popup menu shows the list of available evolutionary game dynamics, which are detailed in the Appendix.These constitute what we call here baseline dynamics, since they do not consider noise.
Parameter µ, called ''probability of random strategy'', or ''noise'', modifies the baseline dynamics by considering a flow of individuals who adopt a random strategy when revising their current one.This selection of random strategies may be based on mutations, errors, or experimentation.Formally, if ẋi = f i (x 1 , x 2 , x 3 ) corresponds to one of the baseline evolutionary dynamics detailed in the Appendix, the evolutionary dynamics considered when µ ̸ = 0 are given by ẋi = Small modifications of this sort can have important qualitative effects on game dynamics (see e.g.[13,14], and [15]).
The user can also easily choose the number of trajectories represented in the simplex (none, few, some or many), the length of the lines showing those trajectories, and the color gradient used to represent the speed of the dynamic.
Parameter maxTimeOut establishes the maximum time that the program is allowed to spend calculating the exact rest points of the dynamics.The default value is 1 s, which is enough in most cases.
Finally, the are two buttons which implement functionality that requires running EvoDyn-3s in Mathematica: • The ''Export figure'' button generates a camera-ready vector pdf file of the phase portrait, including a table showing the components of rest points and a table showing the numerical approximations to the isolated rest points and to their corresponding eigenvalues (for dynamics where the Jacobian is defined).
• The ''Generate full report'' button generates a Mathematica notebook detailing the equations corresponding to the selected dynamics, the phase portrait, and the exact representations of the components of rest points, of the isolated rest points and of the eigenvalues of the Jacobian (if it is defined).The report also includes numerical approximations to the exact values.

The hypercycle system
The hypercycle system is an ordinary differential equation introduced in [16] to model the organization of self-replicating molecules connected in a cyclic, autocatalytic way.The equation can be represented as an instance of the replicator dynamics, so the hypercycle system for a mixture of three macromolecules can be directly modeled with EvoDyn-3s.Fig. 3 shows the values of the parameters required to model this system, and the output thus obtained.It is clear that all interior solutions of this hypercycle system converge to the barycenter of the simplex, as proved by [17] and [18].EvoDyn-3s reports the exact eigenvalues of the Jacobian at the barycenter, i.e.

Zeeman's game
Zeeman's game was introduced in [19] to show that asymptotically stable states under the replicator dynamic do not necessarily correspond to Evolutionarily Stable Strategies (ESS).Fig. 4 shows the values of the parameters required to model Zeeman's game, and the output thus obtained.The interior Nash equilibrium of this game (i.e. the barycenter of the simplex) is not an ESS but is nevertheless asymptotically stable.This is clearly seen by noticing that the real part of the two eigenvalues of the Jacobian at this rest point is negative.EvoDyn-3s reports the exact eigenvalues of the Jacobian at the barycenter, i.e. − 1 3 ± i √ 2 3 , and their numerical approximations −0.333 ± 0.471i.

Impact
The following features of EvoDyn-3s distinguish it from alternate programs such as Dynamo, PDToolbox and egplot, and make it especially useful for research: • The selection of evolutionary dynamics considered in EvoDyn-3s (which is detailed in the Appendix) includes new evolutionary dynamics, based on interactions with partial information.Some of these dynamics are not available on any other program and are studied only in recent papers (e.g.[20,21]).
• Besides decimal approximations, EvoDyn-3s provides the rest points of the selected dynamics and the eigenvalues at isolated rest points (where the Jacobian is defined) using exact arithmetic.This implies, in particular, that if a solution is a rational number that does not admit a finite decimal representation, EvoDyn-3s provides the rational number.Similarly, if a solution is an algebraic number that does not admit a rational representation, EvoDyn-3s provides the algebraic number. 3Exact calculation of rest points is essential for rigorous local stability analysis [20].
• EvoDyn-3s can be run using free software.It does not require expensive software licenses, such as Mathematica or Matlab.EvoDyn-3s can be run with the free Wolfram CDF Player, which offers desktop functionality for Windows, Mac and Linux systems, and it also runs on some mobile platforms (e.g.those running iOS).The only functionality missing when running EvoDyn-3s on free software is the dynamic creation of output files.
• EvoDyn-3s has been designed to be used through a simple and intuitive user interface.It does not require any programming knowledge whatsoever.The user only needs to 3 Beyond degree four, most polynomials do not have roots that can be expressed in terms of radicals.An example of an algebraic number that cannot be expressed in terms of radicals is any of the roots of the polynomial specify the payoff matrix of the game and choose one of the available dynamics using the mouse.
To provide an idea of the interest and potential applications of this type of software, we note that -even with the barrier entry caused by the need to own a Mathematica license and possess some (basic) knowledge of the Mathematica environment -Dynamo has been used as supporting software for cutting-edge research and for generating high-quality graphs in articles published in leading journals in the field, such as Games and Economic Behavior [22,23], Journal of Economic Theory [24], Journal of Theoretical Biology [25], Theoretical Economics [26], Proceedings of the Royal Society of London [27] and Proceedings of the National Academy of Sciences [28].

Conclusions
EvoDyn-3s is open-source software that can generate phase portraits of various evolutionary game dynamics and data to analyze their equilibria.It includes classic and new evolutionary dynamics, it can run on free software and on different platforms, it does not require any programming knowledge and it provides exact results for rest points and their corresponding eigenvalues.Besides being useful to analyze evolutionary dynamics, it is intended to generate camera-ready high-quality graphs for publication.
expected payoff difference π i − π j , as long as this difference is positive.
where the function [a] + is such that [a] + = a if a > 0, and [a] + = 0 if a ≤ 0.
• Brown-von Neumann-Nash [31][32][33][34][35][36].The flow from strategy j to strategy i is proportional to the prevalence x j of strategy j and to the expected payoff difference π • Logit [37][38][39].The flow from strategy j to strategy i is proportional to the prevalence x j of strategy j and to the logit term e

∑
k e π k η −1 − x i Parameter η modulates the impact of payoff differences on the dynamics.As η → 0 the impact of payoff differences gets more acute and the Logit dynamics approach the best response dynamics.
• Single-match imitative logit.This is the mean dynamics of an imitation process in which players obtain a payoff in a random match with another player and occasionally revise their strategy.A revising agent looks at another randomly chosen player (and his obtained payoff) and chooses a strategy with probability proportional to the logit factors of the two obtained payoffs.• Imitate if better in one match [40].This is the mean dynamics of an imitation process in which players obtain a payoff in a random match with another player and occasionally revise their strategy.A revising agent looks at another randomly chosen player.If the payoff obtained by the other player is higher, the strategy of the other player is adopted.• Test another in one match [20,41].This is the mean dynamics of an adaptive process in which players obtain a payoff in a random match with another player and occasionally revise their strategy.A revising agent tests another alternative random strategy in a new random match.If the outcome obtained with the alternative strategy is preferred, that strategy is adopted.where the ''greater than'' function I {a>b} is defined by • Best response to one random player [21,42,43].This is the mean dynamics of an adaptive process in which players occasionally revise their strategy.A revising agent looks at another randomly chosen player and adopts a best response to the strategy of that random player.where BR(j) is the subset of strategies that are a best response to strategy j, i.e., {i : a ij = Max k (a kj )}.

Fig. 2 .
Fig. 2. Overall architecture of EvoDyn-3s.Symbol N in the figure denotes a numerical approximation.

ẋi = x i 3 ∑ j=1 3 ∑ k=1 3 ∑
m=1x j x k x m sign(a ik − a jm ) x m (x j I {a ik >a jm } − x i I {a jm >a ik } )]