Self-Replication of Localized Vegetation Patches in Scarce Environments

Desertification due to climate change and increasing drought periods is a worldwide problem for both ecology and economy. Our ability to understand how vegetation manages to survive and propagate through arid and semiarid ecosystems may be useful in the development of future strategies to prevent desertification, preserve flora—and fauna within—or even make use of scarce resources soils. In this paper, we study a robust phenomena observed in semi-arid ecosystems, by which localized vegetation patches split in a process called self-replication. Localized patches of vegetation are visible in nature at various spatial scales. Even though they have been described in literature, their growth mechanisms remain largely unexplored. Here, we develop an innovative statistical analysis based on real field observations to show that patches may exhibit deformation and splitting. This growth mechanism is opposite to the desertification since it allows to repopulate territories devoid of vegetation. We investigate these aspects by characterizing quantitatively, with a simple mathematical model, a new class of instabilities that lead to the self-replication phenomenon observed.

In arid and semi-arid landscapes around the world, it is common to encounter non-uniform vegetation covers exhibiting large spatial structures, generically called vegetation patterns [1][2][3] . These landscapes are characterised by either water limited resources and/or nutrient-poor territories. In the former case, the potential evapo-transpiration of the plants exceeds the water supply provided by rainfalls. At the level of individual plant, the water scarcity provokes an hydric stress that affects both the plants survivability and growth rate. At the community level, this hydric stress promotes clustering behaviour which induces spatial landscapes fragmentation. It is now generally admitted that this adaptation to hydric stress involves a symmetry-breaking modulational instability leading to the establishment of a stable periodic spatial patterns  .
Vegetation patterns are not always periodic. The spatial distribution of vegetation may consists of isolated or randomly distributed patches or gaps. Such irregular patterns can involve groves within grasslands [27][28][29][30] or spots of bare soil within a grass matrix 31 . They consist of patches which are either isolated or forming clusters. In both cases, such patterns have been interpreted as localized structures 9,[27][28][29][30][31] .
The aperiodic patterning phenomenon is not specific to peculiar soils or plant species. Localized vegetation patches or gaps may develop on soil ranging from sandy and silty to clayey, the nature of vegetation may consist of grasses, shrubs and trees. The extension of a patch can vary from small clumps of grasses (0.5-2 m 2 ) to large groves of mulga (Acacia aneura) trees (100-1000 m 2 ), such as those observed in central Australia 32 . On the other hand, the formation of localized patterns is an important issue not only in plant ecology context and environmental sciences but also it is a multidisciplinary area of research involving physics, chemistry and mathematics 33 .
The fingering instability of planar fronts leading to the formation of labyrinth structures has been reported 55 . Similar phenomenon has been observed in fingering instability of localized structures 56 .
In this article, we investigate the self-replication mechanism in the context of natural vegetation ecosystems. We show that this phenomena is robust as it is observes in a wide range of species and size scales. By analysing satellite images from the semiarid ecosystem of the Catamarca region, Argentina, we show the emergence of characteristic statistical distributions of the vegetation patches and their spatial organisation. We consider a general interaction-redistribution model, where analytical and numerical results show that there exist a critical value of the level of the aridity under which a single circular vegetation patch destabilises, the curvature instability leads to an elliptical deformation followed by patch multiplication. This process continue in time until the system reaches a self-organized vegetation pattern in an hexagonal form. To compare field and numerical observations, we construct an initial condition for the simulations given by randomly distributed patches (which mimic long-range seed spreading), each patch is considered to be in a different stage of the self-replication process (this mimics the different ages of each patch). Under these considerations, we obtain a fair agreement between field and numerical observations, showing how self-replication is one of the mechanisms that mediate the spatial distribution and propagation of the vegetation in scarce environments.
This work is organized as follow: first, we study the spatial distribution and self-replication process in a real ecosystem by remote sensing imagery; secondly, we present a simple mathematical model, which is suitable for the description of vegetation pattern formation in semi-arid ecosystem and exhibits self-replication; thirdly, by numerical simulations and comparison to field observations, we show that self-replication can be an important mechanism for the phytomass repopulation; finally, we present the conclusions and projections of our work.

