{'fraction_non_alphanumeric': 0.02857142857142857, 'fraction_numerical': 0.0044879898862199745, 'mean_word_length': 4.667383127350886, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Efficient Simulations of Individual Based Models for Adaptive Dynamics and the Canonical Equation

We propose a faster algorithm for individual based simulations for adaptive dynamics based on a simple modification to the standard Gillespie Algorithm for simulating stochastic birth-death processes. We provide an analytical explanation that shows that simulations based on the modified algorithm, in the deterministic limit, lead to the same equations of adaptive dynamics as well as same conditions for evolutionary branching as those obtained from the standard Gillespie algorithm. Based on this algorithm, we provide an intuitive and simple interpretation of the canonical equation of adaptive dynamics. With the help of examples we compare the performance of this algorithm to the standard Gillespie algorithm and demonstrate its efficiency. We also study an example using this algorithm to study evolutionary dynamics in a multi-dimensional phenotypic space and study the question of predictability of evolution.


I. INTRODUCTION
Evolution is essentially a birth death process involving mutations and selection. According to the traditional view, evolution optimizes scalar phenotypes like beak size, fecundity, body size etc. The dictum, "Survival of the fittest", leads to the fittest type "winning". Evolution is thus viewed as a dynamical system converging to an equilibrium or a steady state. Such a view fails to take into account interactions among the individuals comprising a population or interactions among various individuals present in a geographical location. The birth and death rates of the individual depend on its phenotypic composition as well as its interaction with the "environment". The "environment" consists of both the "abiotic" environment like temperature, climate etc. and the "biotic" environment which can be, for example, the types of all other individuals with which it interacts, for example, the competition it faces with regards to resources such as food and shelter. Another example of the "biotic" environment is the presence of "prey" or "predator" type of individuals in the population [1]. Studying evolution involves modeling this birth death process in order to address questions about origin of diversity, speciation and predictability of evolution. In order to have such a view, we need to take into account certain ecological mechanisms like the interaction among organisms present in the same geographical location -a phenomenon known as frequency dependence [1][2][3][4]. Adaptive dynamics is a framework that incorporates frequency dependent selection originating from ecological interactions and competition among species [1]. Under this framework, evolution is regarded as a continuous trajectory in the phenotypic space that unravels as rare mutants with a higher "fitness" invade a resident population. The evolutionary dynamics in the above framework is described by the " canonical equation" [5], which, for one * Electronic address: vmadhok@gmail.com dimension, is given by Here γ is a constant of proportionality that scales the rate of evolution. The function, f (x, y), is the measure of fitness of a rare mutant, y, in a monomorphic population with a resident trait x [5]. The derivative of the invasion fitness function, f (x, y), with respect to the mutant phenotype coordinate, y, gives the direction of motion of the evolutionary trajectory. Mutations cause the occurrence of rare types which can possibly invade the resident depending on their relative growth rate with respect to the resident which is given by the selection gradient (RHS of Eq 1).
The derivation of the canonical equation takes into account the stochastic nature of the birth-death process and is based on the ecological processes influencing the birth and death at the level of an individual [5][6][7]. Such a formulation assumes that mutations are small, they occur very rarely and the resident population is taken to be of infinite size [7]. The canonical equation (Eq. 1) describes the dynamics of a monomorphic population in the trait space driven by mutations, selection and invasion.
In order to address questions like the origin of diversity, to study the stochastic effects in evolution and to study the issue of predictability of evolution we need to construct individual based models. Individual based stochastic models based on the deterministic adaptive dynamics are a natural generalization of the latter. In these models, we represent the population as a "cloud" of points in the phenotypic space. Each individual, represented by a point in the cloud, has a given birth and death rate depending on the location of the point in the phenotypic space. Individuals with a higher birth rate have a greater probability to reproduce and the individuals with a higher death rate have a greater probability to die.
Individual-based realizations of the model are typically based on the Gillespie algorithm [8]. Often, simulations of these models take a long time for one to see interesting effects like diversification in the phenotype space. Moreover, once diversification occurs, especially in high dimensional phenotypic spaces, it usually leads to an increase in population of individuals and therefore making simulations even slower. In this work, we accomplish three things. First, we propose a simple modification of the standard Gillespie algorithm and make the individual based simulations more efficient. With the help of numerical simulations, we demonstrate how this leads to faster simulations. Secondly, using this algorithm and under certain assumptions, we give a very intuitive and simple interpretation of the canonical equation of adaptive dynamics. Thirdly, we discuss the issue of the predictability of evolution by considering individuals based models in high dimensional phenotypic spaces and demonstrate the "evolutionary butterfly effect".

II. INDIVIDUAL BASED MODELS TO MACROSCOPIC DETERMINISTIC ADAPTIVE DYNAMICS
We begin with the widely studied logistic model [9], We will call this the "canonical" model, as later on we will also be interested in a general family of such models.
Here N (x, t) is the population density of the phenotype x at time t and K(x) is the carrying capacity of the population consisting entirely of x individuals. The competition between individuals of phenotypes x and y is given by the competition kernel α(x, y) and an individual of phenotype x faces competition with an effective density α(x, y)N (y, t)dy. We set the intrinsic growth parameter r and α(x, x) to be equal to 1.
We consider x and y to be vectors of dimension d describing a multidimensional phenotypic space.
The canonical equations of adaptive dynamics can be derived from the stochastic individual based model in the limit of rare mutations, small mutational effects and infinite population sizes. Under these assumptions, Dieckmann and Law showed that adaptive dynamics is the first order approximation of the mean path averaged over infinitely many realizations of the stochastic simulations obtained from the individual based model. In order to derive the adaptive dynamics of trait x, we assume a monomorphic (Dirac-delta distribution) resident population in trait x with a globally stable equilibrium density given by K(x) independent of the dimension d of x. A key ingredient to obtain the ODEs describing the adaptive dynamics is the invasion fitness function whose gradient determines the direction of selection forces and hence of the deterministic evolutionary trajectory. The invasion fitness function is the difference between the per capita birth and per capita death rates of a rare mutant in a monomorphic population. The invasion fitness, for the logistic map, of a rare mutant y is its per capital rate of growth in a resident population with phenotype x and is given by In the above example for the logistic equation, the selection due to frequency dependence influences the death rates while the birth rates are neutral. In general, birth rates might also depend on interactions between individuals and the fitness function will have a more general form involving competition kernels for both birth and death. From the invasion fitness, we can derive the selection gradient to be Finally, the canonical equation of adaptive dynamics in terms of the selection gradient is given by Equation (5) is a system of d dimensional ODEs that describe a trajectory that may converge to a fixed point, a limit cycle, a quasi periodic orbit or exhibit chaos depending on the nature of the selection gradient s(x) and the initial conditions. Here M (x) describes the mutational process and for simplicity we assume it to be a d × d identity matrix. Individual based models describe the interactions between discrete individuals that possess multi-dimensional adaptive traits. Such a microscopic description is more general and in a sense more fundamental as compared to its deterministic macroscopic approximation like the adaptive dynamics which can be derived from it under suitable assumptions. In order to simulate individual based models, we adapt the Gillespie algorithm [8] for simulating chemical reactions to population biological stochastic processes. Individuals are treated as particles, and birth and death events as chemical reactions. In the deterministic logistic model, the birth rate of phenotype i, ρ i , with trait value x i is 1. When a birth event happens, a new individual is added to the population with a phenotype that is offset from the parent by a small mutation, which happens with probability µ, chosen from a uniform distribution with a small amplitude. Throughout this paper we choose µ = 1.0. Its death rate, δ i is given by N ef f (i)/K(x i ), where N ef f is the effective density experienced by phenotype x. For the individual based model, if there are N individuals, x 1 , ..., x N , at time t then the effective density experienced by individual i is and, therefore, the death rate is given by The individual based simulations are then implemented stochastically in the familiar way by performing one birth or death event [8] that constitutes one computational step that advances the system from time t to t + ∆t. In the above example for the logistic equation, the selection due to frequency dependence influences the death rates while the birth rates are neutral. In general, both the birth as well as the death rates can be affected by frequency dependent interactions among individuals described by competition terms like α birth (x j , x i ) and α death (x j , x i ).

III. A MODIFIED ALGORITHM FOR INDIVIDUAL BASED SIMULATIONS
As mentioned above, individual-based realizations of the model are based on the Gillespie algorithm [8] as described above and consists of assigning each individual i a constant reproduction rate ρ i = 1 and a death rate , as defined by logistic ecological dynamics. The total rate is given by the sum of all individual rates U = i (ρ i + δ i ). A particular "event", birth or death, is chosen at random with probability equal to the rate of this event divided by the total probability rate U . If a birth event is chosen, a new individual is added to the population with a phenotype that is offset from the parent by a small mutation chosen from a uniform distribution with a small amplitude. Modified Gillespie algorithm consists of simply choosing the individual with the maximum death rate and eliminating it when a death event is chosen. That is, we select max(δ 1 , δ 2 , ..., δ n ) and eliminate that individual. If both the birth as well as the death rates are affected by frequency dependent interactions among individuals described by competition terms like α birth (x j , x i ) and α death (x j , x i ), when the birth event is chosen, we select the individual with the maximum birth rate, ρ i , and a new individual is added to the population with a phenotype that is offset from the parent by a small mutation chosen from a uniform distribution with a small amplitude. We treat the death events analogously.
It turns out that this simple modification has remarkable consequences as far as simulations are concerned. We now discuss the consequences of this modification on the efficiency of the algorithm and also on the adaptive dynamics it gener-ates.

A. Analysis of the Algorithm
In order to understand the above algorithm and the reason behind its speed, we first consider a logistic model of the form The above model reduces to Eq. 2 for m = 1. Assuming the resident to be at the equilibrium population density, given by the Dirac delta function, . Then the per capita growth of a mutant y is given by Therefore, the invasion fitness, as defined above is given by From this, we can derive the selection gradient to be which simplifies to Therefore, the above selection gradient is m times the selection gradient given in Eq. 4.
Finally, the canonical equation of adaptive dynamics in terms of the selection gradient is given by which gives the same set of ODEs up to the scaling constant m as obtained previously with the canonical logistic model. We also note that the Eq.13 evolves m times faster than Eq. 5.
In constructing the individual based models for the logistic model in Eq. 8, we note that the birth rates remain constant at ρ i = 1, while the death rate for the individual i is given by Therefore, the vector of these death rates can be expressed as (δ m 1 , δ m 2 , ..., δ m n ) where (δ 1 , δ 2 , ..., δ n ) are the death rates of individuals corresponding to the "canonical" map. Let k be the individual with the highest death rate. This means max(δ m 1 , δ m 2 , ..., δ m n ) = δ m k . When m is sufficiently large, we have δ m k >> δ m i for all i = k. Therefore the probability of choosing k to die is close to 1. But this is exactly what we do in our modified algorithm! We deterministically pick the individual with the highest death rate and eliminate it. Therefore, the modified algorithm can be considered to be an individual based simulation of a birth death process governed by the logistic model given by Eq. 8. The value of m in Eq. 8 should be sufficiently large such that, we satisfy, δ m k >> δ m i for all i = k. We have already shown that the individual based simulations based on this model, in the deterministic limit, give rise to exactly the same equations of adaptive dynamics as the "canonical" model up to the scaling constant m.

B. Conditions for Diversification
We have shown that the individual based simulations based on the modified algorithm gives us the same adaptive dynamics equations as the standard algorithm. ODEs describing the adaptive dynamics are first order differential equations describing the evolution of a monomorphic population in the phenotypic space. In order to study diversification, we need to explore higher orders of the fitness function.
Adaptive diversification has been extensively studied at singular points of the dynamics that also happen to be evolutionarily unstable [1,4,10]. The canonical equation of adaptive dynamics has the form where the functions g i are given as Singular points of the dynamics are the points, x * at which the R.H.S. of Eq. (14) is zero. Singular point, x * is called convergence stable if the dynamics starting at all points sufficiently close to it converge to x * eventually. This occurs when the Jacobian at the singular point, J(x * ) < 0. These points are of great importance in adaptive dynamics as they can be potentially the points at which diversification occurs. For scalar traits, this can be seen by Taylor expanding the invasion fitness function with respect to the mutant y to the second order, we get The first term on the RHS, f (x, x) is zero for all x. Usually, the evolutionary dynamics can be accurately described by the first order term in (16) and for the invasion fitness function in Eq. 10, the second derivative is At equilibrium, when we have ∂f (x,y) ∂y y=x = 0, the first term on the RHS is zero. And the second term coincides with Eq.
17. Therefore, the conditions for potential diversification are identical for both maps. Hence we conclude that to study diversification through the modified algorithm, we can expect to get similar results as we would using the standard Gillespie algorithm.

IV. INTERPRETATION OF THE CANONICAL EQUATION
Canonical equation can be viewed as dynamics that unravel as individuals go through the birth-death process described above. We will give a geometric interpretation of the canonical equation through the birth-death process.
We first consider the time reversal of the dynamics in the logistic model. We first recall that a birth-death process in the logistic model has birth and death rates given by ρ i = 1 and δ i =  =ı α(x i , x j )/K(x i ). If we simply interchange the birth rates with the death rates, i.e., δ i = 1 and ρ i =  =ı α(x i , x j )/K(x i ), this gives us the time reversal of the evolutionary dynamics given by the logistic model. This can also be seen from the expression for the invasion fitness, Eq. 3 where the first term on the RHS represents the per capita birth rate and the second term represents the per capita death rate [5]. Interchanging birth and death rates merely changes the sign of the fitness function and the selection gradient thereby reversing the evolutionary dynamics given in Eq. 5. Therefore, a birth event in the forward time is equivalent to a death event in the time reversed picture and, likewise, a death event in the forward time is equivalent to a birth event in the time reversed picture.
When a birth event is chosen, a new individual is added to the population with a phenotype that is offset from the parent by a small mutation chosen from a uniform distribution with a small amplitude. Therefore, one can view the birth event as the displacement of the population cloud by a very small amount. In the time reversed picture, this birth event becomes a death event and will lead to the displacement of the population cloud in the opposite direction.
Considering the logistic model, we consider a time period ∆τ which is very small compared to the total time in which the dynamics unravels but nonetheless large enough that we have a large number of birth and death events occurring. For infinite population sizes, we will have a large number of birthdeath events happening even when ∆τ becomes very small. We also assume, the number of birth events N b will be equal to the number of death events N d in the time period ∆τ . This is because the population is assumed to be at equilibrium and, the probability of birth is equal to the probability of death event, i.e P (Birth) = P (Death). Therefore, invoking the law of larger numbers, the sequence of birth-death events will be typical sequences with N (birth) = N (death). Consider a monomorphic population cloud at a particular point x in the phenotypic space. In the limit of very small mutations and a sufficiently small ∆τ , a single birth event will cause a displacement of the population cloud given by ∆r i . Since for the logistic model, the birth events are neutral, and therefore i ∆r i = 0. The displacement by a single death event is, ∆q, and the direction of displacement is in the opposite direction to the dying individual k, where δ k = max(δ 1 , δ 2 , ..., δ n ).
. Assuming all individuals except k to be a monomorphic population, δ k = max α(x, x +∆q)/K(x +∆q) ). Thus, we seek to determine a vector displacement ∆q that will maximize δ k . If ∆q is in the direction of the gradient of the function α(x,y)K(x) K(y) with respect to y at the value x. We assume the same direction for ∆q for all the death events occurring in the time interval ∆τ . Therefore, the total net displacement is given by Taking the limits, N death → ∞, |∆q| → 0 and τ → 0 we get the canonical equation of adaptive dynamics, In general and as mentioned above, birth rates might also depend on interactions between individuals and the fitness function will have a more general form involving competition kernels for both birth and death. Let, α birth (x j , x i ) and α death (x j , x i ) be the competition kernels for birth and death respectively. Since we assume an infinite population limit at equilibrium we have N (birth) = N (death). In this case, we will need to add up the displacement by both birth as well as death events. For a single birth event, the displacement is, ∆q b , and the direction of displacement is in the direction of the newly born individual k, where ρ k = max(ρ 1 , ρ 2 , ..., ρ n ).
. Assuming all individuals, except k, as a monomorphic population, ρ k = max α birth (x, x +∆q b )/K(x +∆q b ) ). Thus, we seek to determine a vector displacement ∆q b that will maximize ρ k . If ∆q b is in the direction of the gradient of the function with respect to y at the value x. We assume the same direction for ∆q b for all the birth events occurring in the time interval ∆τ . The calculation for the net displacement, ∆q d , due to death events follows similarly. Therefore, the total net displacement, with N (birth) = N (death), is given by is the per capita growth rate of a mutant y in a resident population. In other words, f (x, y) is the invasion fitness function and we recover canonical equation of adaptive dynamics.
We interpret the above equation as the translation of a population cloud along the direction of instantaneous gradient of the fitness function with respect to the mutant phenotype. It is remarkable that in this derivation, unlike the standard derivation [5] which essentially considers the evolutionary trajectory as a sequence of invasions by the mutant phenotypes, we do not need to explicitly invoke concepts like "invasion fitness" but recover it through the translation of the population cloud undergoing a birth-death process. Equation 21 also suggests a faster way to implement birth death process when we have competition kernels for both birth and death events. At each iteration, we simply choose exactly one birth and one death event. For the birth event, the individual with the maximum birth rate is chosen and a new individual in its vicinity is produced. For the death event, as before, the individual with the maximum death rate is chosen to die. It is clear from the form of the invasion fitness function, Eq. 3, consisting of two separate terms for birth and death that such an algorithm will satisfy the canonical equation of adaptive dynamics. The advantage we get is faster simulations and also the populations in the stochastic simulations does not increase and remains manageable. This is in contrast to the standard algorithm, where the population typically increases after evolutionary branching and the simulations computationally intensive after some time.

