Virtual plant simulations of Africa

EcoSISTEM was designed to scale to much larger areas, supporting many more species. As an illustrative example, here we simulate up to 50,000 plant species over Africa at an 80km grid scale, with a constant background environment of 25°C. When all species are given an equal fitness in the habitat, all 50,000 can co-exist over long time scales of over 100 years (Figure 2A). This can be run on a workstation with 24 threads in just under 5 hours.

We can also explore the behaviour of selective advantage of specialist species over generalists at these scales. When we introduce a specialist species into an African-sized landscape with an existing generalist, the specialist out-competes the generalist and spreads throughout the continent. The larger the selective advantage of the specialist, the faster it is able to invade and colonise across the landscape (Figure 1). These same dynamics can be seen when we introduce a specialist to the full complement of 50,000 species (Figure 1B-D).

SINGLE SPECIES

using EcoSISTEM
using EcoSISTEM.ClimatePref
using EcoSISTEM.Units
using Unitful
using Unitful.DefaultSymbols
using Distances
using StatsBase
using Plots
file = "Africa.tif"
africa = readfile(file, -25°, 50°, -35°, 40°)
active =  Matrix{Bool}(.!isnan.(africa'))

heatmap(active)

# Set up initial parameters for ecosystem
numSpecies = 1; grid = size(africa); req= 10.0kJ; individuals=0; area = 64e6km^2; totalK = 1000.0kJ/km^2

# Set up how much energy each species consumes
energy_vec = SolarRequirement(fill(req, numSpecies))


# Set rates for birth and death
birth = 0.6/year
death = 0.6/year
longevity = 1.0
survival = 0.0
boost = 1.0
# Collect model parameters together
param = EqualPop(birth, death, longevity, survival, boost)

# Create kernel for movement
kernel = fill(GaussianKernel(15.0km, 10e-10), numSpecies)
movement = AlwaysMovement(kernel, Torus())


# Create species list, including their temperature preferences, seed abundance and native status
opts = fill(274.0K, numSpecies)
vars = fill(0.5K, numSpecies)
traits = GaussTrait(opts, vars)
native = fill(true, numSpecies)
# abun = rand(Multinomial(individuals, numSpecies))
abun = fill(div(individuals, numSpecies), numSpecies)
sppl = SpeciesList(numSpecies, traits, abun, energy_vec,
    movement, param, native)
sppl.params.birth

# Create abiotic environment - even grid of one temperature
abenv = simplehabitatAE(274.0K, grid, totalK, area, active)


# Set relationship between species and environment (gaussian)
rel = Gauss{typeof(1.0K)}()

#Create ecosystem
eco = Ecosystem(sppl, abenv, rel)
rand_start = rand(findall(active), 1)[1]
eco.abundances.grid[1, rand_start[1], rand_start[2]] = 100

# EcoSISTEM Parameters
times = 100years; timestep = 1month; record_interval = 1month; repeats = 1
lensim = length(0years:record_interval:times)
abuns = zeros(Int64, numSpecies, prod(grid), lensim)
@time simulate_record!(abuns, eco, times, record_interval, timestep);

abuns = reshape(abuns[1, :, :, 1], grid[1], grid[2], lensim)

anim = @animate for i in 1:lensim
    africa_abun = Float64.(abuns[:, :, i])
    africa_abun[.!(active)] .= NaN
    heatmap(africa_abun, clim = (0, 700_000), background_color = :lightblue, background_color_outside=:white, grid = false, color = cgrad(:algae, scale = :exp), aspect_ratio = 1)
end
gif(anim, "examples/Biodiversity/Africa.gif", fps = 30)

#### SPECIALIST VERSUS GENERALIST ####

specialist_vars = [0.5K, 1.0K, 5.0K, 10.0K, 25.0K, 50.0K]
velocity = zeros(typeof(1.0km/month), length(specialist_vars))
rand_start = rand(findall(active), 1)[1]
for i in eachindex(specialist_vars)
    # Set up initial parameters for ecosystem
    numSpecies = 2; grid = size(africa); req= 10.0kJ; individuals=0; area = 64e6km^2; totalK = 1000.0kJ/km^2

    # Set up how much energy each species consumes
    energy_vec = SolarRequirement(fill(req, numSpecies))


    # Set rates for birth and death
    birth = 0.6/year
    death = 0.6/year
    longevity = 1.0
    survival = 0.1
    boost = 1.0
    # Collect model parameters together
    param = EqualPop(birth, death, longevity, survival, boost)

    # Create kernel for movement
    kernel = fill(GaussianKernel(15.0km, 10e-10), numSpecies)
    movement = AlwaysMovement(kernel, Torus())


    # Create species list, including their temperature preferences, seed abundance and native status
    opts = fill(274.0K, numSpecies)
    vars = [50.0K, specialist_vars[i]]
    traits = GaussTrait(opts, vars)
    native = fill(true, numSpecies)
    # abun = rand(Multinomial(individuals, numSpecies))
    abun = fill(div(individuals, numSpecies), numSpecies)
    sppl = SpeciesList(numSpecies, traits, abun, energy_vec,
        movement, param, native)
    sppl.params.birth

    # Create abiotic environment - even grid of one temperature
    abenv = simplehabitatAE(274.0K, grid, totalK, area, active)

    # Set relationship between species and environment (gaussian)
    rel = Gauss{typeof(1.0K)}()

    #Create ecosystem
    eco = Ecosystem(sppl, abenv, rel)
    eco.abundances.grid[1, rand_start[1], rand_start[2]] = 100

    # EcoSISTEM Parameters
    burnin = 100years; times = 100years; timestep = 1month; record_interval = 1month; repeats = 1
    lensim = length(0years:record_interval:times)
    simulate!(eco, burnin,timestep)
    eco.abundances.grid[2, rand_start[1], rand_start[2]] = 100
    abuns = zeros(Int64, numSpecies, prod(grid), lensim)
    @time simulate_record!(abuns, eco, times, record_interval, timestep);

    abuns = reshape(abuns[:, :, :, 1], numSpecies, grid[1], grid[2], lensim)
    origin = [rand_start[1], rand_start[2]]
    dest = findall(abuns[2, :, :, 1] .> 0)
    inst_velocity = map(1:lensim) do t
        dest = findall(abuns[2, :, :, t] .> 0)
        dists = [euclidean(origin, [dest[i][1], dest[i][2]]) for i in eachindex(dest)] .* getgridsize(eco)
        return maximum(dists)/month
    end
    velocity[i] = mean(inst_velocity)
end

plot(ustrip.(abs.(specialist_vars .- 50.0K)), ustrip.(velocity),
    xlab = "Selective advantage", ylab = "Invasion speed (km/month)",
    label = "", grid = false)

Figure 1: Invasive capacity of a specialist plant species versus a generalist. Selective advantage is the difference in niche width between the specialist and generalist, and invasion speed is calculated as the average distance travelled per month by the specialist.

ONE SPECIALIST VERSUS MANY GENERALISTS

using EcoSISTEM
using EcoSISTEM.ClimatePref
using EcoSISTEM.Units
using Unitful
using Unitful.DefaultSymbols
using JLD2
using Printf
file = "Africa.tif"
africa = readfile(file, -25°, 50°, -35°, 40°)
active =  Matrix{Bool}(.!isnan.(africa'))
# Set up initial parameters for ecosystem
numSpecies = 50_000; grid = size(africa); req= 10.0kJ; individuals=3*10^8; area = 64e6km^2; totalK = 1000.0kJ/km^2

# Set up how much energy each species consumes
energy_vec = SolarRequirement(fill(req, numSpecies))


# Set rates for birth and death
birth = 0.6/year
death = 0.6/year
longevity = 1.0
survival = 0.1
boost = 1.0
# Collect model parameters together
param = EqualPop(birth, death, longevity, survival, boost)

# Create kernel for movement
kernel = fill(GaussianKernel(15.0km, 10e-10), numSpecies)
movement = AlwaysMovement(kernel, Torus())


# Create species list, including their temperature preferences, seed abundance and native status
opts = fill(274.0K, numSpecies)
vars = fill(50.0K, numSpecies)
vars[50_000] = 0.5K
traits = GaussTrait(opts, vars)
native = fill(true, numSpecies)
# abun = rand(Multinomial(individuals, numSpecies))
abun = fill(div(individuals, numSpecies), numSpecies)
sppl = SpeciesList(numSpecies, traits, abun, energy_vec,
    movement, param, native)
sppl.params.birth

# Create abiotic environment - even grid of one temperature
abenv = simplehabitatAE(274.0K, grid, totalK, area, active)


# Set relationship between species and environment (gaussian)
rel = Gauss{typeof(1.0K)}()

#Create ecosystem
eco = Ecosystem(sppl, abenv, rel)
eco.abundances.matrix[50_000, :] .= 0

import EcoSISTEM.simulate!
function simulate!(eco::Ecosystem, times::Unitful.Time, timestep::Unitful.Time, cacheInterval::Unitful.Time, cacheFolder::String, scenario_name::String)
  time_seq = 0s:timestep:times
  counting = 0
  for i in eachindex(time_seq)
      update!(eco, timestep);
      # Save cache of abundances
      if mod(time_seq[i], cacheInterval) == 0year
          @save (joinpath(cacheFolder, scenario_name * (@sprintf "%02d.jld2" uconvert(NoUnits,time_seq[i]/cacheInterval))) abun = eco.abundances.matrix
      end
  end
end

# EcoSISTEM Parameters
burnin = 100years; times = 100years; timestep = 1month; record_interval = 12months;
lensim = length(0years:record_interval:times)
@time simulate!(eco, burnin, timestep)
rand_start = rand(findall(active), 1)[1]
eco.abundances.grid[50_000, rand_start[1], rand_start[2]] = 100
@time simulate!(eco, times, timestep, record_interval, "examples/Biodiversity", "Africa_run");

50,000 SPECIES COEXISTING

using EcoSISTEM
using EcoSISTEM.ClimatePref
using EcoSISTEM.Units
using Unitful
using Unitful.DefaultSymbols
using JLD2
using Printf

file = "Africa.tif"
africa = readfile(file, -25°, 50°, -35°, 40°)
active =  Matrix{Bool}(.!isnan.(africa'))
# Set up initial parameters for ecosystem
numSpecies = 50_000; grid = size(africa); req= 10.0kJ; individuals=3*10^8; area = 64e6km^2; totalK = 1000.0kJ/km^2

# Set up how much energy each species consumes
energy_vec = SolarRequirement(fill(req, numSpecies))


# Set rates for birth and death
birth = 0.6/year
death = 0.6/year
longevity = 1.0
survival = 0.1
boost = 1.0
# Collect model parameters together
param = EqualPop(birth, death, longevity, survival, boost)

# Create kernel for movement
kernel = fill(GaussianKernel(15.0km, 10e-10), numSpecies)
movement = AlwaysMovement(kernel, Torus())


# Create species list, including their temperature preferences, seed abundance and native status
opts = fill(274.0K, numSpecies)
vars = fill(50.0K, numSpecies)
traits = GaussTrait(opts, vars)
native = fill(true, numSpecies)
# abun = rand(Multinomial(individuals, numSpecies))
abun = fill(div(individuals, numSpecies), numSpecies)
sppl = SpeciesList(numSpecies, traits, abun, energy_vec,
    movement, param, native)
sppl.params.birth

# Create abiotic environment - even grid of one temperature
abenv = simplehabitatAE(274.0K, grid, totalK, area, active)


# Set relationship between species and environment (gaussian)
rel = Gauss{typeof(1.0K)}()

#Create ecosystem
eco = Ecosystem(sppl, abenv, rel)

# EcoSISTEM Parameters
burnin = 10years; times = 100years; timestep = 1month; record_interval = 12months;
lensim = length(0years:record_interval:times)
@time simulate!(eco, burnin, timestep)
@time simulate!(eco, times, timestep, record_interval, "examples/Biodiversity", "Africa_run_coexist");

using JLD2
using Plots
using Diversity
abuns = @load "examples/Biodiversity/Africa_run_coexist100.jld2" abun
meta = Metacommunity(abuns)
div = norm_sub_alpha(meta, 0)
sumabuns = reshape(div[!, :diversity], 100, 100)
heatmap(sumabuns,
    background_color = :lightblue,
    background_color_outside=:white,
    grid = false, color = :algae,
    aspect_ratio = 1, layout = (@layout [a b; c d]),
    clim = (0, 50_000), margin = 0.5 * Plots.mm,
    title = "A", titleloc = :left)

abuns = @load "examples/Biodiversity/Africa_run50.jld2" abun
meta = Metacommunity(abuns)
div = norm_sub_alpha(meta, 0)
sumabuns = reshape(div[!, :diversity], 100, 100)
heatmap!(sumabuns,
    background_color = :lightblue,
    background_color_outside=:white,
    grid = false, color = :algae,
    aspect_ratio = 1, subplot = 2,
    clim = (0, 50_000), right_margin = 2.0 * Plots.mm,
    title = "B", titleloc = :left)

abuns = @load "examples/Biodiversity/Africa_run100.jld2" abun
meta = Metacommunity(abuns)
div = norm_sub_alpha(meta, 0)
sumabuns = reshape(div[!, :diversity], 100, 100)
heatmap!(sumabuns,
    background_color = :lightblue,
    background_color_outside=:white,
    grid = false, color = :algae,
    aspect_ratio = 1, subplot = 3,
    clim = (0, 50_000), right_margin = 2.0 * Plots.mm,
    title = "C", titleloc = :left)


abuns = @load "examples/Biodiversity/Africa_run50.jld2" abun
meta = Metacommunity(abuns)
div = norm_sub_rho(meta, 1.0)
sumabuns = reshape(div[!, :diversity], 100, 100)
heatmap!(sumabuns,
    background_color = :lightblue,
    background_color_outside=:white,
    grid = false, color = :algae,
    aspect_ratio = 1, subplot = 4,
     right_margin = 2.0 * Plots.mm,
    title = "D", titleloc = :left, clim = (0, 1))

Figure 2: 100 year simulations of Africa with 50,000 species. (A) Species richness after 100 years of simulation with all species equal. (B) Species richness after 50 years, with one specialist introduced. (C) Species richness after 100 years, with one specialist introduced. (D) Representativeness after 50 years with one specialist introduced (0 is completely unrepresentative of the ecosystem as a whole, 1 is completely representative).