Field observations of self-replication
Andes highlands are semi-arid ecosystems with a low amount of available resources. In particular, the Catamarca region in NW-Argentina (− 23.436253°, − 65.976767° at 3424 m a.s.l.), presents an average annual rainfall that reaches 369 mm (source: grid of climate observations CRU CL 2.0 57 ), with a maximum in January of 71 mm and a minimum in July of 6 mm, temperatures vary from warm in the day to sub-zero in the night. Here, it is well known that Festuca orthophylla which produces tall, evergreen tussocks dominates the landscape over extended areas and periods of time at elevations between 3225 m and 4860 m a.s.l [58][59][60] . This specie is present in a variety of cold climates, adapting to diverse rainfall and soil moisture conditions. Festuca tussock arrange in circular shape compact structures composed by thousands of tightly packed tillers. The size of the tussocks depends on the resources available and weather conditions of their location, for instance, in western Bolivia they can reach 1.6 m high 59,60 . An important characteristic of Festuca is their shallow rooting system, which has been reported to cover an area 6-fold the area of the above ground canopy 59,60 . This quality allows each plant to have access to the resources in a total area equivalent to 6-fold the area of the projected canopy. This root network is the most important mechanism to capture resources in scarce environments, and also allows tussock-tussock competition for resources. This competitive interactions will be important throughout this article in understanding the observed spatial organisation of the tussocks.
This site, was selected in order to have a minimal slope, and no topographic perturbations such as mountains, canyons, rivers or highways. The region studied here covers an area of 109.4 km 2 (384 m × 285 m), the image corresponding to the site was obtained using Google Earth Pro and consisted on a 4800 × 3562 pixels image.
Festuca structures can be easily detected through satellite image analysis for their high light absorption (they appear as dark spots). As mentioned previously Festuca organises in tightly packed structures of circular shape. An example of isolated circular patch is shown in Fig. 2a. However, we have observed that an important number of structure have lost their circular shape, this curvature instability is the mechanism by which a tussock loses its circular shape by growing into an elliptical shape, the death of the tillers in the central part of the elongated structure causes the structure to split into two independent tussocks, we term this process self-replication, different stages of this process are shown by I, II, and III in Fig. 2c. This process is common to a wide range of species and scales, as observed in Fig. 1, where self-replications can be observed for structures in the scale of meters to hundreds of meters.
To address the question on how are Festuca tussocks distributed spatially, we perform a series of studies involving measurements of tussock properties and spatial distribution properties.

