2. Explain example config
Here we will go bit by bit of the config used at the introduction vignette vignette. We will follow the structure presented above.
Content of config_.empty_southamerica.R
Metadata
The first section of the config file contains metadata. It informs the user when and by whom the file has been created, which landscape it was made for and informs about associated publications and implemented processes or exploration intentions of the specific configuration.
# Version: 1.0 # Author: Oskar Hagen # Date: 26.10.2020 # Landscape: SouthAmerica # Publications: R-package gen3sis # Description: Example config used at the introduction vignette and similar to case study global configs in Hagen et al. 2020. # O. Hagen, B. Flück, F. Fopp, J.S. Cabral, F. Hartig, M. Pontarp, T.F. Rangel, L. Pellissier. gen3sis: the general engine for eco-evolutionary simulations on the origins of biodiversity.
Settings
We start our simulation 40 Ma, which corresponds to time step 40 since 1 time-step = 1 myr. We end it a the latest available time-step at the landscape, thus start_time is set to 40 and end_time is set to NA. We limit the maximum total number of species alive in a simulation to 50000 with max_number_of_species and the maximum number of species within one site to 10000 with max_number_of_coexisting_species. Like so, simulations that could possible go off hand and generate too many species are stopped and properly flagged at the sgen3sis summary object. We define which traits we will consider in our simulation with traits_names, in our example, species only have a dispersal ability and a optimum temperature traits. For our model to work, we normalize the environmental values, according to the range informed at environmental_ranges.
Observer
With the observer function, changes over time in any abiotic or biotic information of the virtual world can be recorded by defining the outputs that are saved at specified time steps, and results can be saved and plotted in real-time as the model runs. In our observer function we save the species and plot the richness. Like so, plots are printed as the simulation progresses and also saved inside the plot folder inside output.
Combined plots can be generated inside the observer function as for example:
end_of_timestep_observer = function(data, vars, config){ par(mfrow=c(2,3)) plot_raster_single(data$landscape$environment[,"temp"], data$landscape, "temp", NA) plot_raster_single(data$landscape$environment[,"arid"], data$landscape, "arid", NA) plot_raster_single(data$landscape$environment[,"area"], data$landscape, "area", NA) plot_richness(data$all_species, data$landscape) plot_species_presence(data$all_species[[1]], data$landscape) plot(0,type='n',axes=FALSE,ann=FALSE) mtext("STATUS",1) }
Initialization
The initialization function creates the ancestor species at the start of the simulation. Users can define the number of ancestor from a single to multiple species, their distribution within the paleolandscape and their trait values. In our example, the initial abundance of the 10 ancestor species, as well as the abundance of species on newly colonized sites is 1. We restrict the region of colonization to South America by limiting the range to a spatial extent of c(-95, -24, -68, 13). We then create ten different species, each of which is randomly inhabiting one of the sites of the continent. Thereby all of the sites have an equal probability of being inhabited by a species in the beginning of the simulation. However since the sites vary in their size, the habitation probability per area is not the same for all regions. If you want to correct for that you can consider using an equal area transformation of your coordinates or tweeking sampling functions to account for area correction. We set the optimum temperature trait temp of each population of the initial species to the environmental temperature. This means, that each population is adapted to the local conditions, as we will see later again at the ecology function. Dispersal is set to one, since it is not allowed to evolve, we use it’s maximum value. The function create_ancestor_species is expected to return a list of new species, i.e. the ancestor(s) one(s).
initial_abundance = 1 create_ancestor_species <- function(landscape, config) { range <- c(-95, -24, -68, 13) co <- landscape$coordinates selection <- co[, "x"] >= range[1] & co[, "x"] <= range[2] & co[, "y"] >= range[3] & co[, "y"] <= range[4] new_species <- list() for(i in 1:10){ initial_cells <- rownames(co)[selection] initial_cells <- sample(initial_cells, 1) new_species[[i]] <- create_species(initial_cells, config) #set local adaptation to max optimal temp equals local temp new_species[[i]]$traits[ , "temp"] <- landscape$environment[initial_cells,"temp"] new_species[[i]]$traits[ , "dispersal"] <- 1 plot_species_presence(landscape, species=new_species[[i]]) } return(new_species) }
Core Processes
Dispersal
The dispersal function iterates over all species populations and determines the connectivity between sites and the colonization of new sites in the grid cell. In our example, species dispersal between time-steps is stochastic and similar to all species. It follows a Weibull distribution with shape 1.5 and scale 133. For different dispersal between species, the dispersal trait should be defined. Note that dispersal can be made deterministic by making the function return a fix value.
The histogram of 100 draws from the dispersal function used above.
Speciation
The speciation iterates over every species separately, registers populations’ geographic occupancy (species range), and determines when geographic isolation between population clusters is higher than a user-defined threshold, triggering a lineage splitting event of cladogenesis. The clustering of occupied sites is based on the species’ dispersal capacity and the landscape connection costs. Over time, disconnected clusters gradually accumulate incompatibility, analogous to genetic differentiation. When the divergence between clusters is above the speciation threshold, those clusters become two or more distinct species, and a divergence matrix reset follows. On the other hand, if geographic clusters come into secondary contact before the speciation occurs, they coalesce and incompatibilities are gradually reduced to zero. In our example, speciation takes place after 2 time-steps of isolation and the divergence increase is the same for all species as indicated by get_divergence_threshold. Since our landscape consists of 1 myr time-steps, these 2 time-steps correspond to a span of 2 myr.
Evolution
In the evolution function, clustered populations (exchanging genes) had their trait homogenized and weighted by abundance, meaning that a trait of a population that is doing well in a site, as dictated by the ecology function, will contribute more to the average trait of a cluster. Later, populations mutate based on a normal distribution with standard deviation 0.001, possibly increasing or decreasing species optimum temperature.
# mutate the traits of a species and return the new traits matrix apply_evolution <- function(species, cluster_indices, landscape, config) { trait_evolutionary_power <- 0.001 traits <- species[["traits"]] cells <- rownames(traits) #homogenize trait based on abundance for(cluster_index in unique(cluster_indices)){ cells_cluster <- cells[which(cluster_indices == cluster_index)] mean_abd <- mean(species$abundance[cells_cluster]) weight_abd <- species$abundance[cells_cluster]/mean_abd traits[cells_cluster, "temp"] <- mean(traits[cells_cluster, "temp"]*weight_abd) } #mutations mutation_deltas <-rnorm(length(traits[, "temp"]), mean=0, sd=trait_evolutionary_power) traits[, "temp"] <- traits[, "temp"] + mutation_deltas return(traits) }
Ecology
The ecology function determines the abundance or presence of the populations in occupied sites of each species. The function iterates over all occupied sites and updates the species population abundances or presences on the basis of local environmental values, updated co-occurrence patterns and species traits. The function takes as input the species abundance, species trait, species divergence and clusters, and the landscape values. In our example, we calculate abundances in a site based on how close the species is to the site temperature. We scale the values to avoid small numbers and apply a carrying capacity based on aridity and temperature corrected by the area of a site. If abundances are below 1, species are considered extinct, if total abundance in a site is above the carrying capacity, small abundances are removed progressively and randomly distributed across the present species until total abundance is smaller or equal to the carrying capacity.
apply_ecology <- function(abundance, traits, landscape, config) { abundance_scale = 10 abundance_threshold = 1 #abundance treashold survive <- abundance>=abundance_threshold abundance[!survive] <- 0 abundance <- (( 1-abs( traits[, "temp"] - landscape[, "temp"]))*abundance_scale)*as.numeric(survive) #abundance thhreashold abundance[abundance<abundance_threshold] <- 0 k <- ((landscape[,"area"]*(landscape[,"arid"]+0.1)*(landscape[,"temp"]+0.1))*abundance_scale^2) total_ab <- sum(abundance) subtract <- total_ab-k if (subtract > 0) { # print(paste("should:", k, "is:", total_ab, "DIFF:", round(subtract,0) )) while (total_ab>k){ alive <- abundance>0 loose <- sample(1:length(abundance[alive]),1) abundance[alive][loose] <- abundance[alive][loose]-1 total_ab <- sum(abundance) } #set negative abundances to zero abundance[!alive] <- 0 } return(abundance) }