Linking the Pattern Structures to System Robustness Based on Dynamical Models and Statistical Method

Pattern structures are usually used to describe the spatial and temporal distribution characteristics of individuals. However, the corresponding relationship between the pattern structure and system robustness is not well understood. In this work, we use geostatistical method–semivariogram to study system robustness for different pattern structures based on three dynamical models in different fields. The results show that the structural ratio of different pattern structures including the mixed state of spot and stripe, cold spot, stripe only, and hot spot are more than 75%, which indicated those patterns all have strong spatial dependence and heterogeneity. It was revealed that the systems corresponding to the mixed state of spot and stripe or cold spot are more robust. This article proposed a method to characterize the robustness of the system corresponding to the pattern structure and also provided a feasible approach for the study of “how structures determine their functions.”


INTRODUCTION
Due to some behavior mechanisms of individuals, species present heterogeneous but regular spatial distribution structures in both space and time, which is called as "pattern." These pattern structures exist widely in nature such as the clouds in the sky [1], the patterns on zebra [2], and the ripples on the water [3]. Except for these, Getzin et al. found the gap vegetation pattern-fairy circles in Western Australia [4], the regular stripes vegetation distribution on the hillside of Niger studied by Klausmeier et al. [5], and mussel beds in the intertidal zone show different scales of distribution, namely, large-scale banded distribution at the ecological level and small-scale reticular distribution at individual mussels level [6,7]. There are also thermal convection patterns, spiral wave patterns, and hexagonal patterns observed in the experiment [8,9].
The scientific community has a wide range of interest in the formation mechanism and structural characteristics behind the pattern. Consequently, theories of pattern dynamics have been deeply studied and form a more systematic theoretical research area [10][11][12][13][14][15][16][17][18][19][20][21]. Here, we present some typical works on this topic. In 2001, von Hardenberg et al. pointed out that when the bare state and spot pattern state coexist, and it is also unstable, if exceeding a certain threshold, the system will be completely transformed to a bare state or desertification, that is, the spotted pattern will be used as an early warning signal of desertification [22]. In 2014, Liu et al. revealed that the interaction of self-organization behavior between different scales can improve the robustness, persistence and productivity of mussel ecosystem; in other words, the mosaic patterns with large and small scales imply the ecosystem is more robust [6]. In 2020, Bastiaansen et al. quantified the resilience of ecosystems with spatial patterns by using phase portrait [23]. However, due to the complexity of the system dynamical process and the lack of uniformity on robustness definition, many studies and conclusions are not comprehensive and have certain limitations, even lack of quantitative indicators for the robustness of each pattern structure.
In order to better answer the question "how the pattern structures determine the robustness of the system," we obtain a series of different pattern structures based on three dynamical models in different research fields: vegetation-water coupled model [24], epidemic spatial model [25], and predator-prey model [26], and use geostatistical methods to quantitatively describe and analyze the characteristics of all different pattern structures, so as to find out which type of pattern structures are more robust for the corresponding systems.

CHARACTERIZATION INDEX OF SYSTEM ROBUSTNESS
In ecology, the related concepts of robustness is complicated and imprecise, but most people agree that robustness can be divided into two categories: one is the ability of the system to resist leaving (maintain) the current state after the system is subjected to external disturbance [27]; the other one is the ability of the system to return to the original stable state after suffering disturbance [28]. Therefore, in this study, we will use semi-variogram to analyze the interaction and dependence of each component within the system, and thereby give the quantitative index of the second type of system robustness.
The semivariogram is a mathematical statistical method that can reflect the randomness and structural characteristics of the variable in the spatial distribution and also is the theoretical basis of geostatistics. First, in order to avoid the proportional effect in the study, it is necessary to test the data of the studied variable (Z) and judge whether it is normal distribution or approximate normal distribution; if not, data conversion shall be carried out to make it conform to normal distribution. Afterward, the calculation of the semi-variance function value is carried out. Finally, through simulation, we obtain the best fitted semivariogram model and some important indicators, such as nugget C 0 , sill C 0 + C, range A, and the structural ratio C/(C 0 + C). Figure 1 shows the basic process of this method and the calculation formula of the semivariogram as follows [29]: where γ(h) represents semivariogram, h is the step length, that is the sample point spatial distance, N(h) indicates the total number of data pairs when the sample point distance is h, and Z (x i ) and Z (x i+h ), respectively, are the values of the variable(Z) at the spatial position x i and x i+h . In Figure 1, the nugget C 0 represents the degree of random heterogeneity of variables in the region and the sill (C 0 + C) refers to the maximum variation of the system induced by structural variation and random variation. Then, the structural ratio C/(C 0 + C) indicates the contribution rate of structural factors to spatial heterogeneity, which can reflect the spatial dependence and the spatial heterogeneity of variables [31]. When the structural ratio is less than 25%, it indicates that the variable has weak spatial dependence; if 25% ≤ structural ratio ≤ 75%, it means that the variable has moderate spatial dependent; and the structural ratio > 75% indicates the spatial dependence of the variable is strong [32,33]. Generally speaking, the larger the structural ratio is, the stronger the spatial dependence and spatial heterogeneity of the variables within the system will be. For the species community, the greater the heterogeneity is, the richer the diversity will be. In other words, if the structural ratio of the pattern structure is larger, then it means that the corresponding system is more robust. FIGURE 1 | Flowchart of quantitative system robustness index [30]. The analysis of semivariogram is mainly divided into three phases: normal test of research data, calculation of the semi-variance of variable, and simulation semivariogram to obtain the best fitted model. Structural ratio C/(C 0 + C) can reflect the robustness of the system. Frontiers in Physics | www.frontiersin.org February 2022 | Volume 10 | Article 827929    Table 1 for the other parameter values.

