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.