We recently put together a brief paper describing the work presented here that contains a general overview of the project and some preliminary results and discussion. Below is a summary of the Methods section.
This repository simulates the ecological interactions between species under a stochastic framework of the generalized Lotka-Volterra model undergoing evolution subject to a multi-species Moran process. The generalized Lotka-Volterra (gLV) model is a set of differential equations that is used to express the ecological relationships between a set of S species in an ecological locale. The flexible parameterization of the gLV model allows for it to capture a wide variety of dynamic behaviors important to ecological stability and fluctuations amongst interacting species.
Here, we approach the gLV model under a stochastic framework, modeling entire ecological locale as a multi-species fixed population, Moran process. The stochastic process we consider here is a birth-death process for each individual species, wherein an individual from a single species is chosen to replicate proportional to birth rate and another individual from a different species (in order to treat this system as a continuous-time Markov chain) is chosen to perish proportional to death rate. To incorporate information about inter-species interactions, we introduce the interaction (alternatively called ecological community) matrix from the gLV model and use it to modulate the intrinsic birth and death rates of each of the species. In particular, positive contributions from one species increase a species' birth rate, while negative contributions from one species increase a species' death rate.
We consider an evolutionary perspective by introducing a mutant of a single species at time t = 0. Unlike the standard Moran process in which there is a fixed population size N for a single species comprised of the wild-type and a mutant, here there are additional interacting species that contribute to the relative fitness of the wild-type to the mutant. Because these are mutants of existing species within the ecological locale, we do not generate independent parameters for the mutant but rather perturb the birth and death rates of the wild-type and also perturb the impact of other species in the surrounding environment to the mutant. As such, we create a highly correlated but noisy copy of the wild-type species, generating a realistic mutant.
From this model, we are then able to calculate the fitness of both the wild-type individual in this species of interest as well as the mutant in this same species, yielding a relative fitness. Using the analytical results of the Moran process discussed in class, we are able to calculate the fixation probability of a single mutant individual given the relative fitness and population size. Using the parameters of the model and the continuous-time Markov chain structure of the model, we can then simulate this process with some number of other species in the ecological locale as parameterized by the interaction matrix many, many times, yielding a relatively accurate estimate for the true fixation probability with the additional species. Most importantly, we can then compare the analytical fixation probability from the isolated island approach of the Moran process to the more realistic, ecologically diverse stochastic process discussed above. From this work, we are particularly interested in determining how the fixation probability changes with the addition of more actors and total population size.
The general approach to simulation is to generate and report a base set of parameters and to then simulate the dynamics of such a process according to those parameters. There are three types of parameters to this model: (1) intrinsic birth rate, (2) intrinsic death rate, and (3) interactions with other species.
The intrinsic birth and death rates are specific to each species present in the model, and the two parameters for each wild-type species are drawn from a folded standard Normal distribution (the absolute value of draws from a standard Normal distribution). The birth and death rates for the mutant are then perturbed by taking two independent draws from a Uniform distribution on the interval [-0.1, 0.1], yielding highly similar yet different birth and death rates for the wild-type and mutant from the species of interest. Considering the birth and death rates are draws from a folded standard Normal distribution, we know that approximately 95% of the draws will fall below 2, so we consider the chosen interval to be sizeable and reasonable for perturbation.
Under the May method, the interaction matrix parameters are drawn from a standard Normal distribution, so a total of S (S - 1) parameters are drawn since we set the diagonal elements of the interaction matrix to be 0; the intrinsic birth and death rates handle the self-regulation aspect in this representation of the gLV model. Since these draws are independent and identically distributed, the interaction matrix need not by symmetric. Next, we perturb the interactions of the wild-type species of interest slightly to generate the relevant parameters for the mutant, using the same perturbation method of independent draws from a Uniform distribution on the interval [-0.1, 0.1].
To simulate the gLV model using a stochastic framework as a continuous-time Markov chain, we use the embedded chain transition matrix and transition rates as simulation tools. For each species, we modulate their intrinsic birth and death rates by using entries in the interaction matrix and the populations of those species. This yields adjusted birth and death rates that then factor into the transition rates and then the embedded chain transition matrix. Importantly, we treat the individuals of a single species to each be independent, exponentially distributed random variables such that the overall rate for a single species is linear in the rate for a single individual and proportional to the present number of individuals; mathematically speaking, the minimum of independent Exponential random variables is an Exponential random variable with a rate parameter that is the sum of the comprising Exponentials. As such, we can calculate probabilities for choosing individuals to replicate and die, where we ensure that the individuals that replicate cannot also die such that the Markov chain always transitions to a new state.
We then allow this process to run until a fixation or extinction event occurs in the mutant. The mutant fixes in the population if the wild-type goes extinct. As such, we do not interpret fixation in the context of this stochastic process to mean that the mutant takes over the entire population but rather takes over the entire species sub-population within the ecological locale. From each complete simulation, we determine and store whether or not the mutant fixes or goes extinct, the set of states that the Markov chain occupied through the simulation, and the exponential waiting times for each of the transitions. From this saved information alone, we are able to capture a number of statistics from the simulations, such as the estimated fixation probability (according to the definition given above), the number of times the stochastic model reverts to the classical Moran process (all the other species go extinct so only the wild-type and mutant species of interest), the average number of states occupied until absorption, and the average wait time associated with absorption.
There are two types of figures for this model. There are figures that explain the setup of the model, such as plots of the birth and death rates of all the species and the interaction matrix, both including the mutant. The birth and death rate plots simply show the species-specific birth and death rates for each of the species. The birth and death rates can be easily compared for a single species to ascertain its general fitness (if the birth rate > death rate, then the species is relatively fit; if the birth rate < death rate, then the species is relatively unfit). The species interaction matrix depicts the pairwise interactions between species using a blue-white-red diverging colormap. Thus, positive interactions (which increase fitness of the column species) are shown in red, while negative interactions (which decrease fitness of the column species) are shown in blue. Neutral interactions are shown in white.
We also plot histograms that capture the distributions for simulation length and exponential wait times for a more accurate portrayal of those distributions as opposed to just a measure of central tendency. There are also plots that explain the results of the simulations. Currently, these plots show the trajectories of either wild-type or mutant populations by transition, highlighting cases of fixation or extinction for the desired population. Trajectories that result in eventual extinction of the desired population are colored in red, while those trajectories that result in eventual fixation are colored in blue.