THE CORRESPONDING RELATIONSHIP BETWEEN PATTERN STRUCTURES AND SYSTEM ROBUSTNESS 3.1 Vegetation-Water Dynamical Model
In arid and semi-arid areas, vegetation growth is mainly limited by water resources [5]. In 2004, Gilad et al. considered three feedback mechanisms between vegetation biomass and water resources, infiltration feedback, root augmentation feedback, and soil-water diffusion feedback mechanism (Figure 2), and established a dynamical model for coupling a single vegetation species with water [34,35]. On this basis, scholars made the following assumptions: 1) the lateral extension of vegetation roots is restricted and 2) the infiltration rate of bare soil and vegetation covered area is the same, and then obtained a simplified dimensionless model as follows [24]: In the previous dimensionless model, the biological significances and values of all parameters are shown in Table 1. Next, we study the influence of the shading rate ρ on the robustness of vegetation ecosystems. Fixing other parameters, through numerical simulation, we get a series of vegetation pattern structure, as shown in Figure 3.
For the previous four different vegetation patterns, hot spot, mixed spots and stripes, stripe, and cold spot, we use vegetation density as a research variable and perform a semi-variogram analysis on them. According to Table 2, it is found that the structural ratio C/(C 0 + C) of these four vegetation patterns in descending order are as follows: mixed spots and stripes > cold spot > stripe > hot spot, and the value of the aforementioned four vegetation patterns are both more than 75%, which shows that they all have strong spatial dependence and spatial heterogeneity. However, compared with the other three structures, the vegetation system with mixed spots and stripes has the strongest spatial heterogeneity, which means that the system under this pattern structure will be the most robust. For the vegetation-water dynamical model, the more robust the system is, the lower the possibility of desertification is. The spot-stripe mixed structure is conducive to vegetation diffusion in space, while the hot spot pattern shows that the vegetation presents an isolated high-density distribution, which implies the system is more susceptible to desertification. Based on experiments, Bertolini et al. also found that the vegetation system corresponding to the mixed spot-stripe pattern (labyrinthine-like pattern) is relatively robust [36], and thus, our conclusion is consistent with the previous findings.

Spatiotemporal Dynamical Model for Disease Transmission
The spread of infectious diseases will be affected by the spatial movement of the susceptible or infected, which will eventually lead to spatial patterns of the susceptible and infected individuals

Coefficient of determination.
j Strong spatial dependence.

FIGURE 4 | Schematic diagram of infectious diseases transmission
process. The model satisfies the following conditions: 1) the epidemic cannot be cured after being infected; 2) the susceptible is exposed to the disease repeatedly, leading to a non-linear transmission rate.  [37,38]. In view of the characteristics of acquired immunodeficiency Syndrome (AIDS), hepatitis B virus (HBV), ebola virus, and other infectious diseases that are difficult to be cured after infected, meanwhile considering the complex spatial dynamics of the susceptible and the infected, scholars proposed an epidemic spatial model with non-linear incidence rates [25]. The non-linear incidence rate is caused by twice exposures of susceptible before infection. The transmission mechanism of such infectious disease can be described in Figure 4. The biological significances and values of all parameters in the model are shown in Table 3. The specific mathematical model is as follows: Based on the spatial infectious disease model (3), we studied the impact of the spatial distribution of infected individuals on the spread of infectious diseases in the population. First, fixing other parameters, only changing the transmission rate β, and a series of numerical simulations are carried out on model (3). Finally, we get the pattern structure, as shown in Figure 5.
With the increase in the transmission rate, the spatial distribution of infected individuals present four different pattern structures: hot spot, mixed spots and stripes, stripe, and cold spot ( Figure 5). From Table 4, we find that the   structural ratio C/(C 0 + C) of the four different spatial distributions of the infected are all greater than 75%, and their ratios are arranged in the following order: spot-stripe mixed > cold spot > stripe > and hot spot ( Table 4). Although their structure ratio gap is small, its order is highly consistent with the order of the vegetation patterns, which demonstrates that our conclusion has certain universality.

