From b0e518483d9cf40709f50db5605ebb95a4003682 Mon Sep 17 00:00:00 2001
From: "Ali R. Vahdati" <arv.kavir@gmail.com>
Date: Mon, 29 Aug 2022 15:20:14 +0200
Subject: [PATCH 1/9] Fix the bug in the new UI

---
 src/check_params.jl | 15 ++++++++++++---
 src/load_params.jl  |  2 +-
 src/reproduction.jl | 20 ++++++++++---------
 src/simulation.jl   | 36 +++++++++++++++-------------------
 test/params.jl      | 47 +++++++++++----------------------------------
 5 files changed, 51 insertions(+), 69 deletions(-)

diff --git a/src/check_params.jl b/src/check_params.jl
index f4cff07..d4df8f3 100644
--- a/src/check_params.jl
+++ b/src/check_params.jl
@@ -1,4 +1,4 @@
-function check_yml_params(d)
+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, :bottleneck_function, :reproduction_start_age, :reproduction_end_age]
   model_keys = [:species, :generations, :space, :metric, :periodic, :resources, :interactions, :food_sources, :seed]
   species = d[:species]
@@ -41,9 +41,13 @@ function check_param_shapes(d)
     @assert typeof(dd[:recombination]) <: Real "recombination of species $spname should be type numeric"
     @assert typeof(dd[:initial_energy]) <: Real "Initial energy should be a number"
     # bottleneck should be nothing or array/string
-    @assert typeof(dd[:bottleneck_function]) <: Function || isnothing(dd[:bottleneck_function]) "bottleneck function for species $spname is not a function or `nothing`"
+    @assert typeof(dd[:bottleneck_function]) <: AbstractArray || isnothing(dd[:bottleneck_function]) "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[:migration_threshold]) <: AbstractFloat "migration_threshold of species $spname should be floating point 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)."
   end
   # Space should be 2D
   if !isnothing(d[:space])
@@ -52,7 +56,12 @@ function check_param_shapes(d)
   @assert typeof(d[:space]) <: Tuple || isnothing(d[:space]) "Space should be either a tuple or `nothing`"
   # Resources
   
-  @assert typeof(d[:resources]) <: Function "Resources function not defined"
+  @assert typeof(d[:resources]) <: AbstractArray "Resources function no an array"
+  @assert length(d[:resources]) == d[:generations] + 1 "resources is not as long as generations+1"
+  if !isnothing(d[:space])
+    @assert length(d[:resources][1]) == prod(d[:space]) "Each inner array of resources is not as long as space size" 
+    @assert size(d[:resources][1]) == d[:space] "Inner arrays of resources does not have the same dimension as the space."
+  end
 
   # Interactions
   @assert typeof(d[:interactions]) <: AbstractArray "Interactions should be an array(matrix)"
diff --git a/src/load_params.jl b/src/load_params.jl
index e016071..d84ebdb 100644
--- a/src/load_params.jl
+++ b/src/load_params.jl
@@ -1,6 +1,6 @@
 function load_parameters(paramfile::String)
   include(paramfile)
-  check_yml_params(model_parameters)
+  check_param_names(model_parameters)
   check_param_shapes(model_parameters)
   return model_parameters
 end
diff --git a/src/reproduction.jl b/src/reproduction.jl
index f93d36e..fdc2a22 100644
--- a/src/reproduction.jl
+++ b/src/reproduction.jl
@@ -2,7 +2,7 @@
 """
 Returns a bitarray for sites to be selected from the first (false) and second (true) homologous chromosome.
 """
-function crossing_overs(nsites::Int, ncrossing_overs::Int) 
+function crossing_overs(nsites::Int, ncrossing_overs::Int)
   output = falses(nsites)
   if ncrossing_overs == 0
     return output
@@ -17,7 +17,7 @@ function crossing_overs(nsites::Int, ncrossing_overs::Int)
     if counter < ncrossing_overs && site == breaking_points[counter+1]
       counter += 1
       output[site] = last
-      last = !last 
+      last = !last
     else
       output[site] = last
     end
@@ -45,7 +45,7 @@ function create_gamete(agent, cross_overs, nsites, first::Bool)
 
   q_gamete = agent.q[indices1]
   q_gamete[cross_overs] .= agent.q[indices2][cross_overs]
-  
+
   return epistasisMat_gamete, pleiotropyMat_gamete, q_gamete
 end
 
@@ -68,14 +68,14 @@ function create_one_offspring(ag1::Ind, ag2::Ind, model::ABM)
   gametes2 = create_gamete(ag2, cross_overs2, nsites, rand((true, false)))
 
   sex = rand((true, false))
-  interaction_history = MVector{model.nspecies, Int}(fill(0, model.nspecies))
-  
+  interaction_history = MVector{model.nspecies,Int}(fill(0, model.nspecies))
+
   # Merge gametes
-  nsites2 = nsites*2
-  episMat = MMatrix{nsites2, nsites2}(hcat(gametes1[1], gametes2[1]))
+  nsites2 = nsites * 2
+  episMat = MMatrix{nsites2,nsites2}(hcat(gametes1[1], gametes2[1]))
   pleioMat = MMatrix{model.nphenotypes[species],nsites2}(hcat(gametes1[2], gametes2[2]))
   q = MVector{nsites2}(vcat(gametes1[3], gametes2[3]))
-  abph = get_abiotic_phenotype(species, episMat, pleioMat, q, model) 
+  abph = get_abiotic_phenotype(species, episMat, pleioMat, q, model)
   bph = get_biotic_phenotype(species, episMat, pleioMat, q, model)
   W = abiotic_fitness(abph, species, ag1.pos, model)
   W = adjust_fitness(W, ag1.species, model)
@@ -86,7 +86,9 @@ function create_one_offspring(ag1::Ind, ag2::Ind, model::ABM)
 end
 
 """
-Sexual reproduction for diploid individuals
+Sexual reproduction for diploid individuals.
+
+# TODO: Currently this is assortative mating. Add a parameter to allow nonassortative, and disassortative mating.
 """
 function reproduce!(ag1::Ind, ag2::Ind, model::ABM)
   reproduction_success = 1.0 - phenotypic_distance(ag1, ag2, model)