Spatial distribution analysis.
For studying the properties of the Festuca tussock, we have considered satellite imagery obtained directly from Google Earth Pro. For detecting Festuca patches, we performed image enhancement, which consisted on transforming the image to gray-scale, for the removal of background and improving contrast, a median filter and an adjustment of intensity was applied (Matlab R2105b), the resulting image can be seen in Fig. 2b. By turning the image to binary we could clearly identify the Tussocks. For the analysis objects in contact with the borders of the image were removed as they may not be completely observed, thus, introducing erroneous measurements to the analysis. The spatial resolution of the satellite images is 0.3 m, structures smaller than this size in diameter are not considered, however, this does not affect our spatial distribution analysis as we hypothesise that the spatial distribution of the vegetation will be dominated by the older/bigger tussocks as a consequence of their fully developed shallow root systems.
After the detection of the patches, the boundaries of each object and their properties (area, position, and radius) can be precisely computed. The relation of meters per pixel is extracted directly from the image (0.08 m/ pixel). A total of 3204 structures where detected. The first analysis, corresponds to a detailed characterisation of small distance properties such as nearest neighbour distance (NND), area covered, and equivalent radius of each structure. The equivalent radius is obtaining by comparing the area of each structure with that of a circle of the same area. The nearest neighbour distance is computed as the minimal boundary-to-boundary distance, this will allow us to extract information on the relation between NND and root sphere size, which are useful in understanding the underlying pattern emerging from redistribution and competition for resources. To expose the emergent spatial order in tussocks distribution, spatial Fourier analysis was performed, its circular average allowed us to determine a characteristic wavenumber in the spatial system. For the effect of Fourier analysis, a square sub-figure was selected from the original one to avoid border size effects. Finally, we have computed the Voronoi tessellation for the centres of the structures. This analysis, suitable for studying the regularity of a pattern has been used previously in aerial analysis to address pattern formation in vegetation structures 61 . For both NND and Voronoi cell computation, a subset of structures where selected such that they where far enough from the images border to avoid error induced from non visible structures. As stated previously, the remote sensing image analysis of the Catamarca region in NW-Argentina a total area of 109,44 km 2 (384 m × 285 m) was studied. Structures found in the analysis ranged from an area of 0.09 m 2 (minimum considered) to maximum area of 6.25 m 2 with a mean of 0.95 m 2 with s.d. of 0.79 m 2 . By visual inspection, we have noticed that bigger structures are most probably evolving clusters of structures. The average equivalent radius was found to be 0.50 m with s.d. of 0.22 m as can be seen in Fig. 3a. For calculating the minimal distance between the objects boundaries, only structures which are far enough from the image's border are considered, as the objects not captured in the image could be closer to these structures than the other observable ones. By this consideration, distances between 2837 objects are viable, which vary from 0.25 m to 4.63 m, and averaged a distance of 1.83 m with s.d. of 0.77 m (see Fig. 3b).
Considering the average size and NND, and assuming that the distance between tussock is set under the constrain that the root spheres do not overlap, we can estimate the projected root sphere size of an average tussock as the half of the NND, this results in a root sphere of 1.4 m radius and 6.3 m 2 area, which corresponds to an area of 6.7-fold the area of the average structure. This ratio is in fair agreement with the reported 6-fold ratio reported previously 60 for Festuca in the Bolivian highlands, the 12% difference may suggest that vegetation is not efficiently covering the terrain, leaving fertile terrain unpolulated, this could be also an effect of the death cycle of tussocks.
Although the previous analysis give us insight on the properties of the structures and their distribution, it delivers no information on the spatial organization or the existence of a characteristic wavelength in the system, for this, we perform a Fourier analysis.
For the spatial Fourier analysis, we considered a square sub-figure of the original. This figure contained 1510 structures. The spatial Fourier transform is extensively used by pattern formation community to evaluate the degree of spatial organisation. Pronounced peaks in the 2D Fourier amplitude indicate not only the existence of characteristic lengths in the system but also a preferred spatial direction for the formation of the pattern.
From the spatial Fourier analysis we are unable to detect pronounced peak in the spectrum, indicating that there is no preferred direction for a pattern to form. However, in the circular average of the spectrum, we observe a maximum wavenumber at k max = 2.4 m −1 (see Fig. 4), this is the first sign that the system is arranging in such a way that a characteristic length L c = 2π/k max = 2.6 m emerges. As one would expect this characteristic length is related to the interaction between tussocks, no-overlap (between tussock root spheres) is achieved, in average, if the distances between the centres of the structures is at least 2.4 m according to estimations made in the previous section from the PDFs analysis. Root competition between plants generates a minimal distance between tussocks.  If we assume that the real ecosystem is self-organizing towards a state with a well defined characteristic length, then we should be able to observe some fingerprints of such a process in the macroscale.
One of the problems of detecting such patterns and fingerprints in natural ecosystems is the existence of high scale perturbations such as terrain inhomogeneities, weather conditions, wind, and animal presence among others. All this external noise, alters the interactions between tussocks therefore, mangles their ability to form a regular pattern. The task of finding some type of unitary cell in the tussoks arrangement can be facilitated by the introduction of the Voronoi tessellation 62 . Considering the centre of a structure, the Voronoi cell associated with that structure will correspond to all the points that are closer to its centre than to the center of any other structure (see Fig. 5a). In this sense, the Voronoi tessellation gives us information of the most probable cell arrangement. As it is observed in Fig. 5b, the most probable vertex number is 6, which evidences an underlying tendency of Festuca tussocks to arrange spatially in hexagons. This is reinforced by observing that 6-sided cells seem not to be randomly distributed, rather forming cluster as if the hexagonal pattern was propagating through the system. Important information can also be extracted from the tile area (c.f. Fig. 5c), each tile represent the amount of land that is closer to a certain tussock than to any other, thus, the nutrients present in that portion of soil will be more accessible for the corresponding center tussock. The average tile size is 21.8 m 2 , which corresponds to equivalent radius of 2.63 m, almost twice the radius estimated for the root sphere. We know that when considering multiple structures, the roots spheres determine the minimum distance between them, however, the tile are observed indicates that tussocks are disperse through the terrain and still have space available for increasing the population density. We should remarked that the images considered do not give information of young and smaller tussocks that could be germinating on the suitable terrain.
In the next part of the article we present a simple model for describing pattern formation in the context of vegetation dynamics, the aim will be to theoretically describe the self-phenomena, and to show how important these phenomena is, both in obtaining a qualitative agreements between real field and numerical observations and as a mechanism for repopulating landscapes.

