### Conceptual model

We use a simplistic model to elucidate the effects of variability for the types of environmental gradients and species niche dimensions envisioned for ecotones, specifically for hierarchical competition and the SGH. We use an approach that eliminates the need to know the response of a particular system to specific changes in climate. Instead, we consider the plant’s-eye-view by varying habitat quality as a composite parameter and change it through time with and without increasing variability during a period of improving habitat quality.

The model assumes that interactions among plants are based on the relative size of the focal individual and the summed sizes of its neighbors and that this interaction intensity is a multiplier the habitat quality, which is <1 for competitive effects and >1 for facilitative effects; the range is from 0.5 to 2; [39] for details. The model computes the habitat quality given an initial environmental gradient multiplied by the interaction intensity. This habitat quality is then used to determine the probability of reproduction, the amount of growth, and the probability of mortality.

### Model design

We use a spatially explicit agent based model created in Netlogo [40], Additional file 1. The model is based on a grid of 1000 x 50 cells, wrapped as a cylinder to eliminate edge effects. The plant’s-eye-view of the environment is defined by responses of two species (SpS, for specialist and SpG for generalist, subscript s) on the long axis of the grid set as hierarchical Gaussian curves, E_{ys}, 1-0 across the rows (subscript y) as shown in Fig. 1 (cf. [41]). This approach collapses the various dimensions of niche to a single gradient appropriate for a simple exploratory model. The responses in this model are relative but are based on observations of *Abies lasiocarpa* and *Picea engelmannii*, for SpS, and *Pinus albicaulis*, for SpG, at alpine treeline in western North America (see [42]). *P. albicaulis* is a pioneer, keystone, and foundational species at alpine treeline; it is able to reproduce in high-stress environments without neighbor facilitation, and it grows slowly. *A. lasiocarpa* and *P. engelmanni* often establish next to extant *P. albicaulis* and have somewhat higher growth rates. Information on survivorship is sparse because treeline sites are relatively new (since the Little Ice Age in the Rocky Mountains [43]) but the ages of *P. albicaulis* are much greater than any of the other two species at sites where we have examined tree rings (e.g., [44]). Although the gradient is smooth, the feedback that is simulated creates clumped patterns of species occurrences. Feedback is likely to be more important in the spatial heterogeneity of factors such as microclimate or soils than are purely abiotic processes [45].

Each cell of the grid can be occupied by one individual, and to initialize model runs, all cells are occupied by an individual j, at a random size between 1 and 100, with a probability proportional to E_{ys}. The dynamics of the population are simulated over 400 iterations of recruitment, growth, and death as Monte Carlo processes with treatments implemented in the last 100 iterations. Although the initial size distribution is random rather than realistic, age structure develops during the initial 300 iterations before treatments begin.

The stress gradient version is computed as a logarithmic increase in interaction intensity with the number of neighbors: multiplied by a gradient from 0.5 to 2 across the length of the grid:

$$ \mathrm{d} = \mathrm{number}\ \mathrm{of}\ \mathrm{neighbors} $$

$$ {\mathrm{I}}_{\mathrm{stress}} = .26 + .333\ \ln\ \mathrm{d} \ast \left(2\ \hbox{--}\ 1.5{\mathrm{E}}_{\mathrm{y}}\right) $$

The size gradient version is

$$ \mathrm{d} = {\mathrm{S}}_{\mathrm{n}}/{\mathrm{S}}_{\mathrm{f}} $$

$$ {\mathrm{I}}_{\mathrm{size}}\left|\mathrm{d}<33.333 = 1 + .03\mathrm{d}\right. $$

$$ {\mathrm{I}}_{\mathrm{size}}\left|\mathrm{d}>33.333 = 3.25 - .0375\mathrm{d}\right. $$

where S_{n} is the sum of the sizes of the eight neighbors and S_{f} is the size of the focal individual; for empty cells the intensity for recruitment is calculated with S_{f} = 1 interaction intensities less than 0.5 or greater than 2 are reset to these limits (Fig. 1). A logarithmic form was chosen to give most weight to the first neighbor; cf. [29, 36]. While it would seem that the model would be sensitive to relative size at which competition begins to outweigh facilitation (here at 10 times, with competition reaching its limit at 30), trials with these points set at 3/10 and 33.33/100 produce quantitatively similar results. Although the stress gradient is linear with the environment, spatial patterns develop because of the density of neighbors, especially at the treeline.

Recruitment for SpS and SpG occurs on cells with probability a function of the size of their extant population relative to the size of the grid and the environment of the row:

$$ \mathrm{P}\left({\mathrm{R}}_1\right) = \mathrm{r}\mathrm{N}\ \mathrm{I}\ {\mathrm{E}}_{\mathrm{y}1} $$

$$ \mathrm{P}\left({\mathrm{R}}_2\right) = \mathrm{r}\mathrm{N}\ \mathrm{I}\ {\mathrm{E}}_{\mathrm{y}2} $$

where N is the current population and r is 0.00001 or 0.000008 based on the size of the grid, so that the maxima would be with rN = 1.0 or 0.64 for SpS and SpG, respectively. SpG can establish only on empty cells whereas SpS can replace SpG but with its recruitment rate halved. Recruitment parameters are chosen to reduce dimensionality by setting them to allow replacement recruitment in a hypothetical open, ideal environment. Because of the feedbacks, recruitment in the neighborhood of extant individuals is more common in the treeline area.

SpS has a higher mortality rate than SpG to match its higher reproductive rate:

$$ \mathrm{P}\left({\mathrm{M}}_1\right) = 1 - .5\ \mathrm{I}\ {\mathrm{E}}_{\mathrm{y}1} $$

$$ \mathrm{P}\left({\mathrm{M}}_2\right) = 1 - .4\ \mathrm{I}\ {\mathrm{E}}_{\mathrm{y}2} $$

The maxima are set at half of the maxima for E_{ys} (1.0 and 0.8) based on the logic that once established mortality is relatively rare.

### Model runs

To represent climate change, the value of E_{ys} is reset to E_{y-1s} in iterations 301–400. E_{ys} increases; the climate ameliorates in the plant’s-eye-view. Iterations 1–300 with E_{ys} constant allow all relevant variables to equilibrate.

Climatic variability is created by multiplying the values for all cells in every year by a random number drawn from a normal distribution with a mean of 1 and a standard deviation of 0.25. We use 0.25 based on the standard deviation of standardized tree-ring widths recorded for a > 400 year-old stand of alpine larch at treeline in Glacier National Park, which range from 0.22 to 0.25 (personal communication, Greg Pederson); comparison of the two values has revealed much higher extinction rates for 0.25 when the multiplier doubled over 200 year [11]. We increase the multiplier to 1.5 over 100 year, after which it again becomes constant to represent a new climatic equilibrium. With the multiplier E_{ys} is constrained to not exceed 1.0 or 0.8 for SpS and SpG, respectively.

We present the results of six simulations, but we show the differences between scenarios with constant and with increasing climatic variability for simulations with no feedback, with the stress gradient feedback, and for the size gradient feedback. The metrics of analysis here are the differences in the population sizes of SpS and SpG between the climatic representations averaged from 100 replicate runs, and the proportion of those runs that end in extinction for either species (extinction/200); simulation runs in which extinction occurred during the first 20 iterations, prior to equilibration of population sizes, were discarded and replaced.