Allow variable contribution of abiotic traits

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 * 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
 

docs/src/index.md

@@ -18,6 +18,7 @@ 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.

docs/src/model_description.md

@@ -25,6 +25,7 @@ Each species should have the following parameters in a dictionary objection. The
 * __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. 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.

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
@@ -55,6 +55,13 @@ function check_param_shapes(d)
     @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,21 +71,21 @@ 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])
+  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

@@ -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

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)