GET THE APP
Received Date: October 29, 2011; Accepted Date: November 07, 2011; Published Date: November 09, 2011
Citation: Gaucherel C (2011) Self-Organization of Patchy Landscapes: Hidden Optimization of Ecological Processes. J Ecosys Ecograph 1:105. doi:10.4172/2157-7625.1000105
Copyright: © 2011 Gaucherel C. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and and source are credited.
Visit for more related articles at Journal of Ecosystem & Ecography
One of the main barriers in the understanding of landscape dynamics is the high spatial variability of the surface patterns (vegetation and land cover) to which ecological processes are intimately linked. The aim of this paper is to present some newly found scaling properties for forested landscapes. Furthermore, it is advocated that patchy landscapes can sometimes be self-organized by optimizing some effective functional. A neutral landscape model has been developed in order to test this self-organization hypothesis. This model was built on the basis of a simple function, called the "Hamiltonian" analogically to physical and biological systems. The Hamiltonian is then minimized to optimize the identified landscape interactions. Fully controlled data coming from five different hundred-year runs of a process-based model appeared to be self-similar over five magnitude orders, without being explicitly simulated. The neutral model is able to reproduce the studied observations and to easily model Optimal Patchy Landscapes. The limits to this parsimonious approach that requires only one parameter (the Hamiltonian slope in loglog plot) are also discussed. The links between the effective Hamiltonian and the ecological processes still need to be investigated. Finally, such landscape Hamiltonian function appears to be a fruitful theoretical framework to describe various landscapes and potentially opens the way to a more complete dynamic landscape theory.
Categorical mosaic; Ecological complexity; Self-similar pattern; Multiscale; Rural landscape.
Ecological processes are intimately linked to landscape dynamics observed in the high spatial variability of surface patterns (land cover changes). Human decisions and natural forcing are key drivers in landscape systems [1,2]. They occur over a wide range of spatial scales, from the microscale of the individual landscape unit (patch) trough the mesoscale of farmer property, to the macroscale of regions or states. The goal of this paper is to present some newly found invariance properties for forested landscapes. These properties refer to the invariance of the probability distributions for a forest parameter under a wide change of spatial scales [3-5]. Ecological processes and more specifically the organization of patchy landscapes by land cover geometry and topology involve a certain type of optimization. To check this interpretation of optimal landscape organization by modeling is the second goal of this paper.
The idea that patterns in nature may be obtained by optimality principles based on energy is not new [6,7]. Furthermore, many steadystate patterns of minimal energy dissipation develop self-similar (i.e. fractal) spatial arrangements. More recently, it has been shown that power laws can generally be interpreted as a preferential attachment, which is nothing else than a kind of optimisation . Closer to ecological topics, fractal dimensions have been extensively used as averaged landscape metrics for describing the spatial distribution of various habitats. To decide whether self-similar properties are common in landscape structures is still under debate as proved by the recent papers about scaling laws observed in arid vegetation patterns . Other works recently demonstrated the common self-similar characteristics of landscape heterogeneity when defined by land cover diversity , or the characteristic of disturbance processes .
Yet, there is at present little theory on how to link these properties to generating processes. There would therefore be very helpful to be able to predict the scaling exponents based on ecological interactions in space . If some landscapes exhibit self-similar properties, are rural landscapes also self-organized? A way to test this hypothesis is to try modeling them with parsimonious inputs and an optimization principle, as it has been done with Neutral Landscape Models (NLM, [13,14]. As detailed landscape data are rare over a long time period, I studied a forested landscape evolving under several mechanistic processes such as colonization, competition and regeneration . This reference landscape consisted in a forested sector that covered more than a thousand acres and was made up of 220 patches, allowing a reliable probability distribution analysis for the tree-density of the patch sample. The multiscale study showed how a simple potential force field (based on a Hamiltonian function) is able to capture the selforganization of the reference landscape.
State-of-the-art and justification
Neutral landscape models: Neutral models are more and more used in ecological studies. More specifically, NLM (sensu  were first proposed to verify that (biased) random functions are able, or not, to simulate an observed landscape pattern (not to be confounded with neutral community theory). This approach is a kind of null-hypothesis test. NLM already brought many interesting insights to ecology, namely the definition of landscape analysis indices, critical mosaics for population dynamic studies and the first stage of a general spatial complexity theory . Yet, most of these NLM focused on raster mosaics, while more realistic patchy landscapes (i.e. made of contiguous and uniform patches) are still poorly modeled . In addition to the modeling challenges, a theoretical framework describing patchy landscapes still lacks. I address this question here, as an analogy with physical and biological systems.
The hamiltonian in physics: Nowadays a large scientific community uses Cellular Potts Models (CPM). It began as an extension of the Potts model, itself being an extension of the Ising model, a simple early model of ferromagnetism . In the Ising model of magnetization, individual atoms are characterized by their magnetic moments or spins σ at position with constant interaction energies J with their neighbors. According to Boltzmann theory, the relative probability of any configuration of spins is:
where is the total configuration energy (also called the Hamiltonian), k is the Bolzmann constant and T is the system temperature. The higher the energy of a configuration, the less probable it is. In the absence of any external field, the Ising Hamiltonian is the sum of interactions between all pairs of spins that are nearest-neighbors :
Later on, Potts models were proposed to allow multiple degenerate values of spins (i.e. the energy of a link depends only on whether the neighboring spins are the same or different and not on their particular values) . System dynamics were then studied by the use of the Metropolis algorithm modelling local changes of the individual domain states [19,20]. The Metropolis algorithm leads the system to its equilibrium configuration(s) at temperature T (corresponding to a global energy minimum), where the time evolution is stationary. The T=0 Metropolis algorithm accepts only spin flips which lower the energy, thus bringing the system to a local energy minimum (which is different a priori from the global energy minimum).
This first use of Metropolis methods led to a great expansion of the range of questions that Potts models could address, in particular concerning equilibrium and out-of-equilibrium dynamics. Fruitful analogies have been proposed to model metallic grain growth , liquid soap froths and finally tissue cells . In the latter model, also called Cellular Potts Models (CPM), a continuous tissue is discretized in cell configuration onto a square lattice, for which each site represents a collection of cells of a particular type (cell indices are used here). In this case, the system's Hamiltonian is rewritten:
where (δ (x, y) = 0 if x ≠ y and 1 if x = y) and where the (symmetric) interaction or boundary energy between two sites depends on the cell type at the considered link. The second term of the sum imposes a volume constraint to the cellular system, where v (σ ) is the volume in lattice sites of cell σ, Vt its constant target volume and λvol (τ ) the strength of the constraint. Many improvements of the original CPM model have been developed and successfully compared to observations .
It is worth noting that such a modelling approach is broadly used, even in human sciences. A long time ago, pioneer works emphasized the connection between optimal control theory and classical themes of capital theory . Fundamental questions, the definition of income, how it should be measured or what its relation to wealth is, have been routinely solved using one-dimensional (i.e. non spatial) optimisation by building Hamiltonian functions describing the studied economical system.
Towards environmental hamiltonians
Closer to our ecosystem thoughts, dendritic river networks offer a useful analogy for patchy landscapes . In these hydrographical systems, the observed power law distributions of discharge mass are closely linked to energy dissipation within the river basin . The total cumulative area draining into a river link is used in these works as a surrogate variable for discharge and presents a self-similar behaviour observed in real river basins. This scaling behavior often measured in natural drainage networks reflects a preferential spatial aggregation (i.e. a topology) leading to dendritic patterns. These authors used the exceedence probability of the drainage area characterizing the selfsimilar topology of the hydrographical network to build a Hamiltonian for their system :
where Ai denotes the total contributing drainage area at a point i of the network, and the coefficient linking the local geometry to the global topology of the network. The value of γ used in their studies is 0.5, which comes from simple hydrodynamic considerations. In these works, the authors suggest that minimization of the Hamiltonian should be the leading principle for finding stable river networks (Optimal Channel networks or OCN), with a remarkable success: the observed excess probability of drainage areas, after minimization, is a decreasing power-law with exponent 0.45, in agreement with real river networks.
Following the seminal work of , I adapted geomorphologists' hypotheses. I made here two basic assumptions: i) Landscapes (land cover spatial distribution) can self-organize at large scales when patches interact. Such interactions could have various natures depending on the landscape type. These interactions could either be dispersion and colonization processes or crop rotation systems and regional management for instances ; ii) such landscapes are supposed to globally optimize, through a principle of minimum "cost" expenditure, the mosaic as a whole. This cost, similar to energy in physical systems, refers to the total (natural and/or human) cost of spatial arrangements, which are supposed to be cheaper (less constrained) at large scales. This second principle superimposes to local processes such as dispersion at local scale, competition, etc. Here, we note differences with the hydrographical analogy, dealing with a purely physical system [26,6]: It is only postulated that the optimizing principle is global: it models how species or land covers tended to fill space, and no additional local optimization process is needed
In order to test this hypothesis, I proposed to model Optimal Patchy Landscapes (OPL) using this NLM. One way to proceed is to define a potential force field, which may be minimized in case of self-organized patchy landscapes. Again, this potential can be seen as a summarizing function of the landscape and is hereafter called the "Hamiltonian" function. The hydrographical network configuration (i.e. geometrical) self-similarity looks very close to the composition (i.e. attributive) self-similarity sometimes observed within patchy landscapes [9,10]. Patchy forested landscapes are best characterized using the tree-density spatial distribution of every patch . We found that the tree-density differences between pair-patches better highlights patch interactions in space, while the cumulative tree-density difference distribution between pair-patches and its probability to exceed a defined threshold should be able to capture the scaling features of the landscape [6,26].
The studied site is located in the southeast of France, on the western slope of the Mont Ventou and the elevation is between ~900 and 1600 m. It consists in a forest sector managed by the French national forest service (ONF), with a mixture of species (Abies alba, P. nigra nigricans, Pinus sylvestris, P. uncinata and Fagus silvatica). Most of the real site is subject to active forestry management through thinning and natural or artificial regeneration, but further data are reduced only to natural evolutions. With a total surface area equal to 11.4 km2, the site is made of 220 patches. The boundary coordinates of each patch i designed for the management purposes are delimited areas of roughly uniform stand densities di. The density is defined as the number of stems per hectare. As exhaustive data were not available over such a wide area and over such a long time period, I used the forest dynamics model CAPSIS  to simulate Fir trees and the four other species undergoing natural dynamics in a spatially explicit scheme. About 30 growth and yield models [28,29], mainly including stand-level models, diameterclass models and individual-based models have been integrated into the CAPSIS platform since 1999. The "VentouG" model, developed into the CAPSIS platform, performed a typical simulation hereafter called the reference (Figure 1a). VentouG is a semi-spatialized model (close to "gap-models") that mixes diameter class growth models in a spatialized plot to study the natural evolution of heterogeneous forests evolving towards mixed stands [28,15].
Figure 1: (a) Evolution of spatial distributions of Fir tree densities (stems/ha) in the Mont Ventoux landscape. This multispecies site of about 1100 ha is called the reference landscape (only the Fir specie is displayed). Only year 2000 (observed), years 2050 and 2100 (simulated with the CAPSIS based process model) are shown, with the same coordinates and the same logarithmic colour scale. (b) The same landscape with hypothetical random tree densities is plotted for comparison and for visualizing patch boundaries.
VentouG simulations played the role of observations. This study analyzed the five species dynamics to illustrate how generic the approach is. In addition, I chose to focus on the Fir species (Abies alba) for a more detailed analysis. The dynamical distribution of Fir area is interesting: it is subject to colonization by seeds that come from important residual stands and is sensitive to elevation. We simulated here the evolving tree-densities of every patch according to processes of: i) colonization (exponential kernel dispersion, without external sources to the landscape); ii) regeneration (successive stage growth, with regeneration heights comprised between 0.3 and 1.3 m); and iii) competition (mortality induced by local surroundings) between the stems of the five species [13,15]. All species mentioned above are simulated with the same processes, but with species-dependent calibration: they disperse (at a maximum distance between 1 to 2 km ), they have different fructification power and they compete with intra- as well as interspecies interactions. These ecological processes were explicitly applied on every stem, but no additional organizing processes such as large scale optimization were simulated. In general, simulations have shown a high sensitivity to dispersal, fructification and survival, especially at the seedling stage [28,29]. The run covered a period of 100 years with a time-step of 5 years and started with the real composition observed in year 2000. A fully random patchy landscape at the same site is shown for comparison (Figure 1b).
The landscape hamiltonian
I intended to characterize the landscape observations in terms of Hamiltonian functions. In order to detect a self-similar property of the reference landscape s, I first defined the Hamiltonian H function best characterizing the studied landscape . H is defined using the exceedence probability PE of tree-density differences Δd between pairs of patches for the whole landscape. In case of self-similar landscape, this probability sets:
where τ-1 is the fractal dimension (the slope in loglog plot) of the PE probability curve. Working here with exceedence probability allowed to characterize possible scaling behaviors as other cumulative functions would. The tree-density difference intended to capture ecological patch interactions.
In the empirical studies on river networks, minimization of an Hamiltonian given by equation (4), with γ = 0.5, showed that the exceedence probability distribution for areas obeyed a power law:with τ-1 approximately equal to γ. By analogy with these works, we propose to define our landscape Hamiltonian as:
where Hij is the interaction strength (by analogy with physics, the potential or local force) associated to the tree-density difference Δdij between patches i and j. For simplicity, this study supposed a direct proportionality (equality) between the left and the right part of both equations.
Finally, the total Hamiltonian value summarizing the whole landscape s leads to:
Here again, the HT value is by analogy with physics, the energy or the total "cost" of the landscape. Landscape energy has no clear meaning, but such function is a useful way to synthesize (to gather) all processes occurring in a landscape and leading to the tree-density pattern.
Practically, the H function used here is built from the fit made over fifty percentiles of the exceedence probability of observed tree-density differences. All tree-density differences between pairs of patches at every time step (2000 to 2100 with five a year step) for the reference landscape were computed separately. When computed separately at each time step, the H function shape remained relatively stable along time, but the information was too sparse in this case to cover many magnitude orders. Tree-density differences were then merged in order to extract an averaged H function that expressed the statistical behavior of the dynamic landscape. When H exhibits a self-similar shape, its associated γ exponent given by equation 1 and 2 was computed by fitting the histogram made of the 510510 observed density differences. The fit was obtained by a mean-squares procedure over the 20 to 200 main bins of the histogram. The final fit given by fifty bins was conserved, the others helping to compute the slope's uncertainty. There were too ways for using the computed averaged H function: i) it was possible to compute the HT value at every time step of the reference simulation and to characterize its regimes in time. For this purpose, I normalized all HT values by the biggest value found in the species data; ii) it is possible to use this Hamiltonian to model a new landscape with similar statistical properties.
Starting from the landscape Hamiltonian, it is possible to analyze the observed landscape or to simulate new (optimized) landscapes. It has been demonstrated that it is easy to simulate landscapes using the Hamiltonian function by minimizing its total sum . The minimization stage could be done with the T=0 Gibbs algorithm (i.e. Gibbs process by means of the Metropolis-Hastings algorithm) developed for physical systems studies [20,21,31] and recently applied in ecology too [13,32]. In an adaptation to patchy landscape interactions, it is expressed as follows: a) compute the total landscape value HT (st) at step t using equation (7); b) randomly choose a patch within the landscape and randomly change its tree-density di (bounded or not). Which consequently modifies all the tree-density differences Δdij of pair-patches ij in which this selected patch i is involved; c) compute the new landscape value HT (st +1 ) ; d) and finally compare the latter to the previous HT value. HT (st +1 ) ≤ HT (st) , then keep this new configuration, otherwise keep the previous configuration, and reiterate at stage b); e) accept the convergence after n successless iterations (i.e. without any HT decrease).
It has been demonstrated that this algorithm, also used in Bayesian framework, rapidly converges. No simulated-annealing has been performed in this preliminary study although it is possible and sometimes recommended.
This kind of optimization is a neutral model (NLM) in the sense that is requires very little information, namely the H function and the initial conditions (giving an averaged parameter value), and no ecological processes. In case of a self-similar landscape H possesses only one parameter: the γ exponent. Following hydrographical studies, this exponent connects local processes with an overall topological structure of the landscape (i.e. patch neighborhoods) such that: if 0 < γ << 1, we get a structured aggregation for which fine scales are dominant; if γ >> 1, the mosaic shows patch densities for which broad scale arrangements (ecological links) become dominant.
Starting from the observations (in fact, the CAPSIS mechanistic simulations), it is possible to compute the exceedence probability distribution, then helping to define the landscape Hamiltonian, finally allowing to analyze in turn the original observed landscape (arrow A) or to compute new optimized simulations (arrow B). The flow chart describes the successive stages of the calculation (Figure 2). A landscape analysis (A) consists in computing the evolution of the total HT value of the observed landscape along to its dynamics (here over one century). Simulation study (B) was illustrated here by starting with the observed landscape of Mont Ventoux at year 2000, leading to a fully self-organized landscape at a posterior date. Simulation was made for the five species simultaneously and fir dynamics is analyzed in more detail. The dominant self-similar property was achieved after minimizing the HT value of the evolving landscape. Iterations of the simulation were not supposed to correspond to any real time step. This NLM did not take into account patch elevations.
The exceedence probability of tree-density differences between pairs of patches of the Ventoux landscape as a function of treedensity differences exhibited a strong self-similarity (Figure 3). This directly comes from the ecological processes that we modeled (there is no minimization of H at this stage). Furthermore, the self-similarity covered up to nine (five in log10) magnitude orders for all of the studied species (not shown). This extended scale range was accessible as a large landscape (220 patches, up to 12 km) and all time steps were merged for the analysis. It leaded to slopes equal to – γ = – 0.66; – 0.80; – 0.73; – 0.77; – 0.66 ± 0.02 for Abies alba, P. nigra nigricans, Pinus sylvestris, P. uncinata and Fagus silvatica respectively. Determination coefficients of these fits are roughly equal to r2 = 0.97 ± 0.013, all significant (p <10-5). Each uncertainty associated to the slope was the standard error computed on ten other fits with various numbers of histogram bins for each species. Figure 2a showed the result of this fit for the fifty bin case (Figure 3a circles) on an original diagram gathering 5000 bins (dots). As a reference, the same calculation was performed on the fully random mosaic: it did not show any self-similar behavior (Figure 3 inset).
Figure 3: (a) Exceedence probability (noted PE) of the tree-density differences against tree-density differences (Δd) presented in a loglog plot. 5000 dots were used in the original plot, while circles correspond to the fit with 50 regular bins. The same plot for the random landscape of figure 1b is shown for comparison, while its shape is highlighted with a linear fit (line, inset). (b) Hamiltonian function (noted H) extracted from the exceedence probability of the reference landscape. Here, dots correspond to open circles in figure 2. The slope γ of this loglog plot is highlighted with a mean-square linear fit.
Simply inverting this exceedence probability defined the Hamiltonian H function, in order to favor frequent tree-density differences (Figure 3b). The self-similar behavior detected here indeed indicated that low Fir density differences were more frequent in this landscape. H conserved the exceedence probability slope γ, while local H values (i.e. associated to each pair-patch) are dimensionless. Taking into account very high tree-density differences (extreme points at left in Figure 3a) or not did not change the fit. The four other species also exhibited strong self-similar behavior, with Hamiltonian slopes mentioned above and slightly higher than the one of fir species (Figure 4a).
Figure 4: Hamiltonian function computed from the exceedence probability of the four other tree species involved into the reference landscape: P. nigra nigricans (squares, dashed line); Pinus sylvestris (dots, doted-dashed line); P. uncinata (stars, doted line); Fagus silvatica (circles, plain line). Their respective slopes of this loglog curve are highlighted with a mean-square linear fit.
Landscape analysis and simulation
The evolution of the total HT value of Fir landscape was quite variable over one century (Figure 5A). Higher HT values corresponded to lower tree-density differences and thus to a more homogeneous landscape. There was a clear regime shift at year 2065, after which the fir Ventoux landscape tended to monotonously behave towards a more and more homogeneous landscape. Other species show more monotonous regimes, except maybe for Fagus silvatica starting to be limited by available elevations and competing species at the end of the simulation (Figure 5A). Such dynamics are also sensitive to past history that can vary from one species to another.
Figure 5: Final computation steps described by the Fig. 2 flow chart. (A) Normalized total HT values (i.e. the sum of Hamiltonian values over the whole landscape) against the time steps (in years) of the reference landscape. There is a regime shift around year 2065. (B) Spatial distributions of Fir tree densities (stems/ha) of the Mont Ventoux landscape simulated by the Neutral Landscape Model on the basis of the Hamiltonian function of figure 3a. The inset (b) shows the resulting Hamiltonian function after simulation (similar to figure 2b) in loglog plot.
It was possible to model a landscape self-organization by using the previous H function deduced from the observed fir Ventoux landscape. The total HT minimization procedure led to a new patchy landscape (Figure 5B). The final landscape was perfectly self-organized after only 300 iterations (calculations took a few seconds on a standard PC). The new Hamiltonian function extracted from the fir tree-density difference probabilities (equation (2)) was rigorously self-similar (Figure 5B inset) and exhibited the same exponent as the one used in the Hamiltonian. Fir species filled almost all landscape patches as the NLM did not apply any elevation constraint, contrarily to the reference simulation.
Optimized patchy landscapes
The five landscapes studied here appeared to be self-similar over more than nine (five in log10) magnitude orders (Figure 4a). Data came from a spatially explicit process-based simulation of the Mount Ventoux site over a hundred years , thus extending data and allowing a complete control on the landscape generation mechanisms. The landscape evolution studied combined local processes such as colonization, competition and regeneration [15,30], without any additional organizing process. The landscape dynamic was characterised with the tree-densities of the 220 patches, and then the tree-density differences between all of the pair-patches to capture spatial interactions. The five mono-species landscapes exhibited selfsimilar behaviours, thus illustrating how generic the analysis might be. Their self-similar coefficients ranged from γ = 0.66 to 0.8, mainly due to dispersal and seedling stage survival differences between species. The (simple) way we used to estimate these coefficient values is not important, as we did not try to infer processes underlying observed patterns. This main result strongly suggested that the combination of ecological processes behaves as a self organized landscape.
This observation fitted some others made in previous works [5,12]. For example, rural landscapes often hide a self-organization such as selfsimilar heterogeneity observed over temperate and tropical landscapes . Furthermore, self-organization was interpreted in the present work as global optimization, as we minimized the landscape's Hamiltonian. Indeed, starting from the observed tree-density distribution in space, it was possible to reconstruct a self-similar landscape using a potential force field summarising ecological interactions (Figure 5B). The landscape's self-similar behaviour was the result of a global minimization of the sum of local landscape interactions, as was done in previous studies . The self-similar behavior detected in the forested landscape was based on the observation that low tree density differences were more frequent, which fitted well with the concept of preferential attachment already mentioned to advocate optimizing patterns .
In addition to this interpretation of landscape generation, the optimization process opens a way to an original quantification of landscape features. This study showed how to characterize possible landscape regimes in time on the basis of the integrative Hamiltonian function (Figure 4a). The reference fir landscape studied here clearly showed a regime shift around year 2065 that was not explicitly modelled (Figure 5A). This shift could be the response of the colonization process to the spatial structure encountered in the landscape. The patch's surroundings indeed favoured space exploration in a first stage, and tree-density smoothing (homogenization) in a second stage when the patches were already reached. The other species' dynamics were rather different, more monotonous and smoother, due to different past histories and different calibrated processes (dispersion and mortality).
From a conceptual point of view, one needs to discuss the apparently remarkable fact that optimizing landscapes with respect to an Hamiltonian with parameter γ, always leads to a power-law for the excess probability distribution for density differences with parameter τ −1 = γ . The T=0 Metropolis algorithm we used drives the system into a local minimum of the Hamiltonian, so that there may not be a large departure from the initial condition (already showing an exceedence probability distribution with exponent γ). Yet, the global minimum of this Hamiltonian leads to the trivial and uniform distribution. Detailed studies of optimized river networks [33,34], analytically found the global minimum for the Hamiltonian given by equation (4): for all γ between 0 and 1, and different from 0 and 1, they found τ-1 = 0.05. This was confirmed by simulated annealing simulations, where the global minimum can effectively be found . However, by using T=0 Metropolis algorithm with this Hamiltonian and γ =0.5, one visits the local minima with τ-1 = 0.45 This value is in agreement with those measured in observed river networks. Such study and others led to the concept of “feasible optimality” according to which nature is “unable” to reach the true ground state (the global optimum) when complexity is involved. This could also be the case in patchy landscapes, given that one first checks how one can relate the Hamiltonian to ecological processes.
The present NLM was in agreement with the parsimonious constraint of neutral models [13,14]. The strength of this parsimonious argument could also be a pitfall as the core of the method resides solely in the Hamiltonian H definition. If H is badly defined, the NLM fails or even shows irrelevant landscape features. H is obviously highly dependent on the ecological parameter chosen to describe the landscape. In this forestry study, I focused on tree-density differences closely connected to ecological processes driving landscape dynamics. In other studies, cluster sizes of the vegetation cover were also used [9,11]. Agricultural landscapes for example are strongly dependent on the crop rotation system of farmers and would need other descriptive parameters. Finally, notice that NLM iterations were not directly linked to real time steps, even though this could be investigated if necessary . Similarly, in a first stage here, I supposed a stationary over time H function, associated to an optimisation procedure towards steady-state or equilibrium (Figure 5B). This approximation could (and probably should) obviously be modified in a second stage of the description work of this landscape. This may be done either by forcing (driving) such Hamiltonian-based model or by defining dynamical rules of a new type.
Present computations were also underlying spatial assumptions. First, this work considered in a conservative approach all the existing pair-patches of the landscape (Figure 3), while different choices could be justified as well. For example, H could sum local interactions only coming from neighboring patches, with various (Delaunay or distance based) neighborhoods, as seen in other works [9,25]. Such choice depends on ecological processes advocated for landscape generation. Indeed I argue that the forest's pattern under study involves different processes that could have long distance actions due to seed and pollen dispersion processes. Secondly, the patch's definition may have played a role in such a landscape analysis. This weakness is common to every study in this field: studies focusing on vegetation cluster sizes also have to identify thresholds in order to define vegetation patches [9,12]. Patches were here defined by the sylvicultural management of Mont Ventoux landscape. I checked that slightly modified patch geometry (shapes and inter-distances) as well as topology (neighborhoods) did not change the present results. This is due to the fact that all pairpatches were taken into account for computations. Hence, the present method dedicated to Optimal Patchy Landscapes (OPL) seemed very robust in this sense. Thirdly, it is possible, and even probable, that additional processes influence the ecological pattern described here. We may modify or weight the Hamiltonian definition in order to take into account for the landscape topography, the soils, etc.
A dynamic landscape theory
There is presently little theory interpreting (not necessarily forested) landscape patterns and dynamics. Landscape ecology has been working on that point for a long time and is the origin of many powerful concepts such as heterogeneity . Yet, there is still a lack of theoretical framework aiming at mechanistically interpreting landscape generating processes. No patchy landscape, the most common structure in rural areas, has been to my knowledge described with mathematical equations yet. In raster mode some composition analyses (i.e. averaged land cover percentages) were explored, while most of the landscape's complexity resides in its configuration (its patch arrangements) and its dynamics . Note that composition and configuration properties are highly dependent.
Recent works intended to explore landscape properties, but they often focused on "continuous" (raster) landscapes and did not succeed in giving clear explanations of observed scaling behaviour [9,12]. In this study, I investigated the same question with an ecologically based Hamiltonian function and its optimisation, thus following some pioneer works done on similar soap foams and tissue cells [21,23]. These authors exported a powerful model comprising almost all cells/units interactions into a Hamiltonian function and optional constraints to describe the evolution of the system. The Hamiltonian function is a way to summarise landscape properties within a simple equation (equation 2). By construction, the H function is intimately linked to various ecological processes such as colonization and competition involved in landscape generation. In this study, a higher HT value corresponded to lower tree-density differences and to more homogeneous landscapes. Such homogenization could happen with "older" landscapes, stabilizing through a dominant colonisation process: colonizing species progressively fill the space and saturate it.
To put in a nutshell, this Hamiltonian function appears to be a fruitful theoretical framework to describe various landscapes. However it is still phenomenological (constructed by analogy with other problems) and its derivation from ecological processes is not known. In physical systems, this exponent connects local processes with the overall topological structure of the landscape (i.e. close and far away neighborhoods). Moreover, it is probable that the minimization procedure using a T=0 Metropolis algorithm only makes the landscape relaxed into a local or "metastable" state with respect to the Hamiltonian. It would be interesting to apply the present approach to observed landscapes: we could see how often they are self-organized in terms of various ecological parameters and what the observed γ exponent values are. Yet, it will certainly be difficult to gather the large amount of spatial and temporal data needed for this purpose. Such an analytical tool also gave the opportunity to compare landscapes (with a similar H function) or to study in depth complex landscape dynamics [3,35]. Following works on hydrographical systems , this theoretical framework potentially opens the way to a more complete dynamic landscape theory.
I warmly thank Fabien Campillo and Estelle Pitard for their comments on earlier drafts of this paper.
Make the best use of Scientific Research and information from our 700 + peer reviewed, Open Access Journals