Bernard Hugueny, Howard Cornell, and Susan Harrison. 2007. Metacommunity models predict the local–regional species richness relationship in a natural system. Ecology 88:1696–1706.

Appendix B. Stochastic simulations.

Using Monte-Carlo simulations, we allowed species to differ in their colonization probabilities (m_LL,₁ … m_LL,n; or m_MI,₁ … m_MI,_n), extinction probabilities without competition (e_1,1… e_1,n), and competitive abilities. To vary competitive abilities, we introduced asymmetrical competition by assuming

(B.1)

where n is regional richness, E_i,j is the total extinction probability of species i in patch j, w_k,j= 1 if species k is present in patch j and zero if it is not, and c_i,k is the pairwise competition coefficient which indicates the increase in extinction probability of species i due to the presence of species k. If all the c_i,k are equal to c and all the e_1,i equal to e₁ , Eq. B.1 reduces to Eq. 1. We defined sensitivity to competition of species i (sc_i ) as the increase in extinction probability of species i if all the n-1 other species are present, such that E_i,j = e_1,i + sc_i when all w_k,j= 1. If sc_i differs for each species, there is a competitive hierarchy where the species with the lowest sc_i is the best competitor. Each sc_i value was then randomly partitioned into n-1 c_i,k values such that

.

(B.2)

Note that Monte-Carlo simulations use probabilities and not rates, thus constraining the values of e_1,i and m_MI,_i to fall between 0 and 1. For the same reason e_1,i+ sc_i (extinction rate when all the species occur together) should fall between 0 and 1. When the value of e_1,i is known, the value of sc_i cannot exceed 1- e_1,i and hence sc_i /(1- e_1,i ) should fall between 0 and 1. Another constraint is that, for the Levins-like model, e_1,i cannot exceed m_LL,_i if the metapopulation is to persist. Thus, e_1,i /m_LL,_i must also fall between 0 and 1. Values falling within these bounds were generated by using the distribution Y(m,v) defined by Y=exp(Z)/1+exp(Z), with Z being normally distributed with mean = log[m/(1-m)] and standard deviation= v. By choosing m in the (0,1) interval, Y(m,v) generates values bounded by 0 and 1 with the expected median given by m.

For example, the values of e_1,i were generated as follows. First a value x, was drawn at random from a uniform (0,1) distribution to determine the expected median m of the e_1,i. Then n values of e_1,i were drawn at random from a Y (x, v) distribution. The m_MI,_i or e_1,i /m_LL,_i and sc_i /(1- e_1,i ) were generated similarly by drawing values, u and q, from a uniform (0,1) distribution. Note that e_1,i must be drawn first before the values of m_LL,_i and sc_i can be obtained.

Each simulation is thus based on four parameters (x, u, q, v). The three first parameters are used to generate metacommunities that differ in the average features of their constituent species, more precisely in the median of e_1,i (=x), m_MI,_I or e_1,i /m_LL,_i (=u), and sc_i /(1- e_1,i ) (=q). The parameter v determines dispersion around the median and takes the same value for the three distributions. As v increases, extinction and colonization probabilities become more dissimilar between species and competition becomes more asymmetrical.

Each simulation began with 100 patches and an initial number of species n_t=0. Simulations were initiated with each of the n_t=0 species occupying the fraction of patches it would occupy at equilibrium without competition. Simulations were run for 100 time steps. For each time step and patch, the probability of extinction of each species was computed from Eq. B.1. If a species i did not occur in a patch, its colonization probability was computed as m_{MI, i} (mainland-island) or y_i m_LL,_i (Levins-like), where y_i is the fraction of patches occupied by species i at that time step. The fate of a species in a patch in each time step (extinction or colonization) was determined using a standard Monte-Carlo procedure (Hilborn and Mangel 1997). For example, if E_i,j is the extinction probability of species i in patch j at time t then extinction occurs at time t+1 if E_i,j is higher than a value drawn at random within a (0,1) uniform distribution. The environment was assumed to be constant, thus none of the parameters in the simulations changed through time.

For the Levins-like model, the number of species remaining in the system at t=100 time steps (n_t=100) was determined. The simulation was then run for 50 additional time steps and the number of persisting species at t=150 (n_t=150) was determined. If n_t=150< n_t=100, we assumed that the system had not reached quasi-equilibrium. Otherwise, the mean species richness per patch was computed for the last 50 time steps and divided by n_t=100 to obtain y_observed, the mean frequency of patch occupancy. We then averaged e_1,i , m_LL,i, and c_i,k over all species remaining in the simulation at equilibrium. We used these values in Eqs. 7 and 9 to compute y_predicted, the mean frequency of patch occupancy when all species are identical.

For the mainland-island model, the simulation procedure was simplified because there was no regional extinction, thus n_t=150=n_t=100=n_t=0. As before, the mean species richness per patch was computed for the last 50 time steps and divided by n_t=100 to obtain y_observed. We again averaged e_1,i, , m_MI,i, and c_i,k for each simulation and used these values in Eqs. 7 and 9 to compute y_predicted,.

For each model (Levins-like or mainland-island) and for a given value of v, we ran 20 simulations with n_t=100= 3 , 5 with n_t=100= 8, 5 with n_t=100= 10, 5 with n_t=100=12, and 5 with n_t=100= 14. We thus obtained 40 values of y_observed to compare to 40 corresponding values of y_predicted (Fig. 5).