Our work on the colour ringing project of Rock Pipits along the coast of Gower, and funded by the
British Ecological Society,
was featured on the Gower Ornithological Society’s Annual report!
So far we have mapped several territores, with many re-sights of marked birds. We are currently using the data collected from these observations
to build a theoretical model to better understand the landscape-scale consequences of local individual behaviour in territorial species!
If you want to know more about this species, the Royal Society for the Protection of Birds (RSPB) has an incredible wealth of information on
rock pipits and other species.
This project is supported by the Gower Ringing Group and the
British Trust of Ornithology (BTO), to which this project contributes data.
Lastly, if during your weekend walks around Gower you see any rock pipits wearing colour rings, please take note of them (or even better, a picture!)
and let us know. Thank you for your support! This is a citizen science project.
In this tutorial, Marco explains how to implement the main routines to replicate the effects of the paper on Noise-induced effects in collective dynamics and inferring local interactions from data
Part 1 Gillespie algorithm
The code contains a function for exactly simulating runs from the master equation through the Gillespie algorithm. The Gillespie algorithm is divided in four main parts:
1.- Find the velocity (or propensity) for every reaction. This is done by multiplying the reaction rate (which is conditioned on the reagents encountering each other) with the rate that the reagents do encounter each other. The rate that two reagents encounter each other can be found through the mass action assumption. Usually this is done simply by multiplying the concentration of the reagents, but this is an approximation that holds only for big population sizes because reagents do not interact with themselves, and the concentration of the species with and without one molecule can be quite different in small population sizes and omega. Hence it is possible to use the stochastic law of mass action, equal for one reaction to $ k \Omega \prod \frac{x_i!}{(x_i-1)! \Omega^{s_i}} $ (for details see David Schnoerr et al (2017) J. Phys. A: Math. Theor. 50, 093001). Since there are some factorials, it is hard to compute in practice, and has been rearranged in the code implementation.
2.- Find the waiting time before the next reaction occurs by sampling from an exponential distribution with mean rate equal to the sum of the reaction velocities.
3.- Find the reaction that happens by choosing from the possible reactions with a probability proportional to the reaction velocity.
4.- Update the states based on the reaction that happened.
The function has not been optimized but works for simulating every type of system starting from the microscopic reactions. It performs only one update and must be inserted inside a for loop for producing the entire simulation. As input, it requires:
1.- The current number of molecules/individuals in a specific state through a vector of size equal to the number of possible states and each entry the number of individuals in that state.
2.- A vector containing the reaction rates.
3.- A matrix indicating the stoichiometry of the reagents, with the number of rows equal to the number of reactions, and number of columns equal to the number of possible states, and entries equal to the number of species of that state that must encounter to produce the specific reaction.
4.- A matrix indicating the stoichiometry of the products, with the number of rows equal to the number of reactions, and number of columns equal to the number of possible states, and entries equal to the number of species of that state that are produced from that specific reaction.
5.- A parameter omega indicating the spatial scale of the system.
As output, it gives a vector containing the waiting time, the number of the reaction the occurred, and the new number of individuals/molecules in each state.
It could be a good test to modify the code to simulate clonal logistic growth by using these reactions: A -> 0 with rate 0.05, A -> 2A with rate 0.2, 2A -> A with rate 0.05. You can start multiple replicates by having only one individual and set omega to 50.
Part 2 Replication of the results of the paper “Noise-induced effects in collective dynamics and inferring local interactions from data”
Jhawar J, Guttal V. (2020) Phil. Trans. R. Soc. B 375: 20190381. http://dx.doi.org/10.1098/rstb.2019.0381
The replication of the results consists of two parts, one applied to the voter model, and one applied to the higher order interactions model (see the paper for details). The only difference between the two is that the higher order interactions also has a reaction of the type 2A + B -> 3A, which leads to a non-linear per capita rate of change in the ordinary differential equation and consensus formation with multi-stability in the limit of infinite population size (not present in the voter model). For each model it is shown examples of the dynamics through 20 replicates of the time series obtained iterating the Gillespie algorithm and the probability density of the population being in certain states. Interestingly, also the voter model shows consensus and multi-stability! This is because of the constructive effect of intrinsic noise. To better understand this result, it is possible to fit a Langevin equation to the time series we have simulated and find the functional form of the deterministic (drift) and stochastic (diffusion) terms. Despite the Langevin equation being an approximation, it explicitly defines a term for stochasticity, enabling us to investigate noise, quite cool.
To fit the data, we first have to find the autocorrelation function for figuring out the best time scale (see paper for details).
\[ACF(\tau) = \frac{\langle (M(t) - \langle M(t) \rangle ) (M(t+\tau)- \langle M(t) \rangle) \rangle}{\langle (M(t) - \langle M(t) \rangle)^{2} \rangle }\]
When I see an equation my first though is “how do I calculate that?”. In my opinion Tidyverse has some cool operations which can be used to translate the equation into a specific value. For example, the averages in the ACF are calculated over different things, which can be captured by grouping.
The last part is to find the deterministic and stochastic component, which produce respectively the first and second moments (see paper for details). The first moment is the average “jump” from a specific population state. The second moment is the variance of such “jumps”. In this case there are only two possible individual states, only one degree of freedom to describe the state of the entire population, but in more complex cases a “jump” could be a vector in a multidimensional state space.
Comparing the two models we see that the deterministic component “pushes” the population towards non-consensus in the voter model (equal number of individuals for both options), while there are two stable equilibria separated by and unstable equilibrium in the higher interaction model (equilibria are found when the “push” is equal to zero, and equilibria are stable when the “push” goes towards the equilibrium). For both models the noise is stronger when the population is undecided because there is more uncertainty in the reactions that can occur (the reactions have similar velocities), which is the definition of intrinsic noise. Hence, the voter model is pushed away from its non-consensus equilibrium condition, leading to noise induced consensus. Noise constructs order!
Marco has kindly share his code with us here.
This Friday 12th of January, we resume our activities of the modelling coffee, continuing a tutorial presentation about stochastic calculus,
Gentle introduction to stochastic calculus for population dynamics, by Gui Araujo.
The topics of this week are differential equation models, such as Fokker-Planck and Langevin equations, and a discussion on the concept of
noise and the differences between intrinsic and extrinsic noise.
This week Lucie travelled to Belfast to present her most recent findings on the effects of disturbances on meta food webs at the 2023 Annual Meeting of the British Ecological Society!
An excellent representation of the current research in our lab. It was very well received and Lucie got many interesting questions and comments.
Great news! Abby’s paper on the effects of microbial community structure on coral larval settlement was featured as a Spotlight by Trends in Microbiology!
The original paper: Linking differences in microbial network structure with changes in coral larval settlement can be downloaded here!
