GrowthDynamics.jl
This package implements various models of growing and evolving populations, both with and without spatial structure.
A number of useful observables are defined.
Quick Start
Let's give an overview of the components acting together.
States
A state is a structure comprised of a lattice (see Lattices), each entry representing a member of the population, its value representing its genotype, and various metadata. In unstructured populations the lattice is simply a dummy.
The metadata store information about the current state of the various genotypes, like number of individuals present, their fitness, and so on. Furthermore, a phylogenetic tree is recorded during simulation, enabling access to observables like most-recent common ancestors, or tracking of lineages.
julia> import GrowthDynamics.Populations: spheref
julia> import GrowthDynamics.Lattices: HexagonalLattice
julia> using GrowthDynamics.LatticeDynamics
julia> using GrowthDynamics.Observables
julia> state, _ = spheref(HexagonalLattice, 128, f=1/10, g1=0, g2=1)
(HexagonalLattice{Int64, Matrix{Int64}} 1 genotypes 1789 population, CartesianIndex{2}[CartesianIndex(60, 43), CartesianIndex(62, 43), CartesianIndex(63, 43), CartesianIndex(64, 43), CartesianIndex(65, 43), CartesianIndex(66, 43), CartesianIndex(67, 43), CartesianIndex(68, 43), CartesianIndex(70, 43), CartesianIndex(56, 44) … CartesianIndex(66, 86), CartesianIndex(67, 86), CartesianIndex(68, 86), CartesianIndex(69, 86), CartesianIndex(70, 86), CartesianIndex(71, 86), CartesianIndex(73, 86), CartesianIndex(63, 87), CartesianIndex(65, 87), CartesianIndex(67, 87)])
This prepares a state on a two-dimensional hexagonal lattice of size 128^2
that is unoccupied (genotype 0
is per definition understood as unoccupied.) except for a centered disk of genotype 1
that comprises ~1/10
of the total population.
Every mutation event introduces a new genotype. They are by default labeled consecutively by integers, but custom labels are possible.
Dynamics
Now evolve the population.
julia> eden_with_density!(state; T=128^2, # timesteps mu=1e0, # mutation rate per genome (not site!) d=1/100, # death rate fitness=(s,g_old,g_new)->1.0 + 0.1*randn() # function to assign fitness to new genotypes )
julia> show(state)
HexagonalLattice{Int64, Matrix{Int64}} 6944 genotypes 11384 population
We can plot (done using Makie.jl) the distribution of fitness values to check if it conforms to expectation
Observables
julia> # number of polymorphisms and diversity npolymorphisms(state), mean_pairwise(state)
(11598, 51.29354973680559)
julia> # alleles with frequency larger 0.01 first(sort(allele_spectrum(state, threshold=0.01), :fpop), 6)
6×5 DataFrame Row │ position npop fpop depth samples │ Int64 Int64 Float64 Int64 Int64 ─────┼───────────────────────────────────────────── 1 │ 777851975 114 0.0100141 11384 114 2 │ 619484943 114 0.0100141 11384 114 3 │ 626696516 114 0.0100141 11384 114 4 │ 343118016 115 0.0101019 11384 115 5 │ 157836211 115 0.0101019 11384 115 6 │ 296528449 115 0.0101019 11384 115