V. EXAMPLES
In this section, we illustrate the efficiency of the modified algorithm with the help of numerical examples. We will also address the issue of unpredictability in evolution.

A. Evolution of a single trait
Our first example compares the performance of the standard and modified Gillespie algorithm for a single trait. For this, we use the "canonical" logistic model given in Eq. 2. We consider a Gaussian competition kernel given by The magnitude of σ α tells us about the strength of the competition.
For the carrying capacity, we assume a Gaussian form as well, K 0 refers to the maximum carrying capacity and x 0 is the phenotype with the maximum carrying capacity. The condition for diversification [1], ∂ 2 f (x,y) This can be interpreted as a competition between the relative effect of the competition and carrying capacity at the singular point. There has been great progress theoretically as well as a large amount of emperical evidence for symatric speciation . Fig 1 gives an example of the comparison between the standard implementation of the stochastic birth-death process and our modified procedure. The details of parameter values for the individual based simulations is given in the caption. We see the modified algorithm is far more efficient than the standard Gillespie algorithm. While the latter takes roughly 8 × 10 4 iterations for the first evolutionary branching to occur, the former shows multiple evolutionary branching events over 2 × 10 4 iterations.
B. Evolutionary chaotic dynamics in a multi dimensional phenotypic space and the question of predictability of evolution We assume that the complexity of the epistatic interactions between phenotypic components is determined by the competition kernel α(x, y) : R d × R d → R, which is, in general, a complicated non-linear function. In particular, when we have the competition kernel of the "canonical" form [33], we get For the carrying capacity function, we assume it to be of the form K(x) = exp(− i x 4 i /4) . Then (5) can be written as The above set of ODEs describe the adaptive dynamics of the monomorphic population in the phenotypic space. In [34,35], it was shown that the dynamical system described by Eq. (27) exhibits chaos with almost unit probability for sufficiently high dimensional phenotypic space, i.e. d ≥ 50.
In particular, we are interested when (5) results in chaotic behavior. The above model can also be transformed to give the adaptive dynamics that coincides with the known chaotic systems like the Lorenz model. We just need to fix the coefficients a i , b j and c k of the competition kernel in the above model. For example, specific choice of coefficients will lead to Lorenz model and other popular models in the literature to explore the stochastic effects inherent in the birth death simulations when the corresponding equations for the adaptive dynamics are chaotic. Assuming the adaptive dynamics to be governed by the equations of the Lorenz model, we then proceed to study and construct the corresponding individual based model and simulate it the modified version of the standard algorithm.
Starting two simulations from exactly the same initial conditions and same mutation rates we see that the evolutionary trajectories follow completely different paths after sometime 2. Therefore the intuition about the evolutionary trajectories occupying two completely different regions of phenotypic space is correct. And this occurs purely due to a combination of stochastic effects due to mutations and the underlying  (Fig 1(a)) and the modified algorithm on the right (Fig 1(b)). Parameter values for both (a) and (b): σ k = 2.0, a Gaussian competition kernel with σ α = 0.75. K 0 = 200. Initial population of 200 individuals was initialized about x = 1.0 with a small variance. The mutation size after the birth event has a magnitude σ µ = 0.1 and occurs with a probability µ = 1.0. We see the modified algorithm is far more efficient than the standard Gillespie algorithm. While the latter takes roughly 8 × 10 4 iterations for the first evolutionary branching to occur, the former shows multiple evolutionary branching events over 2 × 10 4 iterations.
chaotic adaptive dynamics. Firstly, this tells us, contrary to the traditional view, that evolution is not a directional optimization process that goes to equilibration [34]. It might be a chaotic trajectory that is ergodic and sensitive to the initial conditions as well as the to the nature of very small mutational changes for exactly the same initial conditions. Even more fascinating is the sensitive dependence of a chaotic evolutionary trajectory on very small mutations, or the "butterfly effect" which casts a serious doubt over our ability to predict evolution over long time scales. Fig 2 demonstrates the "evolutionary butterfly effect" for two instances of individual based simulations starting from exactly the same initial conditions when the underlying adaptive dynamics is chaotic. This is contrasted with individual based simulations of a cyclic adaptive dynamics which remains predictable throughout its course.

