Simulation description
The model was implemented in the rectangular two-dimensional grid with
periodic boundary condition (Fig. 1). In simulations, the per capita
ecological difference among species only contributed to the
probabilities of immigration success; while, the per capita speciation,
birth and death rates among individuals in local community were
equivalent regardless of their identities.
The scenario was specified first by the standard parameters of neutral
model; number of sites, speciation rate ν with point mutation
mode of speciation and niche conservatism, number of individuals in
regional species pool JM , dispersal rate m , and number of
individuals in local community J . Then, a grid of local
communities with coordinate (x , y ) and environmental
gradient E was specified in matrix landscape , and an
initial condition of regional species pool; individuals code, species
code, guilds code, their environmental optimum and tolerance was
specified in matrix pool.t0 (Fig. 2a).
Based on a number of sites, a weights matrix for nearest neighboring
communities with spatial weights 1/8 was generated in a grid onto torus
(matrix nb.mat ). The individuals of regional species pool were
assigned to initial locations at random (matrix LC.t0 ). At each
time step, all individuals were removed and replaced with probabilityν by new species (matrix new.sp ). The matrix of remained
individuals (matrix LC.t0.without.new.sp ) was standardized for
each site by row sum (matrix LC.t0.RA ), and multiplied by
matrix nb.mat to generate a set of probabilities of immigration
success of individuals to each local community under neutral dynamics
(matrix I.RA ). The matrix I.RA was multiplied element
wise by a matrix of habitat associations of individuals to sites (matrixHA.siteBYind ), in that assigned weights depended on the value
of probability density function for normal distribution given the mean
as environmental optimum of individuals and the standard deviation as
tolerance of individuals for the local environment of sites, and
standardized for each site by row maximum. The element wise
multiplication was standardized for each site by row sum to generate a
set of probabilities of immigration success of individuals to each local
community under environmental filtering (matrix EF ). The
parents of each local community were chosen with replacement by a
weighted lottery using a corresponding probability vector in m x
matrix EF + (1-m ) x matrix LC.t0.RA as the sum
of immigrant (left term) and local birth (right term) parents (matrixLC.t1.without.new.sp ). Finally, the matrixLC.t1.without.new.sp was combined with the matrixnew.sp ; then, the extinct lineages of individuals were removed
(matrix LC.t1 ). Also the matrix pool.t0 was adjusted
for regional species pool (matrix pool.t1 ). The ancestors of
new species were identified by the ancestor individuals codes and number
of time steps in matrix LC.t1 , and identified by these codes in
matrix pool.t1 for individuals and species levels; while, the
decendants from a common ancestor individual were identified by same
code, and each lineage was grouped in a single column in matrixLC.t1 and row in matrix pool.t1 . The matrixpool.t1 and matrix LC.t1 were the input data for next
time step (Fig. 2b).
I used a speciation rate ν = 0.001 and a number of individuals in
local community J = 16 individuals. The small system was
simulated in a grid of 20×20 local communities using a JM = 6400
individuals for two lower dispersal rates m = (0.01, 0.09). The
large system was simulated in a grid of 40×40 local communities usingJM = 25600 individuals for the highest dispersal rate m =
0.81. In small system, the environmental gradient E presented
three structures; random, one wave and 16 humps. In large system, the
environmental gradient was quadrupled and presented three structure;
random, two waves and 64 humps. In all cases, the environmental gradient
was bounded between 0 and 1, and repeated continuously across opposite
sides of system (Fig 1). The initial condition of regional species pool
was varied by the number of functional groups g , habitat
associations and population sizes of guilds. In small system, the
environmental optimum of guilds was assigned following a uniform
distribution between 0 and 1, and the tolerance of guilds between 0 and
10. The population sizes of guilds were generated splitting JM atg -1 points. These points were g -1 integers sampled from a
uniform distribution between 1 and JM -1 and sorted by their
values. The simulations in each combination of dispersal rate and
environmental structure were started from the same set of five levels ofg = (1, 8, 40, 160, 500), five habitat associations and five
population sizes of guilds except for the case of g = 1 (i.e. the
simulations were started from 25 habitat associations for which only one
population size of a guild was possible). In large system, the initial
condition of regional species pool was quadrupled. The simulations in
each environmental structure were started from the same set of five
levels of g and five habitat associations of guilds; however,
only one case of population sizes of guilds was explored because the
computation in large system was intensive and time-consuming. So that
825 scenarios, 750 = 2x3x4x5x5 + 2x3x1x25 were simulated in the small
system with two lower dispersal rates and 75 = 1x3x5x5 in the large
system with the highest dispersal rate.
In the neutral model with nearest neighboring communities and ν =
0.001, the approximate scales on which individuals were expected to
diffuse before speciation were three for m = 0.01, nine form = 0.09, and 28 for m = 0.81 (Cencini et al., 2012). The
model was explored using a fixed speciation rate, and a range of system
sizes, dispersal rates, environmental structures and initial conditions
of regional species pool. The model communities were a grid of 10×10
local communities in the center of system, and approximated from an area
encompassing independent biogeographic units to an area packed in a
biogeographic unit with open boundary conditions (i.e. the individuals
diffused across borders of system were expected to originate new species
at the opposite sides of model communities), and presented the three
environmental structures; four humps, linear and random (Fig. 1).
Arguing from the approximate scales, the model communities in large
system with m = 0.81 was packed nine times more compactly within
a biogeographic unit than small system with m = 0.01 in each
environmental structure. The real ecosystem may be huge, and the
extinction may be balanced by minimum speciation rate (Bell, 2003), and
ecologist may study often an area packed much more compactly within a
biogeographic unit. At present; however, the computation in larger
system with lower speciation rate is not possible. The data analyses
were performed on R version 3.6.3 (R Core Team, 2019) using the packages
spdep (Bivand et al., 2019), vegan (Oksanen et al., 2019), PCNM
(Legendre, Borcard, Blanchet, & Dray 2012), randomizr (Coppock, Cooper,
& Fultz 2019), reshape2 (Wickham, 2017), dplyr (Wickham, François,
Henry, Müller, & RStudio 2019), entropart (Marcon & Hérault, 2019),
packfor (Dray, Legendre, & Blanchet 2016), pgirmess (Giraudoux,
Antonietti, Beale, Pleydell, & Treglia 2018), and ape (Paradis et al.,
2019). R sources are available in
https://github.com/takayukiyunoki/spatialIBM.git.