# A meta-community simulation model
# Base R code. Does not require special packages.
# Simulates the dynamics of multiple species in a patchy landscape
# requires two input parameters: 
# [1] sigma, which is a measures of niche breadth(sigma > 0) (equation 2 in the text).
# [2] z, which is the coefficient of exponential decay with distance in the number of arriving propagules into a given patch (equation 4 in the text; smaller values increase arrival rates from distant patches).
# Returns a data frame which includes simulation paramters, landscape characteristics, and species richness.

model = function(sigma, z)
{
	# set run parameters
	Jtotal = 500	# total number of sites 
	Nspecies = 250	# inital number of species
	Tmax = 1000		# number of time steps
	d = 0.25		# mortality rate (fraction of individuals lost in each time step
	Niter = 100		# number of models iterations per parameter combination (a new landscape and a new community are generated each iteration)
	Erichness = c()	# initialize the vector of landscape scale patch richness
	Srichness = c()	# initialize the vector of species richness at the landscape scale
	results = as.data.frame(array(data = 0, dim = c(Niter, 7)))
	colnames(results) = c("Npatches", "sigma", "Nspecies", "Erichness", "Srichness", "z", "Shannon")
	
	for (i in 1:Niter)
	{
		# Generate the landscape
		Npatches = sample(seq(2, 500), 1)	# pick a random number of patches
		Jpatches = rmultinom(1, Jtotal, rep(1, Npatches))	# number of sites per patch
		Jpatches = refill(Jpatches)	# correct patch sizes in cases of zero sites. Uses helper function 'refill' (code below)
		Epatches = sample(seq(1, 500), Npatches, replace = FALSE)	# Assign a unique E value (environmental condition) for each patch.
		XY = matrix(runif(Npatches * 2, 0, 100), nrow = Npatches, ncol = 2) # generate random patch locations
		Ddist = as.matrix(dist(XY, diag = T, upper = T))	# calculate inter-patch distances
		D = 1 / exp(z * Ddist)	# calculate propagule arrival ratios based on an exponential dispersal kernel. z is the distance-decay coefficient.

		# Calculate heterogeneity metrics
		Erichness = c(Erichness, sum(unique(Epatches) > 0))		# calculate patch richness (heterogeneity measure)
		Shannon = -sum(Jpatches / Jtotal * log(Jpatches / Jtotal))	# calculate Shannon's diversity index (another heterogeneity measure)
		
		# Generate species data
		Especies = sample(Epatches, Nspecies, replace = TRUE)		# Assign a niche center for each species
		abundance = array(data = 0, dim = c(Npatches, Nspecies))	# initialize the species abundance table (patch by species array)
		fitness = array(data = 0.0, dim = c(Npatches, Nspecies))	# initialize the array of habitat suitability per species per patch
		for (p in 1:Npatches)
		{
			fitness[p,] = exp(-(Especies - Epatches[p])^2 / (2 * sigma^2))	# calculate habitat suitability per species per patch, based on manuscript equation [2]
			abundance[p, ] = rmultinom(1, Jpatches[p], fitness[p, ])		# determine initial species abundances in patches
		}
		
		# Simulate meta-community dynamics
		for (t in 1:Tmax)
		{
			abundance0 = abundance # keep the abundance levels from the beginning of the current time step
			for (p in 1:Npatches)
			{
				survivors = floor((1 - d) * abundance[p, ])	# abundance of surviving individuals (post mortality; first right-hand term in manuscript equation [1])
				available = Jpatches[p] - sum(survivors)		# number of available sites
				R = (fitness[p, ] * colSums(D[p, ] * abundance0)) / sum(fitness[p, ] * colSums(D[p, ] * abundance0)) # probability of establishment per species (rightmost term, manuscript equation [1])
				recruits = 0
				if (max(R) > 0){recruits =  rmultinom(1, available, R)}	# calculate the number of establishing individuals per species
				abundance[p, ] = survivors + recruits		# update species abundance distribution per patch
			}
		}
		Srichness = c(Srichness, sum(colSums(abundance) > 0, na.rm = TRUE))	# calculate the number of species at the landscape scale
		results[i, ] = c(Npatches, sigma, Nspecies, Erichness[i], Srichness[i], z, Shannon)
	}
	return(results)
}


# Helper function to correct the formation of patches with size 0 by the multinomial procedure in the main model

refill = function(r)
{
	r = sort(r, decreasing = TRUE)
	l = length(r)
	Nneeded = sum(r==0)
	if ((Nneeded == 0) | (max(r) == 1)) return(r)
	r[1] = r[1] - 1
	r[l] = r[l] + 1
	r = refill(r)
	return(r)
}