Localised labyrinthine patterns in ecosystems

Self-organisation is a ubiquitous phenomenon in ecosystems. These systems can experience transitions from a uniform cover towards the formation of vegetation patterns as a result of symmetry-breaking instability. They can be either periodic or localised in space. Localised vegetation patterns consist of more or less circular spots or patches that can be either isolated or randomly distributed in space. We report on a striking patterning phenomenon consisting of localised vegetation labyrinths. This intriguing pattern is visible in satellite photographs taken in many territories of Africa and Australia. They consist of labyrinths which is spatially irregular pattern surrounded by either a homogeneous cover or a bare soil. The phenomenon is not specific to particular plants or soils. They are observed on strictly homogenous environmental conditions on flat landscapes, but they are also visible on hills. The spatial size of localized labyrinth ranges typically from a few hundred meters to ten kilometres. A simple modelling approach based on the interplay between short-range and long-range interactions governing plant communities or on the water dynamics explains the observations reported here.

The appearance of order and structures that involve nonequilibrium exchanges of energy and/or matter have been widely observed in many natural systems including fluid mechanics, optics, biology, ecology, and medicine [1][2][3][4][5][6] . Vegetation populations and vegetation patterns belong to this field of research. Being often undetectable at the soil level, large-scale vegetation patterns have been first observed thanks to the usability of aerial photographs in the early fifties 7 . They appear as extended bands of vegetation alternating periodically with vegetated areas and unvegetated bands. These large-scale botanical organisations have been reported in many semi-arid and arid ecosystems of Africa, Australia, America, and Asia. It is now widely admitted that the origin of these large scale botanical organisations is attributed to a nonequilibrium symmetry-breaking instability leading to the establishment of stable periodic spatial patterns. Extended and periodic vegetation pattern arising in semi-arid and arid ecosystems has been the subject of numerous studies and is by now fairly well-understood issue (see recent overviews [8][9][10] and references therein).
Vegetation patterns are not always periodic and extended in space. They can be spatially localised and aperiodic consisting of isolated or randomly distributed patches on bare soil [11][12][13] or gaps embedded in a uniform vegetation cover 14,15 . They are generated in a regime where the homogeneous cover coexists with periodic vegetation patterns. The interaction between well-separated patches is always repulsive 16,17 while for gaps the interaction alternates between attractive and repulsive depending on the distance separating gaps 14,17 . The localised patches has a more or less circular shape. However, for a moderate aridity condition, the circular shape can exhibit deformation followed by splitting of a single into two new patches. Newer patches will in their turn exhibit deformation and self-replication 18-20 until the system reaches a periodic distribution of patches that occupies the whole space available in landscapes 19,20 . This process leading to spotted periodic patterns can be seen as warnings of ecosystem degradation and may lead to outcome of vegetation recovery. Besides patches self-replication, circular spots can exhibit deformation leading to the formation of arcs and spirals like in isotropic and uniform environmental conditions 21 . The vegetation spirals are not waves since they do not rotate 21 .
In this work, we unveil a new type of vegetation pattern consisting of a localised labyrinth embedded either in a homogeneous cover or surrounded by bare soil. This phenomenon is observed in Africa and Australia by remote sensing imagery. An example of such a botanical self-organisation phenomenon is shown in Fig. 1. They consist of either an irregular distribution of vegetation surrounded by a uniform cover (see Fig. 1a, b), or by a bare state (Fig. 1c, d). We show that localised labyrinths embedded in a uniform cover can be stable even if the environment is isotropic and their formation does not depend on the topography. However, when a localised labyrinth is surrounded by a bare state, they can expand or shrink. We analyse this phenomenon by using three www.nature.com/scientificreports/ well-known self-organisation vegetation models, which support localised labyrinths. We show that localised labyrinths are permanent structures, and they can be observed worldwide involving a range of species and spatial scales. We interpret this phenomenon as a spatial compromise between the extended labyrinth that occupies the whole space available and stable homogeneous states. More precisely, the mechanism leading to their formation is attributed to the pinning-depinning transition that takes place in a parameter space where models exhibit bistability between extended disordered pattern and homogeneous cover.

