julia package to simulate clonal expansions in normal tissues, the simulation models stem cells undergoing (biased) drift. This was used in the following paper: Measuring the distribution of fitness effects in somatic evolution by combining clonal dynamics with dN/dS ratios.
] add StemCellModels
The first step is to setup a
StemCellModel object which sets all the parameters of the simulation.
Ncells = 10^5 SM = StemCellModel(Ncells, Δmut = 0.0, #"bias" of new driver mutations (ie selection coefficient) μp = 0.1, #passenger mutation rate per effective division μd = 0.1, # driver mutation rate per effective division tend = 80.0, # time in years to simulate for r = 0.5, #proportion of stem cells that result in symmetric division λ = 1 # rate of stem cell division (per year) )
Note that the effective of rate of division is given by
rλ, so it's important to consider this product. I would suggest fixing one of these numbers and varying the other. For example
r = 0.5, λ = 1 will result in the same dynamics as
r = 1, λ = 0.5. This is because we only consider stem cells in this model and do not include differentiated cells.
Then we can simulate a population of stem cells using the
sim = runsimulation(SM, progress = true, #include progress bar onedriver = true, #if this is set to true (default) then cells can only accumulate a single driver mutations restart = true # restart simulation if population dies out )
Note on the
onedriver option, setting this to
true means cells can only accumulate a single driver mutation. This is the regime in which the theoretical model described in the paper is valid. This assumption also seems to be mostly born out in the data but may be violated in some cases. If this is set to
false then cells will gain multiple driver mutations but their fitness will be the same as if they had a single driver mutation.
You can access the output using the fields in results object. calling
sim will also print out a summary.
sim.mutationfrequencies_d # frequencies of driver mutations sim.mutationfrequencies_p # frequencies of passenger mutations sim.mutationsize_d #number of cells each mutations is present int sim.mutationsize_p sim.hitchikers #binary array denoting whether each mutation is a hitchiker or not