A new method to infer vegetation boundary movement from 'snapshot' data

Global change may induce shifts in plant community distributions at multiple spatial scales. At the ecosystem scale, such shifts may result in movement of ecotones or vegetation boundaries. Most indicators for ecosystem change require timeseries data, but here a new method is proposed enabling inference of vegetation boundary movement from one ‘snapshot’ (e.g. an aerial photograph or satellite image) in time. The method compares the average spatial position of frontrunners 
of both communities along the vegetation boundary. Mathematical analyses and simulation modeling show that the average frontrunner position of retreating communities is always farther away from a so-called optimal vegetation boundary as compared to that of the expanding community. This feature does not depend on assumptions about plant dispersal or competition characteristics. The method is tested with snapshot data of a northern hardwood-boreal forest mountain ecotone in Vermont, a forest-mire ecotone in New Zealand and a subalpine treeline-tundra ecotone in Montana. The 
direction of vegetation boundary movement is accurately predicted for these case studies, but we also discuss potential caveats. With the availability of snapshot data rapidly increasing, the method may provide an easy tool to assess vegetation boundary movement and hence ecosystem responses to changing environmental conditions.

A new method to infer vegetation boundary movement from 'snapshot' data Maarten B. Eppinga, Carolyn A. Pucko, Mara Baudena, Brian Beckage and Jane Molofsky M. B. Eppinga (m.b.eppinga@uu.nl) and M. Baudena,Dept of Environmental Science,Copernicus Inst. of Sustainable Development,Utrecht Univ.,the Netherlands. MBE also at: Dept of Plant Biology,Univ. of Vermont,63 Carrigan Drive,Burlington,VT 05405,B. Beckage and J. Molofsky,Dept of Plant Biology,Univ. of Vermont,63 Carrigan Drive,Burlington,VT 05405,USA . Global change may induce shifts in plant community distributions at multiple spatial scales. At the ecosystem scale, such shifts may result in movement of ecotones or vegetation boundaries. Most indicators for ecosystem change require timeseries data, but here a new method is proposed enabling inference of vegetation boundary movement from one 'snapshot' (e.g. an aerial photograph or satellite image) in time. The method compares the average spatial position of frontrunners of both communities along the vegetation boundary. Mathematical analyses and simulation modeling show that the average frontrunner position of retreating communities is always farther away from a so-called optimal vegetation boundary as compared to that of the expanding community. This feature does not depend on assumptions about plant dispersal or competition characteristics. The method is tested with snapshot data of a northern hardwood-boreal forest mountain ecotone in Vermont, a forest-mire ecotone in New Zealand and a subalpine treeline-tundra ecotone in Montana. The direction of vegetation boundary movement is accurately predicted for these case studies, but we also discuss potential caveats. With the availability of snapshot data rapidly increasing, the method may provide an easy tool to assess vegetation boundary movement and hence ecosystem responses to changing environmental conditions. A growing body of research reports changes in species distributions, attributed to recent global change (Araújo andNew 2007, Thuiller et al. 2008). On the global and regional scale, there is a general trend of species' ranges moving toward the poles and to higher elevations (Hickling et al. 2006). Also, at the ecosystem scale global change may induce changes in species distributions (Brown et al. 1997). These latter changes are clearly visible in ecosystems that contain spatial vegetation boundaries, because boundaries may shift if the changed conditions favor one community over the other (Eppinga et al. 2009a). Examples include moving ecotones (Allen and Breshears 1998) and spatial spread of better-adapted invasive species (Kriticos et al. 2003).
The most straightforward method to monitor ecosystemscale vegetation boundary movement would be the analysis of a timeseries of images (from Earth Observation data, i.e. satellite-derived or aerial photographs) (Cohen and Lara 2003). Application of this method, however, requires the availability of data from the past, which becomes increasingly scarce as the length of the time series increases or for more remote regions of the globe. Another complication is that seasonal or meteorological variation between images may lead to an inconsistent classification of vegetation (Adams et al. 1995, Weiers et al. 2004, Barbier et al. 2006. This problem becomes more prominent when the number of images in the timeseries increases (Adams et al. 1995, Weiers et al. 2004, Barbier et al. 2006. Therefore, other approaches are needed as well. Some more recent approaches aim to predict future boundary movement by analyzing current properties of the vegetation boundary. These approaches rely on techniques developed in statistical physics. Percolation theory (Bak 1996, Pascual andGuichard 2005), for example, enables quantification of competitive pressure in an area (Milne et al. 1996, Gastner et al. 2009). For application of this method, however, it is necessary to make assumptions about the size of individual interaction neighborhoods for the plant communities involved in order to predict whether communities are likely to decline in cover or be able to expand into new areas (Milne et al. 1996). Another method has shown that the advancement of an expanding plant community front can be accurately predicted by using the physical analogy of a substance propagating into an unstable medium (O'Malley et al. 2009a, b). Although this method does not require specific assumptions about the size of individual interaction neighborhoods (O'Malley et al. 2009b), it focuses on describing the properties of an expanding community front, rather than identifying differences in spatial structure between the retreating front and the expanding front. In this study, we aim to develop an approach that: 1) enables predictions of future vegetation boundary movement from a single snapshot in time (Milne et al. 1996) and 2) does not require specific assumptions about the sizes of plant communities' interaction neighborhoods (O'Malley et al. 2009b), which are often unknown and difficult to quantify (Hastings et al. 2005). In addition, the development of the method includes a statistical testing procedure to formalize predictions on vegetation boundary movement based on empirical snapshot data.
We develop the method by analytical analyses and by analyzing simulations using spatially explicit competition models. A number of spatially explicit competition models have been developed to mimic the formation and movement of vegetation boundaries (Molofsky et al. 2001, Eppinga et al. 2006, Eppstein and Molofsky 2007. Testing the predictions of competition models with empirical data, however, is difficult and not often performed (Eppstein andMolofsky 2007, Tilman 2007). One potential problem is that the formation and movement of vegetation boundaries may be driven by edaphic factors (Crawford andGosz 1982, Wiens et al. 1985) or by scale-dependent spatial feedbacks between vegetation and its abiotic environment (Rietkerk et al. 2004) rather than competitive interactions. However, in a variety of ecosystems, in particular those under relatively benign conditions, the formation and movement of vegetation boundaries may be primarily driven by competitive interactions (Clarke and Hannon 1971, Snow and Vince 1984, Streever and Genders 1997, Peltzer 2001. In this study we focus on this latter type of vegetation boundary. Although not as restrictive as other methods, our use of snapshot data to predict future boundary movement requires a number of assumptions too. It is assumed that the current vegetation distribution is reflecting competitive interactions between communities over a longer time period (i.e. decades). Also, it is assumed that vegetation boundary movement is relatively slow as compared to fluctuations in environmental and meteorological conditions.
To meet these assumptions when applying the method to real ecosystems, we select snapshots of vegetation boundaries that: 1) consist of two discrete communities that are characterized by plants with a relatively long lifespan (mostly trees and shrubs); 2) have been described in the literature, meaning that the direction of vegetation boundary movement is known; 3) are considered to be primarily driven by competitive interactions. First, we consider the upward movement of the northern hardwood-boreal forest ecotone in Vermont (Beckage et al. 2008), where successful tree establishment is limited by competition (Kupfer and Cairns 1996, Pucko et al. unpubl.). Second, we consider a stable forest-mire ecotone in New Zealand (investigated in Agnew et al. 1993, using a snapshot from Google Earth). Both types of plant community are able to modify their habitat in a way that limits growth of the other community (Eppinga et al. 2009a). Third, we consider the advance of trees in a subalpine treeline-tundra ecotone in Montana (Zeng and Malanson 2006), where tree establishment in the tundra is suppressed by the presence of vascular plants (Moir et al. 1999).