Field observations of localised labyrinths
Localised labyrinths observed in nature are large-scale self-organisation patterns. They are satellite images from Africa and Australia obtained by the use of Google Earth software. The landscape of Central Cameroon (zone of forest-savanna mosaic 22 ), shown in Fig. 1a, displays contrasted phases of bare and densely vegetated areas with well-defined scale and symmetry surrounded by more or less uniform woodland. The climate in the zone where we observe the localised labyrinth is humid, with annual averaged precipitation of 1800 mm 23 . The annual averaged of potential evapotranspiration is between 1500 and 1600 mm 24 . The localised labyrinth we observe in Western Australia (see Fig. 1b) consists of localised woodland embedded in the shrubland of Mulga Bush (Acacia Aneura) 25 . In this zone the climate is arid, where the mean annual precipitation is 250 mm 26 and the mean annual potential evapotranspiration is between 1200 and 1300 mm 27 . Besides, the localised labyrinth can be surrounded by bare zones as shown in landscapes of Southwest Niger in a brush-grass Savanna zone 28 (Fig. 1c, d). In this region the climate is semi-arid, the mean annual rainfall is 605 mm, in between June and September 29 , and the annual mean potential evapotranspiration near this zone is 1900 mm 30 . All the climate data is summarized in Table 1 in Methods section. Sparsely populated or bare areas alternate with dense vegetation irregular bands or patches made of microchloa Indica. The field observations suggest that localised labyrinthine structures are formed both in a flat landscape and with topographic variation (see Fig. 2). By their spatial regularity, by their spatial scales ranging from a few hundred meters to ten kilometres, as well as by the composition of their vegetation (tree, shrubs, herbs, and grasses), localised labyrinthine patterns are permanent structures, and they  www.nature.com/scientificreports/ www.nature.com/scientificreports/ can be observed even in non-arid climates. They have neither been observed nor reported. Understanding their formation and maintenance is an important ecological issue. The mechanism underlying the emergence of the localised labyrinth can be captured by using self-organisation mathematical models that can explain vegetation pattern formation within a unified conceptual framework. In this respect, two approaches will be used. The first is based on the relationship between the plants' aerialsubterranean structures, the facilitative and competitive feedbacks which act at the community level, and the plants' spatial propagation by seed dispersion 31,32 . The second approach incorporates explicitly water transport by below ground diffusion and/or above ground run-off [33][34][35] . These models are in reasonable agreement with the field observations [36][37][38] .

