## Abstract

Semiarid ecosystems (including arid, semiarid and dry-subhumid ecosystems) span more than 40% of extant habitats and contain a similar percentage of the human population. Theoretical models and palaeoclimatic data predict a grim future, with rapid shifts towards a desert state, with accelerated diversity losses and ecological collapses. These shifts are a consequence of the special nonlinearities resulting from ecological facilitation. Here, we investigate a simple model of semiarid ecosystems identifying the so-called ghost, which appears after a catastrophic transition from a vegetated to a desert state once a critical rate of soil degradation is overcome. The ghost involves a slowdown of transients towards the desert state, making the ecosystem seem stable even though vegetation extinction is inevitable. We use this model to show how to exploit the ecological ghosts to avoid collapse. Doing so involves the restoration of small fractions of desert areas with vegetation capable of maintaining a stable community once the catastrophic shift condition has been achieved. This intervention method is successfully tested under the presence of demographic stochastic fluctuations.

## 1. Introduction

A major consequence of the habitat-degradation processes associated with the Anthropocene is the accelerated loss of habitats and worldwide species extinctions [1–4]. Moreover, the nature of this decay is likely to be catastrophic, and evidence from field, experimental and palaeoclimatic data along with modelling efforts strongly indicate that such decay will end up in collapse [5–7]. Such ecological shifts are a consequence of the multistable nature of ecosystems, which is especially relevant in semiarid ecosystems. Indeed, recent field studies of vegetation spatial patterns in drylands have identified multiple ecosystem states [8].

Because of the global relevance of drylands [9,10], dedicated efforts have been made towards forecasting the so-called green-desert transition from a vegetated ecosystem to a desert area [6,11–13]. These warning signals are associated with the presence of increasing fluctuations that are known to diverge close to critical points [7] as well as with changes in spatial patterns [14]. The definition and quantification of these warning signals have been a very active research area over the last decade [15]. Less attention has been given to another type of nonlinear phenomenon associated with transient behaviour when a parameter *μ* surpasses a critical value *μ*_{c}. It involves the presence of a special class of very long transients, also called delayed transition or ghost [16–20].

The importance of transients in ecological systems has been recently discussed [21,22]. Importantly, transients have been characterized in experimental systems with *Tribolium* sp. [23] as well as in field data for the crab *Cancer magister* [24]. In a ghost-related transient, the system will appear stationary, thus creating the false impression that the system is stable, despite its assured collapse. In particular, ghosts follow a universal scaling law [17,25] and the length of transients, *τ*, scales following an inverse square-root law:
1.1with *μ* and *μ*_{c} being the control parameter and the bifurcation value, respectively. This scaling law has been identified in physical systems [16,26], in mathematical metapopulation models for autocatalytic species [27], low-dimension [28–30] and high-dimensional replicator models exhibiting cooperation [31,32], and single-species discrete models with Allee effects [33]. In this sense, the presence of very long transients can force us to redefine the concept of steady state [21,22,34] and provide essential clues to the persistence and organization of populations over long timescales.

If ghosts appear in green-desert transitions, the take-home message is simple but leads to a grim picture: a given ecosystem might be on its way towards collapse despite its apparent stability. In this paper, we show that ghosts are indeed expected in the dynamics of semiarid ecosystem models. However, a positive message is that the ghost can be exploited to ensure ecosystem's persistence. We introduce a simple intervention method based on the recovery of small fractions of desert areas with vegetation as a way to keep dryland dynamics near the phase space regions where delays occur, thus allowing the ecosystem to remain stable once the bifurcation towards desertification has occurred.

## 2. Mathematical model

The starting point of our study will be the model developed by Kéfi and co-workers as a coarse-grained approximation of semiarid ecosystems [11]. This model captures the nonlinear interactions due to facilitation and the relevant alternative states found in semi-arid ecosystems (vegetated area, fertile soil with the potential of hosting vegetation and desert areas). In this context, we will not explicitly take into account spatial structure [35–37], an explicit coupling with soil moisture [38,39], or some sources of heterogeneity [40], which can play a relevant role. Here, we aim to understand the properties and management of long transients close to the green-desert shift found in a minimal model capturing the main mechanisms of facilitation.

