Influence of technological progress and renewability on the sustainability of ecosystem engineers populations

Overpopulation and environmental degradation due to inadequate resource-use are outcomes of human's ecosystem engineering that has profoundly modified the world's landscape. Despite the age-old concern that unchecked population and economic growth may be unsustainable, the prospect of societal collapse remains contentious today. Contrasting with the usual approach to modeling human-nature interactions, which are based on the Lotka-Volterra predator-prey model with humans as the predators and nature as the prey, here we address this issue using a discrete-time population dynamics model of ecosystem engineers. The growth of the population of engineers is modeled by the Beverton-Holt equation with a density-dependent carrying capacity that is proportional to the number of usable habitats. These habitats (e.g., farms) are the products of the work of the individuals on the virgin habitats (e.g., native forests), hence the denomination engineers of ecosystems to those agents. The human-made habitats decay into degraded habitats, which eventually regenerate into virgin habitats. For slow regeneration resources, we find that the dynamics is dominated by cycles of prosperity and collapse, in which the population reaches vanishing small densities. However, increase of the efficiency of the engineers to explore the resources eliminates the dangerous cyclical patterns of feast and famine and leads to a stable equilibrium that balances population growth and resource availability. This finding supports the viewpoint of growth optimists that technological progress may avoid collapse.


I. INTRODUCTION
The gloomy consequences of rapid population growth and environmental degradation have been a major public concern, probably first expressed explicitly in the work of Thomas Malthus [1], who believed that informed public policies could mitigate or even avoid altogether the happening of those circumstances. Following this line, in the 1970s the Club of Rome commissioned a major enterprise aiming at modeling the interactions between the earth and human systems using massive computer simulations capable of producing detailed near-future scenarios [2,3]. Although efforts in that direction are still ongoing [4], this approach is somewhat abstruse as the complexity of the simulated models rivals that of the realworld, thus making it difficult to untangle the effects of the model parameters. An alternative approach to address such concerns, which we adopt here, is the study of minimal models that link the population dynamics and the resource dynamics. Works in this line typically build on Lotka-Volterra predator-prey models [5], where man is seen as the predator and the resource base as the prey [6,7].
Yet the predator-prey approach to human-nature dynamics misses the important fact that humans do not primarily prey on nature but modify it. In fact, humans are the ultimate ecosystem engineers, who have shaped the world into a more hospitable place for themselves, most often with unlooked-for long term effects [8]. One such effect is the collapse of human societies likely caused by habitat destruction and overpopulation [9,10], hence the public concern on this matter. Accordingly, here we argue that the baseline population dynamics model that is more suitable to address human-nature interactions is the ecosystem engineers model [11,12], rather than the predator-prey model. The feature of the population dynamics of ecosystem engineers that makes this approach well suited to model those interactions is that the growth of the population is determined by the availability of usable habitats, which in turn are created by the engineers through the modification, and consequent destruction, of virgin habitats. From a mathematical perspective, this feedback results in a density-dependent carrying capacity that produces a very rich dynamic behavior both in the continuous-time [12] and the discrete-time formulations [13].
More pointedly, the Gurney and Lawton model for the population dynamics of ecosystem engineers describes the time evolution of the density of engineers and of the habitats, which can exist in three different forms: virgin, usable and degraded. The transition of a habitat from the virgin form to the usable form occurs due to the work of the engineers. The usable habitats then degrade and eventually recover to become virgin habitats again. Virgin and degraded habitats are unsuitable for the growth of the population. To expedite the numerical analysis, here we use a discrete-time approach where the logistic growth equation of the original model is replaced by the Beverton-Holt equation [14]. The weak nonlinearity of this equation prevents the large fluctuations on the population size produced by the Ricker equation [15], which may result in a spurious chaotic dynamics [13].
We find that the survival of the population depends on two independent factors, viz., the competence of the engineers to transform natural resource base (e.g., native forests) into useful goods (e.g., farms) and the population growth rate per generation. Survival is guaranteed provided that the product of these two factors is greater than the decay rate per generation of the usable habitats. However, survival is not without adversities in the case the regeneration rate per generation of the degraded habitats is too slow. In this case, the dynamics exhibits cycles of prosperity and collapse, in which the population reaches vanishing small densities. Extinction does not occur because the fixed-point associated to it is unstable. The regime of cycles is favored by large population growth rates.
Most interestingly, we find that increase of the efficiency of the engineers to explore the resource base, thus modeling a high-tech society scenario where a few individuals can extract and produce the goods needed for survival, eliminates the dangerous cyclical patterns of feast and famine and leads to a smooth adjustment between population and resources. This finding supports the viewpoint of growth optimists that technological progress may avoid collapse.
In the case the resources are non-renewable, the longterm outcome is a doomsday scenario in which the ecosystem is reduced to degraded habitats only. Since the population is not immediately affected by the disappearance of the virgin habitats, we can interpret the usable habitats as a kind of wealth, i.e., surpluses that humans accumulate and that may extend their survival well beyond the point where the natural resources are exhausted. This observation offers a connection between ecosystem engineers models and economic models of human-nature interactions [7].
The rest of the paper is organized as follows. In Section II, we introduce the discrete time version of Gurney and Lawton model of ecosystem engineers and derive the recursion equations for the density of engineers as well as for the fractions of the three forms of habitats. In that section we also derive and begin the analysis of the local stability of the fixed point solutions of those equations. Section III presents the results of the numerical study of the eigenvalues of the Jacobian matrices that determine the local stability of the fixed points as well as the results of the numerical iteration of the recursion equations. Finally, Section IV is reserved for our concluding remarks.