Mathematical modelling of ecosystems
The absence of the first principles for biological systems in general, and in particular for vegetation populations where phenomena are interconnected makes their mathematical modelling complex. The theory of vegetation pattern formation rests on the self-organisation hypothesis and symmetry-breaking instability that provoke the fragmentation of the uniform cover. The symmetry-breaking instability takes place even if the environment is isotropic 31,33,35 . This instability may be an advection-induced transition that requires the pre-existence of the environment anisotropy due to the topography of the landscape 34,39,40 . Generally speaking, this transition requires at least two feedback mechanisms having a short-range activation and a long-range inhibition. In this respect, we consider three different vegetation models that are experimentally relevant systems: (i) the generic interaction redistribution model describing vegetation pattern formation which incorporates explicitly the facilitation, competition and seed dispersion nonlocal interactions (ii) the local nonvariational partial differential model described by a nonvariational Swift-Hohenberg type of model equation, and (iii) the reaction-diffusion system that incorporate explicetely water transport.
The interaction-redistribution approach. The integrodifferential model. This approach consists of considering a well-known logistic equation with nonlocal plant-to-plant interactions. Three types of interactions are considered: the facilitative M f (r, t) , the competitive M c (r, t) , and the seed dispersion M d (r, t) nonlocal interactions. To simplify further the mathematical modelling, we consider that the seed dispersion obeys a diffusive process M d (r, t) ≈ ∇ 2 b(r, t) , with D the diffusion coefficient, b the biomass density, and ∇ 2 = ∂ 2 /∂x 2 + ∂ 2 /∂y 2 is the Laplace operator acting in the (x,y) plane. The interaction-redistribution reads where i = f , c . ξ i represents the strength of the interaction, N i is a normalisation constant. We assume that their Kernels φ i (r, t) are exponential functions with L i the range of their interactions. The facilitative interaction M f (r, t) favouring vegetation development. They involve the accumulation of nutrients in the neighbourhood of plants, the reciprocal sheltering of neighbouring plants against climatic harshness which improves the water budget in the soil. The range of the facilitative interaction L f operates on the crown size. The competitive interaction operates over a length L c and involves the below-ground structures, i.e., the rhizosphere. In nutrient-poor or/and in water-limited territories, lateral spreading may extend beyond the radius of the crown. This extension of roots relative to their crown size is necessary for the survival and the development of the plant in order to extract enough nutrients and/or water from the soil. When incorporating these nonlocal interactions in the paradigmatic logistic equation, the spatiotemporal evolution of the normalised biomass density b(r, t) in isotropic environmental conditions reads 14 The normalisation is performed with respect to the total amount of biomass supported by the system. The first two terms in the logistic equation with nonlocal interaction Eq. (2) describe the biomass gains and losses, respectively. The third term models seed dispersion. The aridity parameter µ accounts for the biomass loss and gain ratio, which depends on water availability and nutrients soil distribution, topography, etc. The homogeneous cover solutions of Eq. (2) are: b o = 0 which corresponds to the state totally devoid of vegetation, and the homogeneous cover solutions satisfy the equation The homogeneous cover state with higher biomass density is stable and the other is unstable. These solutions are connected by a saddle-node or a tipping point whose coordinates are given by b sn = (� − 1)/�, µ sn = e �−1 /� . The linear stability analysis of vegetated cover ( b s ) with respect to small fluctuations of the from b(r, t) = b s + δb exp{σ t + ik · r} with δb small, yields the dispersion relation Given the spatial isotropy, the growth rate σ (k) is a real quantity. This eigenvalue may become positive for a finite band of unstable modes which triggered the spontaneous amplification of spatial fluctuations towards the formation of periodic structures with a well-defined wavelength. At the symmetry-breaking instability the value of the critical wavenumber k c marking the appearance of a band of unstable modes, and hence the symmetry-breaking www.nature.com/scientificreports/ instability, can be evaluated by two conditions: σ (k c ) = 0 and ∂σ/∂k| k c = 0 . These conditions yield the most unstable mode This critical wavenumber determines the wavelength of the periodic vegetation pattern 2π/k c that emerges from the symmetry-breaking instability. Replacing k c in the condition σ (k c ) = 0 , we can then calculate the critical biomass density b c . The corresponding critical aridity parameter µ c is provided explicitly by the homogeneous steady states Eq. (3).
Local model: a nonvariational Swift-Hohenberg model. The integrodifferential equation (2) can be reduced by means of a multiple-scale analysis to a simple partial differential equation, in the form of nonvariational Swift-Hohenberg equation. This reduction has been performed in the neighbourhood of the critical point associated with the nascent bistability 14,32 . The coordinates of the critical point are: the biomass density b c = 0 , the cooperativity parameter c = 1 , and the aridity parameter µ c = 1 . These coordinates are obtained from Eq. (3) by satisfying the double condition ∂µ/∂b s = 0 and ∂ 2 µ/∂b 2 s = 0 . To apply a multiple-scale analysis it is necessary to define a small parameter that measures the distance from criticality and expand b, µ , and in the Taylor series around their critical values. The symmetry-breaking instability should be close to that critical point. To fulfil this condition, we must consider a small diffusion coefficient in order to include the symmetry-breaking instability in the description of the dynamics of the biomass density. This reduction is valid in the double limit of nascent bistability and close to the symmetry-breaking instability. In this double limit, the time-space evolution of biomass density obeys a non-variational Swift-Hohenberg model 14 where η and κ are, respectively, the deviations of the aridity and cooperativity parameters from their values at the critical point. The linear and nonlinear diffusion coefficients ν , γ , and α depend on the shape of kernels 17 . In addition to the bare state u = 0 , the homogeneous covers obey where the two homogeneous solutions u ± are connected through the saddle-node bifurcation u sn = κ/2, η sn = κ 2 /4 , with κ > 0 . The solution u − is always unstable even in the presence of small spatial fluctuations. The linear stability analysis of vegetated cover ( u + ) with respect to small spatial fluctuations, yields the dispersion relation Imposing ∂σ/∂k| k c = 0 and σ (k c ) = 0 , the critical mode can be determined where u c satisfies 4αu 2 c (2u c − κ) = (2γ u c − ν) 2 . The corresponding aridity parameter η c can be calculated from Eq. (7).
The reaction-diffusion approach. The second approach explicitly adds the water transport by below ground diffusion. The coupling between the water dynamics and the plant biomass involves positive feedbacks that tend to enhance water availability. Negative feedbacks allow for an increase in water consumption caused by vegetation growth, which inhibits further biomass growth.
The modelling considers the coupled evolution of biomass density b(r, t) and groundwater density w(r, t) . In its dimensionless form, this model reads 33 The first term in the first equation describes plant growth at a constant rate ( γ /ω ) that grows linearly with w for dry soil. The quadratic nonlinearity −b 2 accounts for saturation imposed by poor nutrients soil. The term proportional to θ accounts for mortality, grazing or herbivores. The mechanisms of dispersion are modelled by a simple diffusion process. The groundwater evolves due to a precipitation input p. The term (1 − ρb)w in the second equation accounts for the evaporation and drainage, that decreases with the presence of vegetation. The term w 2 b models the water uptake by the plants due to the transpiration process. The groundwater movement follows the Darcy's law in unsaturated conditions; that is, the water flux is proportional to the gradient of the water matric potential 41 . The matric potential is equal to w, under the assumption that the hydraulic diffusivity