Description of the method
The method is explained in three subsections. We first analyze vegetation boundary movement analytically, yielding the basic premise of our method that does not depend on specific competition or dispersal characteristics. In this first subsection, we consider a relatively simple model for vegetation spread, namely unidirectional spread in the Eden model (Barabási and Stanley 1995), which enables analytical treatment. Although there is a vast literature on the spatio-temporal dynamics of advancing fronts (Durrett and Levin 1994, Barabási and Stanley 1995, Van Saarloos 2003, these studies tend to focus on the movement rate of the advancing community, whereas the movement rates of both the advancing and retreating fronts are crucial for the basic premise of our method.
In the second subsection, we will relax some of the simplifying assumptions that are made in the unidirectional Eden model, and test whether the basic premise of our method still holds. In this section, a more complicated model is used that can no longer be solved analytically. Hence we resort to numerical analyses in this subsection.
Finally, the third subsection describes the statistical procedures that can be used to test hypotheses on vegetation boundary movement with Earth Observation derived vegetation boundary maps.

Analytical analysis of vegetation boundary movement
In the following we consider vegetation boundaries in model lattices, in which each cell can be occupied by one individual of a particular plant community. A spatial vegetation boundary is the border between two areas that are each dominated by a different plant community. Within the vegetation boundary zone, individuals of both communities occur. If the two communities are entirely separated spatially, each community consists of a single homogeneous patch (Fig. 1a). This particular spatial configuration maximizes the conditional probability that a random site within the interaction neighborhood of an individual of community i is occupied by another individual of community i. In other words, the perceived density of a community (calculated as the mean frequency of conspecifics occurring within the interaction neighborhood, following Milne et al. 1996) is maximized in this way. We therefore use the term 'optimal boundary position' for the position where the boundary between vegetation communities occurs for this spatial configuration (Fig. 1a).
Thus, for any spatial configuration, the optimal boundary position can be calculated based on the densities of both communities in the entire lattice: if both communities occupy 50% of the cells (as in Fig. 1b), the optimal vegetation boundary occurs exactly in the middle of the lattice when all individuals of one community would be moved to one side of the lattice. Thus, Fig. 1a shows the result of this ordering process for the spatial configuration shown in Fig. 1b, yielding the optimal boundary position for this spatial configuration. A more formal definition of the optimal vegetation boundary position is provided in Supplementary material Appendix 1.
In practice, communities are usually not entirely separated, because patches of both communities are filled with 'islands' of individuals (i.e. cells) of the competitor community (Fig. 1b). This means that there is a vegetation boundary zone in which individuals of both communities are present (exemplified for a single row of a vegetation boundary in Fig. 1c). As a result, the mean position of each community front (calculated as the mean of the positions of the frontrunner for each row, following O' Malley et al. 2009a, b) lies beyond the optimal boundary position (Fig. 1b).
We now examine how the distance between each community front and the optimal boundary develops over time. We assume that community C 1 is taking over community C 2 . That is, community C 1 expands from left to right, creating a vegetation boundary zone (Fig. 1c). In the following, we will refer to the horizontal distance of the vegetation boundary zone as the width of the boundary, and to the vertical distance as its length. To enable analytical treatment of the model, we first consider a simplified version of the well-known contact process model (Durrett and Levin 1994), namely unidirectional spread in the Eden model (Barabási and Stanley 1995). More specifically, for n cells surrounding a C 1 individual there is a probability p per timestep (which can vary between 0 and 1) that a C 2 individual gets replaced by offspring of this C 1 individual. The only transition that can thus occur in the model is C 2 → C 1 , meaning that community C 1 can expand into the C 2 community as if it were open space. Further, it is important to note that we do not consider the influence of propagule pressure. More specifically, the colonization probability of a C 2 cell does not increase when multiple C 1 individuals establish within the neighborhood; the probability of being replaced by a C 1 individual is either 0 (no C 1 individuals in the interaction neighborhood) or p (the number of C 1 individuals in the interaction neighborhood is  1). This simplification enables us to focus on frontrunner dynamics when calculating community expansion. To further enable analytical treatment (i.e. remove spatial correlation effects between rows, Snyder and Nisbet 2000), we consider the colonization process only to occur in the horizontal direction. These assumptions oversimplifying competition will be relaxed in the next subsection (Numerical simulations of vegetation boundary movement).
Thus, each timestep there is a probability p that a new frontrunner colonizes the cell that is n cells in front of the current frontrunner. Averaged over the entire length of the vegetation boundary, this means that these events will cause an average advance per timestep of the C 1 community front of pn cells. For each row of the vegetation boundary, we know that if a colonization of n cells ahead does not occur (probability 1 2 p), there is still a probability p that the cell that is n 2 1 cells ahead of the frontrunner gets colonized. The average advance per timestep of the community fronts due to these events is then (1 2 p) p (n 2 1) cells. Summation of all possible expansion events in this manner yields the average advance per timestep of the C 1 front: In which A 1 indicates the advance rate of the C 1 community front. Note that the text above described the cases for j  0 (advance of n cells) and j  1 (advance of n 2 1 cells) as expressed in Eq. 1. The expansion process always leads to  The actual position of a front is determined by the mean position of the frontrunners. Note that the actual front position is always farther into the competitor community than the optimal boundary position. (c) Processes governing vegetation boundary movement. The dotted cells indicate sites that are colonized by the expanding community (in black, referred to as C 1 in the main text) in the current timestep. On the right side of the boundary, new frontrunners establish due to colonization, which moves the expanding community front to the right at a rate A 1 . Another effect of this movement is that islands of the retreating community (in white, referred to as C 2 in the main text) are becoming part of the boundary zone, at a rate F 2 . On the left side of the boundary, replacement of the frontrunners of the retreating community moves the position of this community front to the right as well, at a rate R 2 . The width of the vegetation boundary zone is comprised by the total number of C 1 individuals in the boundary zone (N 1 ) and the number of C 2 individuals within this zone (N 2 ).
in which N t 2 is the distance of the C 1 community front to the optimal boundary at time t. Figure 2a illustrates that the analytical description of C 2 island formation (Eq. 2, 3) corresponds well with the results of numerical simulations of the unidirectional Eden model. The simulations show that the C 2 community develops toward an equilibrium number of C 2 islands per row, N 2 . As noted above, N 2 is the number of islands occurring in the C 1 community (Fig. 1b) and this is thus the equilibrium distance of the C 1 community to the optimal vegetation boundary. N 2 is obtained by solving Eq. 3 to equilibrium (meaning N 2 t 1 1  N t 2 ), yielding: Thus, in the above, we have derived the development and final equilibrium of the distance to the optimal boundary for the expanding C 1 community (Eq. 2, 3; Fig. 2a). The two community fronts would have an equal distance to the optimal boundary if both had the same number of individuals within the vegetation boundary zone. We know the number of C 2 individuals within the boundary zone at equilibrium (N 2 as defined in Eq. 4), and this number should be equaled by the number of C 1 individuals within the boundary zone. This would yield a total width of the boundary zone of 2 2 2 .
We will now show that a vegetation boundary of this particular width can never be stable, independent of the size of the interaction neighborhood. First, we calculate the density of C 2 individuals when they become part of the boundary, namely at the C 1 community front (on the right side in Fig. 1c). This density is given by the number of C 2 islands formed, F 2 , divided by the total number of both C 1 and C 2 individuals that become part of the vegetation boundary, which is governed by A 1 : In which e 0 2 indicates the density of C 2 individuals at time t  0, specifying one particular point in time at which a set of C 2 islands is formed, which we will refer to as a cohort of C 2 islands. Over time, the density of this C 2 cohort will decline, because more and more C 2 individuals become replaced by C 1 individuals. Because each timestep a fraction p of the remainder of the cohort will be disappearing, the density of this cohort after m timesteps is described by: Note that the dynamics of each cohort are the same, and that the C 1 community front advances at a constant rate (Eq. 1). Therefore, we can infer the time since establishment of each cohort from its distance to the current position of the C 1 community front: the further away (i.e. to the left in Fig. 1c) a C 2 individual is from the C 1 community front, the longer it has been part of the vegetation boundary. Thus, we can define an average distance x between a cohort formed an equilibrium width of the vegetation boundary zone. This means that the advance of the C 1 community front on one side is compensated by an equally large retreat of the C 2 front on the other side (retreat of the C 2 front is indicated with R 2 in Fig. 1c). C 2 individuals become part of the vegetation boundary zone at the C 1 community front. When n  1, individuals of C 1 can establish more than one cell ahead of the current frontrunner position, which leads to 'islands' of C 2 individuals that are then part of the boundary zone (Fig. 1c). Note that an island refers to a single cell in the lattice. We introduce F 2 to indicate the rate of island formation, meaning the number of C 2 individuals per timestep that become part of the vegetation boundary. Note that if there was only one colonization event per row per timestep, F 2 would simply be given by A 1 2 1. Thus, we can use Eq. 1 as a starting point, but need to correct for the possibility of multiple colonization events per row per timestep. With an increasing number of C 1 colonization events, the number of C 2 islands formed decreases. This yields the following equation: In Eq. 2, j indicates the number of cells that the new frontrunner advances beyond the current C 1 front position (note that the colonizer needs to advance at least two cells to create a C 2 island in between), and k is the number of colonization events per row. The first two factors indicate the probability of k cells getting colonized in the current timestep. The third factor indicates how many islands will be formed due to these events. However, k colonization events within j cells may be achieved via multiple different spatial configurations. This is accounted for by the fourth factor in Eq. 2, which is the multiplicity factor from a twostate model as used in statistical physics (Schroeder 2000), but correcting for the fact that the frontrunner advancing j cells is limiting the degrees of freedom. In Eq. 2, correction for multiple colonization effects per row is achieved by considering k  1, and by the third factor in Eq. 2, which describes the number of islands that are formed. A more detailed explanation of the derivation of Eq. 2 is given in Supplementary material Appendix 2. We introduce N 2 to indicate the average number of C 2 islands per row of the vegetation boundary, meaning that these (and only these) islands determine the distance of the competing C 1 community to the optimal vegetation boundary. The dynamics of N 2 are governed by the formation of C 2 islands and by their disappearance. The formation of C 2 islands is given by F 2 as derived above. Because every C 2 individual within the vegetation boundary has a probability p per timestep to disappear, the number of C 2 individuals that disappears each timestep is pN 2 . Both processes can thus be written as: Equation 8 shows the decline in density of C 2 individuals along the width of the vegetation boundary. We can now insert the width of the vegetation boundary for the case that both communities would have an equal distance to the optimal boundary, and calculate the resulting equilibrium C 2 density, D 2 , herein: Note all terms in Eq. 9 have been expressed analytically as functions of the colonization probability p and size of the interaction neighborhood n in previous Eq. 1-8. Figure 2b illustrates that D 2 within the vegetation boundary zone can never (i.e. not for any combination of n  1 and p) equal or exceed 1/2. This means that the distance to the optimal boundary for the retreating community will always become larger than that of the expanding community. Thus, this analysis reveals that the expansion process of a plant community is reflected by a shorter distance to the optimal boundary as compared to the retreating plant community, independent of the assumed size of the interaction neighborhood (n) or probability of establishment (p). We will now examine whether this observation holds if some of the assumptions used above are relaxed. Therefore, we will consider a model system in which: 1) both communities can reproduce and thus replace each other 2) both horizontal and vertical spread are included 3) the types of interaction neighborhoods are more realistic and 4) the types and sizes of interaction neighborhood can be set different for the two communities. As noted in the beginning of the section, this more complicated model system can no longer be solved analytically. We will therefore continue our theoretical analyses using numerical simulations for a range of parameter values and different functional model forms.