VI. CONCLUSION
Individual based models are essential in order to study various aspects of evolution like diversification, predictability, adaptive radiation and frequency dependence. It has been observed that these models follow the adaptive dynamics trajectories remarkably well, even when the underlying adaptive dynamics is chaotic [33]. In this work, we derive the canonical equation of adaptive dynamics in a single step from a modified algorithm for simulating individual based models. We therefore provide a very simple and intuitive way of under-standing the correspondence between adaptive dynamics and individual based models. Using a simple modification to the Gillespie algorithm, we were able to simulate individual based birth-death models with significant efficiency. With the help of examples, we demonstrate the efficiency of the modified algorithm. We also demonstrate how the modified algorithm can be interpreted as simulating the individual based model associated with a generalized logistic equation as given by Eq. 8 and how this yields exactly the same equations for adaptive dynamics as the standard logistic model. We also discuss how the conditions for evolutionary branching for both algorithms remain same. Lastly, we discussed how to simulate the birth death process when both birth and death events have a competition kernel associated with them while, at the same time, keeping the total population size bounded despite of evolutionary branching. Our analysis gives an intuitive interpretation of the canonical equation of adaptive dynamics as the translation of the population cloud in the phenotypic space due to birth-death events. Such an approach has interesting connections to the stochastic quantization in quantum field theory and the birth death events can be viewed as creation and annhilation operators of a quantum field [36]. A complete analysis of these connections is the focus of our future work.
Our simulations also enables us to address important questions regarding the nature of evolutionary dynamics. In [34] it was shown that phenotypic properties can combine in complicated ways and as a result evolutionary trajectories can exhibit chaos. It shows that, contrary to the traditional view, evolution is not a directional optimization process that goes to equilibration. Evolution might be a complicated trajectory in When the adaptive dynamics is the Lorenz attractor (red). Individual based models (black and blue) show the "butterfly effect". Starting from exactly same initial conditions (the green dot), they evolve to different locations on the attractor (the black and the blue dot). (c) Example of another chaotic adaptive dynamics (red) and correspondence with the individual based models (black and blue).
the phenotypic space that is continuously twisting, turning and folding upon itself on a chaotic attractor. Our numerical simulations involving individual based simulations confirm that chaotic trajectories of evolution are sensitive to very small mutational changes. Even when these trajectories start with exactly the same initial conditions, they can end up in completely different regions of the phenotypic space. Therefore, if Gould's tape of life were to be run again, with exactly same initial conditions and dynamics, we can end up with com-pletely different biodiversity on earth. Moreover, the size of the mutations to cause this unpredictability of evolution need not be big as the underlying chaotic dynamics amplify small errors arising due to mutations. We hope our work is useful for simulations involving birth death processes in high dimensional phenotypic spaces as well as the study of evolution of language, culture, religions and agent based models in economics.