Results
Localised labyrinthine vegetation pattern. In our analysis, we focus on the simplest vegetation model that has been derived from the interaction-redistribution approach, namely the non-variational Swift-Hohenberg Eq. (6) described above. This model is appropriate to describe the space-time dynamics of the biomass under resource-limited landscapes such as nutrient limitation or water deprivation. In this case, the average biomass density is low comparing the carrying capacity closed-packing density of unstressed vegetation. The simulated stationary localised vegetation labyrinth is shown in Fig. 3a. Moreover, to confirm the field observation and to show that this phenomenon is model-independent, we conducted numerical simulations of the other two models, the integrodifferential (Eq. (2)) based on the facilitative, competitive, and seed dispersion interactions; and the reaction-diffusion type that explicitly incorporates water transport (Eq. (11)). The results are shown in Fig. 3b and 3c. The parameters used to simulate the different localised labyrinths are listed in Tables 2, 3, and 4 in the "Methods" section. The localised labyrinth consists of one spatially disordered state surrounded by a qualitatively different state. Note that the localised labyrinthine patterns shown in Fig. 3 do not have a round shape. The fact that this shape is not round is attributed to the presence of defects in the disordered pattern since they modify the interface energy. Investigations of fronts propagation between labyrinths and homogeneous states mediated by defects are missing in the literature. The interface separating these two states is stationary leading to a fixed size of a localised labyrinth. It neither grows and invades the uniform cover nor shrinks. The stabilization of localised labyrinth is attributed to the interface pinning phenomenon 42,43 . This phenomenon is characterized by an interface that connects a homogeneous state and a periodic one, which is motionless on a finite region of parameters, pinning range. This pinning effect occurs due to the competition between a global energy symmetry breaking between states that favors the interface propagate in one direction and the spatial modulations that block the interface by introducing potential barriers 42 .
To determine the stability domain of the localised labyrinth, we establish the bifurcation diagram shown in Fig. 4a, where we plot the biomass density as a function of the aridity parameter η . The aridity refers not only to water scarcity but can be also attributed to the nutrient-poor soil. When the aridity is low obviously the uniform vegetated state is stable (blue line) and the bare state (broken line) is unstable. When the aridity parameter is further increased, the homogeneous cover becomes unstable with respect to small fluctuations. Above this symmetry-breaking instability, several branches of solutions emerge sub-critically for η < η c . Example of vegetation patterns that appears follows the well-known sequence made sparse vegetation spots that can be either periodic or localised in space (see i, Fig. 4a), banded vegetation (see ii, Fig. 4a) or a periodic distribution of localised patches setting on the bare state (see iii, Fig. 4a).
An extended labyrinthine pattern can be generated subcritically as indicated by the red line in the bifurcation diagram (see Fig. 4a). The situation which interests us requires that this extended labyrinth exhibits a coexistence with the uniform vegetated state. The coexistence between these two qualitatively different states is the prerequisite condition for the formation of a stable localised labyrinth. However, this condition is necessary but not sufficient, the interface separating these two states exhibits a pinning phenomenon 42 . Indeed, as shown in the inset of Fig. 4a, there exists a finite range of the aridity parameter often called the pinning zone η − p < η < η + p , where localised labyrinthine patterns are stable. Examples of localised labyrinth obtained by numerical simulations for fixed values of the control parameters are shown in Fig. 4a (iv, v, vi). The motionless interface is not necessarily circular, and contains bands perpendicular to it and circular patches. Similar bifurcation diagram is obtained from the integrodifferential model (see Fig. 4b).
Finally, we discuss the situation where the aridity is not homogenous due to the topography. For this purpose, we choose a top hat-like shape for the aridity parameter as shown in Fig. 5a. In this case, numerical simulations of the integrodifferential model Eq. (2) show a stable localised labyrinthine pattern (see Fig. 5a). Note that the localised labyrinthine structures surrounded by bare soil shown in Fig. 1c, d are unstable since the interface propagates. The interface can not be pinned in the absence of spatial oscillations around the bare state. Oscillations around this state are unphysical since the biomass density is a positive defined quantity. However, when the aridity parameter possesses an inverted top hat-like shape, it is possible to pin the interface (see Fig. 5b). In this case, the localised labyrinthine pattern is surrounded by a mosaic extended state, and the mechanism of stabilization is rather due to the inhomogeneity of the aridity parameter.
Deppining mechanism. The spatial location of the localised labyrinth immersed in the bulk of the stable uniform vegetated state depends on the initial condition considered. When ecosystems operate out of the pin- Table 1. Mean annual precipitation, potential evapotranspiration, and aridity index of the regions where localised labyrinthine patterns are observed. For more details on the meteorological data see the references given in the text. www.nature.com/scientificreports/ ning zone, the interface separating the labyrinth and the homogeneous cover propagates due to the depinning transition (see Fig. 6a, b). In this case, depending on the aridity level, the interface propagates from one stable state to another The transition is different when moving the aridity parameter slowly or abruptly. In the second type of variation, when η < η − p , the homogeneous cover invades the system, while when η > η + p , the localised labyrinth survives but it is now embedded by a periodic distribution of gaps (see Fig. 6b).