Model for vegetation dynamics
Pattern formation in vegetated environments has been extensively studied both experimentally and theoretically. It is well known that competition for resources, such as, water and nutrients can lead to spatial self-organization. This behaviour favours the formation of a wide range of patterns that depend on the characteristics of both the environment and underground spatial distribution of roots.
Several models describing vegetation patterns and self-organization in arid and semiarid landscapes have been proposed during last two decades. They can be classified into three types. The first approach, often called generic interaction-redistribution models, are based on the relationship between the structure of individual plants and the facilitation-competition interactions existing within plant communities [4][5][6][7][8][9][10][11][12][13] . The second approach is based on the reaction-diffusion type of models. These take into account the influence of water transport by below ground diffusion and/or above ground run-off [14][15][16][17][18][19][20][21][22] . The third approach focuses on the role of environmental randomness as a source of noise-induced symmetry breaking transitions 23,24,26 . In particular, the formation of localised structures in vegetation, also called localised vegetation patches has been studied in the case of poor resources, isotropic and homogeneous environments. A particular approach is to consider a logistic equation with a non-local term for describing the spatiotemporal evolution of the normalized biomass b(r, t) (the normalization is made with respect to the total amount of biomass supported by the system considered), this equation reads 31 where r and t are the spatial coordinates and time, respectively. The factors k 1 and k 2 account for facilitation and competition mechanisms of the plant-to-plant feedbacks, respectively. The third term corresponds to seed dispersion mechanisms, D is the rate at which plants diffuse, and the kernels Φ in and Φ out weight the incoming and outgoing seed fluxes. The plant-to-plant interactions are considered to be of the form f c f c f c , , , and the interaction strengths φ f,c can be affected by both intrinsic and extrinsic factors 31 . The spatial extension of the plant-to-plant feedbacks are given by the kernels φ f,c , which in the case of the facilitation depend on the overground canopy which can provide a shelter for other plants to grow, conversely, the competition kernel depends on the root sphere size which depletes ground resources, preventing other vegetation to grow. In real vegetation the age, canopy size, and root sphere are related, as older and bigger plants require for a higher amount of nutrients, thus, increasing the root growth, this dependence of the root size on the above ground plant size is known as the allometric factor. If we consider that the interaction kernels correspond to Gaussian fields, which do not depend on allometric factors, via a weak gradient approximation, equation (1) can be reduced to a local, non-variational partial differential equation for the phytomass density ρ(r, t) 27 , which reads This equation contains three positive defined control parameters: η that account for the decrease-to-growth rate ratio; κ is the facilitation-to-competition susceptibility ratio; Δ is proportional to the square root of the facilitation-to-competition range ratio. The parameters Γ and α are the nonlinear diffusion coefficients. The real order-parameter equation (4) constitutes the simplest model of spatial dynamics in which competitive interactions between individuals occur locally. An important feature of this equation is the presence of nonlinear diffusion terms ρ ρ ∇ 2 and ρ ρ ∇ 4 , that render it non-gradient or nonvariational. These nonvariational terms are imputable to the dispersion process, if the dispersion is negligible then equation (4) is similar to the variational Swift-Hohenberg that is regularly derived in spatially extended systems. In that case, the coefficients of ρ ∇ 2 and ρ ∇ 4 are both independent of the biomass density.
The where at each point of the territory, the vegetation production and death are exactly balanced. They should be real and positive. Two situations must be distinguished according to the sign of κ. When κ ≤ 0, only the homogeneous steady state ρ s+ , defines the biomass density, for η < 0. It decreases monotonously with μ and vanishes at η = 0. When κ > 0, the physical part of homogeneous branch of solution extends up to the limit point ρ L = κ/2 and η L = κ 2 /4. In the range 0 < η < η L , the biomass density exhibits a bistable behaviour: the stable homogeneous branches of solutions ρ s− and ρ s+ coexist with the intermediate unstable branch ρ s 0 as shown in Fig. 6. The upper homogeneous state ρ + s undergoes a modulational or spatial instability (Turing instability) characterised by an intrinsic wavelength which measure the distance between two maxima or minima of the plant distribution. The threshold associated with the modulational instability is solution of the following cubic equation There exist more than one threshold associated with the modulational instability. In the following, we focus on parameter regime where the uniform plant distribution exhibit bistability (κ > 0) and a portion of this state becomes unstable with respect to the Turing bifurcation (for η > η L ) as shown in Fig. 6. In this parameter range, any small fluctuation around the uniform plant distribution ρ s+ will trigger spontaneously the evolution of the system towards a stationary, spatially periodic distribution of the biomass density which will invade the whole territory. A detailed nonlinear analysis of two-dimensional periodic vegetation patterns such as stripes (often called tiger bush), and hexagons consisting of either sparsely populated or bare areas alternate with dense vegetations patches have been previously reported 5 .
It should be emphasised that the model presented is an approximation in which allometic factors are not considered, thus, the relation between the age and size of the plants with their interaction kernel is not considered. Moreover, the life-death cycle of the plants is not included in the theoretical construction of the model, this implies that a stable vegetation structure will persist in time without loss of phytomass. Now we will show the mechanism by which a localized patch can lose stability to generate multiple localized patches.