The original model [11] considers a spatial lattice model in which each site can be in one of three alternative states: vegetated (*V*), non-vegetated but fertile (*S*) and destroyed (*D*). Each state can transition to another with some probability, e.g. . In this model, these transition probabilities are given by
Because of the nature of semiarid ecosystems, the transition probabilities are driven not only by the nearest vegetated neighbours (*ρ*_{V}) but also by non-local processes that depend on the fraction of global vegetation (Δ_{V}). Here, *r* is the spontaneous regeneration rate of a site due to, e.g., the arrival of seeds by the wind. The parameters *f* and *d* stand for the facilitation and the soil degradation rates, respectively. Regarding the spread of the vegetation, the dispersion of the seeds is modelled by the constant *δ* (ratio of vegetation due to the whole vegetated sites) and by the development (*b*) and death (*c*) of the seeds. Finally, the rate of vegetation decay is given by a probability *m*.

The microscopic transition rules provide the rates of growth and death that can be incorporated into a time-continuous (mean-field) modelling framework. The resulting equations are expressed in terms of balances between different processes:
Ignoring spatial correlations, we have Δ_{V} = *ρ*_{V} = *V*. This assumption leads us to the following system of differential equations:

The state variables of the system (*D*, *V*, *S*) correspond to the ratio of the area that is occupied by each different state, with *D* + *V* + *S* = 1. Thus, we can reduce the three-variable system to a two-variable one by means of the linear relation *D* = 1 − (*S* + *V*). The two-variable system is now
2.1and
2.2Given that we consider the whole semiarid ecosystems area, no income terms are included. In our case, the laws and conditions that affect the ecosystem depend on the variables in the ecosystem. For this reason, there is no spontaneous recovery of the degraded state to the non-occupied state (i.e. *r* = 0).

## 3. Results

### 3.1. Catastrophic shifts in semiarid ecosystems

The model given by equations (2.1) and (2.2) displays bistability: for a low soil degradation rate, *d*, two possible exclusive stable states are present: an ecosystem with vegetation and a completely desert ecosystem. Actually, this system has four fixed points. The first fixed point is (*V** = 0, *S** = 0), which corresponds, in the case of being stable, to the extinction of vegetation and the destruction of the fertile soil, with a dominance of the desert state. The extinction state is locally stable (see electronic supplementary material, Section S1), and the system displays an abrupt transition. This catastrophic shift leads to a collapse of the vegetated ecosystem once *d* > *d*_{c}, and the system falls into the desert state. These results are summarized in figure 1, where the fraction of vegetated habitat is plotted against *d*. The increase in *d* involves a monotonous decrease of vegetation. As mentioned, at a given critical value *d* = *d*_{c}, an abrupt extinction takes place due to a saddle-node bifurcation. In the bifurcation diagram, note that the two branches of the fixed points (one being stable and the other unstable) approach each other as *d* increases. Specifically, under the parameter values used in our analyses (i.e. *b* = 0.3, *c* = 0.15, *m* = 0.1 and *f* = 0.9), the bifurcation value would be *d*_{c} ≈ 0.22 (see Section S2 for an accurate computation of *d*_{c}).

The phenomenon that we are investigating here takes place right after the bifurcation. When *d* ≳ *d*_{c}, the system has a single steady state: the desert state. However, for a value of *d* close enough to *d*_{c}, the time to extinction *T*_{e} rapidly diverges. This time divergence always occurs when bifurcation values are approached [17] and leads to a paradoxical situation: the system appears to be in a steady state, but it suddenly collapses after a long time period. This phenomenon is also depicted in figure 1 by means of the spatial representation of the ecosystem. The snapshots in figure 1 display the vegetation pattern using four values of *d*. For case (1), the ecosystem is fully vegetated, while for cases (2) and (3), the fraction of vegetation is intermediate. However, case (3) is obtained with a value of *d* past the bifurcation point, but the system remains vegetated for a long time, eventually decaying into a completely desert state, represented by the sphere in case (4).

### 3.2. Dynamics near the green-desert transition