Conclusions
In this paper we have reported for the first time evidence of localised labyrinthine vegetation patterns observed on satellite images from Africa and Australia. We have shown that these localised structures are robustly consisting of either an irregular distribution of vegetation surrounded by a uniform cover or on the contrary surrounded by a bare state. We have shown that the formation of localized labyrinthine patterns is not specific to particular plants or soils. We have found localised labyrinths in ecosystems on flat landscapes and hills. Three relevant models which undergo localised vegetation labyrinthine patterns have been considered; (i) vegetation interaction-redistribution model of vegetation dynamics, which can generate patterns even under strictly homogeneous and isotropic environmental conditions. It is grounded on a spatially explicit formulation of the balance between facilitation and competition. Ecosystems experience transitions towards landscape fragmentation of landscapes (ii) the nonvariational Swift-Hohenberg model that can be derived from the model (i) in the long-wavelength pattern forming regime, and (iii) reaction-diffusion model that incorporates explicitly water transport. We have shown that all these models despite their mathematical structure support the phenomenon  www.nature.com/scientificreports/ of the localised labyrinth. We have established their bifurcation diagram and identified a parameter region, where we have observed a coexistence between a homogeneous cover and an extended labyrinthine structure which are both linearly stable. Within it, there exist a pinning zone of parameters where localised labyrinthine vegetation patterns have been generated as a stable pattern. Note however that localised labyrinth is determined by the initial condition, while their maximum peak biomass remains constant for a fixed value of the system parameters. This phenomenon results from front pinning between qualitatively different coexisting vegetation states. Outside of the pinning region, we have shown that the localised labyrinth either shrink and leads to the formation of regular distribution of circular spots or expand leading to the formation of an extended labyrinth. Finally, we have investigated the formation of localised labyrinth on a hill by considering an inhomogenous aridity parameter. This forcing acts as a trapping potential for the labyrinthine pattern. Owing to its general character, robust localised labyrinthine structures observed and predicted in our analysis should be observed in other systems of various fields of natural science such as fluid mechanics, optics, and medicine. We have documented for the first time the phenomenon of localised vegetation labyrinth by remote observations, using the Google Earth computer program, and numerical simulations of three different theoretical models

