Disentangling the Formation of Contrasting Tree-Line Physiognomies Combining Model Selection and Bayesian Parameterization for Simulation Models

Alpine tree-line ecotones are characterized by marked changes at small spatial scales that may result in a variety of physiognomies. A set of alternative individual-based models was tested with data from four contrasting Pinus uncinata ecotones in the central Spanish Pyrenees to reveal the minimal subset of processes required for tree-line formation. A Bayesian approach combined with Markov chain Monte Carlo methods was employed to obtain the posterior distribution of model parameters, allowing the use of model selection procedures. The main features of real tree lines emerged only in models considering nonlinear responses in individual rates of growth or mortality with respect to the altitudinal gradient. Variation in tree-line physiognomy reflected mainly changes in the relative importance of these nonlinear responses, while other processes, such as dispersal limitation and facilitation, played a secondary role. Different nonlinear responses also determined the presence or absence of krummholz, in agreement with recent findings highlighting a different response of diffuse and abrupt or krummholz tree lines to climate change. The method presented here can be widely applied in individual-based simulation models and will turn model selection and evaluation in this type of models into a more transparent, effective, and efficient exercise.


Introduction
The key to understanding and prediction in science lies in the elucidation of mechanisms underlying observed patterns (Levin 1992;. For example, a key interest in ecology is the establishment of relationships between processes acting at the level of individuals and patterns emerging at the scale of populations or communities (Levin 1992;. One common approach to solving this problem involves the assessment of the ability of mechanistic models to reproduce observed real-world patterns. Nevertheless, the extent to which emergent ecological patterns can be used to disentangle processes acting at lower scales remains disputed, as does the feasibility of using patterns to deduce how mechanisms operating at different scales are linked (Watt 1947;Cale et al. 1989;Levin 1992;Grimm et al. 1996Wiegand et al. 2003;McIntire and Fajardo 2009). Difficulties arise especially because different processes can result in the same patterns (Levin 1992;Wood 1997;Kendall et al. 1999). In this way, linking pattern and process becomes a problem of model selection, that is, of identifying the most likely mechanism behind a given pattern from several competing hypotheses (Chamberlain 1897;Hilborn and Mangel 1997).
Alpine tree-line ecotones are particularly suited to the pursuit of these questions, given that subtle variations in the external environment and biotic processes along the altitudinal gradient result in the emergence of apparent spatial patterns (Slatyer and Noble 1992). The transition from subalpine forest to alpine communities is characterized by changes in tree density and height structure (Holtmeier 2009). Contrasting tree-line physiognomies can be defined according to the smoothness or sharpness of these changes, that is, whether they are gradual or abrupt, along with the presence or absence of distinctive growth forms such as krummholz (stunted individuals or prostrate twisted wood; Butler et al. 2009;Holtmeier 2009). Differences in physiognomy seem to be related to differences in tree-line dynamics. For instance, a recent meta-analysis showed that abrupt tree lines with krummholz are less likely than gradual tree lines to advance in response to climate change (Harsch et al. 2009). Thus, understanding the causes behind the formation of different tree-line physiognomies is a key aspect for predicting the advance or retreat of mountain forest limits. An upward shift of forests, by increasing the terrestrial carbon sink and reducing the availability of habitat for endemic or endangered alpine plants, would have important implications (Grace et al. 2002). Nevertheless, and despite the ubiquity of tree-line ecotones, disentangling the causes behind the formation of different patterns remains an elusive goal (Wardle 1981;Brown 1994;Körner 1998;Kullman and Ö berg 2009).
Theoretical work on tree lines has focused mainly on the emergence of abrupt or diffuse boundaries along smooth abiotic gradients as a result of positive feedback mechanisms promoting habitat amelioration through individual interactions (Wilson and Agnew 1992;Malanson 1997;Alftine and Malanson 2004). This kind of approach ignores changes in tree-line features other than tree cover and simplifies the description of other processes important for tree-line formation. For example, it is usually assumed that responses in physiological and demographic processes are linearly related to changes in environmental factors along the gradient, whereas nonlinear, threshold-type responses are expected in these harsh environments (Tranquillini 1979;Körner 1998). Another line of thought has focused on dispersal limitation (Dullinger et al. 2004;Albert et al. 2008), especially because it determines the ability of plant populations to deal with climate change (Clark et al. 1998;Jeltsch et al. 2008). Nevertheless, its importance to the formation of tree lines remains disputed, given the short distances to source trees and the frequent observation of seedling populations above the tree line (Körner 1998;Batllori et al. 2010). On the other hand, previous work has also indicated that, at least qualitatively, different tree-line physiognomies can be explained in terms of different responses of individual growth and mortality along the tree line (Wiegand et al. 2006). However, a quantitative assessment of the predictions made by different models with actual data gathered in the field is lacking. This is a common challenge in detailed (individual-based) simulation models (e.g., Wiegand et al. 2003Wiegand et al. , 2004aWiegand et al. , 2004b. Here, we take advantage of intensive sampling conducted at four Pyrenean tree lines that differ in their physiognomies to test models designed to represent alternative hypotheses on the functioning of tree lines. We extended for this purpose a semimechanistic individual-based model (IBM) that considers basic demographic processes and interactions such as competition and facilitation (Wiegand et al. 2006). We designed a set of model experiments to assess how changes in growth, mortality, and recruitment along the altitudinal gradient interact to generate real tree lines. The structure and most model rules were strongly constrained by the known biology of the study tree lines, and growth and mortality along the altitudinal gradient were described with flexible functional forms that allow nonlinear responses. The technical task is to find for each of the competing models parameterizations that generate outputs that are consistent with our detailed data and to select a best model, considering both model fit and complexity. However, model fitting is quite challenging in the case of computationally demanding models like IBMs (e.g., Wiegand et al. 2003Wiegand et al. , 2004aDeAngelis and Mooij 2005;, and currently there is no established approach that would parallel model selection for statistical models (e.g., Burnham and Anderson 2002). To overcome these difficulties, we propose the application of Bayesian and Markov chain Monte Carlo techniques for IBM parameterization (Gelman et al. 2003;Robert and Casella 2004) and of informationtheoretical approaches for model selection (Burnham and Anderson 2002).
The overall aim of our study is to find out whether different tree-line physiognomies reflect a local adjustment in a set of prevailing, common mechanisms (i.e., the same model is selected for all four study sites), the array of processes underlying tree-line formation changes among sites (different models are selected for different sites), or local heterogeneities and site idiosyncrasies dominate over prevailing mechanisms (none of the models produces an output compatible with the data). More specifically, we ask three questions that reflect gaps in the current understanding of tree-line formation: (1) What is the role of nonlinear responses in growth and/or mortality for the emergence of different tree-line types? (2) What is the role of dispersal limitation in shaping tree-line patterns at local scales? (3) Which processes caused krummholz at some sites and not at others, and do tree lines with and without krummholz show different dynamics? In this way, we conducted a strong inferential test on the validity and adequacy of the different models (Platt 1964;Hilborn and Mangel 1997;Wiegand et al. 2003), which provides a cue as to the minimal subset of processes important for the emergence of actual tree-line patterns.

Pyrenean Tree-Line Species
The studied ecotones are dominated by Pinus uncinata Ramond ex DC, which reaches its southwestern distribution limit in the Iberian Peninsula. In the Spanish Pyrenees, this species forms most of the alpine tree-line ecotones on any substrate and at any exposure (Ninot et al. 2007). Pinus uncinata is a long-lived, slow-growing, and shade-intolerant conifer. Its small-winged seeds are primarily dispersed by wind early in spring (Cantegrel 1983). This pine forms dense forests between ∼1,700 m and ∼2,300 m a.s.l. At the higher extreme of this interval, canopy density and tree height diminish suddenly, forming the alpine timberline or forest limit (e.g., where the cover of individuals at least 5 m high drops below 30%-40%). Potential tree-line elevation reaches between 2,200 and 2,450 m in the Pyrenees, depending on continentality, exposure, and land form (Carreras et al. 1996).

Study Sites and Characterization of Tree-Line Ecotones
Four Pyrenean tree line ecotones-Ordesa, Tessó , Capifonts, and Portell-were selected for this study (figs. B1, B2 in the online edition of the American Naturalist;  García-Ruiz 1988) and that early in the twentieth century farms were moved to the bottom of valleys. Aerial photographs of the study sites (1946, 1957, and 1990) and demographic analysis indicate no recent tree-line shifts except at Portell. Indeed, a de-tailed analysis of the eastern Pyrenees comparing aerial photographs (1956,2006) indicates that P. uncinata has expanded its distributional range mainly in the lower part of its distribution, recolonizing previously affected areas (Améztegui et al. 2010). A more detailed description of the characteristics and recent dynamics of Pyrenean tree lines can be found elsewhere Gutiérrez 1999, 2004;Batllori and Gutiérrez 2008;Batllori et al. 2010).
The tree line at Ordesa presents an abrupt decrease in tree height. A morphological sequence from unistemmed trees to multistemmed, shrubby krummholz, forming a belt above the tree limit, is clearly identifiable. Some of the krummholz stems grow several meters in length parallel to the floor, and it is not uncommon that they propagate by layering. However, these same individuals can also grow vertically under benign environmental conditions. Krummholz and seedling individuals showed a weak positive spatial interaction, indicating the potential importance of facilitation (Camarero et al. 2000). Tessó is characterized by a gradual decrease in tree height upslope and the lack of a krummholz belt, with the tree line ending in a sparse, low-density seedling belt. Capifonts and Portell present intermediate tree-line physiognomies, with less extreme changes in tree height along the ecotone than at Ordesa. Although they lack a true krummholz belt, an upper belt of small (i.e., height ! 0.5 m), slow-growing individuals (old seedlings) is evident at both sites. True krummholz individuals and old seedlings are functionally similar. They effectively modify snowdrift accumulation, enhancing the survival of nearby seedlings (Batllori et al. 2009 and references therein), but given their smaller size, this facilitative effect is expected to be less important in the case of old seedlings.
At each site, a rectangular plot parallel to the main slope was delimited ( at Ordesa and Tessó ; 140 m # 30 m 180 at Capifonts; and at Portell), m # 40 m 140 m # 40 m and all individual trees were mapped and measured and their age determined by dendrochronological methods (692 P. uncinata individuals at Ordesa, 259 at Tessó , 790 at Capifonts, and 643 at Portell). Data for each individual tree included X and Y coordinates of the center of each main stem, stem height (h), age (a), and vital state (dead or living). Individuals were classified according to their height and age: seedling ( m and years), , and krummholz or old seedling ( m and h ≥ 4 h ! 0.5 years; hereafter krummholz). This classification a 1 10 served also as the basis for characterizing tree-line physiognomy (see "Summary Statistics and Patterns").

Model Formulation
The choice of a spatially explicit, individual-based model was strongly suggested by our data, which are individual based and spatially explicit. We adopted and extended a previous model by Wiegand et al. (2006) that was designed to determine whether autecological processes alone can explain different types of observed tree-line patterns and to identify factors leading to these types along smooth abiotic gradients. The model included the processes and mechanisms that are expected to be involved in the formation of tree lines. Here, we modified the model to allow for nonlinear responses of growth and mortality along the tree-line gradient and for more realistic representation of dispersal limitation within the recruitment process. In the following, we provide a succinct description of the model (schematized in fig. 1), stressing the changes introduced with respect to the previous version. A complete description of the model is available in appendix C in the online edition of the American Naturalist, including the ODD (overview, design concepts, and details) standard protocol for presenting descriptions of IBMs (Grimm et al. 2006).

Model Output
The basic schedule of the model includes the following processes and environmental constraints: individual growth and mortality, competition, seedling establishment, and krummholz-to-seedling facilitation ( fig. 1). All simulations were run in a -m transect with periodic 60 # 180 and absorbent boundaries on the lateral and longitudinal (i.e., altitude) sides, respectively. Each individual tree was defined by four state variables: position (x, y), height (h), status (alive or dead), and age (a). Each model simulation was run until the coefficient of variation (CV) of simulated adult densities was below an arbitrarily set threshold ( , estimated over a moving window of 50 years) CV ! 0.075 that we considered a steady state. Because of the internal stochasticity of the model, some variability was evident after the steady state was reached. However, this process error was more than one order of magnitude smaller than the 90% Bayesian confidence intervals derived from samples of the posterior distribution of the models. In this way, and given our interest in average model behavior, model outputs were averaged over 10 realizations captured from single simulations separated by 50 years.

Growth and Mortality
Growth was described as a Gompertz function dependent on tree age and was directly parameterized with data on maximum tree heights of individuals from the different stages (see supplementary material in Wiegand et al. 2006). For krummholz individuals, horizontal growth is usually more important than vertical growth because they may form mats. This property entered the model through facilitative effects (see "Seedling Establishment"). Mortality declined exponentially with age (Monserud and Sterba 1999) and reached an asymptotic value for older trees. For seedlings, this relationship was in some cases reduced by the positive effect exerted by krummholz mats (see "Seedling Establishment"). Dead individuals were immediately removed from the grid unless they were adults, in which case they were retained for 45 more years (assumed to represent an approximate time for total stem decomposition). Competition between adults followed a "zone-ofinfluence" (ZOI) approach (Bella 1971;Schwinning and Weiner 1998), parameterized with field data (see supplementary in Wiegand et al. 2006). Direct competition within or between the remaining life stages and growth forms was not considered to exert an important role in the dynamics, except in the case of seedlings at establishment.

Seedling Establishment
Seedling establishment was allowed only in safe sites, defined as (1-m 2 ) cells (of a grid superposed to the study transect) not subject to competition by preemption of space, delimited as a circular zone of influence with radius r zoi around adult trees ( , where h is the height r p 0.2h zoi of the focal tree; Wiegand et al. 2006), or as cells free of dead trees and not occupied by more than one krummholz individual. On the basis of the spatial-point-pattern analyses by Camarero et al. (2000) and field experiments by Batllori et al. (2009), a facilitative effect of krummholz on seedling performance was implemented. Grid cells within a radius of 1.5-2.5 m of a krummholz-occupied cell experienced enhanced probabilities of seedling establishment and survival, depending on krummholz height (i.e., ). This was f p 0-0.3 facil necessary to prevent an unrealistic positive effect exerted by the old ( years) but small-sized individuals pres-a 1 10 ent at these sites (i.e., the old-seedlings class; see "Study Sites and Characterization of Tree-Line Ecotones").

Altitudinal-Gradient Response Functions of Growth and Mortality
In mountain ranges, the abiotic conditions change along the altitudinal gradient, thereby producing the tree line. Our hypothesis is that the change in abiotic conditions primarily affects growth and mortality of individual trees and that this effect modulates growth and survival multiplicatively. Instead of using a fixed formulation to represent these a priori unknown relationships, we derived their shape from the data (i.e., a semimechanistic model; Ellner et al. 1998;Kendall et al. 1999). We used two functions, g g (y) and g m (y), with a flexible, logistic functional form to describe the decrease in the individual growth rate and survival (eqq. [C8] and [C9], respectively), depending on the longitudinal position y on the transect, where p is an index for growth (g) or mortality (m), a p and b p are parameters to be fitted, and . The 0 ≤ y ≤ 180 parameter a p ranges between 0 and 0.5 and influences the steepness of the gradient response, with higher values corresponding to steeper gradients. For we have a p 0 p (i.e., no gradient response), and for g (y) p 1 a p 0.5 p p the gradient drops within 20 m from values of near 1 to near 0. The parameter b p shifts the gradient in the ydirection and ranges between 0 and 180 m. Examples for the g p (y) are presented in figure A1 in the online edition of the American Naturalist. For model versions without gradient response, we set .
Seedling establishment was allowed only in safe sites (see "Seedling Establishment"), and krummholz-mediated facilitation (f facil ) increased the establishment probability locally. We tested whether dispersal limitation, included by allowing spatial variation in the probability of seedling emergence, improved the ability of the model to reproduce observed tree-line patterns with respect to a simpler model assuming an infinite, spatially homogeneous seed bank. To this end, the probability of seedling establishment at each cell was a function of the distribution of distances to adults, of tree fecundity (expressed as an allometric function of tree size), and of the form of the dispersal kernel (Ribbens et al. 1994): where p rep and f facil are parameters that scale overall establishment success and facilitation strength, respectively, s(x, y) is seedling density at the location with coordinates (x, y), F i is the fecundity of tree i, and K(x i , y i ) is the dispersal kernel. Note that we assumed that the process is isotropic (i.e., seedling input does not vary as a function of the direction from the source trees).
We examined the lognormal function as a dispersal kernel (i.e., a lognormal is a convex, zero-at-zero leptokurtic kernel), given that its high performance for wind-dispersed species is well documented (Stoyan and Wagner 2001;Greene and Calogeropoulos 2002;Greene et al. 2004 where S and L are kernel parameters to be fitted and d(x i , y i ) is the distance between the locations (x i , y i ) of the reproductive adult tree i and the target cell location (x, y).
Other dispersal kernels (exponential and 2Dt) did not provide a better fit (results not shown).

Model Analyses
The assessment of alternative hypotheses using information-theoretic approaches of model selection is well established in ecology (Hilborn and Mangel 1997;Burnham and Anderson 2002;Johnson and Omland 2004). These methods select a best model or a best set of models by balancing model fit and model complexity. The approach is straightforward for simple models (e.g., regression models) in which the likelihood function can be treated analytically and uncertainty in model parameters can be handled with classical statistical theory. However, this task is much more intricate in the case of complex models such as IBMs, which are usually computer intensive, have a high dimensionality (i.e., a large number of parameters), and may include nonlinear relationships. In these cases, an explicit assessment of a potentially complexly shaped likelihood is complicated by the curse of dimensionality and potential correlations among parameters that confound assessments of uncertainty (Van Oijen et al. 2005;Hobbs and Hilborn 2006). One approach that is increasingly used in ecology to solve these problems is Bayesian parameter estimation with Markov chain Monte Carlo (MCMC) methods. Ecological studies that use Bayesian parameter estimation have been conducted, for example, for inverse modeling of seed dispersal kernels (Clark et al. 1999;Martínez and González-Taboada 2009) and in forest modeling (Van Oijen et al. 2005;Purves et al. 2008). However, these techniques have been rarely applied to more complicated simulation models, such as IBMs. Yet increasing computational power and recent development of refined techniques of computational statistics suggest that we should also be able to apply methods of model selection and model parameterization based on Bayesian and MCMC techniques for individual-based simulation models.

Summary Statistics and Patterns
Exact likelihood estimation using the raw data is intractable for complex simulation models, and one approach is to use instead summary statistics of the data that represent the maximum amount of information in the simplest possible form (Csilléry et al. 2010). These summary statistics should capture different characteristic features of the system that are relevant to the scientific question asked (Wiegand et al. 2003(Wiegand et al. , 2004a(Wiegand et al. , 2004b and are combined to assess the model fit. The focus on multiple summary statistics is important to avoid the problem that different models can reproduce the same pattern . In our case, we need to condense the observed data on tree position, age, and height into summary statistics that provide a good description of tree-line physiognomy. Tree-line physiognomy is usually defined in terms of its abruptness, characterized by changes in either tree density or tree height with altitude, and by the presence or absence of krummholz belts (Holtmeier 2009). Apart from these features, we also considered variation in seedling density along the tree line, because of the presence of apparent seedling belts at some of the sites and the importance of this stage for tree-line dynamics. We therefore used the density of seedlings, adults, and krummholz individuals and the mean height of all individuals (excluding krummholz) at different lo-cations as summary statistics. Other information, including variation in age distribution and the density of saplings and poles along the gradient, was disregarded, given that these features were redundant with the summary statistics selected.
The four summary statistics (indexed by i) were estimated in the same way from the field and model output. On the basis of previous work (Camarero and Gutiérrez 1999), the transects were divided into 5-m-wide, 20-mlong subplots, from which the mean of the individual obs m ij densities or heights and their standard deviation were obs j i calculated for 20-m altitudinal bands (indexed by j). These quantities were employed to assess model fit.

Assessing Model Fit
We assumed a normal likelihood function to describe deviations from the observed mean densities and height along the tree line for each of the four observed summary statistics:

Model Parameterization and Fitting
Extensive model simulations have to be conducted to find the parameterizations that minimize the likelihood function given in equation (5). This requires the use of effective optimization algorithms. Because we were interested in finding not only the minima of the likelihood function but also measures of uncertainty that consider correlations among the parameter estimates, we followed a Bayesian strategy to determine, given the observed data, the most likely combination of parameter distributions for each model. This approach is naturally linked to MCMC methods, taking full advantage of the output of iterative, stochastic optimization algorithms to derive and propagate parameter uncertainty. Because no analytical expression can be derived for the likelihood of the parameters, we take advantage of the sampling distribution of the observed data D The sampling distribution is the joint probability density function of the sample and is equivalent to the likelihood function of model parameters (e.g., Casella and Berger 2001). We then relied on l(vFD) inverse methods to obtain the posterior distribution of model parameters, using Bayes's formula, is the likelihood function and p(v) is the f(DFv) prior distribution of the parameter set v. Bayes's formula is used to gain knowledge of the posterior distribution of the parameters , beginning with some prior beliefs p(vFD) p(v), by updating with the information contained in the sampling distribution . The posterior distribution f(DFv) can be resolved up to a constant, the marginal p(vFD) distribution of the data. The Bayesian approach is a type of probabilistic inversion because it tells us how to infer the parameters v from the data D in terms of how D is determined by the parameter v (Van Oijen et al. 2005).
No prior information was included, apart from appropriate ranges for model parameters and a set of constants directly estimated from our data to describe some processes (e.g., maximum potential growth, maximum scale of competitive effects, and a baseline recruitment level). The interval of the parameters was determined on the basis of either structural constraints or literature information. Parameters were updated with an adaptive algorithm based on Metropolis steps (Robert and Casella 2004). At each step, new values were proposed for a randomly chosen subset of the parameters by means of a truncated normal distribution centered on previous values and with a scale tuned to attain a target acceptance rate of 0.3 (Gilks et al. 1998). The acceptance of proposed parameter sets was based on the ratio of their posterior distributions, obtained via equation (6). A total of 40,000 Monte Carlo steps were run, which included a burn-in period of 10,000 steps. The last 30,000 steps were then thinned with a lag of 30 steps to avoid serial autocorrelation. The resulting 1,000 steps were used to summarize the posterior distributions of the parameters and to propagate uncertainty in parameter estimation to projections of tree-line patterns. A pseudocode of the model-fitting procedure is available in appendix D in the online edition of the American Naturalist.

Model Selection and Simulation Experiments
We took into consideration eight alternative models that arose from the combination of the presence or absence of a functional form for the environmental gradients of growth inhibition and mortality enhancement (eq. [1]) and for the use of a homogeneous distribution of recruits versus dispersal limitation (eq. [2]). For the sake of compactness, different model structures were designated by trinomials, so GrMoDi corresponds to the most complex structure, which includes both growth (Gr) and mortality (Mo) gradients and dispersal limitation (Di). The absence of any of these components is indicated by "Ab"; thus, model AbAbAb does not include any of the above-mentioned processes. On the basis of posterior deviance distributions, we calculated the Akaike Information Criterion (AIC) for each model to rank alternative models in terms of both their complexity and their performance (Burnham and Anderson 2002).

Model Selection
Despite the wide variation in tree-line physiognomy among the four study sites, we found for each site at least one model that matched the observed tree-line patterns simultaneously ( fig. 2). However, not all models were able to generate tree lines consistent with the data (figs. A2-A6 in the online edition of the American Naturalist). For example, models without a reduction in tree growth along the gradient were not able to produce krummholz. The relative performance of the eight alternative models is illustrated by the posterior distributions of the deviance (fig .  B3 in the online edition of the American Naturalist; lower values indicate a better fit). The median deviance (Dev in table 2) indicated substantial differences in the performance of different models for a given site and that the performance of the same model varied among sites.
The model selection procedure revealed that two models received similar support from the data (i.e., ), DAIC ! 2 except at the Ordesa site. The model that included a gradient response in growth (GrAbAb) was the best model at Ordesa and was among the best models at Capifonts and Portell (table 2). At the Tessó site, which shows a smooth tree line without krummholz, the best model (AbMoAb) included a gradient response in survival. Thus, relatively simple models were able to reproduce observed tree lines, but differences in model performance among sites suggest that the importance of different mechanisms varies among localities. In the following, we explore these differences in more detail.
Ordesa. The Ordesa tree line is characterized by abrupt transitions in adult density and mean height, and its most eye-catching feature is the pronounced seedling-krummholz belt ( fig. 2). Formal model selection based on AIC values favored the model GrAbAb unambiguously, with a weight of 93% (table 2). However, while this model matches the height, krummholz, and adult patterns well quantitatively, we observed some differences in the seedling pattern and a slight overestimation of krummholz density at the highest altitudinal band ( fig. 2).

Tessó.
The tree line at Tessó shows smooth transitions in adult density and mean height and does not present a krummholz belt ( fig. 2). All models were able to produce tree lines without krummholz ( fig. A3), but when compared with Tessó data, two of them (AbMoAb and AbAbAb) received similar support and were selected on the basis of AIC (table 2). The model that ranked first by AIC, AbMoAb, included a gradient response in mortality. However, the simplest model, AbAbAb, received similar support (i.e., ) but did not reproduce well the DAIC ! 2 decline in tree height and adult density ( fig. A3). The model GrAbAb, selected for the other sites, ranked close to the AbMoAb model ( ), but it performed DAIC p 2.6 somewhat more poorly in the seedling pattern ( fig. 2).
Capifonts. The tree line at Capifonts showed an old, krummholz-like seedling belt ( fig. 2). Similar to that at Ordesa, the krummholz belt emerged only for models including a gradient response in growth ( ). The mod-Gr * * els GrAbAb and GrMoAb received similar support (table  2), with GrMoAb yielding a somewhat better match in krummholz and seedling patterns, but at the cost of two more parameters.
Portell. The tree line at Portell showed an abrupt transition in adult density, with an old, krummholz-like seedling belt following the adult belt, a smooth height transition, and almost no seedlings ( fig. 2). None of the models reproduced the adult pattern correctly when a match of the low seedling density was demanded ( fig. A5). To find the reasons for this failure, we removed the summary statistics for seedlings from the likelihood function and repeated the analysis. Under those conditions, the models GrAbAb and GrMoAb that were also favored at Capifonts received  (see table 2). Black lines correspond to field data, while the white lines and the shaded areas represent, respectively, the median and the 90% Bayesian confidence intervals derived from samples of the posterior distribution of selected models. Note that in Capifonts the number of altitudinal bands is 9. similar support (table 2). The predictions for seedlings showed that, given the current model structure, the seedling density should be higher at this site ( fig. 2).

Internal Model Functioning and Biological Questions
Analyses of posterior parameter distributions were conducted for a better understanding of the biological con-straints imposed by model structure. For the model selected at each site, we examined the Spearman rank correlation coefficient between (1) the values of the parameters, (2) the partial likelihood of each summary statistics, and (3) the mean of each summary statistic averaged over all the altitudinal bands (table A1 in the online edition of the American Naturalist ). The structure of the Note: The Akaike Information Criteria (AIC) balances model fit and complexity, with a lower AIC value indicating a better model (e.g., a better fit with lower cost in terms of parameters). To ease model comparison, Akaike weights have been also included. Models with lowest AIC and with are highlighted in boldface. n p p number of parameters for each model; Dev p median deviance; w i p DAIC ! 2 Akaike weight for model i.
a Model structures are designated by trinomials; GrMoDi corresponds to the most complex structure, including growth (Gr) and mortality (Mo) gradients and dispersal limitation (Di). The absence of any of these components is indicated by "Ab"; thus, model AbAbAb did not include any of the above-mentioned processes.
model that included nonlinear relationships as well as interdependences between processes resulted in constraints that were common to all or most of the sites. We found high correlations between parameters determining the shape of the gradient response functions (e.g., a g , b g and a m , b m ) or between the two parameters m 0 and e characterizing survival curves. We also found positive correlations between the mortality parameter m 0 and the probability of establishment p rep , which reflects trade-offs between these processes.

Constraints on Gradient Response
Functions. The high correlation between the parameters a g , b g and a m , b m resulted in relatively narrow Bayesian confidence intervals for the gradient response functions ( fig. 3). Figure A7 in the online edition of the American Naturalist shows the projection of the parameter space for the two gradient parameters. In general, the range of the best parameters was much reduced. The gradient responses of Ordesa and Capifonts were similar and abrupt, with the parameter b g having a narrow range and and , respectively, a 1 0.057 a 1 0.066 g g resulting in a thresholdlike response in growth ( fig. A1). In contrast, the parameter a g showed a narrow range at Portell, with 96.6% of all cases ranging below a value of 0.067, thus resulting in a smoother gradient response function for growth (figs. 3, A1). The a m -b m parameter space for Tessó corresponded to a thresholdlike response of increased mortality toward the end of the gradient, although with large variability above altitudes of 60 m ( fig. 3). This is probably due to the lack of a data upslope and also reflects uncertainty in model selection, with a few parameterizations showing a value of , which yields the a p 0 m response function (i.e., the model AbAbAb). g (y) p 1 m Constraints on Mortality Functions and Establishment. The high correlation between the parameters m 0 and e resulted in low uncertainty in the curves that describe mortality in dependence on age ( fig. 4), especially in the case of sites with krummholz. Indeed, survival was lower at these sites than at Tessó . Nevertheless, the effect of facilitation altered this pattern at Ordesa, where seedling survival was enhanced through this process. The posterior distribution of the establishment probability p rep was, in general, narrow and centered at low values between 0.05 and 0.30, and it was especially low at Tessó and Capifonts ( fig. 4). Indeed, p rep showed a strong positive correlation with mean seedling density. Facilitation seemed to play a secondary role in establishment except at Ordesa (see below), given that the posterior distribution of the parameter f facil was indeed not very different from the uniform prior distribution in the remaining sites ( fig. 4).
Krummholz. The strongest condition for emergence of a krummholz belt was a gradient response in growth. This seems to be tautological because krummholz is a stunted growth form (i.e., old individuals with low stature), so a reduction in growth along the gradient without a reduction in mortality will result in krummholz at the upper end of the gradient. However, we found models that were not only able to generate "trivial" krummholz by this general mechanism but also able to quantitatively generate the observed krummholz belt and other characteristic features of the tree line at the same time.
At Ordesa, the partial likelihood of the summary statistic for krummholz density along the transect showed a negative rank correlation with the facilitation parameter f facil (  ; table A1b) and a positive correlation with r p Ϫ0.45 the parameter a g that determines the steepness of the gradient in growth (  ; table A1b). Thus, the likeli-r p 0.39 Sp hood of matching the krummholz pattern at Ordesa increased with increasing facilitation ( fig. 4) and increasing steepness of the gradient response. There was also a weak negative correlation between the two parameters f facil and a g ( ). Thus, the gradient responses of growth r p Ϫ0.33 Sp and facilitation showed a (weak) trade-off. However, no equivalent relationship was found at the other two sites where krummholz-like individuals were present.

Discussion
Tree-line ecotones are characterized by marked spatial changes at small spatial scales that at regional scales result in a variety of physiognomies. Although many factors have been invoked to explain the formation of these contrasting patterns, we have shown here that the main features of real tree lines can emerge from demographic models that consider nonlinear responses in individual rates of growth and mortality with respect to the altitudinal gradient. Our results also indicate that the variation in tree-line physiognomy observed at our sites reflects changes in the relative importance of these nonlinear responses and that other processes, such as dispersal limitation, may play a secondary role. The model selection procedure favored different models for tree lines with and without krummholz. This is consistent with a recent finding that suggests that tree lines with diffuse forms show different functioning than the ones with abrupt or krummholz forms (Harsch et al. 2009). Obtaining these results, which are based on quantitative comparison of model outputs with detailed data, required novel ways of applying techniques of model selection and Bayesian parameterization to individual-based simulation models. In the following, we discuss advantages and problems of the proposed meth-odology, before returning to the ecological insight gained during this process.

Analysis of Individual-Based Models
Uncertainty in model parameters, error propagation, and ad hoc methods of model selection have been the major points of criticism against IBMs (e.g., Wiegand et al. 2003Wiegand et al. , 2004aDeAngelis and Mooij 2005;) and have hindered their wider acceptance. Our approach proved to be very effective in fitting our models to the observed data. Bayesian calibration provided parameter estimates with measures of uncertainty; it considered the correlations among the parameter estimates and yielded estimates of uncertainty in model predictions. Use of the likelihood function, minimized by Bayesian calibration, provided the link for application of established techniques of model selection (Burnham and Anderson 2002). The approach presented here and by others (e.g., Van Oijen et al. 2005;Piou et al. 2009) parallels recent developments in the "new statistics" (Johnson and Omland 2004;Hobbs and Hilborn 2006) that expand statistical regression modeling to a more general class of mathematical models, where the functional forms and definitions of parameters are not chosen for statistical reasons (as in more traditional models used by ecologists for statistical inference; Hobbs and Hilborn 2006) but where variables and parameters explicitly symbolize ecological states and processes. Modified versions of our model could be applied to other systems dominated by environmental-stress gradients, such as marine intertidal zones or desert grasslands. The method of parameterization and model selection presented here can be widely applied in (individual-based) simulation models. Published studies directly amenable to our methods include studies on wolf (Canis lupus; Marucco and McIntire 2010), brown bear (Ursus arctos; Wiegand et al. 2004a, 2004b), and grass-shrub steppes (Paruelo et al. 2008). We expect that these methods may turn model selection and evaluation in this type of model into a more transparent, effective, and efficient exercise.

Summary Statistics and Model Parameterization
Assessment of model fit in IBMs relies on summary statistics, which are aggregated measures that summarize model output for a potentially high number of individuals. This allows for inference that would be intractable if the full data were used. Note that this is also at the heart of Approximate Bayesian Computation (ABC; e.g., Marjoram et al. 2003;Csilléry et al. 2010). Summary statistics should represent the maximum amount of information in the simplest possible form, quantify the features of the system that are relevant for the scientific question asked, and be representative of the system dynamics (DeAngelis and Mooij 2003; Wiegand et al. 2003Wiegand et al. , 2004aCsilléry et al. 2010). Also, to overcome the wellknown problem that different processes can result in the same pattern, it is necessary to consider several independent summary statistics simultaneously (Kendall et al. 1999;Reynolds and Ford 1999;Wiegand et al. 2003;Komuro et al. 2006). Otherwise, the inverse approach may fail because the summary statistics do not contain enough information to constrain the values of model parameters (Wiegand et al. 2004a).
The failure of all models to generate the observed seedling pattern at the Portell site points to potential problems if a summary statistics is not representative. Models that generated the observed low seedling density were unable to match adult densities. However, when allowing our method to find the seedling densities that were most compatible with the other three summary statistics (by removing the seedling summary statistic from the likelihood function), the model GrAbAb, which was also among the best models at Ordesa and Capifonts, was favored. This model predicted that seedling densities at Portell should be much higher than observed. While the densities of krummholz and adult trees, as well as the pattern of mean height along the gradient, represent the cumulative out-come of demographic events and environmental constraints over decades, the observed seedling population is a rather short-term snapshot. Seedling emergence is subject to substantial year-to-year variability (Lloyd and Graumlich 1997;Shiyatov 2003;Camarero and Gutiérrez 2004) and may additionally depend on small-scale heterogeneities not included in the models presented here (Batllori et al. 2009). Note that we also observed some difficulties in generating the pattern of seedling density at the Ordesa site ( fig. 2). In this case, observed seedling densities were higher than those predicted by the models. Thus, the observed seedling densities were probably not representative of equilibrium conditions at Portell and Ordesa sites and were likely to have been caused by year-toyear variability in model parameters regarding seedling emergence and survival. Our modeling exercise thus generated a hypothesis that can be tested and evaluated with additional field investigations. The lesson of this for the general approach is that a critical assessment of the "representativeness" of the summary statistic with respect to the model structure is required. On the other hand, if our model were to explicitly consider the effect of climatic variability on seedling establishment, the seedling summary statistic would probably be more informative.
The fit of complex, spatially explicit, individual-based models is rarely straightforward. The computer-intensive nature of individual-based models and their high dimensionality (i.e., large number of parameters) make commonly employed methods such as the Latin hypercube (McKay et al. 1979;Iman et al. 1981) impractical because of the curse of dimensionality (Bellman 1957). For instance, a preliminary study of our model using a Latin hypercube design did not provide any model fit better than the 1,000 samples from the posterior distributions. A problem is the blind search of this kind of algorithm. Therefore, an effective optimization algorithm is required to find the minima of the likelihood function in the parameter space. This task could be solved with several methods, including hill-climbing (Russell and Norvig 2010), simulated annealing (Kirkpatrick et al. 1983), or evolutionary programming approaches (Fogel et al. 1966;Komuro et al. 2006;Duboz et al. 2010). We adopted a Bayesian approach to take advantage of the natural coupling between the Bayesian approach and MCMC methods, which provides uncertainty estimates in model parameters and predictions that account for correlations between different parameters. This is not possible under alternative approaches proposed for IBMs. Another potential advantage of the Bayesian approach is that it allows the incorporation of evidence from previous experience in a natural way.

Tree-Line Dynamics
Our models were able to fit the detailed data from the four contrasting tree lines well quantitatively. This result shows that tree-line physiognomy at our sites is not determined entirely by local heterogeneities and site idiosyncrasies. However, the array of processes underlying tree-line formation changed among sites with and without krummholz, but in all cases the response of individual growth and mortality along the altitudinal gradient played a leading role in determining tree-line physiognomy. Other processes expected to exert a great influence were not necessary to obtain tree-line patterns compatible with our data, especially in the case of dispersal. Dispersal limitation was included by an explicit dependence of seedling establishment probabilities on the distance to parental plants. In contrast to homogeneous dispersal, more realistic dispersal kernels resulted in better model fits, although they were not favored by the model selection procedure because of the increase in model complexity. Nevertheless, these improvements were not comparable to those resulting from consideration of nonlinear response functions in growth and mortality.
Pinus uncinata is known to have a limited dispersal distance. For example, Camarero et al. (2005) found in a seed-release experiment that most P. uncinata seeds traveled between 4 and 6 m, with a mean distance of 27 m. Our results indicate that dispersal limitation may not play a substantial role in the formation of "equilibrium" treeline patterns, and they support, for example, the view of Körner (1998) that seedling establishment is not determinant in tree-line formation and dynamics but that subsequent processes like the emergence from the warmer boundary layer near the ground or from a shrubby ground cover may be more relevant. Furthermore, Batllori et al. (2010) found that recruitment is not constraining current P. uncinata tree-line dynamics in the Pyrenees. On the other hand, the importance of dispersal is expected to be more apparent at a landscape scale and in a metapopulation context or when predicting the shift of tree-line ecosystems due to climate change (Dullinger et al. 2004;Albert et al. 2008; but see Noble 1993). In this way, although establishment is an essential phase of the life cycle of P. uncinata, its role in shaping the structure of tree-line ecotones turned out to be secondary at the local scale when individual demographic rates are allowed to vary along the altitudinal gradient.
We found that the uncertainty in the predicted gradient response functions was low and showed a characteristic shape at each site ( fig. 3). However, our semimechanistic model does not provide an explanation for the variation of growth and mortality rates along the tree line, because we did not include more detailed mechanisms, such as the ecophysiological processes described in Körner (2003) and Smith et al. (2009). The current model also ignored more detailed temporal changes in the shape of the gradients due to year-to-year climatic variability. Instead, we followed the idea of semimechanistic models (e.g., Ellner et al. 1998;Kendall et al. 1999) and let our model parameterization and selection procedures find the shapes of the gradient functions that were most consistent with the data. This reduced the initial problem of explaining different tree-line patterns to the problem of explaining different response functions for growth and mortality rates along the altitudinal gradient, generated specific hypotheses as to how to improve the model (e.g., including year-to-year variability in seedling emergence and mortality), and provided predictions about the expected shape of this response that could be tested in the field. Tree lines with a pronounced variation in growth along the gradient result in abrupt ecotones with krummholz. On the other hand, nonlinear responses in mortality alone can explain the emergence of smooth gradients. Ordesa and Tessó nearly correspond to these two extremes of a continuum among abrupt tree lines with krummholz and smooth transitions, respectively. The other two study sites, Capifonts and Portell, presented intermediate characteristics. Although the best model at both sites included only a gradient for growth (e.g., as in Ordesa), models simultaneously including nonlinear response functions for mortality and growth received similar support. This reinforces the idea of a continuum between tree lines modulated almost exclusively by the response in individual growth and those regulated only by external causes of mortality.
In a recent meta-analysis, Harsch et al. (2009) found that tree-line sites with a diffuse form were more likely to have advanced in response to climate change than those with an abrupt or krummholz form. According to our results, this can be explained in terms of the different functioning between these two kinds of tree-line physiognomies. If mortality constrains tree-line position, we would expect that a change to more favorable conditions at higher altitudes would cause an advancement of the tree line as a result of increased tree survival. Recruitment and survival at the tree line may be highly sensitive to winter conditions (Gray et al. 2006;Batllori et al. 2009), and a series of favorable years may therefore allow survival of established seedlings at higher altitudes. Interestingly, Harsch et al. (2009) found that tree lines with higher rates of winter warming were more likely to show advance. However, this process may be very stochastic, because a single cold year could destroy the gains made over several warmer winters (Körner 1998). On the other hand, the case of abrupt tree lines with krummholz is quite different. Our model suggests a leading role for growth limitation in determining tree-line structure, so we may expect a slower response under more favorable conditions, involving a conversion of krummholz to flagged krummholz or even larger trees but without immediate establishment of new individuals at higher altitudes. Further, where factors such as topographically induced strong winds are responsible for the presence of stunted growth forms at the tree line, increased temperatures may not promote tree-line migration at higher altitudes at all. Another related result was the greater facilitative effects detected at Ordesa, consisting of an increased recruitment and survival of seedlings mediated by krummholz. The lack of this effect may further slow colonization fronts. These two factors may explain that the probability of detecting tree-line advance is nearly four times higher in smooth tree lines than in abrupt tree lines with krummholz (Harsch et al. 2009).
The strong support found in our study for the notion that the structure of tree-line ecotones is mainly caused by subtle changes in growth and mortality processes along the altitudinal gradient provides several cues for future studies of tree-line dynamics. Our results suggest that understanding and predicting the dynamics of these ecosystems require a characterization of environmental response functions for growth and survival in terms of factors varying along the tree-line ecotone and that different physiognomies can give us a cue as to the processes involved. Environmental factors such as temperature changes along the ecotone, snow-wind interactions, and perhaps in a near future, increased drought stress, are good candidates for characterizing these responses and their temporal variation. Nevertheless, the possibility that human disturbance has influenced the observed tree-line patterns cannot be entirely ruled out. Our results may therefore apply to tree lines that have been subject to certain human impact in the eighteenth and nineteenth centuries. Further research is needed to clarify the relative importance of different external drivers, especially in order to apply our results in the context of expected climate change impacts. Also, dispersal has not emerged as a key process at the scale considered here. In this way, and from a modeling perspective, characterizing the process of dispersal at a regional scale can be enough to predict future tree-line migration, although delayed responses associated with growth limitation cannot be ruled out, especially in the case of tree lines with krummholz. In either case, the current threats to these valuable ecosystems make the development of models capable of predicting the response of alpine tree-line ecotones an urgent need.