The detailed dynamics under different *d* values used in figure 1 are displayed in figure 2 by means of trajectories in the (*V*, *S*) phase space, confined to a simplex, as the sum of all of the variables is constant (i.e. *V* + *S* + *D* = 1). When *d* = 0.001, the attractor involving a fraction of vegetated *V* and of fertile *S* soil is close to the boundary of the simplex, and the basin of attraction of the fixed point (*V* = 0, *S* = 0) is extremely small (see the time series in figure 2(1) and snapshot 1 in figure 1). Typically, for a given value of *d* and for a given initial fraction of vegetated area *V*(0), trajectories in phase space go to the slow manifold (hereafter SM) and then approach the stable fixed point. The SM corresponds to the set of points in the phase space from which the attractor is reached. These dynamics can be clearly observed in the simplexes as well as in the time series of figure 2. As *d* increases, the basin of attraction of the fixed point (0, 0) (shown in red) becomes larger. For this case, the system is less vegetated (see sphere (2) in figure 1). Under this scenario, the dynamics within the basin of attraction of the fixed point involving a vegetated ecosystem makes the orbits to rapidly flow to the SM of this equilibrium point. Once this manifold is reached, there is a slow approach towards this equilibrium (see the thick trajectory in the zoom of figure 2(2).

The two previous scenarios involved a stable green state, provided that *V* (0) and *S*(0) were large enough. Once *d* is slightly above its critical value (figure 2(3)), all of the trajectories end up in the desert state attractor (0, 0). However, the ghost (or saddle-node remnant) appears close to the region where the stable node and the saddle collide, leading to a saddle-node bifurcation [17]. Now, the flows undergo a very slow passage through a bottleneck before reaching the desert state. This is confirmed by the slow decay displayed by the *V*(*t*) time series. The extinction dynamics becomes very fast at high *d* (see time series in figure 2(4)). Indeed, for a large *d*, trajectories usually miss the ‘gap’ on the slow manifold, resulting in shorter transient times towards full desertification.

Another way to see the qualitative changes associated with the emergence of a ghost is provided by the so-called potential functions or quasi-potentials [1]. These surfaces provide information on the stability of the fixed points. Specifically, stable and unstable fixed points appear as valleys and peaks, respectively. The quasi-potential for equations (2.1) and (2.2) in its general form reads (see electronic supplementary material, Section S4 for further details). The different surfaces associated with the previous flows are displayed in the right-hand side of figure 2. Two minima are clearly visible for *d* < *d*_{c}, whereas a drastic change occurs for *d* ≈ *d*_{c}, where the quasi-potential function has a single minimum (figure 2*c*). The previous two-well shape has now been replaced by an elongated, almost flat channel associated with the ghost and the observed bottleneck of the trajectories.

### 3.3. Delayed green-desert transitions

The dependence on *d* of the times to extinction, *T*_{e}, is displayed in figure 3*a*, where we can appreciate the almost stationary behaviour of the green state within the bottleneck, as compared with the steady ecosystem before the critical point. For *d* < *d*_{c} but close to criticality, the green state is stable, with a vegetation cover of *V* ≈ 27% (dashed line in figure 3*a*). If *d* is slightly increased, the system will achieve a seemingly stable fraction of vegetated habitat (of about 25%), which actually still holds at *t* ≈ 3.5 × 10^{5}. However, as the bifurcation has already occurred, the system rapidly collapses after this extremely long delay at *t* ≈ 350 000. Several examples are represented for slightly larger values of *d*. The same behaviour is found, although here, the times to extinction are shorter. Indeed, the dependence of the extinction times on the distance of the bifurcation parameter (*d*) from the bifurcation value (*d*_{c}) is shown to follow the inverse square-root scaling law (figure 3*b*), in agreement with previous works on delayed transitions tied to saddle–node bifurcations [16–18,28,33]. Here, *T*_{e} is plotted against the distance to the bifurcation value using *Θ* = *d* − *d*_{c}, and the power-law dependence is found numerically.

The delaying capacities of the ghost are known to depend on the initial conditions [27]. Usually, for large enough initial conditions, the ghost captures the orbits and the delays appear. This phenomenon is illustrated in electronic supplementary material, figure S3, which displays extinction times for values of *d* close to the bifurcation (by using 10^{3} × 10^{3} different initial conditions). Here, two well-differentiated regions are found. One region where times are extremely long (*T*_{e} ≈ 5 × 10^{5} in electronic supplementary material, figure S3a) exists at the right-hand side of the place where the saddle and the node collided (white area in the simplexes of electronic supplementary material, figure S3). The second one, indicated in black (electronic supplementary material, figure S3a), shows the initial conditions for which the orbits go to the desert state very rapidly. Indeed, extremely close to the bifurcation value, these two regions indicate initial conditions with extremely long transients or initial conditions for which extinction is almost immediate. Once the bifurcation value is further increased, the times decrease substantially. In both cases, long transients are about 1.5 orders of magnitude below the times of electronic supplementary material, figure S3a. Note that in electronic supplementary material, figure S3, we display time series for values at the right side of the simplex (upper panels with long transients) as well as to the left, where the transients are much shorter.

### 3.4. Controlling catastrophic shifts

The nature of the ghost behaviour allows us to enormously expand the extinction time by using small external interventions with suitable frequencies and magnitudes. In our case study, the interventions would involve restoring vegetation in small desert patches. The restoration parameters are the fraction of replanted area (expressed as the amount of vegetation increased, Δ*V*) and the restoration frequency (hereafter ‘freq’). Note that the time lapse between interventions is assumed to be much shorter than the timescale of the ecosystem dynamics. In this context, our restoration events are assumed to be effectively instantaneous. The intervention strategies can range from very rare (low-frequency) interventions affecting large areas to small but very frequent events. The effect of different values for both parameters Δ*V* and freq leads to different dynamical scenarios (see below).

The response of the ecosystem to a given intervention strategy depends on the actual ecosystem conditions and in particular on the increase of the soil degradation rate beyond the tipping point. To evaluate the response of the ecosystem, we have simulated the external interventions incorporating the external increase of vegetated area, performed at a given frequency (see electronic supplementary material, section S5). Figure 4*a*–*c* summarizes the nature of the intervention and its effects. Periodic perturbations create a cyclic behaviour (figure 4*a*), in which the trajectories are constantly re-injected close to the ghost. The effect of the intervention is to preserve (or return) the system within the region where the ghost has a strong delaying effect (see also electronic supplementary material, figure S4). When the pulses fail in re-injecting the trajectory towards the regions influenced by the ghost—i.e. to the right of the ‘gap’ on the slow manifold—the ecosystem decays to the desert state. On the other hand, when trajectories are efficiently sent to the ghost, the ecosystem remains in an stable green state.

To quantify how optimal the interventions are, we have to take into account that they depend on two parameters: freq and Δ*V*. For this reason, we obtain a surface that optimizes the intervention parameters for a given *Θ*, given by a Pareto front (electronic supplementary material, figure 4c). However, to find an optimal intervention strategy for the persistence of a given semiarid ecosystem, it is necessary to design an strategy that optimizes both parameters simultaneously. In our case, the optimality of an intervention is evaluated with a simple cost function, defined as the fraction of desert areas replanted per unit time ((plants/area) · time^{−1}) that is needed to sustain the system avoiding desertification, with the following equation:

In this study, the cost is assumed to increase if the fraction of replanted area and/or the frequency of planting increases. The changes in the cost as the degradation rate increases are shown in figure 4*e*. Note that the relation between the increase of degradation (*Θ*) and the optimal cost of intervention is almost close to a hyperbolic function (linear in log–log scale with exponent 1.08). The slope of this cost changes depending on the frequency of the interventions and on the area replanted. In the case of fixing the desired area to be replanted at each intervention, the optimal frequency will correspond to the Pareto front surface (electronic supplementary material, figure S5a). For larger amounts of added vegetation, Δ*V*, the increase in the cost function and the optimal intervention frequency are lower. The slope goes from 0.5 for large interventions (Δ*V* = 0.1) to 1.08 for very small replantings (Δ*V* = 10^{−4}); see electronic supplementary material, figure S5b.

As previously mentioned, another possible scenario may be found depending on the magnitude and frequency of the interventions, which may not be enough to preserve the ecosystem but would increase the extinction times (see electronic supplementary material, figure S6). Electronic supplementary material, figure S6a displays the increase in the extinction times with respect to the non-interventions case in the space (Δ*V*, freq) as a function of *Θ*. A particular example of this delayed transition is shown in electronic supplementary material, figure S6b. Here, the orbits can be re-injected multiple times until they escape from the influence of the ghost towards ecosystem collapse. Electronic supplementary material, figure S6c displays several cuts done in the three-dimensional space, as displayed in electronic supplementary material, figure S6a. Generically, these kinds of interventions are less costly than those ensuring full ecosystem preservation, and they could be enough to save time, as a larger intervention or ecosystem change can be performed. In the case of semiarid ecosystems, the livestock policy and an active change in the ecosystem (e.g. increase of the humidity) are changes that may return the system to a lower degradation soil condition.

### 3.5. Interventions under stochasticity

The previous results have focused on the catastrophic transition towards the desert state found in the mathematical model, which is deterministic. Here, we extend our analyses by considering stochastic fluctuations due to the ecosystem's finite size, which are not gathered by the deterministic model. To consider stochasticity, we performed stochastic simulations using the *Gillespie method* [41,42] and implemented the transition probabilities introduced in §2. Despite our main goal being to test the robustness of the interventions considering demographic fluctuations, we note that the impact of stochasticity on delayed transitions remains poorly explored. Also, we must note that, as our model has an absorbing boundary (i.e. the desert state), the system will always achieve this absorbing state under stochasticity. However, this mean time to extinction is extremely long, typically being of order , being the system's size or carrying capacity [43], compared to the timescales we are looking at.

Firstly, we computed how the stationary states of vegetation *V* change as the degradation rate of fertile soil increases under stochasticity. Hereafter, sizes of two different systems will be used: and total sites. The simulations for both sizes indicate that the catastrophic extinction of *V* will take place before the predictions of the deterministic model (figure 5*a*,*c*). As expected, the stochastic critical value of *d*, labelled *d*^{s}_{c}, becomes lower for smaller values of . This is the natural outcome of the so-called noise-induced transitions. In noise-induced transitions, the systems' size (which determines the intensity of the noise) can act as a bifurcation parameter [44]. In our model, noise is decreasing the threshold of the degradation rate of fertile soil causing the catastrophic transition. Interestingly, the dynamics for *d* ≳ *d*^{s}_{c} also display the characteristic long plateau typical of delayed transitions (figure 5*a*.1), which now is found just after the stochastic bifurcation. The delay for the stochastic system was quantified by computing the mean extinction times for *V* for values of *d* above *d*^{s}_{c}. The results, displayed in electronic supplementary material, figure S7 for both values of , indicate that extinction times increase by about two orders of magnitude when comparing *d* = 0.3 to *d* ≳ *d*^{s}_{c}. This delaying behaviour is clearly seen in the time series of electronic supplementary material, figure S7. For values of *d* close to *d*^{s}_{c} (see panels (1) and (4) in electronic supplementary material, figure S7), the stochastic trajectories undergo the plateau, with some of them being extremely long. The delaying behaviour is lost when *d* increases further away from *d*^{s}_{c} (panels (2, 3) and (5, 6) in electronic supplementary material, figure S7).

The same intervention strategy performed in the previous section was tested with stochastic simulations. Figure 5 indicates that the intervention method is able to sustain vegetation after the catastrophic extinction (yellow areas in panels (*b*) and (*d*)). Specifically, we display the probability of survival of the vegetated state, *P*_{S}, for different values of the intervention parameters (freq, Δ*V*), setting *d* = 0.22 > *d*^{s}_{c}. *P*_{S} is computed as the number of replicates surviving at *t* = 10^{7} from a total of 50 independent runs. The region where vegetation is sustained becomes larger for larger values of . The dynamics tied to the interventions are displayed in figure 5 by means of sustained stochastic trajectories (displayed up to *t* = 10^{5} but sustained until *t* = 10^{7}). An example of the stochastic delayed dynamics is displayed in the simplex (*V*, *S*), similarly to the analysis shown in figure 4*a* for the deterministic model. Finally, several other examples of successful and failed interventions are displayed in electronic supplementary material, figure S8. We must note that strong stochasticity (small ) and being far away from the transition involve a much larger frequency of interventions, increasing the intervention costs for sustaining the vegetation. An example of the former can be seen in panels (*b*) and (*d*) of figure 5, as the boundary separating the vegetated state from the desert state is displaced to the upper-right corner of the space (freq, Δ*V*), meaning that the cost (as defined in the previous section) increases. That is, for larger , the vegetation can be sustained by using lower frequencies and lower values of Δ*V* to keep the ecosystem stable at a lower cost.

## 4. Conclusion

In this study, we have introduced a dynamical restoration approach aimed at protecting drylands facing tipping points. This approach is based on a simple but well-established model for semiarid ecosystems [11]. The model exhibits bistable behaviour, specifically a green state (vegetated) and a desert state, separated by a sharp breakpoint involving a saddle–node bifurcation (first-order phase transition). We have shown that an intrinsic property of this model, namely the presence of a ghost, can be exploited to restore endangered semiarid vegetation by means of small, periodic replantings.

Our model considers three variables: vegetated soil, fertile soil and desert. The presence of a ghost has been identified as a function of soil degradation levels. Because of the apparently stationary nature of transients close to the tipping point, a given system might have already crossed the bifurcation but display vegetation cover until a rapid decline occurs. The shift is irreversible: a decrease in degradation *d* will not restore previous vegetated states once the desert state has been achieved. As a consequence, prevention strategies should help to avoid collapse, even when the tipping point has already being crossed. Such a goal has been shown to be achievable. We have quantified how these interventions (in terms of newly vegetated areas and frequency of application) might be useful depending on how far the system is placed beyond the bifurcation. Moreover, this intervention strategy has been tested successfully under demographic stochastic fluctuations, under which delayed transitions are also found. As expected, stochasticity involves a different bifurcation value, which results in lower critical degradation rates of fertile soil at decreasing systems' sizes. The lower soil degradation rates make the ecosystem more fragile under stochasticity. However, the intervention methodology also has been shown to work under demographic noise. Here, it is important to highlight two aspects. The first is that delayed transitions are extremely local phenomena. That is, the ghost effect rapidly disappears as *d* increases beyond *d*_{c}. Another important point is that our approach neglects spatial correlations, which are known to be crucial in ecological systems [45–47]. It is known that space usually enlarges transients [21,22]. However, the impact of space in delayed transitions remains unknown for ecological systems as well as for other spatially extended systems exhibiting saddle-node bifurcations.

Another important point is that the success of the interventions might depend on the spatial patterns involved in each habitat. Depending on the local or regional climatic conditions (e.g., moisture, slope and degradation) the vegetation can be homogeneous or organized in gaps, stripes and spots [48–50]. That is, for the same fraction of vegetated area introduced during the interventions, differences in the dynamics could arise depending on where vegetation is replanted and thus on the specific spatial patterns of the existing vegetation at the moment of applying the interventions. These points should be addressed in future research. Finally, a practical intervention will require the right choices, such as regarding the specific biological components used to restore local patches. Restoration ecology has already considered different strategies to speed the recovery of degraded land [51], including explicit consideration of alternative states [52,53]. Drylands have also been studied under the restoration perspective [54], and several strategies have been presented, involving different restoration approaches [55,56]. In our suggested scenario, plants play a central role as islands of fertility [57] along with the rich communities inhabiting soil crusts [57–59]. Both plants and soil crust interact and could be used as potential restoration agents, but we should not discard other, perhaps complementary stategies, such as synthetic biology approaches [60,61] or resource augmentation [55,62] with a local enrichment of soil patches in nutrients and/or moisture. The different methods outlined above should also be considered in potential practical alternatives to our basic intervention approach.

## Data accessibility

All the needed methods and results to replicate this article are part of the manuscript and the electronic supplementary material.

## Author's contributions

R.S. and B.V. designed the mean field model. B.V. and J.S. analysed the mathematical model and processed the data. B.V. performed the stochastic simulations. All authors wrote the manuscript and gave final approval for publication.

## Competing interests

We have no competing interests.

## Funding

This study was supported by an European Research Council Advanced Grant (SYNCOM), by the Botin Foundation, by Banco Santander through its Santander Universities Global Division, by the grant FIS2015-67616-P, the PR01018-EC-H2020-FET-Open MADONNA project and by the Santa Fe Institute. This work has also counted with the support of Secretaria d'Universitats i Recerca del Departament d'Economia i Coneixement de la Generalitat de Catalunya. J.S. has been also partially funded by a ‘Ramón y Cajal’ Fellowship (RYC-2017-22243) and by the CERCA Programme of the Generalitat de Catalunya. The research leading to these results has also received funding from ‘la Caixa’ Foundation as well as from a MINECO grant awarded to the Barcelona Graduate School of Mathematics under the ‘María de Maeztu’ Programme (grant MDM-2014-0445).

## Acknowledgements

The authors thank Sergi Valverde for facilitating the software used to build the spheres in figure 1 as well as Antoni Guillamon, Tomás Alarcón and Ernest Fontich for their useful comments. We want to especially thank J. Tomás Lázaro for his technical assistance. We also thank Nuria Conde and the other members of the Complex Systems Lab for helpful discussions as well as Victor de Lorenzo and Fernando Maestre and his team for useful discussions.

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.4127387.

- Received February 5, 2018.
- Accepted May 29, 2018.

- © 2018 The Author(s)

Published by the Royal Society. All rights reserved.