Methods
Google Earth data. The satellite images (cf. Fig. 1
The elevation profiles in Fig. 2 are obtained from Google Earth. This software uses digital elevation data from the Shuttle Radar Topography Mission at a resolution of 30 m 44,45 . The error, at a 90% confidence level, associated to the absolute height data is less than 6 m for the territories considered here (Africa and Australia) 44 .
Climate data. Localised labyrinthine patterns are observed in Central Cameroon (Fig. 1a), Western Australia (Fig. 1b), and Southwest Niger (Fig. 1c, d). The climate types of these regions are humid, arid, and semiarid, respectively. The climatic classification is based on the aridity index (see Table 1), which is the ratio of mean annual precipitation and potential evapotranspiration 46 . Note that the aridity index is small (big) when the aridity parameter ( η or µ ), defined in the interaction-redistribution approach subsection, is big (small).
Numerical simulations data. Numerical simulations of models under consideration were solved in square grids with Runge-Kuttta 4 time integrator. The spatial derivatives were approximated using finite difference scheme with a three point stencil using periodic boundary conditions. In the integrodifferential simulation, the convolution integrals were solved in Fourier space through DFT algorithms. The detail of the parameters used in the numerical simulations are listed in the Tables below.
Generation of numerical localised labyrinthine patterns. The localised labyrinthine patterns are initialised in a region of parameters where the uniform vegetation cover and the labyrinthine pattern coexist, in particular, in a pinning zone (see Fig. 4). The initial condition consists of a circular patch of labyrinthine pattern in the centre of the simulation box, embedded in a homogenous background (see Fig. 7). After a transient accommodation of the biomass field, the stable localised labyrinth emerges. The dynamics towards equilibrium in the integrodifferential model Eq. (2) is resumed in Fig. 7 and the Supplementary Video S1. Fig. 4 were determined with analytical and direct numerical integration techniques of the governing equations. The blue and black curves account for the vegetated state and the bare one, respectively. The curves are solid when the corresponding state is stable, and broken if unstable. The critical points in which the different states change their stability are determined by linear analysis, detailed in the interaction-redistribution approach subsection.

Computation of the bifurcation diagrams. The bifurcation diagrams in
The red curve is the stable branch of labyrinthine patterns, and it is determined by direct numerical integration of the governing equations (using the algorithm explained above). Starting from a vegetated state with a small amplitude noise perturbation, in the region where the uniform vegetation state is unstable, a stable labyrinthine pattern can emerge (see (ii) in Fig. 4). The stability range of the labyrinth state, that is (i) and (iii) in Figure 5. Localised labyrinthine patterns generated by inhomogenous aridity in the integrodifferential model. The spatially forced pattern can be supported by (a) the vegetated state (top hat-like shape µ parameter), and (b) the bare state (inverted top hat-like shape µ parameter). From numerical simulations, the figure was created using Inkscape 1.0 (https:// inksc ape. org/ relea se/ inksc ape-1. 0/). www.nature.com/scientificreports/ Fig. 4, are found by decreasing/increasing the aridity parameter starting from the labyrinthine pattern (see the black arrows in Fig. 4). The blue triangles account for the stable branch of the localised labyrinthine pattern. The initial condition is a stable localised labyrinth state (cf. state (iv) in Fig. 4). The aridity is decreased until the localised labyrinthine pattern becomes a localised hexagonal pattern, which determines the left boundary of the pinning region ( η − p or µ − p ). On the other hand, the right boundary of the pinning region ( η + p or µ + p ) is determined by increasing the aridity until the localised labyrinthine pattern invades all the system (see Fig. 6).  Table 3 on "Methods" section. The sequence t 1 = 1 to t 6 = 4 · 10 5 accounts for the evolution towards equilibrium of the localised labyrinthine pattern, starting from a circular patch of a labyrinth state embedded in a vegetated background ( t 1 ). The curve in the right shows the evolution of the average biomass density b , that is the double integral of the two dimensional biomass field b divided by the area of the simulation box (see the Integrodifferential model subsection). The stable labyrinthine pattern is reached in t 5 ≈ 10 5 iterations of the RK4 time integrator, when there is no change in b . From numerical simulations, the figure was created using Inkscape 1.0 (https:// inksc ape. org/ relea se/ inksc ape-1. 0/).