Merge pull request #32 from kavir1698/variable_selection

allow time-variable selection coefficients

CHANGELOG.md

@@ -1,6 +1,8 @@
 # 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.
+* Allow differential contribution of abiotic traits to fitness.
 
 # 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

@@ -18,8 +18,10 @@ See [Tutorial](@ref) for running the model, and [Model description](@ref) for a
 * Space is a grid on which individuals can migrate.
 * Possibility to model a complex environment with spatio-temporally varying resources.
 * Allows time varying optimal phenotypic value per site per species.
+* Allow differential contribution of traits to fitness.
 * 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,8 @@ 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).
+* __phenotype\_contribution\_to\_fitness__: The relative contribution of each __abiotic phenotype__ to fitness. This parameter can be `nothing` to denote equal contribution of all abiotic phenotypes, or it can be a vector of numbers with as many elements as there are abiotic phenotypes.
 * __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`).

examples/example1.jl

@@ -29,6 +29,7 @@ using EvoDynamics
 #   :growth_rate => 1.0,
 #   :expression_array => [0.28, 0.46],
 #   :selection_coefficient => 0.1,
+# :phenotype_contribution_to_fitness => nothing,
 #   :mutation_probabilities => [0.9, 0.0, 0.0],
 #   :mutation_magnitudes => [0.05, 0.0, 0.0],
 #   :N => [100],

examples/example2.jl

@@ -48,6 +48,7 @@ using Plots
 #   :growth_rate => 2.8,
 #   :expression_array => [0.28878032859775615, 0.4629421231828499, 0.26092147517051467, 0.952859489607121, 0.9638502824424, 0.05038142018016245, 0.05930756376654234, 0.033459292878885716, 0.32421526342800044, 0.9029235877297073, 0.7670060809312949, 0.12766808941531993, 0.8656895869985795, 0.342191940658253],
 #   :selection_coefficient => 0.1,
+#   :phenotype_contribution_to_fitness => nothing,
 #   :mutation_probabilities => [0.9, 0.0, 0.0],
 #   :mutation_magnitudes => [0.05, 0.0, 0.01],
 #   :N => [1000, 0, 0, 0],
@@ -93,6 +94,7 @@ using Plots
 #   :growth_rate => 1.5,
 #   :expression_array => [0.24923147816626035, 0.7155732641738595, 0.9655184311211502, 0.8638149724268844, 0.5075272565823061, 0.9189652626508431, 0.7897640036022151, 0.17091233765481717],
 #   :selection_coefficient => 0.1,
+#   :phenotype_contribution_to_fitness => nothing,
 #   :mutation_probabilities => [0.9, 0.0, 0.0],
 #   :mutation_magnitudes => [0.05, 0.0, 0.01],
 #   :N => [100, 0, 0, 0],

examples/paramfile1.jl

@@ -19,7 +19,8 @@ species1 = Dict(
   :growth_rate => 1.0,
   :expression_array => [0.28, 0.46],
   :selection_coefficient => 0.1,
-  :mutation_probabilities => [0.9, 0.0, 0.0],
+  :phenotype_contribution_to_fitness => nothing,
+  :mutation_probabiliti => [0.9, 0.0, 0.0],
   :mutation_magnitudes => [0.05, 0.0, 0.0],
   :N => [100],
   :environmental_noise => 0.01,

examples/paramfile2.jl

@@ -39,6 +39,7 @@ species1 = Dict(
   :growth_rate => 2.8,
   :expression_array => [0.28878032859775615, 0.4629421231828499, 0.26092147517051467, 0.952859489607121, 0.9638502824424, 0.05038142018016245, 0.05930756376654234, 0.033459292878885716, 0.32421526342800044, 0.9029235877297073, 0.7670060809312949, 0.12766808941531993, 0.8656895869985795, 0.342191940658253],
   :selection_coefficient => 0.1,
+  :phenotype_contribution_to_fitness => nothing,
   :mutation_probabilities => [0.9, 0.0, 0.0],
   :mutation_magnitudes => [0.05, 0.0, 0.01],
   :N => [1000, 0, 0, 0],
@@ -86,6 +87,7 @@ species2 = Dict(
   :growth_rate => 1.5,
   :expression_array => [0.24923147816626035, 0.7155732641738595, 0.9655184311211502, 0.8638149724268844, 0.5075272565823061, 0.9189652626508431, 0.7897640036022151, 0.17091233765481717],
   :selection_coefficient => 0.1,
+  :phenotype_contribution_to_fitness => nothing,
   :mutation_probabilities => [0.9, 0.0, 0.0],
   :mutation_magnitudes => [0.05, 0.0, 0.01],
   :N => [100, 0, 0, 0],

src/check_params.jl

@@ -1,5 +1,5 @@
 function check_param_names(d)
-  species_keys = [:name, :number_of_genes, :number_of_phenotypes, :abiotic_phenotypes, :biotic_phenotypes, :migration_phenotype, :migration_threshold, :ploidy, :epistasis_matrix, :pleiotropy_matrix, :growth_rate, :expression_array, :selection_coefficient, :mutation_probabilities, :mutation_magnitudes, :N, :vision_radius, :check_fraction, :environmental_noise, :optimal_phenotypes, :age, :recombination, :initial_energy, :bottlenecks, :reproduction_start_age, :reproduction_end_age, :abiotic_variance, :biotic_variance, :mating_scheme]
+  species_keys = [:name, :number_of_genes, :number_of_phenotypes, :abiotic_phenotypes, :biotic_phenotypes, :migration_phenotype, :migration_threshold, :ploidy, :epistasis_matrix, :pleiotropy_matrix, :growth_rate, :expression_array, :selection_coefficient, :mutation_probabilities, :mutation_magnitudes, :N, :vision_radius, :check_fraction, :environmental_noise, :optimal_phenotypes, :age, :recombination, :initial_energy, :bottlenecks, :reproduction_start_age, :reproduction_end_age, :abiotic_variance, :biotic_variance, :mating_scheme, :phenotype_contribution_to_fitness]
   model_keys = [:species, :generations, :space, :metric, :periodic, :resources, :interactions, :food_sources, :seed]
   species = d[:species]
   for sp in species
@@ -42,14 +42,26 @@ 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"
     # Checking the dimensions of optimal phenotypes
     @assert length(dd[:optimal_phenotypes]) == (d[:generations] + 1) "Species $spname: optimal phenotypes should be defined for every generation including step 0 (generations + 1)."
     @assert length(dd[:optimal_phenotypes][1]) == prod(d[:space]) "Species $spname: at each generation, optimal phenotypes should be defined for all sites."
     @assert length(dd[:optimal_phenotypes][1][1]) == length(dd[:abiotic_phenotypes]) "Species $spname : optimal phenotypes at each site should be defined for all abiotic phenotypes. If there is only one abiotic phenotype, define the optimal phenotype in a vector (vector of length 1)."
+    if isnothing(dd[:phenotype_contribution_to_fitness])
+      dd[:phenotype_contribution_to_fitness] = [1.0 for phen in dd[:abiotic_phenotypes]]
+    else
+    @assert typeof(dd[:phenotype_contribution_to_fitness]) <: AbstractArray "Species $spname: `phenotype_contribution_to_fitness` should be either `nothing` or a vector of floating numbers"
+    @assert eltype(dd[:phenotype_contribution_to_fitness]) <: AbstractFloat "Species $spname: elements of `phenotype_contribution_to_fitness` should be floating point numbers."
+    @assert length(dd[:phenotype_contribution_to_fitness]) == length(dd[:abiotic_phenotypes]) "Species $spname: length of `phenotype_contribution_to_fitness` should be as long as abiotic phenotypes."
+    end
     # mating scheme
     @assert typeof(dd[:mating_scheme]) <: Int "Species $spname: mating scheme should be an integer." 
     @assert in(dd[:mating_scheme], [-1, 0, 1]) "Species $spname: mating scheme should be one of the following values: -1, 0, 1."

src/interactions.jl

@@ -71,22 +71,22 @@ phenotypic_distance(ag1::Ind, ag2::Ind, model) = ag1.species == ag2.species ? ph
 """
 Abiotic distance to the optimal values
 """
-function abiotic_distance(phenotype, optimal, variance)
+function abiotic_distance(phenotype, optimal, variance, phenotypic_contributions)
   distance = 0.0
   counter = 0.0
   for i in eachindex(phenotype)
     d = abs(0.5 - cdf(Normal(optimal[i],  variance), phenotype[i]))
-    distance += d
-    counter += 1.0
+    distance += (d * phenotypic_contributions[i])
+    counter += phenotypic_contributions[i]
   end
   distance += distance
   distance /= counter
   return distance
 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])
+  rawW = 1.0 - abiotic_distance(abiotic_phenotype, return_opt_phenotype(species, pos, model), model.abiotic_variances[species], model.phenotype_contribution_to_fitness[species])
+  # W = adjust_abiotic_fitness(rawW, model.selectionCoeffs[species][model.step[1]+1])
   return rawW
 end
 

src/migration.jl

@@ -33,7 +33,7 @@ end
 function check_site(agent, site, model)
   phenotype = agent.abiotic_phenotype
   optimal = return_opt_phenotype(agent.species, site, model)
-  abiotic_distance(phenotype, optimal, model.abiotic_variances[agent.species])
+  abiotic_distance(phenotype, optimal, model.abiotic_variances[agent.species], model.phenotype_contribution_to_fitness[agent.species])
 end
 
 function pick_site(agent, sites, nsites, model)

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}}}}
@@ -57,6 +57,7 @@ struct Params{F<:AbstractFloat,I<:Int,N<:AbstractString}
   biotic_variances::Vector{F}
   abiotic_variances::Vector{F}
   mating_schemes::Vector{Int}
+  phenotype_contribution_to_fitness::Vector{Vector{F}}
 end
 
 """
@@ -173,6 +174,7 @@ function create_properties(dd)
   biotic_variances = [allspecies[i][:biotic_variance] for i in 1:nspecies]
   abiotic_variances = [allspecies[i][:abiotic_variance] for i in 1:nspecies]
   mating_schemes = [allspecies[i][:mating_scheme] for i in 1:nspecies]
+  phenotypic_contributions = [allspecies[i][:phenotype_contribution_to_fitness] for i in 1:nspecies]
 
   # reshape single value matrices to (1,1)
   if length(resources) == 1
@@ -181,7 +183,7 @@ function create_properties(dd)
     dd[:interactions] = reshape(dd[:interactions], 1, 1)
   end
 
-  properties = Params(ngenes, nphenotypes, growthrates, selectionCoeffs, ploidy, optvals, Mdists, Ddists, Ns, Ed, generations, nspecies, Ns, migration_traits, vision_radius, check_fraction, migration_thresholds, step, nnodes, biotic_phenotyps, abiotic_phenotyps, max_ages, names, dd[:food_sources], dd[:interactions], resources, dd[:resources], recombination, initial_energy, bottlenecks, repro_start, repro_end, biotic_variances, abiotic_variances, mating_schemes)
+  properties = Params(ngenes, nphenotypes, growthrates, selectionCoeffs, ploidy, optvals, Mdists, Ddists, Ns, Ed, generations, nspecies, Ns, migration_traits, vision_radius, check_fraction, migration_thresholds, step, nnodes, biotic_phenotyps, abiotic_phenotyps, max_ages, names, dd[:food_sources], dd[:interactions], resources, dd[:resources], recombination, initial_energy, bottlenecks, repro_start, repro_end, biotic_variances, abiotic_variances, mating_schemes, phenotypic_contributions)
 
   return properties, (epistasisMatS, pleiotropyMatS, expressionArraysS)
 end
@@ -270,13 +272,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
 

test/params.jl

@@ -30,6 +30,7 @@ species1 = Dict(
   :growth_rate => 2.0 , # max number of offsprings per mating mean of a Poisson
   :expression_array => [0.28878032859775615, 0.4629421231828499, 0.26092147517051467, 0.952859489607121, 0.9638502824424, 0.05038142018016245, 0.05930756376654234, 0.033459292878885716, 0.32421526342800044, 0.9029235877297073, 0.7670060809312949, 0.12766808941531993, 0.8656895869985795, 0.342191940658253],  # ngenes x ploidy long
   :selection_coefficient => 0.02,
+  :phenotype_contribution_to_fitness => nothing,
   :abiotic_variance => 1.0, # variance of a normal distribution used in determining the phenotypic distance of agents to the optimal environmental phenotypes. The larger the variance, the less important is the distance.
   :biotic_variance => 10.0, # same as above but for determining the biotic phenotypic distance between two individuals (used in any kind of interaction). Larger values mean that all pairs are equally likely to interact, regardless of their phenotype difference. 
   :mutation_probabilities => [0.99, 0.99, 0.99],  # for gene expression array, pleiotropy matrix and epistasis matrix, respectively
@@ -63,6 +64,7 @@ species2 = Dict(
   :growth_rate => 1.1,
   :expression_array => [0.24923147816626035, 0.7155732641738595, 0.9655184311211502, 0.8638149724268844, 0.5075272565823061, 0.9189652626508431, 0.7897640036022151, 0.17091233765481717],
   :selection_coefficient => 0.05,
+  :phenotype_contribution_to_fitness => nothing,
   :abiotic_variance => 1.0,
   :biotic_variance => 10.0,
   :mutation_probabilities => [0.9, 0.0, 0.0],

test/simulation_tests.jl

@@ -103,9 +103,9 @@ end
   @test EvoDynamics.phenotypic_distance(agent, agent2, model) == EvoDynamics.phenotypic_distance(agent2, agent, model)
 
   optphen = EvoDynamics.return_opt_phenotype(agent2.species, agent.pos, model)
-  @test EvoDynamics.abiotic_distance(agent2.abiotic_phenotype, optphen, model.abiotic_variances[agent2.species]) > 0
+  @test EvoDynamics.abiotic_distance(agent2.abiotic_phenotype, optphen, model.abiotic_variances[agent2.species], model.phenotype_contribution_to_fitness[agent2.species]) > 0
   agent2.abiotic_phenotype .= optphen
-  @test EvoDynamics.abiotic_distance(agent2.abiotic_phenotype, optphen, model.abiotic_variances[agent2.species]) == 0
+  @test EvoDynamics.abiotic_distance(agent2.abiotic_phenotype, optphen, model.abiotic_variances[agent2.species], model.phenotype_contribution_to_fitness[agent2.species]) == 0
   @test EvoDynamics.abiotic_fitness(agent2, model) == 1.0
 
   EvoDynamics.interact!(agent, agent2, model)