II. DISCRETE-TIME POPULATION DYNAMICS
The population dynamics of ecosystem engineering was modeled originally by Gurney and Lawton in the 1990s using a continuous-time model [12], while discretetime versions were considered only much more recently [13,16]. We recall that the main advantage of discretetime formulations is that they can easily be extended to incorporate a spatial dependence of the dynamic variables [17], as well as of the model parameters, as neatly demonstrated in the classic studies of the spatial dynamics of host-parasitoid systems [18,19] (see [20,21] for more recent contributions). However, discrete-time nonlinear systems usually exhibit chaotic behavior, which is somewhat problematic from a biological perspective [22] and are not observed in their continuous-time counterparts [23]. Here we show that this drawback can be circumvented by a suitable choice of the equation that determines the growth of the population. In particular, whereas in previous studies the logistic equation of the continuous-time Gurney and Lawton model was replaced by the highly nonlinear Ricker equation [15], here we use the growth equation of Beverton and Holt [14], whose much weaker nonlinearity greatly limits the betweengenerations fluctuations and produces a stauncher representation of the continuous-time model.
We assume that the population of engineers at generation t is composed of E t individuals and that the carrying capacity of the environment depends on the number of units of usable habitats available at generation t, which we denote by H t . Thus, using the Beverton-Holt equation we write the expected number of engineers at generation t + 1 as where r > 1 is the intrinsic growth rate of the population of engineers and the time-dependent carrying capacity is (r − 1) H t . This model describes a system of ecosystem engineers due to the assumption that the usable habitats are produced by engineers working on virgin habitats [12]. More pointedly, assuming that there are V t units of virgin habitats at generation t, a fraction C (E t ) V t of them will be transformed in usable habitats by the engineers at the next generation t + 1. Here C (E t ) is any function that satisfies 0 ≤ C (E t ) ≤ 1 for all E t and C (0) = 0. In other words, if there are no engineers at generation t, then virgin habitats cannot be transformed in usable ones. Hence C (E t ) measures the efficiency of the engineer population to build usable habitats from the raw materials provided by the virgin habitats. Usable habitats decay into degraded habitats that are useless to the engineers, in the sense that they are inhabitable and lack the raw materials needed to build usable habitats. Denoting by δ ∈ [0, 1] the probability that a unit of usable habitat decays and becomes a unit of degraded habitat in one generation, we write the expected number of units of usable habitats at generation t + 1 as In this paper we assume that the decay probability δ is a density independent parameter. We recall that in the Gurney and Lawton model the resources are represented by the virgin habitats, whose probability of change into usable habitats is density dependent, i.e. C = C (E t ). For example, in an island scenario, the virgin habitats can be thought of as the native forests whereas the usable habitats are the lands cleared for crops, whose degradation, due mainly to erosion and soil depletion of nutrients, is more suitably modeled by a constant decay probability δ, rather than by a density-dependent one. (i.e., δ = δ (E t )). However, in an economic context where the habitats represent, for instance, oil wells, the decay probability should be modeled by a density-dependent parameter, of course. To take into account the possibility that the degraded habitats will eventually recover and become virgin habitats again, we introduce the parameter ρ ∈ [0, 1] that gives the probability that one unit of degraded habitat becomes one unit of virgin habitat in one generation. Thus, denoting the units of degraded habitats at generation t by D t we write Thus ρ is the regeneration rate per generation that measures the renewability of the resource base. Finally, recalling that virgin habitats are destroyed by the action of the engineers and created by the restoration of the degraded habitats, we write the recursion equation for the expected number of units of virgin habitats as which completes the description of the feedback interactions between habitats and engineers.
Since habitats can only be transformed from one form to another, the total number of units of habitats Accordingly, we define also the density of engineers e t = E t /T that, differently from the habitat fractions, may take on values greater than 1. In terms of these intensive quantities, the above recursion equations are rewritten as where we have used To complete the model we must specify the densitydependent probability c (e t ), which measures the engineers' efficiency to transform the virgin habitats into usable ones. The function c (e t ) incorporates the cooperation strategies used to build structures (e.g., beaver dams, termite mounds and anthills) to respond to external and internal threats to the population survival [24,25]. For humans, this function incorporates the beneficial effects (from their perspective) of the technological advancements that enable a more efficient harvesting of natural resources [26]. Here we consider the continuous function where α ∈ [0, 1] is the efficiency parameter that measures the competence of the engineers in transforming natural resources into useful goods. The range of variation of α guarantees that c (e t ) ≤ 1 regardless of the value of e t . For α 1 we have a scenario where the efficient exploration of the virgin habitats requires a high density of engineers. In this low-tech scenario, we expect to observe an Allee effect, in which the population quickly goes extinct below a critical density [27]. We note that when e t > 1 all the available virgin habitats are transformed into usable ones in just a single generation. This arbitrary threshold value does not produce any significant effect on our results since we find that e t < 1 at all times in the parameter regions of interest, i.e., in the regions close to the transitions between the fixed point attractors and the limit cycles.

