The base parameter values were set as follows: species pool (i.e. metacommunity) size (J _M) = 1,000,000 individuals; fundamental biodiversity number (θ) = 100; local community size (J ) = 1,000; mortality rate (d ) = 0.1; immigration rate (m ) = 0.1. Based on the base parameter values of this model, the parameter-dependency patterns were studied by changing the four targeted parameters (θ, J , d , and m ). For the fundamental biodiversity number (θ), 100 numbers were randomly selected between 1 and 1,000, without allowing duplicates. For the local community size (J ), 25 numbers were randomly selected for each order of magnitude (10²–10⁵); thus, 100 numbers were selected in total. For both mortality rate (d ) and immigration rate (m ), 100 numbers were determined between 0.01–1.00, in 0.01 increments.

At the beginning of each simulation, local community members were chosen depending on the set number of local community size, and 1,100 time steps, including the initial state, were conducted for each simulation (Fig. 1b). After the simulation, the information of local community dynamics between time step 1 (T1) and time step 99 (T99) were excluded from the following analysis for testing parameter dependency because the randomness of the initial sampling could affect the community dynamics at the beginning phase. The dynamics from T100 to T1100 were targeted for the main analysis, and thus, a 1001 time-series dataset was used for each simulation (Fig. 1b). The Bray–Curtis dissimilarity index (Odum, 1950) and Sørensen dissimilarity index (Sørensen, 1948) were calculated for each simulation.

The preliminary analysis revealed that the temporal distance-decay curve rarely reached one (i.e. completely dissimilar); thus, upper limits exist at less than one in temporal beta diversity under neutral dynamics. Therefore, to describe the form of the simulated temporal distance-decay patterns of both the Bray–Curtis and Sørensen dissimilarity indices, a negative exponential curve was fitted for each dataset (Fig. 1c). The equation for the negative exponential function is as follows:

\(TBD=-\delta\times e^{-\varepsilon\ \times\ td}+\zeta\),

where TBD is the temporal beta diversity value, e is Napier’s constant, td is the temporal distance between a pair of communities, and δ, ε, and ζ are the parameters of the negative exponential function (Table 1). The changes in the curve shape with parameter changes are summarised in Fig. 2. The parameter ζ determines the position of the upper limit of the curve; thus, a larger ζ indicates higher upper limits in the temporal distance-decay of temporal beta diversity (Fig. 2c). Parameter δ determines the relative position of the intercept against parameter ζ (i.e. the upper limit). The ζ – δ value indicates the intercept of the curve; thus, a larger δ indicates a lower position of the intercept (Fig. 2a). Parameter ε determines the curvature of the curve; thus, a larger ε indicates a higher curvature of the curve (Fig. 2b). All estimated parameter values without errors were used in the analysis to check parameter dependency in the present study. Although Martín-Devasa, Martínez‐Santalla, Gómez‐Rodríguez, Crujeiras, & Baselga, (2022) recently developed a method to compare the fit of three different functions on distance-decay curves, judging the best-fitted function for a distance-decay curve was the main use of the method and ecologically interpreting the changes in parameter values was difficult based on their methods. Therefore, the three-parameter negative exponential function was used in this study owing to the differences in purpose.