allow time-variable selection coefficients

CHANGELOG.md

@@ -1,6 +1,7 @@
 # v0.16
 
-Users can modify the sequence of events in the simulations by providing their own stepping functions in the `runmodel` functions.
+* Users can modify the sequence of events in the simulations by providing their own stepping functions in the `runmodel` functions.
+* Allows time-variable selection coefficients per species so that selection can be relaxed or increase during the simulation.
 
 # v0.15
 

Project.toml

@@ -1,7 +1,7 @@
 name = "EvoDynamics"
 uuid = "c8768967-421d-4a02-8523-37736f3dbe06"
 authors = ["Ali R. Vahdati <[email protected]>", "Carlos Melian"]
-version = "0.15.1"
+version = "0.16.0"
 
 [deps]
 Agents = "46ada45e-f475-11e8-01d0-f70cc89e6671"

docs/src/index.md

@@ -20,6 +20,7 @@ See [Tutorial](@ref) for running the model, and [Model description](@ref) for a
 * Allows time varying optimal phenotypic value per site per species.
 * Possibility of killing individuals at certain times and sites.
 * Modeling both haploid and diploid species.
+* Time-variable selection strength.
 * Includes mutation and recombination.
 * Easy data collection.
 * Runs replicates and collects data in a table automatically.

docs/src/model_description.md

@@ -24,7 +24,7 @@ Each species should have the following parameters in a dictionary objection. The
 * __pleiotropy\_matrix__: A _binary matrix_ (0s and 1s) with size _number of phenotypes_ times _l_. The pleiotropy matrix specifies the phenotypes that each locus affects.
 * __expression\_array__: A _vector_ of size _l_ that represent the expression amount of each locus determining its effect size.
 * __growth\_rate__: Mean of a Poisson distribution for number of offsprings per reproduction. This number is fully realized when fitness of a haploid individual is 1, or the distance between the biotic phenotypes of two diploid individuals is 0. Otherwise, actual growth rate is a fraction of this value. 
-* __selection\_coefficient__: A number between 0 and 1 that determines the importance of fitness. 0 is be a model without selection.
+* __selection\_coefficient__: A number between 0 and 1 that determines the importance of fitness. 0 is be a model without selection. If you want time-variable selection coefficient, provide a vector as long as number of generations plus 1 (to account for generation zero).
 * __mutation\_probabilities__: A _vector of three numbers_ each of which specifies the probability for a different type of mutations: mutation probability of the _expression\_array_, _pleiotropy\_matrix_, and _epistasis\_matrix_, respectively.
 * __mutation\_magnitudes__: A _vector of numbers_ with the same size as _mutation\_probabilities_ that determines the magnitude of mutation for each of the three categories. Specifically, the numbers are the variances of normal distributions with mean 0 for expression array and epistasis matrices, and probability of changing a 0 and 1 in in the pleiotropy matrix. When mutating an offspring, it first checks whether there will be a mutation (_mutation\_probabilities_), and if positive, a mutation with a magnitude determined by _mutation\_magnitudes_ occurs.
 * __N__: A _vector of integers_ for the initial number of individuals at each site.

docs/src/tutorial.md

@@ -98,4 +98,8 @@ addprocs(4)
 @everywhere param_file = "params.jl"
 @everywhere EvoDynamics.load_parameters(param_file)
 adata, mdata, models = runmodel(param_file, replicates=10, parallel=true)
-```
+```
+
+## Modifying the sequence of events and specific functions
+
+If you want to change the sequence of events within each time-step, you can do so by providing your own stepping function to the `runmodel` function. With this flexibility, you can also keep the sequence as it is, but change one specific event within the sequence by using a different function for it. A good staring point would be to copy the `agent_step!` function from the package (at `src/simulation.jl`).

src/check_params.jl

@@ -42,7 +42,12 @@ function check_param_shapes(d)
     @assert typeof(dd[:initial_energy]) <: Real "Initial energy should be a number"
     # bottleneck should be nothing or array/string
     @assert typeof(dd[:bottlenecks]) <: AbstractArray || isnothing(dd[:bottlenecks]) "bottleneck function for species $spname is not an `Array` or `nothing`"
-    @assert typeof(dd[:selection_coefficient]) <: AbstractFloat "selection coefficient of species $spname should be floating point number"
+    @assert typeof(dd[:selection_coefficient]) <: AbstractFloat || typeof(dd[:selection_coefficient]) <: AbstractFloat "selection coefficient of species $spname should be either a floating point number or an array of floating point numbers as long as generations plus 1 (for generation zero)."
+    if typeof(dd[:selection_coefficient]) <: AbstractFloat 
+      dd[:selection_coefficient] = [dd[:selection_coefficient] for gen in 0:d[:generations]]
+    else
+      @assert length(dd[:selection_coefficient]) == (length(dd[:generations]) + 1) "Species $spname: selection coefficients should be as long as generations plus 1 (to account for generation zero generation)" 
+    end
     @assert typeof(dd[:migration_threshold]) <: AbstractFloat "migration_threshold of species $spname should be floating point number"
     @assert typeof(dd[:abiotic_variance]) <: AbstractFloat "Species $spname: abiotic variance should be a floating number"
     @assert typeof(dd[:biotic_variance]) <: AbstractFloat "Species $spname: biotic variance should be a floating number"

src/interactions.jl

@@ -86,7 +86,7 @@ end
 
 function abiotic_fitness(abiotic_phenotype, species::Int, pos, model::ABM)
   rawW = 1.0 - abiotic_distance(abiotic_phenotype, return_opt_phenotype(species, pos, model), model.abiotic_variances[species])
-  # W = adjust_abiotic_fitness(rawW, model.selectionCoeffs[species])
+  # W = adjust_abiotic_fitness(rawW, model.selectionCoeffs[species][model.step[1]+1])
   return rawW
 end
 

src/simulation.jl

@@ -25,7 +25,7 @@ struct Params{F<:AbstractFloat,I<:Int,N<:AbstractString}
   ngenes::Vector{I}
   nphenotypes::Vector{I}
   growthrates::Vector{F}
-  selectionCoeffs::Vector{F}
+  selectionCoeffs::Vector{Vector{F}}
   ploidy::Vector{I}
   optvals::Vector{Vector{Matrix{Vector{F}}}}
   mutProbs::Vector{Vector{DiscreteNonParametric{Bool,Float64,Vector{Bool},Vector{Float64}}}}
@@ -270,13 +270,13 @@ end
 
 function adjust_fitness!(agent::Ind, model::ABM)
   W = agent.W < 0 ? 0.0 : agent.W
-  newW = 1.0 - ((1.0 - W) * model.selectionCoeffs[agent.species])
+  newW = 1.0 - ((1.0 - W) * model.selectionCoeffs[agent.species][model.step[1]+1])
   agent.W = newW
 end
 
 function adjust_fitness(W, species, model::ABM)
   W2 = W < 0 ? 0.0 : W
-  newW = 1.0 - ((1.0 - W2) * model.selectionCoeffs[species])
+  newW = 1.0 - ((1.0 - W2) * model.selectionCoeffs[species][model.step[1]+1])
   return newW
 end