Numerical simulations of vegetation boundary movement
The interaction neighborhood of plant communities depends on competition and dispersal processes, which are often difficult to quantify. For example, deriving the dispersal neighborhood to parameterize mathematical models has proven to be a major challenge in modeling invasions of exotic species (Hastings et al. 2005). In the following, however, we will show that also in a more complicated competition model the key result derived above holds regardless of the size of the interaction neighborhood.
In all simulations, the invading front of a competitively stronger community propagates into the area occupied by a If it is assumed that both community fronts have an equal distance to the optimal boundary, then stability of this spatial configuration requires that the density of the retreating community, D 2 , equals 1/2. However, plotting Eq. 9 presented in the main text reveals that D 2 will always be smaller than 1/2 in this situation, meaning that the vegetation boundary will always develop into a situation where the retreating community front is farther away from the optimal boundary than the expanding community front.
m timesteps ago and the current C 1 community front position: x mA  1 Inserting Eq. 7 into Eq. 6, we obtain the spatial proportional density of C 2 as a function of distance to the current C 1 community front position: In which g ij is the competition coefficient determining the negative effect of community j on community i, d j the mortality rate of community j and K i is the carrying capacity of community i. In Eq. 12, the factor G t i determines the growth rate of community i, and it is described by (Eppstein and Molofsky 2007): In which F t i is the relative frequency of community i at timestep t, i.e. the density of the community that is being experienced relative to its carrying capacity. Because the effect of one individual on its surroundings may depend on community type, community densities are weighed by means of the competition coefficients. Thus, we can write for , as used in Eq. 13. Since each cell can only be occupied by a single individual, overshoot needs to be prevented in the model by making the growth rate of a community dependent on interspecific interactions. More specifically, growth of a community is negatively affected by increasing growth potential of competitor communities (Eppstein and Molofsky 2007). In Eq. 13, H t i is the quality of a habitat for colonization by community i. For brevity, we refer to H t i as 'fitness', defined as: In which b i is a density-independent fitness component and a i,j is a coefficient representing (positive or negative) species interaction feedback, which relates the density of community j to the fitness of plant community i (note that the summation includes intraspecific feedback, i  j). Combining Eq. 12, 13 and 14 then leads to the following discrete time integro-difference equation: Where S is the total number of vegetation communities in the system (here we only consider S  2), b is a parameter determining the intensity of disturbance, creating a fraction of cells in the lattice devoid of vegetation in the next timestep, r i is the potential density that community i could reach the next timestep in the absence of disturbance and dynamics of other plant communities, and r MAX is the maximum density of individuals that can be achieved in the lattice in the absence of any kind of stress or disturbance (here r MAX ≡ 1). It is important to note that the perceived density of community i depends on the size of the interaction neighborhood of this community. This is indicated by the double subscripts in Eq. 10, with r t i,i representing the density of plant community i that occurs within the inter action neighborhood of plant community i (i.e. the density of community i as perceived by community i). Different sizes and types of interaction neighborhoods can be considered. These neighborhood-dependent community densities can thus be turned into a spatially explicit equation by defining the following convolution integral: Where r t i,j (x) is the density of community j as perceived by community i at location x at time t, W is the model domain, k i (x 2 x′) is the interaction neighborhood of community i, and f t j is the frequency of community j in the lattice, which indicates presence (1) or absence (0) of individuals of community j for each cell of the lattice. In this study, we define k using the negative exponential, Gaussian and fattailed kernels (Supplementary material Appendix 3), and we also consider four and eight cell interaction neighborhoods.
In Eq. 10, the potential density is described by a logistic growth function that includes intra-and interspecific competition (Eppstein and Molofsky 2007): The non-spatial parameters in most analyses are set to mimic the standard spatially explicit two-species competition model, of which the spatial dynamics are well known (Bolker et al. 2003). Our point here, however, is to see how these well known spatial dynamics are reflected in the distance of the two vegetation communities to the optimal boundary over time. Starting a simulation with both communities entirely separated, there are three possible outcomes (Bolker et al. 2003).
propagates most rapidly (Fig. 3b), meaning that this community increases in cover (and also the vegetation boundary moves in the same direction as the expanding community front). The weaker community, however, does not go extinct, meaning that the communities mix throughout the lattice. Due to this mixing effect, the vegetation boundary will eventually disappear and both communities will be present in the entire lattice (Fig. 3a). Persistence of the weakest community might be inferred from the density across the lattice, because the community densities are relatively constant along the vegetation boundary zone (Fig. 3c). Mixing occurs when both communities experience similar interspecific competition that is lower than the experienced intraspecific competition. Second, the vegetation boundary may be stable, meaning that neither community is expanding into the other community's area (Fig. 3d). In the absence of an environmental gradient, this situation occurs in the model when both communities experience similarly stronger interspecific competition than intraspecific competition. This lack of boundary movement is reflected by the distance to the optimal boundary not being significantly different between communities ( Fig. 3e) and the density profile showing a consistent trend along the vegetation boundary (Fig. 3f ).
Third, one community can expand if it experiences considerably less interspecific competition than the other community (Fig. 3g). This type of boundary movement is reflected by the expanding community consistently having its front closest to the optimal boundary (Fig. 3h). Together with a consistent trend of the density profile along the vegetation boundary (Fig. 3h, i), these observations reflect competitive exclusion of the weaker community. Thus, the simulations show that movement of the vegetation boundary can be detected at any particular moment in time, because the front of the expanding community is significantly closer to the optimal boundary at any moment than the front of the contracting community (Fig. 3h).
The analytical approach in the previous section considered equilibrium distances to the optimal boundary. An important observation in the model simulations is that the retreating community front is also farther away from the optimal boundary when community fronts are still far away of frontrunners of both community types, and considering the mean difference in relative distance to the optimal boundary: More specifically, we used a bootstrap technique (Efron and Tibshirani 1993) to estimate d and its standard deviation. We generated 100 000 bootstrap replicates to estimate the mean of d and its standard error (which is given by the standard deviation of d within the bootstrap replicates, Efron and Tibshirani 1993), using a random permutation function as implemented in Matlab (ver. 7.13, Mathworks 2011;Eppinga et al. 2010).
Whether the mean value of d differed significantly from 0 was determined by the fraction of the bootstrap samples that had an index value of 0 or the opposite sign of d (Fox 2008, Carvalho et al. 2010. For example, if a mean d value was calculated to be 11.0, and only two percent of all bootstrap replicates had an index value that was 0 or negative, the conclusion was that d was significantly positive, with a p-value of 0.02. This procedure can be written in equation form as: In which H() again depicts the Heaviside function, d i the index value of the ith bootstrap replicate, d the mean index value of all bootstrap replicates and n the total number of bootstrap replicates. Based on the assumptions of a onesample t-test, the variance in x → T,1 2 x → T,2 , together with the sample size (i.e. the length of the vegetation boundary in the image) can also be used to calculate the power that a particular image provides. More specifically, for the case in which x → T,1 2 x → T,2  0, power is calculated as: In which erf () is the error function, t crit defines the desired significance of the t-test, t is the effect size (the difference between distances of the community fronts to the optimal boundary), and s 2 the variance. It is important to note that this analysis comprises a retrospective power analysis. This retrospective power analysis can only provide useful additional information to the t-test if the power analysis is done for an a priori determined effect size (Thomas 1997).
Together with the observed variance in an image, the power can then be calculated. Based on the model simulations in this study, we used an effect size of x → T,1 2 x → T,2  0.2, which was approximately the smallest difference between an expanding and a retreating front as observed in the model simulations. Further, we set as a criterion that the 95% confidence intervals around the front positions of the two communities should not overlap, meaning from reaching their equilibrium distance (Fig. 3). Almost immediately after model initialization, the clear pattern of community front distances emerges (Fig. 3). This suggests that the method can also be used in situations in which the vegetation boundary width is still in a transient state. In real ecosystems this situation is likely to occur.
We further verified that these simulation results were consistent for communities with varying levels of competitive strength (determined by the g coefficients), for different functional model forms in which communities also have complex feedback interactions (determined by the a coefficients), and for communities with varying sizes of commonly used interaction neighborhoods (determined by k) (Supplementary material Appendix 3). Although competitive strength, feedbacks and the sizes of the interaction neighborhood of both communities all influence rates of front advance and retreat (Supplementary material Appendix 3), the front of the expanding community is closer to the optimal boundary than the retreating front at any particular moment in time in all situations (Supplementary material Appendix 3).

Statistical analyses
The spatial positions of community frontrunners not only yield a mean community front position, but also its standard deviation. This enables statistical testing of whether the two community fronts differ in their distance to the optimal boundary. Also, it is possible to calculate the statistical power that a particular snapshot provides for this test, which we will now explain.
Along a vegetation boundary with a length of y individuals, the position of the y frontrunners of community i is determined as follows: In which x → f,i indicates the spatial positions of frontrunners, x → O the optimal boundary position (which is the same for all y elements in the vector x → O ) and x → i encompasses the distances of the frontrunners of community i to the optimal boundary. The frequency distribution of x → i is typically positively skewed, but the distribution becomes close to normal when log-transformed (O'Malley et al. 2009b). Here, we use log(x 1 1) transformation, and ensure that negative positions (i.e. frontrunners that are behind the optimal vegetation boundary) have the same weight as positive positions by using the Heaviside function. More specifically, the distances of frontrunners to the optimal boundary are transformed as follows: in which x → T,i represents the vector of transformed frontrunner positions of community i and H() is the Heaviside function. As noted above, frontrunner positions, and hence distances to the optimal boundary within a community front are correlated, and can thus not be considered to be independent observations. This hampers direct employment of conventional statistical tests. To address this problem, we used a randomization technique considering pairs highly aggregated (Fig. 5b). This was similar to the model simulations of stable boundaries, which also showed the relatively low degree of mixing as compared to moving or mixing vegetation boundaries (Fig. 3). Indeed, analysis of the community fronts showed that both communities were close to the optimal boundary, without a significant difference between the fronts (Fig. 5b). Examination of the density profile (Supplementary material Appendix 5) also showed a pattern consistent with that of a stable boundary as observed in model simulations (Fig. 3) and the power analysis revealed that the number of frontrunner observations was sufficient (Table 1). These results were in compliance with previous research describing this vegetation boundary as being stable (Agnew et al. 1993). In this particular case, the boundary is possibly maintained either by fire (Mark and Smith 1975), or by a particular species of the mire community, manuka tree Leptospermum scoparium, a mire species that grows vigorously close to the ecotone, especially after long periods without fire (Agnew et al. 1993). Some scattered stands of Leptospermum also seem to occur on the open part of the mire (Fig. 5a), as also described for the pakihi studied by Agnew et al. (1993). The occurrence of Leptospermum may hamper tree expansion onto the mire (Agnew et al. 1993), a shift that is frequently observed in European boreal mires, probably due to recent climate change (Eppinga et al. 2009a).

Example 3: tree advancement in a subalpine treeline-tundra ecotone (Montana, USA)
Historical records suggest that the subalpine treeline in Glacier National Park (Montana) has advanced into adjacent tundra at higher elevations during the 19th century (Bekker 2005). Zeng and Malanson (2006) provide an image of the leading edge of a treeline that expanded more recently (Fig. 6a). By analyzing this image, however, we reached an opposite conclusion: the treeline front was much farther from the optimal boundary than the tundra front (Fig. 6a). This suggested that the tundra front is expanding into the subalpine tree community (Fig. 6a). A likely cause for the discrepancy between the observed result and our prediction was that only the leading edge of the vegetation boundary was included in the image. The ecotone between the subalpine forest and tundra communities extends further into the forest than was captured in the image with tundra patches present at lower elevations. These tundra patches may also be changing composition toward that of subalpine meadows or, eventually, to subalpine forest (G. P. Malanson pers. comm.). Analyses of model simulations that considered the entire vegetation boundary (Zeng and Malanson 2006) yielded results that were consistent with the subalpine treeline advancing into the tundra community (Fig. 6b, c). However, when a snapshot captured only the leading edge of this vegetation boundary region from these model simulations, the result incorrectly suggested the advance of tundra (Fig. 6d). Hence, the wrong prediction made from the real system image may have been related to an incomplete representation of the vegetation boundary. t crit  2.807. This power test is of particular relevance to images in which no significant difference in distance to the optimal boundary is found between two communities. In this case, the power test can be used to test whether the lack of a significant difference could be due to a limited number of observations (i.e. length of the vegetation boundary in the image).

Application to Earth Observation derived boundary maps
Example 1: upward movement of a northern hardwood-boreal forest ecotone (Vermont, USA) The ecotone on the Vermont mountains separates the northern hardwood community (sugar maple Acer saccharum, American beech Fagus grandifolia, yellow birch Betula alleghaniensis) present on lower elevations from the boreal forest community (red spruce Picea rubens, balsam fir Abies balsamea, montane paper birch Betula papyrifera var. cordifolia) found on higher elevations (Beckage et al. 2008). Since 1962, this ecotone has shifted ~ 100 m upward, coinciding with a regional temperature increase of 1.1°C (Beckage et al. 2008). In other words, the northern hardwood community has been expanding at the expense of the boreal forest community. We analyzed a satellite image from Camels Hump ( Fig. 4; see Supplementary material Appendix 4 for image processing details), one of the mountains that was studied by Beckage et al. (2008), by selecting six rectangular frames oriented so that the vegetation boundaries in these frames were (approximately) perpendicular to the mountain slope (Fig. 4). For all frames, the northern hardwood community was significantly closer to the optimal boundary ( Fig. 4b-f, Table 1), which corroborated field observations of this community expanding. The density profiles of the frames (Supplementary material Appendix 5) looked most similar to that of an advancing boundary, meaning that the boreal forest may become excluded entirely at lower elevations. However, the density profiles were not as clear as in model simulations meaning that the possibility of a mixing boundary cannot be excluded yet. This was in compliance with the results presented in Beckage et al., which show a decrease of the boreal forest at lower elevations though complete exclusion has not yet been observed (Beckage et al. 2008).

Example 2: a stable forest-mire ecotone (New Zealand)
The forest-mire ecotone at the edge of so-called pakihis (shrubby mire openings in forested peatlands) in New Zealand has been identified as a stable boundary (Agnew et al. 1993). Because the original study used a transect approach and did not provide an image of the system, we used Google Earth to obtain an image close to the field site described in Agnew et al. (1993) (Fig. 5a, see Supplementary material Appendix 4 for image processing details). From this image, it was evident that both communities were  4), or rock outcrops higher on the mountain. (c-h) Selecting six frames in which the vegetation boundary is approximately perpendicular to the slope, the northern hardwood front is consistently closer to the optimal boundary than the boreal forest front (see Table 1 for test statistics).
This example highlights the importance of the snapshot encompassing a representative part of the vegetation boundary. This can be tested by examining the density of the two communities across the image. If the image contains a large enough part of the vegetation boundary, the densities of the two communities converge to a constant value near the edges of the image, indicating that the outcome of the analysis is not strongly affected by selecting a smaller or larger image (see Supplementary material Appendix 5 for details).

Discussion
The method developed in this study suggests that future spatial vegetation boundary movement can be predicted based on a single snapshot in time. Previously developed methods based on percolation theory require specific assumptions about the interaction neighborhoods of plant species (Milne et al. 1996, Gastner et al. 2009), which may be unknown and difficult to quantify (Hastings et al. 2005). Other approaches focused on properties of the expanding community (O'Malley et al. 2009a, b), rather than differences in spatial structure between expanding and retreating community fronts. Analytical analyses and model simulations confirmed that the method presented here does not rely on specific assumptions about the inter action neighborhood (Eq. 9; Fig. 2b, 3; Supplementary material Appendix 3) or the specific formulation of the interactions between the vegetation communities (Supplementary material Appendix 3).
The models that were used to derive our method assume that the distribution of vegetation communities is only determined by direct and indirect local competitive interactions. Application of the method to Earth Observation derived vegetation boundary maps showed that the method could accurately predict vegetation boundary movement, suggesting that the method can still be used if the model assumptions are somewhat relaxed. We assumed constant parameter values in most model simulations, but even when environmental fluctuations were included the method still accurately described boundary dynamics over longer periods of time (Supplementary material Appendix 3). Further testing will be needed, however, to assess under which conditions spatial vegetation patterns are sufficiently driven by local species interactions to infer future movement based on a snapshot of that pattern. Other factors that could drive the spatial vegetation patterns include edaphic factors (Wiens et al. 1985), feedbacks that change from positive to negative with increasing spatial scale (Eppinga et al. 2009b), herbivory (Silliman et al. 2005) and long-distance dispersal processes (Shigesada et al. 1995) possibly leading to the formation of satellite populations Subalpine treeline-tundra, frame a (X T,1 : tundra, X T,2 : treeline) 3.6 (0.1) 5.4 (0.1)  0.00001 tundra advances Subalpine treeline-tundra, frame b (X T,1 : tundra, X T,2 : treeline) 5.8 (0.01) 5.4 (0.02) 0.00001 treeline advances Subalpine treeline-tundra, frame c (X T,1 : tundra, X T,2 : treeline) 5.9 (0.03) 4.4 (0.08) 0.00001 treeline advances Subalpine treeline-tundra, frame d (X T,1 : tundra, X T,2 : treeline) 5.1 (0.2) 5.3 (0.5)  0.00001 tundra advances * X T,i is the transformed measure for the distance of the community front to the optimal boundary, as explained in the main text. There is no significant difference between the communities in their distance to the optimal boundary position (see Table 1 for test statistics). Only the optimal boundary position is shown, which overlaps with the community front positions. Figure 6. Snapshots (Earth Observation derived boundary maps and models) of a vegetation boundary between the subalpine treeline (black) and tundra (white) communities in Glacier National Park, Montana (USA) (taken from Zeng and Malanson 2006). (a) In contrast to field observations, the snapshot (derived by Zeng and Malanson (2006) from an ADAR image) suggests advance of the tundra into the alpine treeline. (b-c) Snapshots from model simulations by Zeng and Malanson (2006) suggest that the treeline is advancing at the expense of the tundra, which agrees with field observations. (d) Performing the analysis on a subset of the vegetation boundary shown in (b) yields an opposite conclusion on the vegetation boundary movement. Hence, the erroneous prediction from the Earth Observation derived boundary map may be due to an incomplete representation of the vegetation boundary (see Table 1 for test statistics). (Moody and Mack 1988). The basic premise of our method is less likely to hold when one or several of the above factors drive the observed vegetation pattern. Similar to the approach followed in this study, a promising starting point for future studies would be to test the method in systems in which the vegetation boundary movement has been described in previous studies. This approach would provide a rather concrete test of the importance of local species interactions (which are at the core of competition models) with empirical data. Deriving concrete predictions from competition models that can be tested with empirical data is generally difficult to do and not often performed (Eppstein andMolofsky 2007, Tilman 2007). Further testing of the method would benefit from application to timeseries data, meaning that timeseries data would still be needed in the testing stage. An interesting avenue to pursue further in this stage would be to examine whether the method could also be used to assess the velocity of boundary movement in addition to the movement direction. Once the conditions for successful application of the method have been identified, the challenges involved in the generation of timeseries could be circumvented. Under these conditions, the method presented here would provide a relatively easy way to predict future vegetation boundary movement.
Although consistent vegetation classification is less challenging for a single snapshot as compared to a timeseries of images (Adams et al. 1995, Weiers et al. 2004, Barbier et al. 2006, there are also potential caveats in classifying a single snapshot. The vegetation patterns on the Earth Observation derived boundary maps will depend on the pixel size of the image. It is important that the pixel size is smaller than the typical dispersal distance of the two vegetation communities involved, so that community islands in the vegetation boundary zone can be noticed. In this study, pixel size for the forest-mire and alpine treeline-tundra ecotones was approximately 1 by 1 m and for the northern hardwood-boreal forest ecotone 10 by 10 m. These pixel sizes are small when compared to the typical dispersal kernels of tree species (Clark et al. 1999), which dominate the expanding communities in our case studies.
Another point of attention is that misclassification of vegetation communities can still occur in a single image. Remote sensing techniques, however, can yield probabilities of a pixel belonging to a certain class (Lees and Ritman 1991). Rather than assuming a perfect classification when calculating frontrunner positions, classification probabilities could be used to calculate a weighed frontrunner position. For example, the frontrunner of a community in a particular row of the boundary may be much farther away from the optimal vegetation boundary than the second individual of that community in the same row. If the frontrunner pixel, however, has a very low probability of having been classified correctly, the actual frontrunner position in that row would be corrected to a position closer to the optimal boundary. Note that this procedure would be very similar to the averaged probability used in the theoretical section of this paper, which is an appropriate approach provided that the number of observations is sufficient (Fig. 2a).
Our findings could be especially useful for monitoring and predicting ecosystem responses to global change (Supplementary material Appendix 3). Due to global change, many vegetation boundaries are moving or are expected to move during the coming century (Hickling et al. 2006). For example, mountain ecotones are largely climate dependent, meaning that the projected increase in global temperature may change competitive interactions in a way that leads to accelerated upward movement of these ecotones (Beckage et al. 2008). Also, for some invasive species, there is a lag phase between establishment and expansion (Shigesada et al. 1995, Hastings et al. 2005. One possible cause that can trigger expansion is changing environmental conditions that favor the invader (Kriticos et al. 2003). Of particular concern are the projected increases in extreme events and climatic variability (Meehl et al. 2007) that may cause more frequent ecosystem disturbances and trigger invasive species expansion (Seabloom et al. 2003).
The increasing availability of snapshot data, aerial photographs and satellite imagery (e.g. through Google Earth), offers a great potential for monitoring ecosystem responses to global change. This study suggests that it may even be possible to predict future ecosystem changes from a single snapshot in time. It is our opinion that this method deserves further application to Earth Observation derived vegetation boundary maps, as it provides a concrete test for widely used competition models and it may develop into an easy tool for quickly predicting how community distributions are moving in response to changes in environmental conditions.