Predator-Prey Model With Spatial Diffusion
The predator function response is an indispensable part of describing the predator-prey model. This function can reflect the influence of the competition between predators on the predation efficiency, such as the coyotes, and jackrabbits in the western wilderness of North America and the plankton experiment [39]. In particular, when the predator can capture a large amount of prey per unit time, or when saturation is not considered, a ratio-dependent functional response is obtained [40]. However, with the addition of spatial diffusion, this ratiodependent predator-prey model will produce a heterogeneous spatial distribution structure [26]. As a result, we consider this ratio-dependent predator-prey model with diffusion terms. After introducing dimensionless variables, the model is simplified to the following equations [26] zx The dimensionless variables in model (4) and their values are shown in Table 5. Fixing other parameters and changing the prey consumption rate by predator γ, through numerical simulation, we obtain prey pattern structures, as shown later ( Figure 6).
On the basis of the previous method, we conduct a semivariogram analysis on the previous two spatial distribution structures of prey. The detailed results are shown in Table 6. According to the heterogeneity classification standard, we find that the structural ratio C/(C 0 + C) of these two patterns is significantly higher than 75%. However, by comparison, the structural ratio of the cold spot pattern is larger than that of the stripe pattern, that is, its corresponding system is more robust.  Compared with the stripe structure, the cold spot pattern makes the prey to gather together and is not easy to be captured by the predators, and thus, the system is more robust.
The range A is an important index in geostatistics that can reflect the spatial heterogeneity or spatial dependence scale of regional variables. From the perspective of the range, we find that all the ranges obtained by analyzing the pattern structure above are larger than the sampling interval − − 1. Therefore, it shows that the sampling interval used in the study is credible for unbiased estimation of this area.

CONCLUSION AND DISCUSSION
System robustness and pattern structure are two important characteristics to portray the spatiotemporal complexity of systems. However, the analysis of their corresponding relationship is lack of systematic research results. Focusing on the scientific problem of "which pattern structure implies the robustness of the system, and which pattern structure means the vulnerability of the system," we combined three dynamical models from different fields, vegetation-water dynamical model [24], epidemic spatiotemporal dynamical model [25], and predator-prey spatial evolution dynamical model [26], and used geostatistical methods to analyze a series of different pattern structures produced by them. Finally, we found that the robustness of the system corresponding to different pattern structures is arranged as follows: spots and stripes mixed > , cold spot > , stripe > , and hot spot. The research results may provide some early warning signals for desertification prevention, infectious disease prevention and control, biodiversity conservation, and other related fields.
In addition, we explored the systems' resilience of different vegetation pattern structures (Figure 7), combined with the latest research methods [41]. We once again confirmed that the recovery rates of ecosystems with hot spot are the worst, compared with the other three pattern structures. And it also revealed that the recovery rate of the ecosystem corresponding to the spots and stripes mixed is greater than that of the cold spot. However, the only difference is that the recovery rate of system with stripe is the largest, which may be because the added disturbance increases the local density and finally affects the characteristics of the strip structure. The evaluation of system robustness generally requires a large number of monitoring data or experimental data [42,43]. However, our study is mainly based on the dynamical model to obtain the pattern structures, which can not only dynamically reflect the evolution of the system in time and space but also quantitatively predict the future spatiotemporal distribution structure. At the same time, the semi-variogram analysis method can comprehensively analyze and compare the spatial characteristics of each pattern structure, overcome the complexity of the previous analysis system robustness, and provide a new idea for describing the corresponding relationship between the pattern structure and the system robustness.
It is worth noting that our research ignores the scale dependence and does not combine with real data, while the pattern structure has different scales and can be completely corresponded to the real data. In this case, we can improve and verify our theoretical results based on GIS and big data analysis. Furthermore, the pattern structure can be divided into steady-state and non-steady-state structures, and this study only focuses on the former case. While for the non-steady-state pattern, its existence in the real world is more extensive, and hence, the correspondence between the unsteady structure and the robustness of the system is also an important scientific issue. We hope these questions will be systematically addressed in future research.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
All authors have made great contributions to the writing of study and approved the submitted version. G-QS, YP, LL, CL, YW, and ZW established dynamical modeling. G-QS, YP, CL, and YW participated in the program design and provided valuable comments on the manuscript writing. G-QS, YP, and LL collected and processed the relevant published data. GF, ZJ, and B-LL guided and improved the manuscript.