Self-replication as an extended pattern forming mechanism
When increasing the aridity parameter η, i.e. decreasing the amount of resources available, the structures that appear first are gaps. They consist of spots of spare vegetation. They exist until they lose their stability towards the formation of localized vegetated patches. The region where these localised patches are stable is limited by aridity values η I and η II shown in Fig. 6. When a localised vegetation patch is stressed by decreasing the aridity below η I , the patch exhibits an elliptically deformation followed by its splitting as shown in Fig. 7.
This self-replicating process continues until the system is entirely occupied by spots. Only spots which have available space around them are able to replicate. Because of this, only spots located in the edges can replicate. In the real ecosystem available space can be generated by the death of a plant by natural or external perturbations (animals, fires), this is the reason we can observe self-replication throughout all the territory analysed previously.
For long time evolution, transition from a single patch to a self-organized hexagonal pattern through a self-replication phenomenon is shown in Fig. 7, defects in the biomass distribution are attributed to boundary conditions used to numerically simulate Eq. (4). This hexagonal regularity is not observed in the arrangement of Festuca tussocks observed in Fig. 2b as vegetation in a real ecosystem is not nucleated by a single spot but rather developed from the random seed spreading by wind and animals thus generating multiple tussocks in different locations, each with the possibility of splitting to spread through the terrain.
To study the importance of self replication in the large scale organisation and distribution of the phytomass, we consider the following numerical approach: we construct an initial condition given by a 1000 × 1000 points field, which contains 1849 randomly distributed vegetation patches, constructed by a two-dimensional Poisson point process with rate r = 0.002. This random distribution aims to mimic the natural long-range seed spreading mechanisms, such as, wind, birds, or terrestrial animals which can transport seeds through long distances, these factors are not considered in the local interaction-redistribution model, but could be incorporated in the general non-local logistic equation (1). For including life-cycle factors in our approach, we consider that 185 randomly selected structures are in some stage of the self-replication process, these range from single to fully split patches. To conserve isotropy and homogeneity the direction of the splitting of each patch is also chosen randomly. This artificially constucted field is then considered as the initial condition for the simulation of Eq. (4). The aridity level in the simulation is set bellow the η I threshold, allowing each of the spots to continue the self-replication process, a portion of the resulting simulated field after 5000 iterations (with temporal step dt = 0.03) is shown in Fig. 8a. If we let the system evolve a sufficiently long time, then the system would reach an hexagonal pattern as observed in Fig. 7t 6 . However, when analysing the system in a intermediate state of evolution, we observe that the spatial distribution of the structures approaches qualitatively to the field observations of the Festuca tussocks in the Catamarca region, Argentina.
By computing the NND between the structures edges we obtain a PDF similar to that observed in the remote sensing analysis (see Fig. 8b). The dispersion around an average value of 9.5 A.U. (with a s.d. of 4.7) is generated by the self-replication process which alters the shape, consequently, the distance between structures. Structures at very small distances are not observed as the competitive interaction between spots causes a repelling force that generate a fast distancing between them once the splitting has occurred.
The Voronoi tessellation of the intermediate field shown in Fig. 8c exhibits a distribution where clusters of 6-sided cells are forming, in a system evolving towards and hexagonal organisation. In this stage, the tile vertex count (see Fig. 8d) and the cell area distributions (Fig. 8e) show that the characteristics previously observed in the field analysis are qualitatively well described by the generic interaction-redistribution model considered, under the initial conditions given. The existence of the self-replication mechanism is fundamental in the spatial organization of the vegetation. If only a Poisson point process is considered for the location of the localized patches, without self-replication, similar distributions can be obtained, even for the Voronoi tessellations, however the Fourier transform of such field . The localized patch suffers a curvature instability subsequently accompanied by the emergence of two spots. In turn these spots suffer a similar instability thus generating more spots, which begin to invade the system generating the emergence of a hexagonal pattern with several defects. is homogeneous, showing no peaks, on the contrary, when giving the liberty of self-replication to the localized patches, immediately we observe the emergence of spatial organization, given by peaks in the Fourier transform, these peaks appear distributed in the form of a ring similar to the one exhibited by the field observations, however they differ in that the numerical one shows a wide central peak (see Fig. 8f).
Despite the simplicity of the model considered for the description of the vegetation dynamics, we have shown that the self-replication induced by the diminution of the aridity parameter and the distribution induced by competitive interactions, mediate the spatial organisation of the vegetation in semi-arid ecosystems.

Conclusions
We have studied the self-replication phenomenon in the context of vegetation dynamics. Through remote sensing analysis of the Andean highlands we showed the existence of self-replication in Festuca tussocks, in this process the shrub is affected by a modulational instability that deform the structure from its circular shape into an elliptical shape, process after which the tussocks split into two new structures, we have also observed this process in a variety of species and size scales. By statistical analysis we have encounter characteristic distributions which are signatures of an underlying self-organization process. Though a general interaction-redistribution model that exhibits self-replication of localised structures we have shown that under certain initial conditions, the self replication and competitive interactions are sufficient conditions to exhibit spatial properties as the ones observed in natural ecosystems. These properties of the vegetation dynamical systems are the underlying mechanisms which mediate the extended self-organization of tussocks in arid and semi-arid ecosystems.
Self-replication in vegetation gives new lights on the way plants propagate and populate scarce environments.