A. Fixed-point solutions
In this section we study the fixed-point solutions of the system of recursion equations (5)-(7) that are obtained by setting e t+1 = e t = e * , h t+1 = h t = h * and v t+1 = v t = v * . Next we consider separately the two cases, e * = 0 (zero-engineers fixed point) and e * > 0 (finite-engineers fixed point). The local stability of the fixed point solutions are determined by the eigenvalues of the Jacobian matrix of the system of recursion equations (5), (6), (7), where x t = e t /h t and c (e t ) = dc (e t ) /de t .

Zero-engineers fixed point
In the absence of the engineers (i.e., e * = 0) all usable habitats will degrade and then recover to virgin habitats (provided that ρ > 0), so we have h * = 0 and v * = 1. To circumvent the indetermination of the ratio e * /h * , we write an equation for x t dividing eq. (5) by eq. (6), whereĉ (e t ) = c (e t ) /e t . Thus the fixed point x t+1 = x t = x * > 0 is given by the positive solution of the quadratic equation where we have usedĉ (0) = lim e * →0 c (e * ) /e * = c (0) = α. The quadratic equation (11) has two real roots, one positive and the other negative, for all values of the model parameters. In particular, for δ = α (r − 1) the solution is x * = r − 1 and we can easily show that x * > r − 1 for δ > α (r − 1) and x * < r − 1 for δ < α (r − 1), a result that will proven helpful to verify the local stability of this fixed point.
In fact, the local stability of the zero-engineers fixed point is determined by requiring that the three eigenvalues of the Jacobian matrix have norms less than 1. The eigenvalues of the Jacobian matrix (12) are λ 1 = r/ (1 + x * ), λ 2 = r/ (1 + x * ) 2 − αx * and λ 3 = 1 − ρ. The condition | λ 1 |< 1 is satisfied provided that x * > r − 1, which is guaranteed for δ > α (r − 1), whereas the condition | λ 3 |< 1 is always satisfied since ρ < 1. The proof that the condition | λ 2 |< 1 is always satisfied whenever λ 1 < 1 is very simple. Since α > 0 and x * > 0 we can write Next we need to show that λ 2 > −1. To do so we use eq. (11) to write αx * = λ 1 − 1 + δ so that We note that this stability analysis describes the convergence (or divergence) of the system of recursion equations for e t , h t and v t , eqs. (5), (6) and (7), respectively, to the fixed point e * = h * = 0 and v * = 1 and that equations (10) and (11) were used only to specify the ratio lim t→∞ e t /h t that is required to evaluate the Jacobian (9). We mention this because use of recursion equations for e t , x t and v t instead, yields a different Jacobian for the zero-engineers fixed point, which exhibits the same eigenvalues λ 1 and λ 3 as the Jacobian (12), but a different eigenvalue λ 2 . However, the local stability of the fixed point is again determined solely by the norm of the common eigenvalue λ 1 so this issue has no effect in our findings regarding the zero-engineers fixed point.
In sum, the condition for the stability of the zeroengineers fixed point is δ < α (r − 1), which, as expected, does not depend on the regeneration rate ρ since there is no shortage of virgin habitats (i.e., v t ≈ 1) in the vicinity of this fixed point.

