Configs.md

@@ -8,11 +8,10 @@ Configs
 ### 1.1.1  Global equatorial ancestor (Pantropical diversity, Hagen and Skeels et al 2021)
 
 This function intialises a single species that occupies all non-arid equatorial grid cells.
-Requires 3 traits: dispersal, temp, temp_width
-Specify paramters:  params$initialAbundance
+
+Required traits: dispersal, temp, temp_width
 
 ~~~R
-initial_abundance = params$initialAbundance
 
 create_ancestor_species <- function(landscape, config) {
 
@@ -34,12 +33,11 @@ create_ancestor_species <- function(landscape, config) {
 }
 ~~~
 
-### 1.1.2
-
-Global ancestor (Evolutionary Speed Hypothesis, Skeels et al 2021)
+### 1.1.2 Global ancestor (Evolutionary Speed Hypothesis, Skeels et al 2021)
 
 This function intialises a single species that occupies all global grid cells.
-Requires 4 traits: dispersal, temp, prec, body_size
+
+Required traits: dispersal, temp, prec, body_size
 
 ~~~R
 
@@ -74,7 +72,11 @@ create_ancestor_species <- function(landscape, config) {
 ### 1.2.1 Weibull dispersal kernel (Pantropical diversity, Hagen and Skeels et al 2021)
 
 This function draws the dispersal kernel from a Weibull distriution with shape and scale parameters.
-Specify paramters: params$WeibullShape, params$WeibullScale
+
+Explored paramters: 
+params$WeibullShape = [2,15]
+params$WeibullScale = 222 # requires maximum of 2 decimal degree resolution to
+
 ~~~R
 
 get_dispersal_values <- function(n, species, landscape, config) {
@@ -88,8 +90,12 @@ get_dispersal_values <- function(n, species, landscape, config) {
 ### 1.3.1 Species-level Brownian Motion (Pantropical diversity, Hagen and Skeels et al 2021)
 
 This function evolves a single temperature-niche trait under a random walk (Brownian Motion) model.
-Requires 1 trait: temp
-Specify paramters:   params$Tev
+
+Explored paramters:   
+params$Tev = [0.001, 0.02]
+
+Required trait: temp
+Required environment: temperature (temp)
 
 ~~~R
 apply_evolution <- function(species, cluster_indices, landscape, config) {
@@ -111,10 +117,15 @@ apply_evolution <- function(species, cluster_indices, landscape, config) {
 
 This function evolves a temperature-niche trait under a Ornstein Uhlenbeck model of trait evolution and evolves a body size trait under a Brownian Motion model. The OU model has a random walk element as in the Brownian motion model and also a parameter (psi) that draws the value towards an optimal value (based on the mean temperature value across the clusters geographic range). When psi=0 OU = BM. The OU and BM models each have an independant rate parameter (sigma squared). This function evolves each geographic cluster seperately (instead of each species speerately or each grid cell seperately).
 
-Requires 2 traits: temp and body_size
-Specify paramters:  params$sigma_squared_bs, params$sigma_squared_t, params$psi
+Explored paramters:  
+params$sigma_squared_bs = [0.001, 0.02]
+params$sigma_squared_t = [0.001, 0.015]
+params$psi = 0
 
+Required traits: temp and body_size
+Requires environment: temperature (temp) 
 
+~~~R
 sigma_squared_bs <-  params$sigma_squared_bs
 sigma_squared_t <-  params$sigma_squared_t
 psi <-  params$psi
@@ -144,7 +155,7 @@ apply_evolution <- function(species, cluster_indices, landscape, config) {
   return(traits)
   
 }
-
+~~~
 
 # 1.4 Speciation
 
@@ -152,8 +163,8 @@ apply_evolution <- function(species, cluster_indices, landscape, config) {
 
 This function increases divergence between populations by 0.1 at each timestep and speciation occurs after divergence crosses threshold S. Values < 1 mean populations seperate more slowly than they come back toegther, as populations that come into secondary contact seperate at rate of 1 per time step.
 
-Specify paramters:   params$S
-Requires 2 traits: temp, temp_width
+Explored paramters: 
+params$divergence_threshold = [1.5,3] 
 
 ~~~R
 divergence_threshold = params$S
@@ -163,6 +174,133 @@ get_divergence_factor <- function(species, cluster_indices, landscape, config) {
 }
 ~~~
 
+
+### 1.4.2 Exponential temperature-dependant divergence (Evolutionary Speed Hypothesis, Skeels et al 2021)
+
+This function increases divergence between populations each timestep as an exponential function of the Sum of the mean temperature values across each diverging populations geographic ranges. The exponent is scaled with the parameter, lambda.  Speciation occurs after divergence crosses threshold S. Values < 1 mean populations seperate more slowly than they come back toegther, as populations that come into secondary contact seperate at rate of 1 per time step.
+
+Explored paramters: 
+params$divergence_threshold = [2,10] 
+params$lambda = [2,5]
+
+Required traits: body_size
+Required environment: temperature (temp) 
+
+~~~R
+divergence_threshold =  params$divergence_threshold
+lambda =  params$lambda 
+
+get_divergence_factor <- function(species, cluster_indices, landscape, config) {
+  divergences <- vector("numeric", length(unique(cluster_indices)))
+  
+  divergence_matrix <- matrix(0, ncol=length(unique(cluster_indices)),  nrow=length(unique(cluster_indices)))
+  cluster_indices_order <- unique(cluster_indices)[order(unique(cluster_indices))]
+  # order them, find pairwise temp dist
+  for(i_x in 1:length(cluster_indices_order)){
+    mean_temp_i <- mean(landscape$environment[names(species$abundance[which(cluster_indices %in% cluster_indices_order[i_x])]),"temp"], na.rm=T)
+    for(j_x in 1:length(cluster_indices_order)){
+      mean_temp_j <- mean(landscape$environment[names(species$abundance[which(cluster_indices %in% cluster_indices_order[j_x])]),"temp"], na.rm=T)
+      divergence_matrix[i_x, j_x] <- (sum(mean_temp_i, mean_temp_j)/(2))^lambda
+    }
+  }
+  rownames(divergence_matrix) <- cluster_indices_order
+  colnames(divergence_matrix) <- cluster_indices_order
+  
+  if(length(divergence_matrix) == 1){ divergence_matrix <- as.numeric(divergence_matrix)}else{
+    
+    diag(divergence_matrix) <- 0
+  }
+  
+  return(divergence_matrix)
+}
+~~~
+
+### 1.4.3 Exponential body size-dependant divergence (Evolutionary Speed Hypothesis, Skeels et al 2021)
+
+This function increases divergence between populations each timestep as an exponential function of 1 - the sum of the mean body size values across each diverging populations geographic ranges. The exponent is scaled with the parameter, lambda.  Speciation occurs after divergence crosses threshold S. Values < 1 mean populations seperate more slowly than they come back toegther, as populations that come into secondary contact seperate at rate of 1 per time step.
+
+Explored paramters: 
+params$divergence_threshold = [2,10] 
+params$lambda = [2,5]
+
+Required traits: body_size
+Required environment: temperature (temp) 
+~~~R
+divergence_threshold =  params$divergence_threshold
+lambda =  params$lambda 
+
+get_divergence_factor <- function(species, cluster_indices, landscape, config) {
+  divergences <- vector("numeric", length(unique(cluster_indices)))
+  
+  divergence_matrix <- matrix(0, ncol=length(unique(cluster_indices)),  nrow=length(unique(cluster_indices)))
+  
+  cluster_indices_order <- unique(cluster_indices)[order(unique(cluster_indices))]
+  # order them, find pairwise temp dist
+  for(i_x in 1:length(cluster_indices_order)){
+    # make it 1 minus body size, so high values = low body size, therefore small bpdy size species have higher rates
+    mean_body_size_i <- 1-mean(species$traits[names(species$abundance[which(cluster_indices %in% cluster_indices_order[i_x])]),"body_size"], na.rm=T)
+    for(j_x in 1:length(cluster_indices_order)){
+      mean_body_size_j <- 1-mean(species$traits[names(species$abundance[which(cluster_indices %in% cluster_indices_order[j_x])]),"body_size"], na.rm=T)
+      divergence_matrix[i_x, j_x] <- ((sum(mean_body_size_i, mean_body_size_j)/2))^lambda
+    }
+  }
+  rownames(divergence_matrix) <- cluster_indices_order
+  colnames(divergence_matrix) <- cluster_indices_order
+  
+  if(length(divergence_matrix) == 1){ divergence_matrix <- as.numeric(divergence_matrix)}else{
+    
+    diag(divergence_matrix) <- 0
+  }
+  
+  return(divergence_matrix)
+}
+
+~~~
+
+
+### 1.4.3 Exponential temperature + body size-dependant divergence (Evolutionary Speed Hypothesis, Skeels et al 2021)
+
+This function increases divergence between populations each timestep as an exponential function of the sum of the mean temperature values across each diverging populations geographic ranges and 1 - the sum of mean body size values. The exponent is scaled with the parameter, lambda.  Speciation occurs after divergence crosses threshold S. Values < 1 mean populations seperate more slowly than they come back toegther, as populations that come into secondary contact seperate at rate of 1 per time step.
+
+Explored paramters: 
+params$divergence_threshold = [2,10] 
+params$lambda = [2,5]
+
+Required traits: body_size and temperature
+Required environment: temperature (temp) 
+
+~~~R
+divergence_threshold =  params$divergence_threshold
+lambda =  params$lambda 
+
+get_divergence_factor <- function(species, cluster_indices, landscape, config) {
+  divergences <- vector("numeric", length(unique(cluster_indices)))
+  
+  divergence_matrix <- matrix(0, ncol=length(unique(cluster_indices)),  nrow=length(unique(cluster_indices)))
+  
+  cluster_indices_order <- unique(cluster_indices)[order(unique(cluster_indices))]
+  # order them, find pairwise temp dist
+  for(i_x in 1:length(cluster_indices_order)){
+    # make it 1 minus body size, so high values = low body size, therefore small bpdy size species have higher rates
+    mean_body_size_i <- 1-mean(species$traits[names(species$abundance[which(cluster_indices %in% cluster_indices_order[i_x])]),"body_size"], na.rm=T)
+    for(j_x in 1:length(cluster_indices_order)){
+      mean_body_size_j <- 1-mean(species$traits[names(species$abundance[which(cluster_indices %in% cluster_indices_order[j_x])]),"body_size"], na.rm=T)
+      divergence_matrix[i_x, j_x] <- ((sum(mean_body_size_i, mean_body_size_j)/2))^lambda
+    }
+  }
+  rownames(divergence_matrix) <- cluster_indices_order
+  colnames(divergence_matrix) <- cluster_indices_order
+  
+  if(length(divergence_matrix) == 1){ divergence_matrix <- as.numeric(divergence_matrix)}else{
+    
+    diag(divergence_matrix) <- 0
+  }
+  
+  return(divergence_matrix)
+}
+
+~~~
+
 # 1.5 Ecology
 
 ### 1.5.1 Temperature and aridity thresholds (Pantropical diversity, Hagen and Skeels et al 2021)
@@ -178,4 +316,78 @@ apply_ecology <- function(abundance, traits, landscape, config) {
   abundance[which(landscape[, "prec"] > 0.5)] <- 0
   return(abundance)
 }
+~~~
+
+### 1.5.2 Environmental carrying capacities (Evolutionary Speed Hypothesis, Skeels et al 2021)
+
+This function using the temperature niche optima and the sites aridity value to determine the local abundance of each species. There are zero-sum dynamics operating so increases in one species will lead to decreases in otehr species if the site is at equilibrium.
+
+The size (N) of population i in site j is determined by a gaussian function of resource use efficiency - the distance between the temperature value in the site (Tj), and the population’s temperature optimum (Ti): 
+
+Nij = K * exp(-(Ti - Tj/omega)2 ) 		equation 1
+
+where omega is a parameter that determines the strength of environmental filtering, with small values leading to a sharper decline in abundance as the species temperature niche optimum (Ti) becomes more different from the temperature of the site (Tj). Nij will equal K in the absence of competitors if population i were perfectly adapted to the site. The carrying capacity for each site (K) was independent of temperature but decreased exponentially with the aridity index in each site (Aj) according to the function:
+
+K = Kc * exp(-1*Aj)  			equation 2
+
+Where Kc was a constant [30000]. We model a zero-sum game where sites have finite resources available and this places ecological limits on the maximum number of individuals in a site across populations of all species present (Nj). Therefore, if the total community abundance across all species, Nj, is at or above K the addition of any new species will decrease resources available to all present species. In saturated communities (where Nj >= K), the realised abundance of each species (N̂ij) is apportioned according to the resource use efficiency of each species: 
+
+N̂ij=Nij∗min(Nj,K)/Nj 			equation 3
+
+Local extinction occurs deterministically if N̂ij = 0 or stochastically as a sigmoidal function of N̂ij:
+
+1/(1+exp(-μd * (μt -  N̂ij))))		equation 4
+
+Where μt is the population size threshold below which extirpation in site j becomes more likely and μd is the rate of decay of the function. Extinction of a species occurs when it no longer occupies any sites.
+
+Explored parameters:
+params$K = 30000
+params$omega   [0.01, 0.035]
+params$aridity_cost = 2
+params$inflexion = 1000
+params$decay = 0.1
+
+Rrequired traits: temp, body_size
+Required Environments: tmeperature (temp), aridity (prec)
+~~~R
+
+
+K_opt_max <-  params$K    # environmental filtering parameter (maximum abundance of all individuals in a mesic cell)
+omega <- params$omega     # environmental filtering parameter (higher values equal less restrictive niche)
+aridity_cost <-  params$aridity_cost # environmental filtering parameter (higher values equal less abundance in arid cells)
+
+x0 <- params$inflexion # extinction parameter: inflexion point of extinction probability curve
+decay <- params$decay  # extinction parameter: inflexion point of extinction probability curve
+
+apply_ecology <- function(abundance, traits, landscape, config) {
+  
+  # Realised carrying capacity (reduced in arid cells)
+  K <- K_opt_max * exp(1-aridity_cost*landscape[, "prec"])
+  
+  # difference between landscape and species optimum temperature niche
+  diff_t <- abs(traits[, "temp"] - landscape[, "temp"])
+
+  # potential population size based on resource use efficieny 
+  # omega is strength of environmental filtering
+  # will equal K when tthe species is perfectly adapted
+  # modified from McPeek 2008/2007
+  Nij <- K * exp(-(diff_t/omega)^2)
+  
+  # These lines make carrying capcaity a zero-sum game - species abundance is scaled based on their resourse use efficiency (Nij)
+  Nij_sum <- sum(Nij)
+  Nij_hat <- Nij * (sort(c(Nij_sum, K),partial=1)[1] / Nij_sum)
+  Nij_hat[which(is.na(Nij_hat))] <- 0
+  # now do an extinction filter based on population size
+  prob_extinction <- (1/(1+exp(-decay*(x0 -  Nij_hat))))
+  Nij_hat[which(sapply(prob_extinction, FUN=function(x){sample(c(0,1),1, prob= c(x, 1-(x)))})==0)] <- 0 
+  
+  # absolute precipitation threshold 
+  #if(landscape[, "prec"] > 0.5){Nij <- rep(0, length(Nij))}
+  
+  return(Nij_hat)
+}
+
+
+
+
 ~~~
\ No newline at end of file