diff --git a/src/simulation.jl b/src/simulation.jl
index e5a5ecc..ece89a5 100644
--- a/src/simulation.jl
+++ b/src/simulation.jl
@@ -19,13 +19,13 @@ mutable struct Ind{B<:AbstractFloat,C<:AbstractArray,D<:AbstractArray,E<:Abstrac
   isalive::Bool
 end
 
-struct Params{F<:AbstractFloat,I<:Int,N<:AbstractString,X<:Function, Fun, FunVec}
+struct Params{F<:AbstractFloat,I<:Int,N<:AbstractString}
   ngenes::Vector{I}
   nphenotypes::Vector{I}
   growthrates::Vector{F}
   selectionCoeffs::Vector{F}
   ploidy::Vector{I}
-  optvals::Vector{X}
+  optvals::Vector{Vector{Matrix{Vector{F}}}}
   mutProbs::Vector{Vector{DiscreteNonParametric{Bool,Float64,Vector{Bool},Vector{Float64}}}}
   mutMagnitudes::Vector{Vector{UnivariateDistribution{S} where S<:ValueSupport}}
   N::Vector{Vector{I}}
@@ -46,10 +46,10 @@ struct Params{F<:AbstractFloat,I<:Int,N<:AbstractString,X<:Function, Fun, FunVec
   food_sources::Matrix{F}
   interactions::Matrix{F}
   resources::Matrix{I}
-  resources_org::Fun
+  resources_org::Vector{Matrix{I}}
   recombination::Vector{Poisson{F}}
   initial_energy::Vector{F}
-  bottlenecks::FunVec
+  bottlenecks::Vector{Vector{Matrix{F}}}
   repro_start::Vector{Int}
   repro_end::Vector{Int}
 end
@@ -166,11 +166,11 @@ function create_properties(dd)
   bottlenecks = [allspecies[i][:bottleneck_function] for i in 1:nspecies]
   repro_start = [allspecies[i][:reproduction_start_age] for i in 1:nspecies]
   repro_end = [allspecies[i][:reproduction_end_age] for i in 1:nspecies]
-  resources = dd[:resources](0)
+  resources = dd[:resources][1]
 
   # reshape single value matrices to (1,1)
   if length(resources) == 1
-    resources = reshape(resources, 1,1)
+    resources = reshape(resources, 1, 1)
     dd[:food_sources] = reshape(dd[:food_sources], 1, 1)
     dd[:interactions] = reshape(dd[:interactions], 1, 1)
   end
@@ -181,17 +181,17 @@ function create_properties(dd)
 end
 
 function return_opt_phenotype(species::Int, site::Int, model::ABM)
-  ss = vertex2coord(site, model)
-  model.optvals[species](ss, model)
+  # ss = vertex2coord(site, model)
+  model.optvals[species][model.step[1]+1][site]
 end
 
 function return_opt_phenotype(species::Int, site::Tuple{Int,Int}, model::ABM)
-  model.optvals[species]( site, model)
+  model.optvals[species][model.step[1]+1][site...]
 end
 
 function model_step!(model::ABM)
   model.step[1] += 1
-  model.resources .= model.resources_org(model.step[1])
+  model.resources .= model.resources_org[model.step[1]+1]
 end
 
 function agent_step!(agent::Ind, model::ABM)
@@ -206,13 +206,13 @@ function agent_step!(agent::Ind, model::ABM)
   burn_energy!(agent)
   # consume basic energy if agent can
   consume_food!(agent, model)
-  # interact with other species
-  interact!(agent, model)
   # Kill the agent if it doesn't have energy
   if agent.isalive && agent.energy < 0
     remove_agent!(agent, model)
     return
   end
+  # interact with other species
+  interact!(agent, model)
   # survive!(agent, model)
   if !agent.isalive
     return
@@ -221,23 +221,19 @@ function agent_step!(agent::Ind, model::ABM)
   migrate!(agent, model)
   # reproduction for the haploid
   reproduce!(agent, model)
-  
+
   if agent.isalive && agent.age ≥ model.max_ages[agent.species]
     remove_agent!(agent, model)
   end
   # bottleneck
-  if agent.isalive && model.bottlenecks[agent.species](agent, model)
-    if EvoDynamics.bottleneck(agent, model)
+  bn_prob = model.bottlenecks[agent.species][model.step[1]+1][agent.pos...]
+  if agent.isalive && bn_prob > 0.0
+    if rand() < bn_prob
       remove_agent!(agent, model)
     end
   end
 end
 
-function bottleneck(agent, model)
-  # @info "Package bottleneck function"
-  return false
-end
-
 function remove_agent!(agent, model)
   agent.isalive = false
   kill_agent!(agent, model)
diff --git a/test/params.jl b/test/params.jl
index 84f9a1f..ea47efb 100644
--- a/test/params.jl
+++ b/test/params.jl
@@ -1,30 +1,5 @@
-# NB keep the order of sections 1, 2, and 3.
-
-### 1. functions
-##----------------------------------------------------------------
-
-function bn(agent::AbstractAgent, model::ABM)
-  return false
-end
-
-function bn1(agent::AbstractAgent, model::ABM)
-  return false
-end
-
-function optphens1(site::Tuple{Int,Int}, model::ABM)
-  return [1.5]
-end
-
-function optphens2(site::Tuple{Int,Int}, model::ABM)
-  phenotypic_values = [[0.7614758101208934 2.3911353663016586; 2.2091361414343313 1.7858661540288683; 0.7920974352892358 0.7630236717263839; 2.587205421882114 2.311211631439866],
-    [0.6437641305315445 2.097629262577973; 0.9954545836033715 1.1314248848669466; 2.469792530385348 1.490299526620522; 1.6158867433451882 0.2566056477862022]]
-  agent_site = (site[2] - 1) * size(model.space)[2] + (site[1])
-  if model.step[1] > 7
-    return phenotypic_values[1][agent_site, :]
-  else
-    return phenotypic_values[2][agent_site, :]
-  end
-end
+generations = 14
+space = (2, 2)
 
 function env_resources(time::Int)
   if time < 7
@@ -34,7 +9,7 @@ function env_resources(time::Int)
   end
 end
 
-## 2. Species properties
+## 1. Species properties
 ##----------------------------------------------------------------
 
 species1 = Dict(
@@ -60,13 +35,13 @@ species1 = Dict(
   :N => [1000, 0, 0, 0],  # number of individuals per site at time 0
   :environmental_noise => 0.0,  # variance of a normal distribution with mean 0
   # each row is the optimal phenotypes for each site for all abiotic traits. There are as many element as number of sites times number of abiotic traits. The first N elements are for the first abiotic trait, where N is the number of sites, and so on.
-  :optimal_phenotypes => optphens1,
+  :optimal_phenotypes => [fill([1.5 for p in 1:1], space...) for t in 0:generations], # optimal phenotypes per site and time for each phenotype
   :age => 5,  # max age
   :reproduction_start_age => 1,
   :reproduction_end_age => 5,
   :recombination => 1, # Mean of a Poisson distributions for number of crossing overs
   :initial_energy => 0, # A parameter for parental care of infants. Values more than 0 indicate that newly born individuals can survive for a number of times without requiring food from the environment/other species. The consumption rate (i.e. how many generations this initial energy suffices) is determined by the sum of the corresponding rows in "food sources"
-  :bottleneck_function => bn  # The function has two argument: agent and model. It returns true or false for death or survival, respectively.
+  :bottleneck_function => [fill(0.0, space...) for t in 0:generations]  # an array of arrays with probablity of external death at each site and generation.
 )
 
 species2 = Dict(
@@ -89,27 +64,27 @@ species2 = Dict(
   :mutation_magnitudes => [0.05, 0.0, 0.01],
   :N => [1000, 0, 0, 0],
   :environmental_noise => 0.01,
-  :optimal_phenotypes => optphens2,
+  :optimal_phenotypes => [fill([1.5 for i in 1:2], space...) for t in 0:generations],
   :age => 3,
   :reproduction_start_age => 1,
   :reproduction_end_age => 3,
   :recombination => 1,
   :initial_energy => 0,
-  :bottleneck_function => bn1
+  :bottleneck_function => [fill(0.0, space...) for t in 0:generations]
 )
 
-## 3. Model parameters
+## 2. Model parameters
 ##----------------------------------------------------------------
 
 #NB this dict should be called model_parameters
 
 model_parameters = Dict(
   :species => [species1, species2],
-  :generations => 14,  # number of simulation steps
-  :space => (2, 2),
+  :generations => generations,  # number of simulation steps
+  :space => space,
   :metric => "chebyshev",  # how many neighbors a space site has. "chebyshev" metric means that the r-neighborhood of a position are all positions within the hypercube having side length of 2*floor(r) and being centered in the origin position. "euclidean" metric means that the r-neighborhood of a position are all positions whose cartesian indices have Euclidean distance ≤ r from the cartesian index of the given position.
   :periodic => false,  # whether boundaries of the space are connected
-  :resources => env_resources,  # a function that returns a vector of integers for the available resources per site at a given time. It accepts only the model object as the input.
+  :resources => [env_resources(x) for x in 0:generations],  # a function that returns a vector of integers for the available resources per site at a given time. It accepts only the model object as the input.
   :interactions => [0.1 0.0; 0.0 -1.0],  # How individuals from different species interact. value  is strength of interaction (between 0 and 1). Sign is the direction of interaction where positive means similar individuals interact more strongly and negative is dissimilar ones tend to interact more. 
   :food_sources => [1.0 0.0; 0.5 0.0],  # What each species feeds on (consumption rate). Has priority over interactions. Non-zero diagonal means the food resource is from the environment. It will be read from rows (species order) to columns (species order).
   :seed => 3

From c5beaa3d8c5045c0f01bd59b90a4fba6afd8f33e Mon Sep 17 00:00:00 2001
From: "Ali R. Vahdati" <arv.kavir@gmail.com>
Date: Mon, 29 Aug 2022 19:19:13 +0200
Subject: [PATCH 2/9] Add feature and fix bugs

---
 src/check_params.jl |  7 ++++++-
 src/interactions.jl | 12 +++++++-----
 src/migration.jl    |  2 +-
 src/reproduction.jl | 10 +++++++---
 src/simulation.jl   | 18 +++++++++++-------
 test/params.jl      | 20 +++++++++++++-------
 6 files changed, 45 insertions(+), 24 deletions(-)

diff --git a/src/check_params.jl b/src/check_params.jl
index d4df8f3..13ba98e 100644
--- a/src/check_params.jl
+++ b/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, :bottleneck_function, :reproduction_start_age, :reproduction_end_age]
+  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, :bottleneck_function, :reproduction_start_age, :reproduction_end_age, :abiotic_variance, :biotic_variance, :mating_scheme]
   model_keys = [:species, :generations, :space, :metric, :periodic, :resources, :interactions, :food_sources, :seed]
   species = d[:species]
   for sp in species
@@ -44,10 +44,15 @@ function check_param_shapes(d)
     @assert typeof(dd[:bottleneck_function]) <: AbstractArray || isnothing(dd[:bottleneck_function]) "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[: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)."
+    # 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."
   end
   # Space should be 2D
   if !isnothing(d[:space])
diff --git a/src/interactions.jl b/src/interactions.jl
index 56c5366..5d10dfa 100644
--- a/src/interactions.jl
+++ b/src/interactions.jl
@@ -28,7 +28,7 @@ get_abiotic_phenotype(ag::Ind, model::ABM) = get_abiotic_phenotype(ag.species, a
 """
 Calculates the average distance between two sets of phenotypes.
 
-Distance is a probability and is calculated by the number of standard deviation that phenotype 2 is different from phenotype 1.
+Distance is a probability and is calculated by the number of standard deviations that phenotype 2 is different from phenotype 1.
 
 Variance can stand for selection coefficient, with a larger variance, less sever the difference.
 """
@@ -50,7 +50,7 @@ function phenotypic_distance(phenotype1, phenotype2, interaction_coeff, variance
   return distance
 end
 
-phenotypic_distance(ag1::Ind, ag2::Ind, model) = phenotypic_distance(ag1.biotic_phenotype, ag2.biotic_phenotype, model.interactions[ag1.species, ag2.species], variance)
+phenotypic_distance(ag1::Ind, ag2::Ind, model) = phenotypic_distance(ag1.biotic_phenotype, ag2.biotic_phenotype, model.interactions[ag1.species, ag2.species], model.biotic_variances[ag1.species])
 
 """
 Abiotic distance to the optimal values
@@ -69,7 +69,7 @@ function abiotic_distance(phenotype, optimal, variance)
 end
 
 function abiotic_fitness(abiotic_phenotype, species::Int, pos, model::ABM)
-  rawW = 1.0 - abiotic_distance(abiotic_phenotype, return_opt_phenotype(species, pos, model), variance)
+  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])
   return rawW
 end
@@ -101,12 +101,14 @@ function interact!(ag1::Ind, ag2::Ind, model::ABM)
 
   if sp1 != sp2 && model.food_sources[sp1, sp2] > 0  # predation
     d = phenotypic_distance(ag1, ag2, model)
-    if rand() < d
+    prob = interaction_power(ag1, ag2, d, model)
+    if rand() < prob
       eat!(ag1, ag2, model)
     end
   elseif sp1 != sp2 && model.food_sources[sp2, sp1] > 0 # predation
     d = phenotypic_distance(ag1, ag2, model)
-    if rand() < d
+    prob = interaction_power(ag1, ag2, d, model)
+    if rand() < prob
       eat!(ag2, ag1, model)
     end
   else # interaction
diff --git a/src/migration.jl b/src/migration.jl
index 67ded4e..e0449a6 100644
--- a/src/migration.jl
+++ b/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, variance)
+  abiotic_distance(phenotype, optimal, model.abiotic_variances[agent.species])
 end
 
 function pick_site(agent, sites, nsites, model)
diff --git a/src/reproduction.jl b/src/reproduction.jl
index fdc2a22..56fae6c 100644
--- a/src/reproduction.jl
+++ b/src/reproduction.jl
@@ -87,11 +87,15 @@ end
 
 """
 Sexual reproduction for diploid individuals.
-
-# TODO: Currently this is assortative mating. Add a parameter to allow nonassortative, and disassortative mating.
+Whether it is assortative or dissortative depends on the sign of the corresponding element in the interaction matrix. 
 """
 function reproduce!(ag1::Ind, ag2::Ind, model::ABM)
-  reproduction_success = 1.0 - phenotypic_distance(ag1, ag2, model)
+  if model.mating_schemes[ag1.species] == 0.0
+    reproduction_success = 1.0
+  else
+    d = phenotypic_distance(ag1.biotic_phenotype, ag2.biotic_phenotype, model.mating_schemes[ag1.species], model.biotic_variances[ag1.species])
+    reproduction_success = interaction_power(ag1, ag2, d, model)
+  end
   growth_rate = model.growthrates[ag1.species]
   nchildren = rand(Poisson(reproduction_success * growth_rate))
   for c in 1:nchildren
diff --git a/src/simulation.jl b/src/simulation.jl
index ece89a5..046c838 100644
--- a/src/simulation.jl
+++ b/src/simulation.jl
@@ -52,10 +52,11 @@ struct Params{F<:AbstractFloat,I<:Int,N<:AbstractString}
   bottlenecks::Vector{Vector{Matrix{F}}}
   repro_start::Vector{Int}
   repro_end::Vector{Int}
+  biotic_variances::Vector{F}
+  abiotic_variances::Vector{F}
+  mating_schemes::Vector{Int}
 end
 
-const variance = 1.0
-
 """
     model_initiation(dd)
 
@@ -100,7 +101,7 @@ function model_initiation(dd)
         if model.ploidy[sp] == 2
           sex = rand((true, false))
         end
-        interaction_history = MVector{model.nspecies,Int}(fill(0, model.nspecies))
+        interaction_history = MVector{model.nspecies,Int}(fill(-1, model.nspecies))
         initial_energy = model.initial_energy[sp]
         add_agent!(model.nodes[pos], model, sp, biotic_ph, abiotic_ph, epistasisMat[sp], pleiotropyMat[sp], expressionArrays[sp], 0, sex, interaction_history, initial_energy, W, true)
       end
@@ -167,6 +168,9 @@ function create_properties(dd)
   repro_start = [allspecies[i][:reproduction_start_age] for i in 1:nspecies]
   repro_end = [allspecies[i][:reproduction_end_age] for i in 1:nspecies]
   resources = dd[:resources][1]
+  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]
 
   # reshape single value matrices to (1,1)
   if length(resources) == 1
@@ -175,7 +179,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)
+  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)
 
   return properties, (epistasisMatS, pleiotropyMatS, expressionArraysS)
 end
@@ -206,13 +210,13 @@ function agent_step!(agent::Ind, model::ABM)
   burn_energy!(agent)
   # consume basic energy if agent can
   consume_food!(agent, model)
+  # interact with other species
+  interact!(agent, model)
   # Kill the agent if it doesn't have energy
   if agent.isalive && agent.energy < 0
     remove_agent!(agent, model)
     return
   end
-  # interact with other species
-  interact!(agent, model)
   # survive!(agent, model)
   if !agent.isalive
     return
@@ -257,7 +261,7 @@ function survive!(agent::Ind, model::ABM)
 end
 
 """
-Kills the agent by chance given its fitness `W`
+Kills the agent by chance given its fitness `W`.
 """
 function abiotic_survive!(agent::Ind, model::ABM)
   if agent.isalive && rand() > agent.W
diff --git a/test/params.jl b/test/params.jl
index ea47efb..037d525 100644
--- a/test/params.jl
+++ b/test/params.jl
@@ -1,11 +1,11 @@
-generations = 14
+generations = 24
 space = (2, 2)
 
 function env_resources(time::Int)
   if time < 7
-    return [200 158; 183 190]
+    return [2000 1580; 1830 1900]
   else
-    return [200 158; 183 190] .- 50
+    return [2000 1580; 1830 1900] .- 50
   end
 end
 
@@ -29,7 +29,9 @@ species1 = Dict(
   :pleiotropy_matrix => Bool[1 1 0 0 0 0 0 1 0 1 0 1 0 1; 1 0 1 1 1 0 1 1 1 0 0 1 1 0; 1 0 1 0 1 0 0 0 0 0 1 0 1 0; 0 0 0 0 1 0 0 1 1 1 1 0 1 0],
   :growth_rate => 0.8 , # 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.2,
+  :selection_coefficient => 0.02,
+  :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.9, 0.0, 0.0],  # for gene expression array, pleiotropy matrix and epistasis matrix, respectively
   :mutation_magnitudes => [0.05, 0.0, 0.01], # same as above
   :N => [1000, 0, 0, 0],  # number of individuals per site at time 0
@@ -39,6 +41,7 @@ species1 = Dict(
   :age => 5,  # max age
   :reproduction_start_age => 1,
   :reproduction_end_age => 5,
+  :mating_scheme => 0, # 0 means number of children between a pair is independent of the phenotype of the pair. 1 means the more similar they are, the more children they will have (assortative mating). -1 means the more dissimilar they are, the more children they will have (dissortative mating).
   :recombination => 1, # Mean of a Poisson distributions for number of crossing overs
   :initial_energy => 0, # A parameter for parental care of infants. Values more than 0 indicate that newly born individuals can survive for a number of times without requiring food from the environment/other species. The consumption rate (i.e. how many generations this initial energy suffices) is determined by the sum of the corresponding rows in "food sources"
   :bottleneck_function => [fill(0.0, space...) for t in 0:generations]  # an array of arrays with probablity of external death at each site and generation.
@@ -59,7 +62,9 @@ species2 = Dict(
   :pleiotropy_matrix => Bool[1 0 0 0 1 1 0 0; 1 0 0 1 1 0 1 1; 0 1 1 1 1 1 0 1; 1 0 1 0 0 1 0 1; 1 1 0 1 1 1 1 1],
   :growth_rate => 1.2,
   :expression_array => [0.24923147816626035, 0.7155732641738595, 0.9655184311211502, 0.8638149724268844, 0.5075272565823061, 0.9189652626508431, 0.7897640036022151, 0.17091233765481717],
-  :selection_coefficient => 0.5,
+  :selection_coefficient => 0.05,
+  :abiotic_variance => 1.0,
+  :biotic_variance => 10.0,
   :mutation_probabilities => [0.9, 0.0, 0.0],
   :mutation_magnitudes => [0.05, 0.0, 0.01],
   :N => [1000, 0, 0, 0],
@@ -68,6 +73,7 @@ species2 = Dict(
   :age => 3,
   :reproduction_start_age => 1,
   :reproduction_end_age => 3,
+  :mating_scheme => 0,
   :recombination => 1,
   :initial_energy => 0,
   :bottleneck_function => [fill(0.0, space...) for t in 0:generations]
@@ -85,7 +91,7 @@ model_parameters = Dict(
   :metric => "chebyshev",  # how many neighbors a space site has. "chebyshev" metric means that the r-neighborhood of a position are all positions within the hypercube having side length of 2*floor(r) and being centered in the origin position. "euclidean" metric means that the r-neighborhood of a position are all positions whose cartesian indices have Euclidean distance ≤ r from the cartesian index of the given position.
   :periodic => false,  # whether boundaries of the space are connected
   :resources => [env_resources(x) for x in 0:generations],  # a function that returns a vector of integers for the available resources per site at a given time. It accepts only the model object as the input.
-  :interactions => [0.1 0.0; 0.0 -1.0],  # How individuals from different species interact. value  is strength of interaction (between 0 and 1). Sign is the direction of interaction where positive means similar individuals interact more strongly and negative is dissimilar ones tend to interact more. 
-  :food_sources => [1.0 0.0; 0.5 0.0],  # What each species feeds on (consumption rate). Has priority over interactions. Non-zero diagonal means the food resource is from the environment. It will be read from rows (species order) to columns (species order).
+  :interactions => [0.1 1.0; 1.0 -1.0],  # How individuals from different species interact. value  is probability of interaction (between 0 and 1). Sign is the direction of interaction where positive means similar individuals interact more strongly and negative is dissimilar ones tend to interact more. If you want full interaction, use 1 or -1 for nondiagonals. Diagonals are only used for reproduction (in case of diploids) and competition/cooperation (in case of same sexes in diploids, and and interaction between haploids of the same species).
+  :food_sources => [1.0 0.0; 1.5 0.0],  # What each species feeds on (consumption rate). Non-zero diagonal means the food resource is from the environment. It will be read from rows (species order) to columns (species order).
   :seed => 3
 )
\ No newline at end of file

From fc2f36b1a21606f647ff6220dcefdd907918d567 Mon Sep 17 00:00:00 2001
From: "Ali R. Vahdati" <arv.kavir@gmail.com>
Date: Mon, 29 Aug 2022 22:58:53 +0200
Subject: [PATCH 3/9] bug fixes

---
 src/simulation.jl |  6 +++---
 test/params.jl    | 18 +++++++++---------
 2 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/src/simulation.jl b/src/simulation.jl
index 046c838..04fd7d5 100644
--- a/src/simulation.jl
+++ b/src/simulation.jl
@@ -103,7 +103,7 @@ function model_initiation(dd)
         end
         interaction_history = MVector{model.nspecies,Int}(fill(-1, model.nspecies))
         initial_energy = model.initial_energy[sp]
-        add_agent!(model.nodes[pos], model, sp, biotic_ph, abiotic_ph, epistasisMat[sp], pleiotropyMat[sp], expressionArrays[sp], 0, sex, interaction_history, initial_energy, W, true)
+        add_agent!(model.nodes[pos], model, sp, biotic_ph, abiotic_ph, epistasisMat[sp], pleiotropyMat[sp], expressionArrays[sp], 1, sex, interaction_history, initial_energy, W, true)
       end
     end
   end
@@ -199,8 +199,6 @@ function model_step!(model::ABM)
 end
 
 function agent_step!(agent::Ind, model::ABM)
-  # update age
-  agent.age += 1
   # abiotic survive
   abiotic_survive!(agent, model)  # the agent first survives then feeds, reproduces, and interacts.
   if !agent.isalive
@@ -236,6 +234,8 @@ function agent_step!(agent::Ind, model::ABM)
       remove_agent!(agent, model)
     end
   end
+  # update age
+  agent.age += 1
 end
 
 function remove_agent!(agent, model)
diff --git a/test/params.jl b/test/params.jl
index 037d525..e4cad33 100644
--- a/test/params.jl
+++ b/test/params.jl
@@ -1,11 +1,11 @@
-generations = 24
+generations = 44
 space = (2, 2)
 
 function env_resources(time::Int)
   if time < 7
-    return [2000 1580; 1830 1900]
+    return [2000 1980; 1830 1900]
   else
-    return [2000 1580; 1830 1900] .- 50
+    return [2000 1980; 1830 1900] .+ 50
   end
 end
 
@@ -27,7 +27,7 @@ species1 = Dict(
   :epistasis_matrix => [1.0 0.43 -0.41 0.38 0.48 -0.43 -0.1 -0.08 -0.09 -0.5 0.41 0.44 -0.21 -0.12; 0.34 1.0 -0.19 -0.36 0.38 -0.28 0.24 -0.22 0.12 0.12 -0.12 -0.39 0.21 0.26; 0.05 0.27 1.0 0.04 0.01 -0.14 0.3 -0.28 0.43 -0.13 0.2 -0.02 0.25 -0.39; -0.12 0.33 -0.48 1.0 -0.4 -0.48 -0.22 -0.36 -0.24 -0.07 -0.12 -0.49 -0.37 0.27; 0.25 0.25 -0.14 0.49 1.0 0.28 -0.34 -0.49 0.45 -0.14 0.26 -0.13 -0.44 -0.17; -0.47 0.19 -0.24 0.41 0.08 1.0 0.11 0.03 0.15 0.49 0.04 0.41 -0.19 0.13; 0.37 0.09 -0.11 0.4 0.42 0.45 1.0 -0.01 -0.47 0.07 0.5 0.44 -0.18 -0.2; -0.32 0.15 0.4 -0.24 -0.21 0.5 0.22 1.0 -0.33 0.48 -0.49 0.07 0.5 -0.07; 0.02 -0.16 0.33 0.48 -0.42 0.39 0.2 -0.11 1.0 0.46 -0.06 0.22 -0.3 0.31; 0.41 -0.18 -0.16 -0.4 0.01 0.04 0.07 0.2 -0.37 1.0 -0.33 0.49 0.05 -0.42; 0.03 0.25 0.14 -0.36 0.28 -0.18 0.09 -0.2 0.46 -0.48 1.0 -0.21 0.41 -0.46; 0.0 0.44 -0.34 -0.42 0.37 -0.04 0.43 -0.25 0.21 0.19 0.29 1.0 -0.02 0.06; 0.46 -0.1 0.14 -0.22 -0.26 0.13 -0.5 -0.41 -0.31 0.0 -0.15 0.29 1.0 0.17; -0.22 0.21 0.46 -0.01 -0.35 -0.11 0.25 -0.03 0.18 -0.38 -0.4 -0.28 0.05 1.0],
   # nphenotype x (ngenes x ploidy)
   :pleiotropy_matrix => Bool[1 1 0 0 0 0 0 1 0 1 0 1 0 1; 1 0 1 1 1 0 1 1 1 0 0 1 1 0; 1 0 1 0 1 0 0 0 0 0 1 0 1 0; 0 0 0 0 1 0 0 1 1 1 1 0 1 0],
-  :growth_rate => 0.8 , # max number of offsprings per mating mean of a Poisson
+  :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,
   :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.
@@ -37,7 +37,7 @@ species1 = Dict(
   :N => [1000, 0, 0, 0],  # number of individuals per site at time 0
   :environmental_noise => 0.0,  # variance of a normal distribution with mean 0
   # each row is the optimal phenotypes for each site for all abiotic traits. There are as many element as number of sites times number of abiotic traits. The first N elements are for the first abiotic trait, where N is the number of sites, and so on.
-  :optimal_phenotypes => [fill([1.5 for p in 1:1], space...) for t in 0:generations], # optimal phenotypes per site and time for each phenotype
+  :optimal_phenotypes => [fill([1.4 for p in 1:1], space...) for t in 0:generations], # optimal phenotypes per site and time for each phenotype
   :age => 5,  # max age
   :reproduction_start_age => 1,
   :reproduction_end_age => 5,
@@ -60,14 +60,14 @@ species2 = Dict(
   :ploidy => 1,
   :epistasis_matrix => [1.0 -0.01 -0.01 -0.34 -0.42 0.18 -0.09 0.13; 0.28 1.0 -0.21 -0.11 0.12 -0.46 0.21 -0.39; 0.0 0.28 1.0 0.1 0.09 0.3 0.1 0.17; 0.0 -0.42 0.05 1.0 -0.11 -0.07 -0.27 0.43; 0.31 0.47 -0.42 -0.34 1.0 -0.4 0.16 0.11; -0.19 0.29 -0.33 -0.49 0.17 1.0 -0.1 -0.28; -0.43 0.38 -0.06 0.39 0.21 -0.5 1.0 -0.08; 0.26 0.27 -0.44 -0.08 0.47 -0.27 -0.27 1.0],
   :pleiotropy_matrix => Bool[1 0 0 0 1 1 0 0; 1 0 0 1 1 0 1 1; 0 1 1 1 1 1 0 1; 1 0 1 0 0 1 0 1; 1 1 0 1 1 1 1 1],
-  :growth_rate => 1.2,
+  :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,
   :abiotic_variance => 1.0,
   :biotic_variance => 10.0,
   :mutation_probabilities => [0.9, 0.0, 0.0],
   :mutation_magnitudes => [0.05, 0.0, 0.01],
-  :N => [1000, 0, 0, 0],
+  :N => [100, 0, 0, 0],
   :environmental_noise => 0.01,
   :optimal_phenotypes => [fill([1.5 for i in 1:2], space...) for t in 0:generations],
   :age => 3,
@@ -91,7 +91,7 @@ model_parameters = Dict(
   :metric => "chebyshev",  # how many neighbors a space site has. "chebyshev" metric means that the r-neighborhood of a position are all positions within the hypercube having side length of 2*floor(r) and being centered in the origin position. "euclidean" metric means that the r-neighborhood of a position are all positions whose cartesian indices have Euclidean distance ≤ r from the cartesian index of the given position.
   :periodic => false,  # whether boundaries of the space are connected
   :resources => [env_resources(x) for x in 0:generations],  # a function that returns a vector of integers for the available resources per site at a given time. It accepts only the model object as the input.
-  :interactions => [0.1 1.0; 1.0 -1.0],  # How individuals from different species interact. value  is probability of interaction (between 0 and 1). Sign is the direction of interaction where positive means similar individuals interact more strongly and negative is dissimilar ones tend to interact more. If you want full interaction, use 1 or -1 for nondiagonals. Diagonals are only used for reproduction (in case of diploids) and competition/cooperation (in case of same sexes in diploids, and and interaction between haploids of the same species).
-  :food_sources => [1.0 0.0; 1.5 0.0],  # What each species feeds on (consumption rate). Non-zero diagonal means the food resource is from the environment. It will be read from rows (species order) to columns (species order).
+  :interactions => [0.1 1.0; 1.0 -1.0],  # How individuals from different species interact. value  is probability of interaction (between 0 and 1). Sign is the direction of interaction where positive means similar individuals interact more strongly and negative is dissimilar ones tend to interact more. If you want full interaction, use 1 or -1 for nondiagonals. Diagonals are used for competition/cooperation.
+  :food_sources => [1.0 0.0; 1.7 0.0],  # What each species feeds on (consumption rate). Non-zero diagonal means the food resource is from the environment. It will be read from rows (species order) to columns (species order).
   :seed => 3
 )
\ No newline at end of file

From 1caa432c04d0950c636d5d69c8e9b3816b755212 Mon Sep 17 00:00:00 2001
From: "Ali R. Vahdati" <arv.kavir@gmail.com>
Date: Tue, 30 Aug 2022 13:27:40 +0200
Subject: [PATCH 4/9] Reorder operations and bug fixes

---
 Project.toml        |  2 +-
 src/interactions.jl |  4 +++-
 src/reproduction.jl |  6 ++++--
 src/simulation.jl   | 10 ++++++----
 test/params.jl      |  8 ++++----
 5 files changed, 18 insertions(+), 12 deletions(-)

diff --git a/Project.toml b/Project.toml
index d5f1da7..20b45a8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "EvoDynamics"
 uuid = "c8768967-421d-4a02-8523-37736f3dbe06"
 authors = ["Ali R. Vahdati <arv.kavir@gmail.com>", "Carlos Melian"]
-version = "0.14.0"
+version = "0.15.0"
 
 [deps]
 Agents = "46ada45e-f475-11e8-01d0-f70cc89e6671"
diff --git a/src/interactions.jl b/src/interactions.jl
index 5d10dfa..86d4e47 100644
--- a/src/interactions.jl
+++ b/src/interactions.jl
@@ -113,7 +113,9 @@ function interact!(ag1::Ind, ag2::Ind, model::ABM)
     end
   else # interaction
     if sp1 == sp2 && ag1.sex != ag2.sex && model.ploidy[sp1] == 2 && in_reproduction_age(ag1, model) && in_reproduction_age(ag2, model)  # reproduce
-      reproduce!(ag1::Ind, ag2::Ind, model::ABM)
+      # reproduce!(ag1, ag2, model)
+      ag1.mate = ag2.id
+      ag1.time_met_other_sex = model.step[1]
     else
       ix_value1 = model.interactions[sp1, sp2]
       ix_value2 = model.interactions[sp2, sp1]
diff --git a/src/reproduction.jl b/src/reproduction.jl
index 56fae6c..2608fd3 100644
--- a/src/reproduction.jl
+++ b/src/reproduction.jl
@@ -81,7 +81,7 @@ function create_one_offspring(ag1::Ind, ag2::Ind, model::ABM)
   W = adjust_fitness(W, ag1.species, model)
   initial_energy = model.initial_energy[species]
 
-  offspring = add_agent!(ag1.pos, model, ag1.species, bph, abph, episMat, pleioMat, q, 0, sex, interaction_history, initial_energy, W, true)
+  offspring = add_agent!(ag1.pos, model, ag1.species, bph, abph, episMat, pleioMat, q, 1, sex, interaction_history, initial_energy, W, true, 0, -1)
   return offspring
 end
 
@@ -117,9 +117,11 @@ function reproduce!(agent::Ind, model::ABM)
     interaction_history = deepcopy(agent.interaction_history)
     interaction_history[1:end] .= 0
     for c in 1:nchildren
-      offspring = add_agent!(agent.pos, model, agent.species, agent.biotic_phenotype, agent.abiotic_phenotype, agent.epistasisMat, agent.pleiotropyMat, agent.q, 0, agent.sex, interaction_history, model.initial_energy[agent.species], agent.W, true)
+      offspring = add_agent!(agent.pos, model, agent.species, agent.biotic_phenotype, agent.abiotic_phenotype, agent.epistasisMat, agent.pleiotropyMat, agent.q, 1, agent.sex, interaction_history, model.initial_energy[agent.species], agent.W, true, 0, -1)
       mutate!(offspring, model)
     end
+  elseif model.ploidy[agent.species] == 2 && agent.interaction_history[agent.species] == model.step[1] && agent.time_met_other_sex == model.step[1]
+    reproduce!(agent, model[agent.mate], model)
   end
 end
 
diff --git a/src/simulation.jl b/src/simulation.jl
index 04fd7d5..ca61fd3 100644
--- a/src/simulation.jl
+++ b/src/simulation.jl
@@ -17,6 +17,8 @@ mutable struct Ind{B<:AbstractFloat,C<:AbstractArray,D<:AbstractArray,E<:Abstrac
   energy::B  # determines whether need to feed (energy=0) or not.
   W::B  # survival probability
   isalive::Bool
+  mate::Int  # id of mate. Only relevant for diploids
+  time_met_other_sex::Int
 end
 
 struct Params{F<:AbstractFloat,I<:Int,N<:AbstractString}
@@ -103,7 +105,7 @@ function model_initiation(dd)
         end
         interaction_history = MVector{model.nspecies,Int}(fill(-1, model.nspecies))
         initial_energy = model.initial_energy[sp]
-        add_agent!(model.nodes[pos], model, sp, biotic_ph, abiotic_ph, epistasisMat[sp], pleiotropyMat[sp], expressionArrays[sp], 1, sex, interaction_history, initial_energy, W, true)
+        add_agent!(model.nodes[pos], model, sp, biotic_ph, abiotic_ph, epistasisMat[sp], pleiotropyMat[sp], expressionArrays[sp], 1, sex, interaction_history, initial_energy, W, true, 0, -1)
       end
     end
   end
@@ -208,7 +210,7 @@ function agent_step!(agent::Ind, model::ABM)
   burn_energy!(agent)
   # consume basic energy if agent can
   consume_food!(agent, model)
-  # interact with other species
+  # interact with other species: predation, cooperation, competition.
   interact!(agent, model)
   # Kill the agent if it doesn't have energy
   if agent.isalive && agent.energy < 0
@@ -219,10 +221,10 @@ function agent_step!(agent::Ind, model::ABM)
   if !agent.isalive
     return
   end
+  # reproduction for both haploids and diploids
+  reproduce!(agent, model)
   # migrate
   migrate!(agent, model)
-  # reproduction for the haploid
-  reproduce!(agent, model)
 
   if agent.isalive && agent.age ≥ model.max_ages[agent.species]
     remove_agent!(agent, model)
diff --git a/test/params.jl b/test/params.jl
index e4cad33..446b1d8 100644
--- a/test/params.jl
+++ b/test/params.jl
@@ -1,4 +1,4 @@
-generations = 44
+generations = 5
 space = (2, 2)
 
 function env_resources(time::Int)
@@ -38,9 +38,9 @@ species1 = Dict(
   :environmental_noise => 0.0,  # variance of a normal distribution with mean 0
   # each row is the optimal phenotypes for each site for all abiotic traits. There are as many element as number of sites times number of abiotic traits. The first N elements are for the first abiotic trait, where N is the number of sites, and so on.
   :optimal_phenotypes => [fill([1.4 for p in 1:1], space...) for t in 0:generations], # optimal phenotypes per site and time for each phenotype
-  :age => 5,  # max age
+  :age => 4,  # max age
   :reproduction_start_age => 1,
-  :reproduction_end_age => 5,
+  :reproduction_end_age => 4,
   :mating_scheme => 0, # 0 means number of children between a pair is independent of the phenotype of the pair. 1 means the more similar they are, the more children they will have (assortative mating). -1 means the more dissimilar they are, the more children they will have (dissortative mating).
   :recombination => 1, # Mean of a Poisson distributions for number of crossing overs
   :initial_energy => 0, # A parameter for parental care of infants. Values more than 0 indicate that newly born individuals can survive for a number of times without requiring food from the environment/other species. The consumption rate (i.e. how many generations this initial energy suffices) is determined by the sum of the corresponding rows in "food sources"
@@ -91,7 +91,7 @@ model_parameters = Dict(
   :metric => "chebyshev",  # how many neighbors a space site has. "chebyshev" metric means that the r-neighborhood of a position are all positions within the hypercube having side length of 2*floor(r) and being centered in the origin position. "euclidean" metric means that the r-neighborhood of a position are all positions whose cartesian indices have Euclidean distance ≤ r from the cartesian index of the given position.
   :periodic => false,  # whether boundaries of the space are connected
   :resources => [env_resources(x) for x in 0:generations],  # a function that returns a vector of integers for the available resources per site at a given time. It accepts only the model object as the input.
-  :interactions => [0.1 1.0; 1.0 -1.0],  # How individuals from different species interact. value  is probability of interaction (between 0 and 1). Sign is the direction of interaction where positive means similar individuals interact more strongly and negative is dissimilar ones tend to interact more. If you want full interaction, use 1 or -1 for nondiagonals. Diagonals are used for competition/cooperation.
+  :interactions => [0.1 1.0; 1.0 1.0],  # How individuals from different species interact. value  is probability of interaction (between 0 and 1). Sign is the direction of interaction where positive means similar individuals interact more strongly and negative is dissimilar ones tend to interact more. If you want full interaction, use 1 or -1 for nondiagonals. Diagonals are used for competition/cooperation.
   :food_sources => [1.0 0.0; 1.7 0.0],  # What each species feeds on (consumption rate). Non-zero diagonal means the food resource is from the environment. It will be read from rows (species order) to columns (species order).
   :seed => 3
 )
\ No newline at end of file

From b84758fd73037664bd21ab8b6c28226f1e94240a Mon Sep 17 00:00:00 2001
From: "Ali R. Vahdati" <arv.kavir@gmail.com>
Date: Tue, 30 Aug 2022 17:45:53 +0200
Subject: [PATCH 5/9] update random number and examples

---
 examples/example1.jl   |  40 ++++++-------
 examples/paramfile1.jl |  31 ++++------
 examples/paramfile2.jl | 131 ++++++++++++++++++++++++-----------------
 src/interactions.jl    |  12 ++--
 src/migration.jl       |   2 +-
 src/reproduction.jl    |  14 ++---
 src/simulation.jl      |  23 ++++----
 test/params.jl         |   2 +-
 8 files changed, 132 insertions(+), 123 deletions(-)

diff --git a/examples/example1.jl b/examples/example1.jl
index 334963b..bdc2635 100644
--- a/examples/example1.jl
+++ b/examples/example1.jl
@@ -8,18 +8,10 @@ using EvoDynamics
 # A simple one-species model with no spatial structure. Model parameters are in a .jl file as follows:
 
 # ```julia
-# ## 1. Functions
-# function bn(agent::AbstractAgent, model::ABM)
-#   return false
-# end
+# generations = 20
+# space = (1, 1)
 
-# function optphens(site::Tuple{Int,Int}, model::ABM)
-#   return [1.5]
-# end
-
-# env_resources(time::Int) = [200]
-
-# ## 2. Species parameters
+# ## 1. Species parameters
 
 # species1 = Dict(
 #   :name => "a",
@@ -30,43 +22,47 @@ using EvoDynamics
 #   :migration_phenotype => 0,
 #   :migration_threshold => 3.5,
 #   :vision_radius => 0,
-#   :check_fraction => 0,
+#   :check_fraction => 0.0,
 #   :ploidy => 1,
 #   :epistasis_matrix => [1.0 0.0; 0.0 1.0],
 #   :pleiotropy_matrix => Bool[1 0; 0 1],
 #   :growth_rate => 1.0,
 #   :expression_array => [0.28, 0.46],
-#   :selection_coefficient => 0.5,
+#   :selection_coefficient => 0.1,
 #   :mutation_probabilities => [0.9, 0.0, 0.0],
 #   :mutation_magnitudes => [0.05, 0.0, 0.0],
 #   :N => [100],
 #   :environmental_noise => 0.01,
-#   :optimal_phenotypes => optphens,
-#   :bottleneck_function => bn,
+#   :optimal_phenotypes => [fill([1.5 for p in 1:1], space...) for t in 0:generations],
+#   :bottleneck_function => [fill(0.0, space...) for t in 0:generations],
 #   :age => 2,
 #   :recombination => 0,
-#   :initial_energy => 0
+#   :initial_energy => 0,
 #   :reproduction_start_age => 1,
-#   :reproduction_end_age => 1,
+#   :reproduction_end_age => 2,
+#   :abiotic_variance => 1.0,
+#   :biotic_variance => 1.0,
+#   :mating_scheme => 1
 # )
 
-# ## 3. Model parameters
+# ## 2. Model parameters
 
 # #NB this dict should be called model_parameters
 # model_parameters = Dict(
 #   :species => [species1],
-#   :generations => 5,
-#   :space => nothing,
+#   :generations => generations,
+#   :space => space,
 #   :metric => "chebyshev",
 #   :periodic => false,
-#   :resources => env_resources,
+#   :resources => [fill(200, 1, 1) for i in 0:generations],
 #   :interactions => [-0.1],
 #   :food_sources => [1.0],
 #   :seed => nothing
 # )
 # ```
 
-param_file = "../../examples/paramfile1.jl" #hide
+param_file = "../../examples/paramfile1.jl"
+
 agentdata, modeldata, model = runmodel(param_file);
 
 modeldata
\ No newline at end of file
diff --git a/examples/paramfile1.jl b/examples/paramfile1.jl
index 8b429a4..f0ec348 100644
--- a/examples/paramfile1.jl
+++ b/examples/paramfile1.jl
@@ -1,15 +1,7 @@
-## 1. Functions
-function bn(agent::AbstractAgent, model::ABM)
- return false
-end
+generations = 20
+space = (1, 1)
 
-function optphens(site::Tuple{Int, Int}, model::ABM)
-  return [1.5]
-end
-
-env_resources(time::Int) = [200]
-
-## 2. Species parameters
+## 1. Species parameters
 
 species1 = Dict(
   :name => "a",
@@ -26,30 +18,33 @@ species1 = Dict(
   :pleiotropy_matrix => Bool[1 0; 0 1],
   :growth_rate => 1.0,
   :expression_array => [0.28, 0.46],
-  :selection_coefficient => 0.5,
+  :selection_coefficient => 0.1,
   :mutation_probabilities => [0.9, 0.0, 0.0],
   :mutation_magnitudes => [0.05, 0.0, 0.0],
   :N => [100],
   :environmental_noise => 0.01,
-  :optimal_phenotypes => optphens,
-  :bottleneck_function => bn,
+  :optimal_phenotypes => [fill([1.5 for p in 1:1], space...) for t in 0:generations],
+  :bottleneck_function => [fill(0.0, space...) for t in 0:generations],
   :age => 2,
   :recombination => 0,
   :initial_energy => 0,
   :reproduction_start_age => 1,
   :reproduction_end_age => 2,
+  :abiotic_variance => 1.0,
+  :biotic_variance => 1.0,
+  :mating_scheme => 1
 )
 
-## 3. Model parameters
+## 2. Model parameters
 
 #NB this dict should be called model_parameters
 model_parameters = Dict(
   :species => [species1],
-  :generations => 5,
-  :space => nothing,
+  :generations => generations,
+  :space => space,
   :metric => "chebyshev",
   :periodic => false ,
-  :resources => env_resources,
+  :resources => [fill(200, 1,1) for i in 0:generations],
   :interactions => [-0.1],
   :food_sources => [1.0],
   :seed => nothing
diff --git a/examples/paramfile2.jl b/examples/paramfile2.jl
index 17a222b..87daf01 100644
--- a/examples/paramfile2.jl
+++ b/examples/paramfile2.jl
@@ -1,34 +1,7 @@
-## 1. Functions
+generations = 14
+space = (2, 2)
 
-function bn(agent::AbstractAgent, model::ABM)
- return false
-end
-
-function optphens2(site::Tuple{Int,Int}, model::ABM)
-  opt1 = [0.8214794136627462, 0.49335401288317304, 0.3435556249656141, 0.25050033180075226]
-  opt2 = opt1 .+ 0.1
-  ss = EvoDynamics.coord2vertex(site, model)
-  if model.step[1] > 6
-    return [opt1[ss]]
-  else
-    return [opt2[ss]]
-  end
-end
-
-function optphens3(site::Tuple{Int,Int}, model::ABM)
-  phenotypic_values = [[0.7614758101208934 2.3911353663016586; 2.2091361414343313 1.7858661540288683; 0.7920974352892358 0.7630236717263839; 2.587205421882114 2.311211631439866],
-    [0.6437641305315445 2.097629262577973; 0.9954545836033715 1.1314248848669466; 2.469792530385348 1.490299526620522; 1.6158867433451882 0.2566056477862022]]
-  agent_site = (site[2] - 1) * size(model.space)[2] + (site[1])
-  if model.step[1] > 7
-    return phenotypic_values[1][agent_site, :]
-  else
-    return phenotypic_values[2][agent_site, :]
-  end
-end
-
-env_resources2(time::Int) = [200 220; 183 190]
-
-## 2. Species parameters
+## 1. Species parameters
 
 species1 = Dict(
   :name => "a",
@@ -37,26 +10,49 @@ species1 = Dict(
   :abiotic_phenotypes => [1],
   :biotic_phenotypes => [2, 3],
   :migration_phenotype => 4,
-  :migration_threshold => 3.5,
+  :migration_threshold => 1.5,
   :vision_radius => 1,
   :check_fraction => 0.5,
   :ploidy => 2,
-  :epistasis_matrix => [1.0 0.43 -0.41 0.38 0.48 -0.43 -0.1 -0.08 -0.09 -0.5 0.41 0.44 -0.21 -0.12; 0.34 1.0 -0.19 -0.36 0.38 -0.28 0.24 -0.22 0.12 0.12 -0.12 -0.39 0.21 0.26; 0.05 0.27 1.0 0.04 0.01 -0.14 0.3 -0.28 0.43 -0.13 0.2 -0.02 0.25 -0.39; -0.12 0.33 -0.48 1.0 -0.4 -0.48 -0.22 -0.36 -0.24 -0.07 -0.12 -0.49 -0.37 0.27; 0.25 0.25 -0.14 0.49 1.0 0.28 -0.34 -0.49 0.45 -0.14 0.26 -0.13 -0.44 -0.17; -0.47 0.19 -0.24 0.41 0.08 1.0 0.11 0.03 0.15 0.49 0.04 0.41 -0.19 0.13; 0.37 0.09 -0.11 0.4 0.42 0.45 1.0 -0.01 -0.47 0.07 0.5 0.44 -0.18 -0.2; -0.32 0.15 0.4 -0.24 -0.21 0.5 0.22 1.0 -0.33 0.48 -0.49 0.07 0.5 -0.07; 0.02 -0.16 0.33 0.48 -0.42 0.39 0.2 -0.11 1.0 0.46 -0.06 0.22 -0.3 0.31; 0.41 -0.18 -0.16 -0.4 0.01 0.04 0.07 0.2 -0.37 1.0 -0.33 0.49 0.05 -0.42; 0.03 0.25 0.14 -0.36 0.28 -0.18 0.09 -0.2 0.46 -0.48 1.0 -0.21 0.41 -0.46; 0.0 0.44 -0.34 -0.42 0.37 -0.04 0.43 -0.25 0.21 0.19 0.29 1.0 -0.02 0.06; 0.46 -0.1 0.14 -0.22 -0.26 0.13 -0.5 -0.41 -0.31 0.0 -0.15 0.29 1.0 0.17; -0.22 0.21 0.46 -0.01 -0.35 -0.11 0.25 -0.03 0.18 -0.38 -0.4 -0.28 0.05 1.0],
-  :pleiotropy_matrix => Bool[1 1 0 0 0 0 0 1 0 1 0 1 0 1; 1 0 1 1 1 0 1 1 1 0 0 1 1 0; 1 0 1 0 1 0 0 0 0 0 1 0 1 0; 0 0 0 0 1 0 0 1 1 1 1 0 1 0],
-  :growth_rate => 0.8,
+  :epistasis_matrix => [
+    1.0 0.43 -0.41 0.38 0.48 -0.43 -0.1 -0.08 -0.09 -0.5 0.41 0.44 -0.21 -0.12;
+    0.34 1.0 -0.19 -0.36 0.38 -0.28 0.24 -0.22 0.12 0.12 -0.12 -0.39 0.21 0.26;
+    0.05 0.27 1.0 0.04 0.01 -0.14 0.3 -0.28 0.43 -0.13 0.2 -0.02 0.25 -0.39;
+    -0.12 0.33 -0.48 1.0 -0.4 -0.48 -0.22 -0.36 -0.24 -0.07 -0.12 -0.49 -0.37 0.27;
+    0.25 0.25 -0.14 0.49 1.0 0.28 -0.34 -0.49 0.45 -0.14 0.26 -0.13 -0.44 -0.17;
+    -0.47 0.19 -0.24 0.41 0.08 1.0 0.11 0.03 0.15 0.49 0.04 0.41 -0.19 0.13;
+    0.37 0.09 -0.11 0.4 0.42 0.45 1.0 -0.01 -0.47 0.07 0.5 0.44 -0.18 -0.2;
+    -0.32 0.15 0.4 -0.24 -0.21 0.5 0.22 1.0 -0.33 0.48 -0.49 0.07 0.5 -0.07;
+    0.02 -0.16 0.33 0.48 -0.42 0.39 0.2 -0.11 1.0 0.46 -0.06 0.22 -0.3 0.31;
+    0.41 -0.18 -0.16 -0.4 0.01 0.04 0.07 0.2 -0.37 1.0 -0.33 0.49 0.05 -0.42;
+    0.03 0.25 0.14 -0.36 0.28 -0.18 0.09 -0.2 0.46 -0.48 1.0 -0.21 0.41 -0.46;
+    0.0 0.44 -0.34 -0.42 0.37 -0.04 0.43 -0.25 0.21 0.19 0.29 1.0 -0.02 0.06;
+    0.46 -0.1 0.14 -0.22 -0.26 0.13 -0.5 -0.41 -0.31 0.0 -0.15 0.29 1.0 0.17;
+    -0.22 0.21 0.46 -0.01 -0.35 -0.11 0.25 -0.03 0.18 -0.38 -0.4 -0.28 0.05 1.0],
+
+  :pleiotropy_matrix => Bool[
+    1 1 0 0 0 0 0 1 0 1 0 1 0 1;
+    1 0 1 1 1 0 1 1 1 0 0 1 1 0;
+    1 0 1 0 1 0 0 0 0 0 1 0 1 0;
+    0 0 0 0 1 0 0 1 1 1 1 0 1 0],
+
+  :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.5,
+  :selection_coefficient => 0.1,
   :mutation_probabilities => [0.9, 0.0, 0.0],
   :mutation_magnitudes => [0.05, 0.0, 0.01],
   :N => [1000, 0, 0, 0],
   :environmental_noise => 0.01,
-  :optimal_phenotypes => optphens2,
-  :bottleneck_function => bn,
+  :optimal_phenotypes => [fill([1.4 for p in 1:1], space...) for t in 0:generations],
+  :bottleneck_function => [fill(0.0, space...) for t in 0:generations],
   :age => 5,
   :recombination => 1,
-  :initial_energy => 0,
-  :reproduction_start_age => 1,
+  :initial_energy => 1,
+  :reproduction_start_age => 2,
   :reproduction_end_age => 5,
+  :abiotic_variance => 1.0,
+  :biotic_variance => 1.0,
+  :mating_scheme => 0
 )
 
 species2 = Dict(
@@ -66,39 +62,64 @@ species2 = Dict(
   :abiotic_phenotypes => [1,2],
   :biotic_phenotypes => [3, 4],
   :migration_phenotype => 5,
-  :migration_threshold => 3.4,
+  :migration_threshold => 1.4,
   :vision_radius => 1,
   :check_fraction => 0.5,
   :ploidy => 1,
-  :epistasis_matrix => [1.0 -0.01 -0.01 -0.34 -0.42 0.18 -0.09 0.13; 0.28 1.0 -0.21 -0.11 0.12 -0.46 0.21 -0.39; 0.0 0.28 1.0 0.1 0.09 0.3 0.1 0.17; 0.0 -0.42 0.05 1.0 -0.11 -0.07 -0.27 0.43; 0.31 0.47 -0.42 -0.34 1.0 -0.4 0.16 0.11; -0.19 0.29 -0.33 -0.49 0.17 1.0 -0.1 -0.28; -0.43 0.38 -0.06 0.39 0.21 -0.5 1.0 -0.08; 0.26 0.27 -0.44 -0.08 0.47 -0.27 -0.27 1.0],
-  :pleiotropy_matrix => Bool[1 0 0 0 1 1 0 0; 1 0 0 1 1 0 1 1; 0 1 1 1 1 1 0 1; 1 0 1 0 0 1 0 1; 1 1 0 1 1 1 1 1],
-  :growth_rate => 1.2,
+  :epistasis_matrix => [
+    1.0 -0.01 -0.01 -0.34 -0.42 0.18 -0.09 0.13;
+    0.28 1.0 -0.21 -0.11 0.12 -0.46 0.21 -0.39;
+    0.0 0.28 1.0 0.1 0.09 0.3 0.1 0.17;
+    0.0 -0.42 0.05 1.0 -0.11 -0.07 -0.27 0.43;
+    0.31 0.47 -0.42 -0.34 1.0 -0.4 0.16 0.11;
+    -0.19 0.29 -0.33 -0.49 0.17 1.0 -0.1 -0.28;
+    -0.43 0.38 -0.06 0.39 0.21 -0.5 1.0 -0.08;
+    0.26 0.27 -0.44 -0.08 0.47 -0.27 -0.27 1.0],
+
+  :pleiotropy_matrix => Bool[
+    1 0 0 0 1 1 0 0;
+    1 0 0 1 1 0 1 1;
+    0 1 1 1 1 1 0 1;
+    1 0 1 0 0 1 0 1;
+    1 1 0 1 1 1 1 1],
+
+  :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.5,
+  :selection_coefficient => 0.1,
   :mutation_probabilities => [0.9, 0.0, 0.0],
   :mutation_magnitudes => [0.05, 0.0, 0.01],
-  :N => [1000, 0, 0, 0],
+  :N => [100, 0, 0, 0],
   :environmental_noise => 0.01,
-  :optimal_phenotypes => optphens3,
-  :bottleneck_function => bn,
+  :optimal_phenotypes => [fill([1.4, 1.6], space...) for t in 0:generations],
+  :bottleneck_function => [fill(0.0, space...) for t in 0:generations],
   :age => 3,
-  :recombination => 1,
+  :recombination => 0,
   :initial_energy => 0,
   :reproduction_start_age => 1,
   :reproduction_end_age => 3,
+  :abiotic_variance => 1.0,
+  :biotic_variance => 1.0,
+  :mating_scheme => 0
 )
 
-## 3. Model parameters
+## 2. Model parameters
 
-#NB this dict should be called model_parameters
+#NB this dictionary should be called model_parameters
 model_parameters = Dict(
   :species => [species1, species2],
-  :generations => 14,
-  :space => (2, 2),
+  :generations => generations,
+  :space => space,
   :metric => "chebyshev",
   :periodic => false,
-  :resources => env_resources2,
-  :interactions => [0.1 0.0; 0.0 -1.0],
-  :food_sources => [1.0 0.5; 0.0 0.0], 
+  :resources => [
+    [2000 2200;
+    1830 1900] for i in 0:generations],
+
+  :interactions => [0.0 1.0;
+                    1.0 0.0],
+
+  :food_sources => [1.0 0.0; 
+                    2.5 0.0],
+
   :seed => nothing
 )
\ No newline at end of file
diff --git a/src/interactions.jl b/src/interactions.jl
index 86d4e47..5b15dfe 100644
--- a/src/interactions.jl
+++ b/src/interactions.jl
@@ -3,7 +3,7 @@ function get_biotic_phenotype(species, epistasisMat, pleiotropyMat, expressionAr
   abp = model.biotic_phenotypes[species]
   x = pleiotropyMat[abp, abp] * (epistasisMat[abp, abp] * expressionArray[abp])  # phenotypic values
   d = model.E[species]
-  z = x .+ rand(d)
+  z = x .+ rand(model.rng, d)
   return z
 end
 
@@ -16,7 +16,7 @@ function get_abiotic_phenotype(species, epistasisMat, pleiotropyMat, expressionA
   if d.σ == 0 && d.μ == 0
     randd = 0.0
   else
-    randd = rand(d)
+    randd = rand(model.rng, d)
   end
   z = x .+ randd
   return z
@@ -102,13 +102,13 @@ function interact!(ag1::Ind, ag2::Ind, model::ABM)
   if sp1 != sp2 && model.food_sources[sp1, sp2] > 0  # predation
     d = phenotypic_distance(ag1, ag2, model)
     prob = interaction_power(ag1, ag2, d, model)
-    if rand() < prob
+    if rand(model.rng) < prob
       eat!(ag1, ag2, model)
     end
   elseif sp1 != sp2 && model.food_sources[sp2, sp1] > 0 # predation
-    d = phenotypic_distance(ag1, ag2, model)
-    prob = interaction_power(ag1, ag2, d, model)
-    if rand() < prob
+    d = phenotypic_distance(ag2, ag1, model)
+    prob = interaction_power(ag2, ag1, d, model)
+    if rand(model.rng) < prob
       eat!(ag2, ag1, model)
     end
   else # interaction
diff --git a/src/migration.jl b/src/migration.jl
index e0449a6..eeac3f8 100644
--- a/src/migration.jl
+++ b/src/migration.jl
@@ -8,7 +8,7 @@ function migrate!(agent::Ind, model::ABM)
 
   sites = collect(nearby_positions(agent, model, model.vision_radius[agent.species]))
   if model.check_fraction[agent.species] == 0
-    destination = rand(sites)
+    destination = rand(model.rng, sites)
   else
     nsites = length(sites)
     nsites_selected = ceil(Int, model.check_fraction[agent.species] * nsites)
diff --git a/src/reproduction.jl b/src/reproduction.jl
index 2608fd3..241697d 100644
--- a/src/reproduction.jl
+++ b/src/reproduction.jl
@@ -59,15 +59,15 @@ function create_one_offspring(ag1::Ind, ag2::Ind, model::ABM)
   if model.recombination[species] == 0
     nco1, nco2 = 0, 0
   else
-    nco1, nco2 = rand(model.recombination[species], 2)
+    nco1, nco2 = rand(model.rng, model.recombination[species], 2)
   end
 
   cross_overs1 = crossing_overs(nsites, nco1)
   cross_overs2 = crossing_overs(nsites, nco2)
-  gametes1 = create_gamete(ag1, cross_overs1, nsites, rand((true, false)))
-  gametes2 = create_gamete(ag2, cross_overs2, nsites, rand((true, false)))
+  gametes1 = create_gamete(ag1, cross_overs1, nsites, rand(model.rng, (true, false)))
+  gametes2 = create_gamete(ag2, cross_overs2, nsites, rand(model.rng, (true, false)))
 
-  sex = rand((true, false))
+  sex = rand(model.rng, (true, false))
   interaction_history = MVector{model.nspecies,Int}(fill(0, model.nspecies))
 
   # Merge gametes
@@ -97,7 +97,7 @@ function reproduce!(ag1::Ind, ag2::Ind, model::ABM)
     reproduction_success = interaction_power(ag1, ag2, d, model)
   end
   growth_rate = model.growthrates[ag1.species]
-  nchildren = rand(Poisson(reproduction_success * growth_rate))
+  nchildren = rand(model.rng, Poisson(reproduction_success * growth_rate))
   for c in 1:nchildren
     offspring = create_one_offspring(ag1, ag2, model)
     mutate!(offspring, model)
@@ -111,8 +111,8 @@ function reproduce!(agent::Ind, model::ABM)
   if model.ploidy[agent.species] == 1 && in_reproduction_age(agent, model)
     growth_rate = model.growthrates[agent.species]
     # W = agent.W >= 0.0 ? agent.W : 0.0
-    # nchildren = rand(Poisson(growth_rate * W))
-    nchildren = rand(Poisson(growth_rate))
+    # nchildren = rand(model.rng, Poisson(growth_rate * W))
+    nchildren = rand(model.rng, Poisson(growth_rate))
     nchildren == 0 && return
     interaction_history = deepcopy(agent.interaction_history)
     interaction_history[1:end] .= 0
diff --git a/src/simulation.jl b/src/simulation.jl
index ca61fd3..ca7f75a 100644
--- a/src/simulation.jl
+++ b/src/simulation.jl
@@ -101,7 +101,7 @@ function model_initiation(dd)
         W = adjust_fitness(W, sp, model)
         sex = false
         if model.ploidy[sp] == 2
-          sex = rand((true, false))
+          sex = rand(model.rng, (true, false))
         end
         interaction_history = MVector{model.nspecies,Int}(fill(-1, model.nspecies))
         initial_energy = model.initial_energy[sp]
@@ -202,7 +202,7 @@ end
 
 function agent_step!(agent::Ind, model::ABM)
   # abiotic survive
-  abiotic_survive!(agent, model)  # the agent first survives then feeds, reproduces, and interacts.
+  abiotic_survive!(agent, model)
   if !agent.isalive
     return
   end
@@ -217,7 +217,6 @@ function agent_step!(agent::Ind, model::ABM)
     remove_agent!(agent, model)
     return
   end
-  # survive!(agent, model)
   if !agent.isalive
     return
   end
@@ -232,7 +231,7 @@ function agent_step!(agent::Ind, model::ABM)
   # bottleneck
   bn_prob = model.bottlenecks[agent.species][model.step[1]+1][agent.pos...]
   if agent.isalive && bn_prob > 0.0
-    if rand() < bn_prob
+    if rand(model.rng) < bn_prob
       remove_agent!(agent, model)
     end
   end
@@ -257,8 +256,6 @@ function survive!(agent::Ind, model::ABM)
     remove_agent!(agent, model)
   elseif agent.age ≥ model.max_ages[agent.species]
     remove_agent!(agent, model)
-    # elseif rand() > agent.W 
-    #   remove_agent!(agent, model)
   end
 end
 
@@ -266,7 +263,7 @@ end
 Kills the agent by chance given its fitness `W`.
 """
 function abiotic_survive!(agent::Ind, model::ABM)
-  if agent.isalive && rand() > agent.W
+  if agent.isalive && rand(model.rng) > agent.W
     remove_agent!(agent, model)
   end
 end
@@ -299,19 +296,19 @@ end
 function mutate!(agent::Ind, model::ABM)
   mutated = false
   # mutate gene expression
-  if rand(model.mutProbs[agent.species][1])
-    agent.q .+= rand(model.mutMagnitudes[agent.species][1], model.ngenes[agent.species] * model.ploidy[agent.species])
+  if rand(model.rng, model.mutProbs[agent.species][1])
+    agent.q .+= rand(model.rng, model.mutMagnitudes[agent.species][1], model.ngenes[agent.species] * model.ploidy[agent.species])
     mutated = true
   end
   # mutate pleiotropy matrix
-  if rand(model.mutProbs[agent.species][2])
-    randnumbers = rand(model.mutMagnitudes[agent.species][2], size(agent.pleiotropyMat))
+  if rand(model.rng, model.mutProbs[agent.species][2])
+    randnumbers = rand(model.rng, model.rng, model.mutMagnitudes[agent.species][2], size(agent.pleiotropyMat))
     agent.pleiotropyMat[randnumbers] .= .!agent.pleiotropyMat[randnumbers]
     mutated = true
   end
   # mutate epistasis matrix
-  if rand(model.mutProbs[agent.species][3])
-    agent.epistasisMat .+= rand(model.mutMagnitudes[agent.species][3], size(agent.epistasisMat))
+  if rand(model.rng, model.mutProbs[agent.species][3])
+    agent.epistasisMat .+= rand(model.rng, model.mutMagnitudes[agent.species][3], size(agent.epistasisMat))
     mutated = true
   end
   # update biotic and abiotic phenotypes and W
diff --git a/test/params.jl b/test/params.jl
index 446b1d8..e7a726b 100644
--- a/test/params.jl
+++ b/test/params.jl
@@ -43,7 +43,7 @@ species1 = Dict(
   :reproduction_end_age => 4,
   :mating_scheme => 0, # 0 means number of children between a pair is independent of the phenotype of the pair. 1 means the more similar they are, the more children they will have (assortative mating). -1 means the more dissimilar they are, the more children they will have (dissortative mating).
   :recombination => 1, # Mean of a Poisson distributions for number of crossing overs
-  :initial_energy => 0, # A parameter for parental care of infants. Values more than 0 indicate that newly born individuals can survive for a number of times without requiring food from the environment/other species. The consumption rate (i.e. how many generations this initial energy suffices) is determined by the sum of the corresponding rows in "food sources"
+  :initial_energy => 0, # A parameter for parental care of infants. Values more than 0 indicate that newly born individuals can survive for a number of times without requiring food from the environment/other species. The consumption rate (i.e. how many generations this initial energy suffices) is determined by the sum of the corresponding rows in "food sources". Note that having initial energy larger than zero can lead to infinite population growth because agents without food can reproduce and their offsprings also reproduce without food. Use with care, for example, when start age of reproduction is larger than 1.
   :bottleneck_function => [fill(0.0, space...) for t in 0:generations]  # an array of arrays with probablity of external death at each site and generation.
 )
 

From 0745fda8657a221a23f989db97712c1580fb9424 Mon Sep 17 00:00:00 2001
From: "Ali R. Vahdati" <arv.kavir@gmail.com>
Date: Tue, 30 Aug 2022 18:30:18 +0200
Subject: [PATCH 6/9] fix tests

---
 test/api_tests.jl | 1 -
 test/params.jl    | 2 +-
 2 files changed, 1 insertion(+), 2 deletions(-)

diff --git a/test/api_tests.jl b/test/api_tests.jl
index cd0b917..89fb30e 100644
--- a/test/api_tests.jl
+++ b/test/api_tests.jl
@@ -1,4 +1,3 @@
-#TODO: more tests
 @testset "Model" begin
   adata, mdata, models = runmodel(param_file)
   @test size(mdata, 1) == 15
diff --git a/test/params.jl b/test/params.jl
index e7a726b..5dfde6e 100644
--- a/test/params.jl
+++ b/test/params.jl
@@ -1,4 +1,4 @@
-generations = 5
+generations = 14
 space = (2, 2)
 
 function env_resources(time::Int)

From 38011923cd75e5907412611c87a257ebfbef65a2 Mon Sep 17 00:00:00 2001
From: "Ali R. Vahdati" <arv.kavir@gmail.com>
Date: Wed, 31 Aug 2022 15:50:47 +0200
Subject: [PATCH 7/9] Fix leaking seed

---
 src/interactions.jl | 2 +-
 src/migration.jl    | 2 +-
 src/reproduction.jl | 8 ++++----
 3 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/interactions.jl b/src/interactions.jl
index 5b15dfe..9ecc00e 100644
--- a/src/interactions.jl
+++ b/src/interactions.jl
@@ -158,7 +158,7 @@ function target_species_ids(agent, model::ABM)
   elseif foundlen < nspecies
     return allids
   else
-    species_ids = sample(allids, nspecies, replace=false)
+    species_ids = sample(model.rng, allids, nspecies, replace=false)
     return species_ids  
   end
 end
diff --git a/src/migration.jl b/src/migration.jl
index eeac3f8..fb33a51 100644
--- a/src/migration.jl
+++ b/src/migration.jl
@@ -12,7 +12,7 @@ function migrate!(agent::Ind, model::ABM)
   else
     nsites = length(sites)
     nsites_selected = ceil(Int, model.check_fraction[agent.species] * nsites)
-    sites_checked = EvoDynamics.sample(sites, nsites_selected, replace=false)
+    sites_checked = sample(model.rng, sites, nsites_selected, replace=false)
     # check sites and move to the best one
     distance, site_index = pick_site(agent, sites_checked, nsites_selected, model)
     current_place = check_site(agent, agent.pos, model)
diff --git a/src/reproduction.jl b/src/reproduction.jl
index 241697d..ced914b 100644
--- a/src/reproduction.jl
+++ b/src/reproduction.jl
@@ -2,7 +2,7 @@
 """
 Returns a bitarray for sites to be selected from the first (false) and second (true) homologous chromosome.
 """
-function crossing_overs(nsites::Int, ncrossing_overs::Int)
+function crossing_overs(nsites::Int, ncrossing_overs::Int, RNG)
   output = falses(nsites)
   if ncrossing_overs == 0
     return output
@@ -10,7 +10,7 @@ function crossing_overs(nsites::Int, ncrossing_overs::Int)
     output[1:2:end] .= true
     return output
   end
-  breaking_points = sample(1:(nsites-1), ncrossing_overs, replace=false, ordered=true)
+  breaking_points = sample(RNG, 1:(nsites-1), ncrossing_overs, replace=false, ordered=true)
   last = true
   counter = 0
   for site in 1:nsites
@@ -62,8 +62,8 @@ function create_one_offspring(ag1::Ind, ag2::Ind, model::ABM)
     nco1, nco2 = rand(model.rng, model.recombination[species], 2)
   end
 
-  cross_overs1 = crossing_overs(nsites, nco1)
-  cross_overs2 = crossing_overs(nsites, nco2)
+  cross_overs1 = crossing_overs(nsites, nco1, model.rng)
+  cross_overs2 = crossing_overs(nsites, nco2, model.rng)
   gametes1 = create_gamete(ag1, cross_overs1, nsites, rand(model.rng, (true, false)))
   gametes2 = create_gamete(ag2, cross_overs2, nsites, rand(model.rng, (true, false)))
 

From 3dc1e0f8e6c42387cb7ea93a95e75a92fbc409ab Mon Sep 17 00:00:00 2001
From: "Ali R. Vahdati" <arv.kavir@gmail.com>
Date: Wed, 31 Aug 2022 17:50:31 +0200
Subject: [PATCH 8/9] Add more tests and fix bugs

---
 src/EvoDynamics.jl       |   5 ++
 src/api.jl               |   2 -
 src/interactions.jl      |  20 ++++++-
 src/reproduction.jl      |   6 +-
 src/simulation.jl        |   4 +-
 test/params.jl           |   4 +-
 test/simulation_tests.jl | 116 +++++++++++++++++++++++++++++++++++++--
 7 files changed, 140 insertions(+), 17 deletions(-)

diff --git a/src/EvoDynamics.jl b/src/EvoDynamics.jl
index c657289..4bc22a6 100644
--- a/src/EvoDynamics.jl
+++ b/src/EvoDynamics.jl
@@ -9,6 +9,11 @@ using StaticArrays
 using YAML
 using RandomNumbers
 
+export runmodel, model_initiation, load_parameters
+# export abiotic_survive!, burn_energy!, consume_food!, interact!, remove_agent!, reproduce!, migrate!, target_species_ids, eat!, abiotic_distance, abiotic_fitness, interaction_power, phenotypic_distance, get_abiotic_phenotype, get_biotic_phenotype, update_fitness!, vertex2coord, coord2vertex, mutate!, adjust_fitness!, create_gamete, crossing_overs, return_opt_phenotype
+# export check_site, pick_site, get_migration_trait
+# export step!
+
 include("simulation.jl")
 include("interactions.jl")
 include("migration.jl")
diff --git a/src/api.jl b/src/api.jl
index c3a5cb3..cfbce7f 100644
--- a/src/api.jl
+++ b/src/api.jl
@@ -1,5 +1,3 @@
-export runmodel
-
 "Returns a tuple whose entries are the mean fitness of each species."
 function mean_fitness_per_species(model::ABM)
   mean_fitness = Array{Float64}(undef, model.nspecies)
diff --git a/src/interactions.jl b/src/interactions.jl
index 9ecc00e..cd5f728 100644
--- a/src/interactions.jl
+++ b/src/interactions.jl
@@ -32,7 +32,7 @@ Distance is a probability and is calculated by the number of standard deviations
 
 Variance can stand for selection coefficient, with a larger variance, less sever the difference.
 """
-function phenotypic_distance(phenotype1, phenotype2, interaction_coeff, variance::Real)
+function phenotypic_distance_different_species(phenotype1, phenotype2, interaction_coeff, variance::Real)
   distance = 0.0
   counter = 0.0
   for ph1 in phenotype1
@@ -50,7 +50,23 @@ function phenotypic_distance(phenotype1, phenotype2, interaction_coeff, variance
   return distance
 end
 
-phenotypic_distance(ag1::Ind, ag2::Ind, model) = phenotypic_distance(ag1.biotic_phenotype, ag2.biotic_phenotype, model.interactions[ag1.species, ag2.species], model.biotic_variances[ag1.species])
+function phenotypic_distance_same_species(phenotype1, phenotype2, interaction_coeff, variance::Real)
+  distance = 0.0
+  counter = 0.0
+  for (ph1, ph2) in zip(phenotype1, phenotype2)
+      d = abs(0.5 - cdf(Normal(ph1, variance), ph2)) # 0.5 is the value when ph1 and ph2 are identical
+      distance += d
+      counter += 1.0
+  end
+  distance += distance  # Double the distance so that it is in range 0 to 1.
+  distance /= counter  # average distance
+  if interaction_coeff < 0
+    distance = 1.0 - distance
+  end
+  return distance
+end
+
+phenotypic_distance(ag1::Ind, ag2::Ind, model) = ag1.species == ag2.species ? phenotypic_distance_same_species(ag1.biotic_phenotype, ag2.biotic_phenotype, model.interactions[ag1.species, ag2.species], model.biotic_variances[ag1.species]) : phenotypic_distance_different_species(ag1.biotic_phenotype, ag2.biotic_phenotype, model.interactions[ag1.species, ag2.species], model.biotic_variances[ag1.species])
 
 """
 Abiotic distance to the optimal values
diff --git a/src/reproduction.jl b/src/reproduction.jl
index ced914b..578e1d1 100644
--- a/src/reproduction.jl
+++ b/src/reproduction.jl
@@ -81,7 +81,7 @@ function create_one_offspring(ag1::Ind, ag2::Ind, model::ABM)
   W = adjust_fitness(W, ag1.species, model)
   initial_energy = model.initial_energy[species]
 
-  offspring = add_agent!(ag1.pos, model, ag1.species, bph, abph, episMat, pleioMat, q, 1, sex, interaction_history, initial_energy, W, true, 0, -1)
+  offspring = add_agent!(ag1.pos, model, ag1.species, deepcopy(bph), deepcopy(abph), deepcopy(episMat), deepcopy(pleioMat), deepcopy(q), 1, sex, deepcopy(interaction_history), initial_energy, W, true, 0, -1)
   return offspring
 end
 
@@ -93,7 +93,7 @@ function reproduce!(ag1::Ind, ag2::Ind, model::ABM)
   if model.mating_schemes[ag1.species] == 0.0
     reproduction_success = 1.0
   else
-    d = phenotypic_distance(ag1.biotic_phenotype, ag2.biotic_phenotype, model.mating_schemes[ag1.species], model.biotic_variances[ag1.species])
+    d = phenotypic_distance(ag1, ag2, model)
     reproduction_success = interaction_power(ag1, ag2, d, model)
   end
   growth_rate = model.growthrates[ag1.species]
@@ -117,7 +117,7 @@ function reproduce!(agent::Ind, model::ABM)
     interaction_history = deepcopy(agent.interaction_history)
     interaction_history[1:end] .= 0
     for c in 1:nchildren
-      offspring = add_agent!(agent.pos, model, agent.species, agent.biotic_phenotype, agent.abiotic_phenotype, agent.epistasisMat, agent.pleiotropyMat, agent.q, 1, agent.sex, interaction_history, model.initial_energy[agent.species], agent.W, true, 0, -1)
+      offspring = add_agent!(agent.pos, model, agent.species, deepcopy(agent.biotic_phenotype), deepcopy(agent.abiotic_phenotype), deepcopy(agent.epistasisMat), deepcopy(agent.pleiotropyMat), deepcopy(agent.q), 1, agent.sex, interaction_history, model.initial_energy[agent.species], agent.W, true, 0, -1)
       mutate!(offspring, model)
     end
   elseif model.ploidy[agent.species] == 2 && agent.interaction_history[agent.species] == model.step[1] && agent.time_met_other_sex == model.step[1]
diff --git a/src/simulation.jl b/src/simulation.jl
index ca7f75a..d1ec617 100644
--- a/src/simulation.jl
+++ b/src/simulation.jl
@@ -105,7 +105,7 @@ function model_initiation(dd)
         end
         interaction_history = MVector{model.nspecies,Int}(fill(-1, model.nspecies))
         initial_energy = model.initial_energy[sp]
-        add_agent!(model.nodes[pos], model, sp, biotic_ph, abiotic_ph, epistasisMat[sp], pleiotropyMat[sp], expressionArrays[sp], 1, sex, interaction_history, initial_energy, W, true, 0, -1)
+        add_agent!(model.nodes[pos], model, sp, deepcopy(biotic_ph), deepcopy(abiotic_ph), deepcopy(epistasisMat[sp]), deepcopy(pleiotropyMat[sp]), deepcopy(expressionArrays[sp]), 1, sex, deepcopy(interaction_history), initial_energy, W, true, 0, -1)
       end
     end
   end
@@ -302,7 +302,7 @@ function mutate!(agent::Ind, model::ABM)
   end
   # mutate pleiotropy matrix
   if rand(model.rng, model.mutProbs[agent.species][2])
-    randnumbers = rand(model.rng, model.rng, model.mutMagnitudes[agent.species][2], size(agent.pleiotropyMat))
+    randnumbers = rand(model.rng, model.mutMagnitudes[agent.species][2], size(agent.pleiotropyMat))
     agent.pleiotropyMat[randnumbers] .= .!agent.pleiotropyMat[randnumbers]
     mutated = true
   end
diff --git a/test/params.jl b/test/params.jl
index 5dfde6e..9d5f3e6 100644
--- a/test/params.jl
+++ b/test/params.jl
@@ -32,8 +32,8 @@ species1 = Dict(
   :selection_coefficient => 0.02,
   :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.9, 0.0, 0.0],  # for gene expression array, pleiotropy matrix and epistasis matrix, respectively
-  :mutation_magnitudes => [0.05, 0.0, 0.01], # same as above
+  :mutation_probabilities => [0.99, 0.99, 0.99],  # for gene expression array, pleiotropy matrix and epistasis matrix, respectively
+  :mutation_magnitudes => [0.01, 0.1, 0.1], # same as above
   :N => [1000, 0, 0, 0],  # number of individuals per site at time 0
   :environmental_noise => 0.0,  # variance of a normal distribution with mean 0
   # each row is the optimal phenotypes for each site for all abiotic traits. There are as many element as number of sites times number of abiotic traits. The first N elements are for the first abiotic trait, where N is the number of sites, and so on.
diff --git a/test/simulation_tests.jl b/test/simulation_tests.jl
index 1df1654..55cde5d 100644
--- a/test/simulation_tests.jl
+++ b/test/simulation_tests.jl
@@ -1,6 +1,33 @@
+@testset "Basic functions" begin
+  dd = EvoDynamics.load_parameters(param_file)
+  model = EvoDynamics.model_initiation(dd)
+  EvoDynamics.remove_agent!(model[1], model)
+  @test !haskey(model.agents, 1) == true
+
+  EvoDynamics.consume_food!(model[2], model)
+  @test model.resources == [1999 1980; 1830 1900]
+
+  @test EvoDynamics.coord2vertex((2, 1), model) == 2
+  @test EvoDynamics.coord2vertex((1, 2), model) == 3
+  @test EvoDynamics.coord2vertex((2, 2), model) == 4
+  @test EvoDynamics.vertex2coord(2, model) == (2, 1)
+  @test EvoDynamics.vertex2coord(3, model) == (1, 2)
+  @test EvoDynamics.vertex2coord(4, model) == (2, 2)
+
+  @test model[2].W == model[3].W
+  model[2].q[1] -= 0.05
+  EvoDynamics.update_fitness!(model[2], model)
+  @test model[2].W != model[3].W
+end
+
 @testset "Mutation" begin
   dd = EvoDynamics.load_parameters(param_file)
   model = EvoDynamics.model_initiation(dd)
+
+  @test EvoDynamics.nagents(model) == 1100
+  @test size(model.space) == (2, 2)
+  @test length(EvoDynamics.ids_in_position((1, 1), model)) == EvoDynamics.nagents(model)
+
   ag1 = model.agents[1]
   ag2 = model.agents[2]
 
@@ -18,14 +45,91 @@
   EvoDynamics.update_fitness!(ag1, model)
   @test ag1.W != Worg
 
+  # no leakage. using deepcopy for creating new agents
+  @test model[3].q == model[4].q
+  @test model[3].pleiotropyMat == model[4].pleiotropyMat
+  @test model[3].epistasisMat == model[4].epistasisMat
+  EvoDynamics.mutate!(model[3], model)
+  @test model[3].q != model[4].q
+  @test model[3].pleiotropyMat != model[4].pleiotropyMat
+  @test model[3].epistasisMat != model[4].epistasisMat
+  @test model[5].q == model[4].q
+  @test model[5].pleiotropyMat == model[4].pleiotropyMat
+  @test model[5].epistasisMat == model[4].epistasisMat
 end
 
-# #TODO
-# @testset "Migration" begin
+@testset "Seeding" begin
+  dd = EvoDynamics.load_parameters(param_file)
+  dd[:seed] = 1
+  model = model_initiation(dd)
+  EvoDynamics.step!(model, EvoDynamics.agent_step!, 1)
+
+  dd2 = EvoDynamics.load_parameters(param_file)
+  dd2[:seed] = 2
+  model2 = EvoDynamics.model_initiation(dd2)
+  EvoDynamics.step!(model2, EvoDynamics.agent_step!, 1)
 
-# end
+  dd3 = EvoDynamics.load_parameters(param_file)
+  dd3[:seed] = 1
+  model3 = EvoDynamics.model_initiation(dd3)
+  EvoDynamics.step!(model3, EvoDynamics.agent_step!, 1)
 
-# #TODO
-# @testset "Reproduction" begin
+  @test collect(keys(model.agents)) != collect(keys(model2.agents))
+  @test collect(keys(model.agents)) == collect(keys(model3.agents))
+end
 
-# end
+@testset "Interactions" begin
+  dd = EvoDynamics.load_parameters(param_file)
+  model = EvoDynamics.model_initiation(dd)
+
+  agent = model[1]
+  bph = EvoDynamics.get_biotic_phenotype(agent, model)
+  aph = EvoDynamics.get_abiotic_phenotype(agent, model)
+
+  @test length(bph) == length(dd[:species][1][:biotic_phenotypes])
+  @test length(aph) == length(dd[:species][1][:abiotic_phenotypes])
+
+  # identical phenotypes in the same species results in the same distance:
+  agent2 = model[1001]
+  d = EvoDynamics.phenotypic_distance(agent, agent2, model)
+  dsame = EvoDynamics.phenotypic_distance(agent, agent, model)
+  dsame2 = EvoDynamics.phenotypic_distance(agent2, agent2, model)
+  @test dsame == dsame2 == 0
+
+  # identical phenotypes in two different species results in non-zero distance:
+  agent.biotic_phenotype .= agent2.biotic_phenotype
+  d2 = EvoDynamics.phenotypic_distance(agent, agent2, model)
+  @test d2 > 0.0
+  @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
+  agent2.abiotic_phenotype .= optphen
+  @test EvoDynamics.abiotic_distance(agent2.abiotic_phenotype, optphen, model.abiotic_variances[agent2.species]) == 0
+  @test EvoDynamics.abiotic_fitness(agent2, model) == 1.0
+
+  EvoDynamics.interact!(agent, agent2, model)
+  @test agent.isalive == false
+
+  @test length(EvoDynamics.target_species_ids(agent, model)) == 2
+  @test length(EvoDynamics.target_species_ids(agent2, model)) == 2
+
+  # reproduction
+  nagentsbefore = EvoDynamics.nagents(model)
+  EvoDynamics.reproduce!(agent2, model)
+  EvoDynamics.reproduce!(agent2, model)
+  EvoDynamics.reproduce!(agent2, model)
+  nagentsafter = EvoDynamics.nagents(model)
+  @test nagentsbefore < nagentsafter
+
+  agent3 = model[10]
+  agent4 = model[11]
+  agent3.sex = true
+  agent4.sex = false
+
+  nagentsbefore = EvoDynamics.nagents(model)
+  EvoDynamics.reproduce!(agent3, agent4, model)
+  EvoDynamics.reproduce!(agent3, agent4, model)
+  nagentsafter = EvoDynamics.nagents(model)
+  @test nagentsbefore < nagentsafter
+end

From 30a0515704e575542dbc01e35f62c90aebf1dbb1 Mon Sep 17 00:00:00 2001
From: "Ali R. Vahdati" <arv.kavir@gmail.com>
Date: Thu, 1 Sep 2022 16:15:57 +0200
Subject: [PATCH 9/9] Update the docs

---
 docs/src/index.md             |  14 ++--
 docs/src/model_description.md |  57 ++++++++------
 examples/example1.jl          |   2 +-
 examples/example2.jl          | 138 ++++++++++++++++++++--------------
 examples/paramfile1.jl        |   2 +-
 examples/paramfile2.jl        |  12 +--
 src/check_params.jl           |   4 +-
 src/reproduction.jl           |   2 +-
 src/simulation.jl             |   2 +-
 test/params.jl                |  10 +--
 10 files changed, 138 insertions(+), 105 deletions(-)

diff --git a/docs/src/index.md b/docs/src/index.md
index 9b1a84d..da1cffa 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -1,6 +1,8 @@
 # EvoDynamics.jl Documentation
 
-EvoDynamics.jl tries to bridge the gap in studying biological systems at small and large scales. Some studies only focus on single genes affecting single phenotypes, some studies only analyze gene interactions, some focus on populations, and some on species interactions. EvoDynamics.jl is a framework to study the effect of interactions between all these levels. It includes explicit pleiotropy, epistasis, selection acting on multiple phenotypes, different phenotypes affecting fitness at different amounts, arbitrary spatial structure, migration, and interacting species.
+EvoDynamics.jl is a framework to bridge the gap in studying biological systems at small and large scales. Specifically, it allows modeling biological dynamics at evolutionary and ecological scale. Genotypes affect phenotypes, and they interact with ecology (i.e. environment and other species). The flow of causal effect runs in both directions. Such changes affect adaptation, population demography, and species coexistence. Such eco-evolutionary feedbacks occur when the speed of evolution is fast; populations need to constantly adapt to changing conditions. This happens in antagonistic interactions between species (e.g. predator-prey or parasitism), competitive interactions, and mutualistic interactions.  
+
+Some studies only focus on single genes affecting single phenotypes, some studies only analyze gene interactions, some focus on populations, and some on species interactions. EvoDynamics.jl is a framework to study the effect of interactions between all these levels. It includes an explicit genotype-phenotype architecture (pleiotropy and epistasis), selection acting on multiple phenotypes, different phenotypes affecting fitness at different amounts, arbitrary spatial structure, migration, and interacting species.
 
 Figure below shows different biological levels controlled by the model.
 
@@ -11,13 +13,13 @@ See [Tutorial](@ref) for running the model, and [Model description](@ref) for a
 ## Features
 
 * Possibility to model complex food webs with various asymmetrical interactions.
-* Individuals interact given their phenotypes.
+* Individuals interact with one another and with the environment given their phenotypes.
 * Connecting genome structure to phenotypes to populations.
 * Space is a grid on which individuals can migrate.
-* Possibility to model a complex environment with various amounts of resources at each site.
-* Allows time varying optimal phenotypic value per site.
-* Possibility of killing specific individuals at certain times and sites.
-* Can model both haploid and diploid species.
+* Possibility to model a complex environment with spatio-temporally varying resources.
+* 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.
 * Includes mutation and recombination.
 * Easy data collection.
 * Runs replicates and collects data in a table automatically.
diff --git a/docs/src/model_description.md b/docs/src/model_description.md
index a18ddbe..894d664 100644
--- a/docs/src/model_description.md
+++ b/docs/src/model_description.md
@@ -2,47 +2,57 @@
 
 ## Parameters
 
-### Species specific parameters
+All parameters of a model are written in a julia (`.jl`) file. Parameters are divided into two sections: species specific and model specific. Parameters for each species and for the model are written in separate dictionary (`Dict`) objects. The keys of the dictionary should be Symbols (use colon : before the names).
 
-Each species should have the following parameters. The order that you write these parameters does not matter.
+An easy way to parameterize a model is to copy an existing set of parameters from the [Examples](@ref) and modify them.
 
-* __name__: a name for the species.
+### 1. Species specific parameters
+
+Each species should have the following parameters in a dictionary objection. The order that you write these parameters does not matter.
+
+* __name__: A `String` as the name of the species.
 * __number\_of\_genes__: An _integer_ for number of genes that the species has.
-* __ploidy__: Either 1 for haploid or 2 for diploid genomes. Diploids can have recombination.
+* __ploidy__: Either 1 for haploid or 2 for diploid genomes. Diploids may undergo recombination.
 * __number\_of\_phenotypes__: An _integer_ for the number of phenotypes that the species has.
 * __abiotic\_phenotypes__: An _array of integers_ (e.g. "[1,2]") specifying abiotic phenotypes among all phenotypes. Abiotic phenotypes determine how the species interacts with the environment.
 * __biotic\_phenotypes__: An _array of integers_ (e.g. "[3]") specifying biotic phenotypes among all phenotypes. Biotic phenotypes determine how the species interacts with other individuals from the same or different species.
-* __migration\_phenotype__: An _integer_ specifying the phenotype that determines migration trait. If the species does not migrate, put 0.
+* __migration\_phenotype__: An _integer_ specifying the phenotype that corresponds to the migration trait. If the species does not migrate, put 0.
 * __migration\_threshold__: The phenotypic value of the migration phenotype after which an individual can migrate.
 * __vision\_radius__: A _number_ determining the radius of neighboring sites that the agent can see before migration.
-* __check\_fraction__: A _number_ between 0 and 1 showing the fraction of the visible sites to the agent that it can check and decide whether to migrate to.
+* __check\_fraction__: A _floating number_ between 0 and 1 showing the fraction of the visible sites to the agent that it can check and decide whether to migrate to.
 * __epistasis\_matrix__: An epistasis matrix is of size $l \times l$, where _l_ is the product of _number of genes_ and _ploidy_. Epistasis matrix specifies the direction (positive or negative) and size of effect of one locus on other loci. For example, if at row 1 and column 2 is a value 0.2, it means that locus 1 affects locus 2 by increasing the effect of locus 2 (because its positive) with 20% of the effect of locus 1.
 * __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 the maximum mean when fitness of a haploid individual is 1, or the distance between the biotic phenotypes of two diploid individuals is 0.
-* __selection\_coefficient__: A number between 0 and 1 that determines the importance of fitness. 0 would be a model without selection.
-* __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.
+* __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.
+* __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.
 * __environmental\_noise__: A number for the variance of a normal distribution with mean 0 that will be added to the phenotypes.
-* __optimal\_phenotypes__: Name of a function that returns a _vector of integers_ for optimal phenotypes of the species at a specific site.  The vector should be as long as there are `abiotic phenotypes`. The function should have two arguments: `site `of type `Tuple{Int, Int}`and `model `of type `ABM `from the `Agents.jl `package. The definition of the function should be in a file referred to in the `functions file` parameter (see blow). This function can be time dependent.
+* __optimal\_phenotypes__: A _vector of matrices_. The vector holds optimal phenotypes at each time step. Each matrix, is the optimal phenotypes for each site in space. The optimal phenotype at each site is a vector as long as `abiotic_phenotypes`. The vector is as long as model steps (generations) plus 1 (for time zero).
 * __age__: An _integer_ for maximum age of individuals of this species.
 * __reproduction\_start\_age__: The age at which individuals can reproduce.
 * __reproduction\_end\_age__: The age after which individuals cannot reproduce.
 * __recombination__: Mean of a Poisson distributions for number of crossing overs per sexual reproduction.
-* __initial\_energy__: A parameter for parental care of infants. Values more than 0 indicate that newly born individuals can survive for a number of times without requiring food from the environment/other species. The consumption rate (i.e. how many generations this initial energy suffices) is determined by the sum of the corresponding rows in "food sources" model parameter.
-* __bottleneck\_function__: Name of a function in `functions file` that accepts two arguments: `agent` and `model`. It returns true or false, where true means the agent is killed. This function models killing agents due to a process above the dynamics of the model, for example, hunting or environmental disaster. Since it accepts both the agent and the model, there are information about the position of the agent, its properties, and the time step of the simulation to take into account.
+* __initial\_energy__: A parameter for parental care of infants. Values more than 0 indicate that newly born individuals can survive for a number of times without requiring food from the environment/other species. The consumption rate is always constant for all species at one unit per time step. Use with care. Having initial energy larger than zero can lead to infinite population growth because agents without food can reproduce and their offsprings also reproduce without food. For example, this happens when start age of reproduction is 1 and initial energy is larger than 0.
+* __bottlenecks__: An vector of matrices each of which has the same size as space and contains the probability of death at each site. The vector is as long as model steps (generations) plus 1 (for time zero). Use this if you want to impose a killing event in the model at specific times and places. Otherwise, just use 0.0 for all the values in matrices.
+* __mating\_scheme__: One of the following: -1, 0, 1. Only used in sexual reproduction. Zero means randomly paired couple are likely to have as many offsprings as any other pair. Their phenotypes does not affect their reproduction success. -1 is disassortative mating, where the more different the phenotypes of the pair, the more children they have. 1 is the opposite, assortative mating.
+* __abiotic\_variance__: 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__: Same as _abiotic\_variance_ 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.
+
+### 2. Model parameters
 
-### Model parameters
+This dictionary should be named "model_parameters".
 
+* __species__: A vector containing the name of the species _dictionaries_.
 * __generations__: An _integer_ for the number of steps the model will run.
-* __space__: A _vector of two integers_ that determine the size of a grid for space.
+* __space__: A _tuple of two integers_ (e.g. `(3,2)`) that determines the size space. If you do not want any spatial structure, use `(1,1)`.
 * __metric__: Either "chebyshev" or "euclidian". Determines how many neighbors a space site has. "chebyshev" metric means that the r-neighborhood of a position are all positions within the hypercube having side length of 2*floor(r) and being centered in the origin position. "euclidean" metric means that the r-neighborhood of a position are all positions whose cartesian indices have Euclidean distance ≤ r from the cartesian index of the given position.
 * __periodic__: _Boolean__ (true or false) to determine whether boundaries of the space are connected or not.
-* __resources__: A function that returns a _vector of integers_ determining available resources (e.g. vegetation) per site. The function takes the `time_step` as its only argument. The function is defined in the functions file (parameter blow).
-* __interactions__: A species-species interaction _matrix of numbers_ determining how individuals from different species interact. Each value  is strength of interaction (between 0 and 1). Sign (+/-) is the direction of interaction where positive means similar individuals interact more strongly and negative is dissimilar ones tend to interact more.
-* __food\_sources__: A species-species food _matrix of numbers_ determining what each species feeds on (consumption rate). Non-zero diagonal means the food resource is from the environment. Off-diagonals mean an species (in the rows) feeds on another species (in the columns). Numbers can be zero or any positive number. The magnitude of the number determines how many generations can an individual live off of given one unit of the food source. For example, if a diagonal is 2, it means that the species will eat one unit of the environmental resources and that is enough for it to live two steps.
-* __seed__: Either an _integer_ or _Null_ for random number generator seed.
+* __resources__: A _vector of matrices_ where each matrix is the size of space and determines the amount of resources per site. The vector is as long as model steps (generations) plus 1 (for time zero).
+* __interactions__: A species-species interaction _matrix of floating numbers_ determining how individuals from different species interact. Each value  is strength of interaction (between 0 and 1). Sign (+/-) is the direction of interaction where positive means similar individuals interact more strongly and negative is dissimilar ones tend to interact more.
+* __food\_sources__: A species-species food _matrix of floating numbers_ determining what each species feeds on (consumption rate). Non-zero diagonal means the food resource is from the environment. Off-diagonals mean an species (in the rows) feeds on another species (in the columns). Numbers can be zero or any positive number. The magnitude of the number determines how much energy the predator gains by eating the prey. All species use 1 unit of energy per time step.
+* __seed__: Either an _integer_ or `nothing` for a random number.
 
 ## Simulation outline
 
@@ -53,10 +63,9 @@ The simulations are fully agent-based, meaning that agents do not receive any mo
 3. Burn energy.
 4. Eat from the environment if the species can.
 5. Interact with other individuals.
-   1. If meeting another individual of the same species but with different sex, try reproduction.
-6. If haploid, reproduce.
+6. Reproduce.
 7. Survive. Agents die if they are too old, or do not have enough energy, or by chance given its fitness.
-8. Go through the bottleneck function if it should be activated at the current time.
+8. Kill agents will a probability given in `bottlenecks`.
 
 ### Migration
 
@@ -84,13 +93,13 @@ $ | 0.5 - \text{cdf}(\text{Normal}(\text{ph1},  1), \text{ph2})| $
 
 , where cdf is cumulative density function, Normal is a normal distribution with mean ph1 (phenotypic value 1) and variance 1, and ph2 is phenotypic value 2.
 
-If the two individuals are not predator-prey, and they are from the same species but different sexes and they are both in reproduction age, then they try to reproduce. Otherwise, if they are from different species and the two species interact with each other (the corresponding values in the _interactions_ matrix is non-zero) they interact with each other. If the individual 1 to 2 has a positive value in the _interactions_ matrix, then it increases the fitness of individual 2.  If the element is zero, it does not increase or decrease the fitness of individual 2. And it the element is negative, it decreases the fitness of individual 2. Similarly, individual 2 can affect the fitness of individual 1. The two interactions do not need to be symmetric. The magnitude of change in the fitness due to interaction is determined by the phenotypic distance between the two individuals and interaction value in the _interactions_ matrix.
+If the two individuals are not predator-prey, and they are from the same species but different sexes and they are both in reproduction age, then they are marked to reproduce in the next stage. Otherwise, if they are from different species and the two species interact with each other (the corresponding values in the _interactions_ matrix is non-zero) they interact with each other. If the individual 1 to 2 has a positive value in the _interactions_ matrix, then it increases the fitness of individual 2. If the element is zero, it does not increase or decrease the fitness of individual 2. And it the element is negative, it decreases the fitness of individual 2. Similarly, individual 2 can affect the fitness of individual 1. The two interactions do not need to be symmetric. The magnitude of change in the fitness due to interaction is determined by the phenotypic distance between the two individuals and interaction value in the _interactions_ matrix.
 
 ### Reproduction
 
 If an individual is haploid and is in reproduction age, it reproduces an identical offspring to itself except that the expression array, pleiotropy matrix, and the epistasis matrix of the offspring will mutate given the probabilities and magnitudes in the _mutation probabilities_ and _mutation magnitudes_ matrices. The number of offsprings is a random number from a Poisson distribution with a mean equal to the species' growth rate times the fitness of the individual.
 
-If two diploid individuals are mated to reproduce, their reproduction success is proportional to their phenotypic similarity. Number of offsprings is a random number from a Poisson distribution with mean equal to the reproduction success of the two individuals times the growth rate of the species.
+If two diploid individuals are mated to reproduce, their reproduction success is proportional to their phenotypic distance (depending on `mating\_scheme`). Number of offsprings is a random number from a Poisson distribution with mean equal to the reproduction success of the two individuals times the growth rate of the species.
 
 Each offspring of diploid individuals inherits a gamete from each parent. A gamete is a half of the expression array, pleiotropy matrix, and epistasis matrix. If recombination is allowed (> 0), then the gametes undergo crossing over. The number of crossing overs is a random number from a Poisson distribution with mean equal to the _recombination_ parameter.
 
diff --git a/examples/example1.jl b/examples/example1.jl
index bdc2635..74ad741 100644
--- a/examples/example1.jl
+++ b/examples/example1.jl
@@ -34,7 +34,7 @@ using EvoDynamics
 #   :N => [100],
 #   :environmental_noise => 0.01,
 #   :optimal_phenotypes => [fill([1.5 for p in 1:1], space...) for t in 0:generations],
-#   :bottleneck_function => [fill(0.0, space...) for t in 0:generations],
+#   :bottlenecks => [fill(0.0, space...) for t in 0:generations],
 #   :age => 2,
 #   :recombination => 0,
 #   :initial_energy => 0,
diff --git a/examples/example2.jl b/examples/example2.jl
index 2fa4622..cfde78a 100644
--- a/examples/example2.jl
+++ b/examples/example2.jl
@@ -4,41 +4,15 @@
 # Here we create a predator prey model of two species, one haploid and one diploid.
 
 using EvoDynamics
+using Plots
 
 # Model parameters are in a .jl file as follows:
 
 # ```julia
-# ## 1. Functions
-
-# function bn(agent::AbstractAgent, model::ABM)
-#   return false
-# end
-
-# function optphens2(site::Tuple{Int,Int}, model::ABM)
-#   opt1 = [0.8214794136627462, 0.49335401288317304, 0.3435556249656141, 0.25050033180075226]
-#   opt2 = opt1 .+ 0.1
-#   ss = EvoDynamics.coord2vertex(site, model)
-#   if model.step[1] > 6
-#     return [opt1[ss]]
-#   else
-#     return [opt2[ss]]
-#   end
-# end
-
-# function optphens3(site::Tuple{Int,Int}, model::ABM)
-#   phenotypic_values = [[0.7614758101208934 2.3911353663016586; 2.2091361414343313 1.7858661540288683; 0.7920974352892358 0.7630236717263839; 2.587205421882114 2.311211631439866],
-#     [0.6437641305315445 2.097629262577973; 0.9954545836033715 1.1314248848669466; 2.469792530385348 1.490299526620522; 1.6158867433451882 0.2566056477862022]]
-#   agent_site = (site[2] - 1) * size(model.space)[2] + (site[1])
-#   if model.step[1] > 7
-#     return phenotypic_values[1][agent_site, :]
-#   else
-#     return phenotypic_values[2][agent_site, :]
-#   end
-# end
-
-# env_resources2(time::Int) = [200 220; 183 190]
-
-# ## 2. Species parameters
+# generations = 14
+# space = (2, 2)
+
+# ## 1. Species parameters
 
 # species1 = Dict(
 #   :name => "a",
@@ -47,26 +21,47 @@ using EvoDynamics
 #   :abiotic_phenotypes => [1],
 #   :biotic_phenotypes => [2, 3],
 #   :migration_phenotype => 4,
-#   :migration_threshold => 3.5,
+#   :migration_threshold => 1.5,
 #   :vision_radius => 1,
 #   :check_fraction => 0.5,
 #   :ploidy => 2,
-#   :epistasis_matrix => [1.0 0.43 -0.41 0.38 0.48 -0.43 -0.1 -0.08 -0.09 -0.5 0.41 0.44 -0.21 -0.12; 0.34 1.0 -0.19 -0.36 0.38 -0.28 0.24 -0.22 0.12 0.12 -0.12 -0.39 0.21 0.26; 0.05 0.27 1.0 0.04 0.01 -0.14 0.3 -0.28 0.43 -0.13 0.2 -0.02 0.25 -0.39; -0.12 0.33 -0.48 1.0 -0.4 -0.48 -0.22 -0.36 -0.24 -0.07 -0.12 -0.49 -0.37 0.27; 0.25 0.25 -0.14 0.49 1.0 0.28 -0.34 -0.49 0.45 -0.14 0.26 -0.13 -0.44 -0.17; -0.47 0.19 -0.24 0.41 0.08 1.0 0.11 0.03 0.15 0.49 0.04 0.41 -0.19 0.13; 0.37 0.09 -0.11 0.4 0.42 0.45 1.0 -0.01 -0.47 0.07 0.5 0.44 -0.18 -0.2; -0.32 0.15 0.4 -0.24 -0.21 0.5 0.22 1.0 -0.33 0.48 -0.49 0.07 0.5 -0.07; 0.02 -0.16 0.33 0.48 -0.42 0.39 0.2 -0.11 1.0 0.46 -0.06 0.22 -0.3 0.31; 0.41 -0.18 -0.16 -0.4 0.01 0.04 0.07 0.2 -0.37 1.0 -0.33 0.49 0.05 -0.42; 0.03 0.25 0.14 -0.36 0.28 -0.18 0.09 -0.2 0.46 -0.48 1.0 -0.21 0.41 -0.46; 0.0 0.44 -0.34 -0.42 0.37 -0.04 0.43 -0.25 0.21 0.19 0.29 1.0 -0.02 0.06; 0.46 -0.1 0.14 -0.22 -0.26 0.13 -0.5 -0.41 -0.31 0.0 -0.15 0.29 1.0 0.17; -0.22 0.21 0.46 -0.01 -0.35 -0.11 0.25 -0.03 0.18 -0.38 -0.4 -0.28 0.05 1.0],
-#   :pleiotropy_matrix => Bool[1 1 0 0 0 0 0 1 0 1 0 1 0 1; 1 0 1 1 1 0 1 1 1 0 0 1 1 0; 1 0 1 0 1 0 0 0 0 0 1 0 1 0; 0 0 0 0 1 0 0 1 1 1 1 0 1 0],
-#   :growth_rate => 0.8,
+#   :epistasis_matrix => [
+#     1.0 0.43 -0.41 0.38 0.48 -0.43 -0.1 -0.08 -0.09 -0.5 0.41 0.44 -0.21 -0.12
+#     0.34 1.0 -0.19 -0.36 0.38 -0.28 0.24 -0.22 0.12 0.12 -0.12 -0.39 0.21 0.26
+#     0.05 0.27 1.0 0.04 0.01 -0.14 0.3 -0.28 0.43 -0.13 0.2 -0.02 0.25 -0.39
+#     -0.12 0.33 -0.48 1.0 -0.4 -0.48 -0.22 -0.36 -0.24 -0.07 -0.12 -0.49 -0.37 0.27
+#     0.25 0.25 -0.14 0.49 1.0 0.28 -0.34 -0.49 0.45 -0.14 0.26 -0.13 -0.44 -0.17
+#     -0.47 0.19 -0.24 0.41 0.08 1.0 0.11 0.03 0.15 0.49 0.04 0.41 -0.19 0.13
+#     0.37 0.09 -0.11 0.4 0.42 0.45 1.0 -0.01 -0.47 0.07 0.5 0.44 -0.18 -0.2
+#     -0.32 0.15 0.4 -0.24 -0.21 0.5 0.22 1.0 -0.33 0.48 -0.49 0.07 0.5 -0.07
+#     0.02 -0.16 0.33 0.48 -0.42 0.39 0.2 -0.11 1.0 0.46 -0.06 0.22 -0.3 0.31
+#     0.41 -0.18 -0.16 -0.4 0.01 0.04 0.07 0.2 -0.37 1.0 -0.33 0.49 0.05 -0.42
+#     0.03 0.25 0.14 -0.36 0.28 -0.18 0.09 -0.2 0.46 -0.48 1.0 -0.21 0.41 -0.46
+#     0.0 0.44 -0.34 -0.42 0.37 -0.04 0.43 -0.25 0.21 0.19 0.29 1.0 -0.02 0.06
+#     0.46 -0.1 0.14 -0.22 -0.26 0.13 -0.5 -0.41 -0.31 0.0 -0.15 0.29 1.0 0.17
+#     -0.22 0.21 0.46 -0.01 -0.35 -0.11 0.25 -0.03 0.18 -0.38 -0.4 -0.28 0.05 1.0], 
+#   :pleiotropy_matrix => Bool[
+#     1 1 0 0 0 0 0 1 0 1 0 1 0 1
+#     1 0 1 1 1 0 1 1 1 0 0 1 1 0
+#     1 0 1 0 1 0 0 0 0 0 1 0 1 0
+#     0 0 0 0 1 0 0 1 1 1 1 0 1 0],
+#   :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.5,
+#   :selection_coefficient => 0.1,
 #   :mutation_probabilities => [0.9, 0.0, 0.0],
 #   :mutation_magnitudes => [0.05, 0.0, 0.01],
 #   :N => [1000, 0, 0, 0],
 #   :environmental_noise => 0.01,
-#   :optimal_phenotypes => optphens2,
-#   :bottleneck_function => bn,
+#   :optimal_phenotypes => [fill([1.4 for p in 1:1], space...) for t in 0:generations],
+#   :bottlenecks => [fill(0.0, space...) for t in 0:generations],
 #   :age => 5,
 #   :recombination => 1,
-#   :initial_energy => 0,
-#   :reproduction_start_age => 1,
+#   :initial_energy => 1,
+#   :reproduction_start_age => 2,
 #   :reproduction_end_age => 5,
+#   :abiotic_variance => 1.0,
+#   :biotic_variance => 1.0,
+#   :mating_scheme => 0
 # )
 
 # species2 = Dict(
@@ -76,45 +71,72 @@ using EvoDynamics
 #   :abiotic_phenotypes => [1, 2],
 #   :biotic_phenotypes => [3, 4],
 #   :migration_phenotype => 5,
-#   :migration_threshold => 3.4,
+#   :migration_threshold => 1.4,
 #   :vision_radius => 1,
 #   :check_fraction => 0.5,
 #   :ploidy => 1,
-#   :epistasis_matrix => [1.0 -0.01 -0.01 -0.34 -0.42 0.18 -0.09 0.13; 0.28 1.0 -0.21 -0.11 0.12 -0.46 0.21 -0.39; 0.0 0.28 1.0 0.1 0.09 0.3 0.1 0.17; 0.0 -0.42 0.05 1.0 -0.11 -0.07 -0.27 0.43; 0.31 0.47 -0.42 -0.34 1.0 -0.4 0.16 0.11; -0.19 0.29 -0.33 -0.49 0.17 1.0 -0.1 -0.28; -0.43 0.38 -0.06 0.39 0.21 -0.5 1.0 -0.08; 0.26 0.27 -0.44 -0.08 0.47 -0.27 -0.27 1.0],
-#   :pleiotropy_matrix => Bool[1 0 0 0 1 1 0 0; 1 0 0 1 1 0 1 1; 0 1 1 1 1 1 0 1; 1 0 1 0 0 1 0 1; 1 1 0 1 1 1 1 1],
-#   :growth_rate => 1.2,
+#   :epistasis_matrix => [
+#     1.0 -0.01 -0.01 -0.34 -0.42 0.18 -0.09 0.13
+#     0.28 1.0 -0.21 -0.11 0.12 -0.46 0.21 -0.39
+#     0.0 0.28 1.0 0.1 0.09 0.3 0.1 0.17
+#     0.0 -0.42 0.05 1.0 -0.11 -0.07 -0.27 0.43
+#     0.31 0.47 -0.42 -0.34 1.0 -0.4 0.16 0.11
+#     -0.19 0.29 -0.33 -0.49 0.17 1.0 -0.1 -0.28
+#     -0.43 0.38 -0.06 0.39 0.21 -0.5 1.0 -0.08
+#     0.26 0.27 -0.44 -0.08 0.47 -0.27 -0.27 1.0],
+#   :pleiotropy_matrix => Bool[
+#     1 0 0 0 1 1 0 0
+#     1 0 0 1 1 0 1 1
+#     0 1 1 1 1 1 0 1
+#     1 0 1 0 0 1 0 1
+#     1 1 0 1 1 1 1 1],
+#   :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.5,
+#   :selection_coefficient => 0.1,
 #   :mutation_probabilities => [0.9, 0.0, 0.0],
 #   :mutation_magnitudes => [0.05, 0.0, 0.01],
-#   :N => [1000, 0, 0, 0],
+#   :N => [100, 0, 0, 0],
 #   :environmental_noise => 0.01,
-#   :optimal_phenotypes => optphens3,
-#   :bottleneck_function => bn,
+#   :optimal_phenotypes => [fill([1.4, 1.6], space...) for t in 0:generations],
+#   :bottlenecks => [fill(0.0, space...) for t in 0:generations],
 #   :age => 3,
-#   :recombination => 1,
+#   :recombination => 0,
 #   :initial_energy => 0,
 #   :reproduction_start_age => 1,
 #   :reproduction_end_age => 3,
+#   :abiotic_variance => 1.0,
+#   :biotic_variance => 2.0,
+#   :mating_scheme => 0
 # )
 
-# ## 3. Model parameters
+# ## 2. Model parameters
 
-# #NB this dict should be called model_parameters
+# #NB this dictionary should be called model_parameters
 # model_parameters = Dict(
 #   :species => [species1, species2],
-#   :generations => 14,
-#   :space => (2, 2),
+#   :generations => generations,
+#   :space => space,
 #   :metric => "chebyshev",
 #   :periodic => false,
-#   :resources => env_resources2,
-#   :interactions => [0.1 0.0; 0.0 -1.0],
-#   :food_sources => [1.0 0.5; 0.0 0.0],
-#   :seed => nothing
+
+#   :resources => [
+#     [2000 2200
+#      1830 1900] for i in 0:generations],
+
+#   :interactions => [0.0 1.0
+#                     1.0 0.0],
+
+#   :food_sources => [1.0 0.0
+#                     1.0 0.0], 
+
+#   :seed => 32
 # )
 # ```
 
 param_file = "../../examples/paramfile2.jl"  #hide
 agentdata, modeldata, model = runmodel(param_file);
 
-modeldata
\ No newline at end of file
+modeldata
+
+plot(1:size(modeldata, 1), getindex.(modeldata[:, 3], 1), xlabel="Time", ylabel="N", label="Prey")
+plot!(1:size(modeldata, 1), getindex.(modeldata[:, 3], 2), label="Predator")
\ No newline at end of file
diff --git a/examples/paramfile1.jl b/examples/paramfile1.jl
index f0ec348..b37a56e 100644
--- a/examples/paramfile1.jl
+++ b/examples/paramfile1.jl
@@ -24,7 +24,7 @@ species1 = Dict(
   :N => [100],
   :environmental_noise => 0.01,
   :optimal_phenotypes => [fill([1.5 for p in 1:1], space...) for t in 0:generations],
-  :bottleneck_function => [fill(0.0, space...) for t in 0:generations],
+  :bottlenecks => [fill(0.0, space...) for t in 0:generations],
   :age => 2,
   :recombination => 0,
   :initial_energy => 0,
diff --git a/examples/paramfile2.jl b/examples/paramfile2.jl
index 87daf01..c8354a6 100644
--- a/examples/paramfile2.jl
+++ b/examples/paramfile2.jl
@@ -1,4 +1,4 @@
-generations = 14
+generations = 30
 space = (2, 2)
 
 ## 1. Species parameters
@@ -44,7 +44,7 @@ species1 = Dict(
   :N => [1000, 0, 0, 0],
   :environmental_noise => 0.01,
   :optimal_phenotypes => [fill([1.4 for p in 1:1], space...) for t in 0:generations],
-  :bottleneck_function => [fill(0.0, space...) for t in 0:generations],
+  :bottlenecks => [fill(0.0, space...) for t in 0:generations],
   :age => 5,
   :recombination => 1,
   :initial_energy => 1,
@@ -91,14 +91,14 @@ species2 = Dict(
   :N => [100, 0, 0, 0],
   :environmental_noise => 0.01,
   :optimal_phenotypes => [fill([1.4, 1.6], space...) for t in 0:generations],
-  :bottleneck_function => [fill(0.0, space...) for t in 0:generations],
+  :bottlenecks => [fill(0.0, space...) for t in 0:generations],
   :age => 3,
   :recombination => 0,
   :initial_energy => 0,
   :reproduction_start_age => 1,
   :reproduction_end_age => 3,
   :abiotic_variance => 1.0,
-  :biotic_variance => 1.0,
+  :biotic_variance => 2.0,
   :mating_scheme => 0
 )
 
@@ -119,7 +119,7 @@ model_parameters = Dict(
                     1.0 0.0],
 
   :food_sources => [1.0 0.0; 
-                    2.5 0.0],
+                    1.0 0.0],
 
-  :seed => nothing
+  :seed => 32
 )
\ No newline at end of file
diff --git a/src/check_params.jl b/src/check_params.jl
index 13ba98e..4e62348 100644
--- a/src/check_params.jl
+++ b/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, :bottleneck_function, :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]
   model_keys = [:species, :generations, :space, :metric, :periodic, :resources, :interactions, :food_sources, :seed]
   species = d[:species]
   for sp in species
@@ -41,7 +41,7 @@ function check_param_shapes(d)
     @assert typeof(dd[:recombination]) <: Real "recombination of species $spname should be type numeric"
     @assert typeof(dd[:initial_energy]) <: Real "Initial energy should be a number"
     # bottleneck should be nothing or array/string
-    @assert typeof(dd[:bottleneck_function]) <: AbstractArray || isnothing(dd[:bottleneck_function]) "bottleneck function for species $spname is not an `Array` or `nothing`"
+    @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[: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"
diff --git a/src/reproduction.jl b/src/reproduction.jl
index 578e1d1..47ed3bf 100644
--- a/src/reproduction.jl
+++ b/src/reproduction.jl
@@ -87,7 +87,7 @@ end
 
 """
 Sexual reproduction for diploid individuals.
-Whether it is assortative or dissortative depends on the sign of the corresponding element in the interaction matrix. 
+Whether it is assortative or disassortative depends on the sign of the corresponding element in the interaction matrix. 
 """
 function reproduce!(ag1::Ind, ag2::Ind, model::ABM)
   if model.mating_schemes[ag1.species] == 0.0
diff --git a/src/simulation.jl b/src/simulation.jl
index d1ec617..7e00fe1 100644
--- a/src/simulation.jl
+++ b/src/simulation.jl
@@ -166,7 +166,7 @@ function create_properties(dd)
   names = Dict(i => allspecies[i][:name] for i in 1:nspecies)
   recombination = [Poisson(allspecies[i][:recombination]) for i in 1:nspecies]
   initial_energy = [AbstractFloat(allspecies[i][:initial_energy]) for i in 1:nspecies]
-  bottlenecks = [allspecies[i][:bottleneck_function] for i in 1:nspecies]
+  bottlenecks = [allspecies[i][:bottlenecks] for i in 1:nspecies]
   repro_start = [allspecies[i][:reproduction_start_age] for i in 1:nspecies]
   repro_end = [allspecies[i][:reproduction_end_age] for i in 1:nspecies]
   resources = dd[:resources][1]
diff --git a/test/params.jl b/test/params.jl
index 9d5f3e6..09e90e7 100644
--- a/test/params.jl
+++ b/test/params.jl
@@ -41,10 +41,10 @@ species1 = Dict(
   :age => 4,  # max age
   :reproduction_start_age => 1,
   :reproduction_end_age => 4,
-  :mating_scheme => 0, # 0 means number of children between a pair is independent of the phenotype of the pair. 1 means the more similar they are, the more children they will have (assortative mating). -1 means the more dissimilar they are, the more children they will have (dissortative mating).
+  :mating_scheme => 0, # 0 means number of children between a pair is independent of the phenotype of the pair. 1 means the more similar they are, the more children they will have (assortative mating). -1 means the more dissimilar they are, the more children they will have (disassortative mating).
   :recombination => 1, # Mean of a Poisson distributions for number of crossing overs
-  :initial_energy => 0, # A parameter for parental care of infants. Values more than 0 indicate that newly born individuals can survive for a number of times without requiring food from the environment/other species. The consumption rate (i.e. how many generations this initial energy suffices) is determined by the sum of the corresponding rows in "food sources". Note that having initial energy larger than zero can lead to infinite population growth because agents without food can reproduce and their offsprings also reproduce without food. Use with care, for example, when start age of reproduction is larger than 1.
-  :bottleneck_function => [fill(0.0, space...) for t in 0:generations]  # an array of arrays with probablity of external death at each site and generation.
+  :initial_energy => 0, # A parameter for parental care of infants. Values more than 0 indicate that newly born individuals can survive for a number of times without requiring food from the environment/other species. Note that having initial energy larger than zero can lead to infinite population growth because agents without food can reproduce and their offsprings also reproduce without food. Use with care, for example, when start age of reproduction is larger than 1.
+  :bottlenecks => [fill(0.0, space...) for t in 0:generations]  # an array of matrices with probablity of external death at each site and generation.
 )
 
 species2 = Dict(
@@ -76,7 +76,7 @@ species2 = Dict(
   :mating_scheme => 0,
   :recombination => 1,
   :initial_energy => 0,
-  :bottleneck_function => [fill(0.0, space...) for t in 0:generations]
+  :bottlenecks => [fill(0.0, space...) for t in 0:generations]
 )
 
 ## 2. Model parameters
@@ -90,7 +90,7 @@ model_parameters = Dict(
   :space => space,
   :metric => "chebyshev",  # how many neighbors a space site has. "chebyshev" metric means that the r-neighborhood of a position are all positions within the hypercube having side length of 2*floor(r) and being centered in the origin position. "euclidean" metric means that the r-neighborhood of a position are all positions whose cartesian indices have Euclidean distance ≤ r from the cartesian index of the given position.
   :periodic => false,  # whether boundaries of the space are connected
-  :resources => [env_resources(x) for x in 0:generations],  # a function that returns a vector of integers for the available resources per site at a given time. It accepts only the model object as the input.
+  :resources => [env_resources(x) for x in 0:generations],
   :interactions => [0.1 1.0; 1.0 1.0],  # How individuals from different species interact. value  is probability of interaction (between 0 and 1). Sign is the direction of interaction where positive means similar individuals interact more strongly and negative is dissimilar ones tend to interact more. If you want full interaction, use 1 or -1 for nondiagonals. Diagonals are used for competition/cooperation.
   :food_sources => [1.0 0.0; 1.7 0.0],  # What each species feeds on (consumption rate). Non-zero diagonal means the food resource is from the environment. It will be read from rows (species order) to columns (species order).
   :seed => 3