Finite-engineers fixed point
Since e * > 0 we have h * = e * / (r − 1) and v * = 1 − e * γ/ (r − 1) with e * given by the solution of the equation where we have introduced the notation γ = (1 + δ/ρ). Let us assume that e * < 1 so we can use eq. (8) to rewrite eq. (13) as the quadratic equation The local stability of the two solutions of this equation is determined by the the Jacobian matrix (9) that is rewrit-ten as (15) The main difficulty of the analysis of the eigenvalues of this Jacobian matrix is that two eigenvalues are complex, leading to the enthralling dynamic behavior exhibited by the recursion equations (5), (6), (7) with the prescription (8). However, this problem can easily be solved numerically as follows. First we determine the critical points of the characteristic cubic polynomial, which allows us to obtain one real root of the characteristic equation with high numerical accuracy using the bisection method. (We recall that the critical points of a function are those values of the variable where the slope of the function is zero.) Then this root is factored out and the remaining roots of the characteristic equation are found by solving a quadratic equation. Actually, in the case of complex roots, we do no need to solve this quadratic equation since the local stability is determined by the norm of those roots that is given by the square root of the ratio between the constant and the quadratic coefficient. The fixed point solution is locally stable provided that the norms of all three roots of the characteristic polynomial are less than 1. We find that the smaller root of the fixed-point equation (14) is always locally unstable. Now we consider the case e * > 1. Equation (13) yields that holds for δ < (r − 2) / (1 + 1/ρ). In this case, the Jacobian matrix (9) becomes Hence λ 1 = 1/r < 1 is always a root of the characteristic polynomial and the other two roots are given by the quadratic equation g (λ) = λ 2 −λ (1 − ρ − δ)+ρδ = 0. In the case of complex roots, we have | λ | 2 = ρδ < 1. In the case of real roots, we use that g(1) > 0 and g(−1) > 0 and that the critical point of g(λ) is between −1 and 1 to guarantee that the absolute value of those roots are less than 1. Hence the fixed point (16) is always locally stable.

