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