III. RESULTS
The local stability of the fixed points e * = 0 and 0 < e * < 1 is determined by the requirement that the norms of the eigenvalues of the Jacobian matrices (12) and (15), respectively, are less than 1. Accordingly, Fig.  1 shows the regions of stability of those fixed points. There is a large region in the space of parameters (δ, α)  where the fixed points are unstable and the solutions of the recursion equations exhibit periodic behavior. The drop-shaped region delimits the regions of stability of the finite-engineers fixed point. To better understand the results shown in Fig. 1, we first examine whether eq. (14) admits a physical solution with a vanishingly small density of engineers. In fact, for e * 1 we have , (18) so that the density of engineers vanishes continuously as the line α = δ/ (r − 1) is approached from the region of large α, provided the denominator in eq. (18) is positive.
i.e., αγ − (1 − α) (r − 1) > 0. In Fig. 1 we indicate with a red bullet the point δ = δ c ≈ 0.202 and α = α c ≈ 0.067 at which both the numerator and the denominator in eq. (18) vanish, meaning that e * = 0 is a double root of the quadratic equation (14). Thus, there is a continuous transition between the finite and the zero-engineers fixed point at α = δ/ (r − 1) for δ > δ c . Next, we focus on the fixed-point solutions at the transition line α = δ/ (r − 1). In this case, eq. (14) has the solution e * = 0 and provided that α > (r − 1) / ((r − 1) − γ), which guarantees that e * < 1. Otherwise, we have e * = (r − 1) / (δ + γ) > 1 (see Section II A 2). Figure 2 shows the finite-engineers fixed points at the transition, where the stable and unstable segments are flagged with different colors. The figure reveals that the fixed point e * > 0 is stable for small δ, as expected, and that it is also stable close to the point (δ c , α c ). This indicates that this point does not belong to the drop-shaped curve shown in Fig. 1, which delimits the regions of stability of the finiteengineers fixed point. This facet will become much more evident when we consider larger values of the recovery probability ρ.
Another interesting feature of Fig. 1 is the appearance of a region of bistability of the two fixed points for small values of δ and α < δ/3. The dependence on the initial density of engineers, as well as on the initial habitats fractions, that characterizes the bistability regime results in an Allee effect where the environment cannot sustain a low density of engineers, as the change of the virgin habitats into usable ones can only be achieved by highdensity populations. Figure 3 illustrates the time dependence of the density of engineers e t and of the fraction of virgins habitats v t for three representative regions of Fig. 1, where the zero-engineers fixed point (e * = 0) is unstable, i.e, for δ < α (r − 1), so the results do not depend on the initial conditions. The time dependence of the fraction of usable habitats h t is similar to that of e t . The recursion equations (5)-(7) were iterated using quadruple precision. The oscillatory convergence to the finite-engineers fixed point (e * > 0) is the effect of the complex eigenvalues of the Jacobian matrix (15). These oscillations may lead to extremely low values of the density engineers (e.g., e t on the order of 10 −10 ), hence the need of a high precision numerical scheme to iterate the recursion equations.  , α) where the fixed points are locally stable for ρ = 0.005 and r = 4 (red curve), r = 3 (green curve) and r = 2 (magenta curve).
The conventions for the labelling of the regions are the same as in Fig. 1, but we have omitted the regions where the fixed point e * > 1 appears. The zero-engineers fixed point is stable below the dashed blue line α = δ/ (r − 1).
Since e * = 0 is unstable, the population eventually recovers from these near-extinction events as illustrated in the middle and right panels of Fig. 3. The nature of the transition between the periodic solutions and the zero-engineers fixed point is clarified in Fig. 4 where we fix the decay rate per generation to δ = 0.15 and show e t and v t for different values of α in distinct panels. The left and middle panels, which show results for values of α very close to the stability boundary α = 0.05, reveal what happens: the period of the oscillations increases as that boundary is approached and diverges when the zero-engineers fixed point becomes stable, so that e * = 0 and v * = 1 for all times. The right panel in Fig. 4 reminds us that the periodic solutions are not necessarily associated to near-extinction events and that they can reproduce smooth oscillations with the usual phase shift of one-quarter of a cycle between engineers and virgin habitats, which characterizes the predator-prey models [28].
Next we consider the effects of varying the parameters r and ρ that were kept fixed in the above analysis. Figure 5 shows the regions of stability of the fixed points for three values of the intrinsic growth rate r. To keep the size of the region where the zero-engineers fixed point is stable we use the scaled variable δ/ (r − 1) as the inde-pendent variable. The results indicate that increase of r increases the size of the regions of limit cycles and decreases the size of the regions of bistability, so variation of r produces quantitative effects only. This is in stark contrast with the effects of varying the regeneration rate per generation ρ shown in Fig. 6. In fact, increase of ρ decreases the size of the region of limit cycles, which eventually disappears altogether for large ρ, since there is no region left in the space of parameters where the fixed points e * > 0 and e * = 0 are simultaneously unstable. For instance, the curve for ρ = 0.3 shown in the right panel of Fig. 6 is given by the condition that the discriminant of the quadratic equation (14) is zero, so for large ρ the finite-engineers fixed point is locally stable wherever it is real. The piecewise curve for ρ = 0.1 is given by the composition of the conditions that e * is real and that the norms of the eigenvalues of the Jacobian (15) are less than 1. In addition, increase of ρ greatly increases the region of bistability as well as the region where e * > 1.
To conclude, we consider now a simple but instructive scenario where the natural resources are not renewable, i.e., ρ = 0. Figure 7 shows the time evolution of the population towards the doomsday fixed point v * = e * = h * = 0 and d * = 1 where only the degraded habitats are left. We have not considered this fixed point in Section II A 1 because it exists and is stable solely for ρ = 0 and δ > 0, regardless of the values the other model parameters. The interesting aspect this figure highlights is that the population survives for a long time after the irreversible disappearance of the virgin habitats, thus revealing that the usable habitats introduce a delay on the influence of the virgin habitats on the population that is absent in the usual predator-prey models. However, this delay appears in more sophisticated models of human-nature interactions in which humans accumulate surpluses (i.e., wealth) and use them when natural resources are scarce or unavailable [7].

IV. DISCUSSION
The derivation of the discrete-time recursion equations (5)-(8) assumes non-overlapping generations, which is clearly not the case for human populations. Our justification to use this approach is purely pragmatic since it is much easier and much less error-prone to iterate those recursions than to solve numerically the differential equations of a continuous-time model, particularly in the region of reversible collapses where the density of engineers is very close to zero and the periods of the limit cycles are exceedingly long. However, as the solutions do not exhibit short-time fluctuations and are all smooth in the scale of a few generations we do not expect any significant differences between the discrete and the continuous time approaches. In addition, the discrete-time approach can easily be extended to describe space-dependent problems, resulting in coupled map lattices that can be thoroughly studied numerically. In particular, in this context a dif-ferent type of population collapse can be observed when the expanding colonies of engineers reach the boundaries of their territories, thus ending the exploration of new virgin habitats and resulting in a sharp drop on the population density [16].
In agreement with the findings for the predator-prey type models [6,7], we find that a key parameter to determine the long-term outcome of the population dynamics of ecosystem engineers is the regeneration rate per generation ρ, which determines the renewability of the resource base (see Fig. 6). Cycles of prosperity, collapse and revival occur in the case the degraded habitats have a slow regeneration rate. In the extreme case of non-renewable resources (i.e., ρ = 0), the outcome of the dynamics is the irreversible collapse of the entire ecosystem (see Fig.  7). In the other extreme, when the degraded habitats rapidly recover into virgin habitats, the limit cycles disappear altogether and the outcome of the dynamics is determined by the basins of attraction of the finite and zero-engineers fixed points.
The local stability analysis of the zero-engineers fixed point yields a simple condition for the survival of the population, viz., the product of the parameter that measures the effective population growth rate (i.e., r −1) and the parameter that measures the efficiency of the transformation of virgin into usable habitats (i.e., α) must be greater than the decay rate per generation of the usable habitats, i.e., α (r − 1) > δ. Hence increase of the efficacy of the engineers to explore the virgin habitats or increase of their intrinsic effective growth rate will both guarantee the survival of the population, regardless of the regeneration rate ρ > 0 of the degraded habitats. The characteristics of the steady-states for large r and large α, however, are completely distinct.
On the one hand, increase of r favors limit cycles characterized by long periods of time when e t 1, which we refer to as reversible collapses. These collapses are reversible because they occur in regions of the model parameters where the zero-engineers fixed point is unstable. This scenario is expected since a large value of r produces a large density of engineers, which, according to the prescription (8), can transform all virgin habitats in a single generation, thus producing the cyclical pattern of feast and famine illustrated in the left and middle panels of Fig. 4.
On the other hand, increase of α favors the finiteengineers fixed points and leads to a stable balance between population and resources. In our setup, a small α corresponds to a low-tech society where the exploration of the virgin habitats requires a high density of individuals. Interestingly, our model predicts that technology improvements that allow a small population to explore efficiently the available natural resources, a situation that is modeled by increasing α while keeping all the other parameters fixed, can eliminate the dangerous cycles of prosperity, collapse and revival.
This result is especially notable because it neatly expresses the argument put forward by the growth opti- mists, namely, that technological progress can reduce the impact of resource depletion and avoid environmen-tal and economic collapse [29]. Of course, this happens because our model does not account (among many other things) for the Jevons effect that occurs when technological progress increases the efficiency with which a resource is extracted, but the rate of consumption of that resource rises due to increasing demand [30]. In our model, increasing demand amounts to considering the decay rate per generation dependent on the availability of usable habitats, i.e., δ = δ (h t ), so that the more goods are available, the higher their consumption by the population. The possibility of supporting the growth optimists viewpoints as well as addressing quantitatively the Jevons effect, the analysis of which we postpone for a future contribution, highlights the potential of the ecosystem engineers approach to model both ecological and economic aspects of the human-nature dynamics.