*Article* **Quantitatively Inferring Three Mechanisms from the Spatiotemporal Patterns**

#### **Kang Zhang <sup>1</sup> and Wen-Si Hu <sup>2</sup> and Quan-Xing Liu 1,2,3,\***


Received: 18 November 2019; Accepted: 7 January 2020; Published: 10 January 2020

**Abstract:** Although the diversity of spatial patterns has gained extensive attention on ecosystems, it is still a challenge to discern the underlying ecological processes and mechanisms. Dynamical system models, such partial differential equations (PDEs), are some of the most widely used frameworks to unravel the spatial pattern formation, and to explore the potential ecological processes and mechanisms. Here, comparing the similarity of patterned dynamics among Allen–Cahn (AC) model, Cahn–Hilliard (CH) model, and Cahn–Hilliard with population demographics (CHPD) model, we show that integrated spatiotemporal behaviors of the structure factors, the density-fluctuation scaling, the Lifshitz–Slyozov (LS) scaling, and the saturation status are useful indicators to infer the underlying ecological processes, even though they display the indistinguishable spatial patterns. First, there is a remarkable peak of structure factors of the CH model and CHPD model, but absent in AC model. Second, both CH and CHPD models reveal a hyperuniform behavior with scaling of −2.90 and −2.60, respectively, but AC model displays a random distribution with scaling of −1.91. Third, both AC and CH display uniform LS behaviors with slightly different scaling of 0.37 and 0.32, respectively, but CHPD model has scaling of 0.19 at short-time scales and saturation at long-time scales. In sum, we provide insights into the dynamical indicators/behaviors of spatial patterns, obtained from pure spatial data and spatiotemporal related data, and a potential application to infer ecological processes.

**Keywords:** Allen–Cahn model; Cahn–Hilliard model; spatial patterns; spatial fluctuation; dynamic behaviors

#### **1. Introduction**

The reaction–diffusion systems have been posed to address the spatially extended interactions among species or chemical substances since 1937 [1], when the traveling waves behaviors of the reaction–diffusion equations were discovered and studied. To date, the spatial systems are widely used to explore the spatiotemporal behaviors going beyond the initial frameworks such as existence and stability of spatial solutions in biology, geology, physics, and ecology [2–9]. The solutions to reaction–diffusion equations display diverse behaviors including the formation of traveling waves [10], wave-like phenomena [11], spatial self-organized patterns [12], and coarsening behaviors [13].

In mathematics, the spatial extend models can be classified into the following three categories with the simplest version. Firstly, the model describes the Brownian dispersion and growth-mortality processes on single species [14,15]. Herein, the species abundance can dynamically fluctuate around their equilibria that is limited by carrying capacity. Secondly, we can describe the mass-conserving

systems with the Cahn–Hilliard model if species have a biased movement behavior such as a local density-dependent relationship [13,16–18]. In fact, both of these processes may exist in an ecological system that lead to a third theoretical model, the Cahn–Hilliard model with demographical processes of species [19–21]. These models display similar patterns for the reasonable parameters on the same systems. How to infer the underlying mechanisms from the spatial behaviors are crucial themes for ecologists. Many spatial hallmarks were proposed to characterize their spatiotemporal behaviors including spatial correlation functions [22,23], spatial structure factor [24], LS law [25–27], the exponent of cluster size [28], and the exponent of spatial fluctuations [29], but few studies are hitherto proposed to comprehensively compare the hallmarks on the above-mentioned models.

In this paper, we present a comprehensive comparison of three typical PDE systems in two-dimensional space: the Allen–Cahn model, the mass-conserving Cahn–Hilliard model, and the Cahn–Hilliard model with demographical processes. The central theme is the characterization of the spatiotemporal dynamics of such systems displaying a spatially patterned behavior. A key insight from our spatial analysis, provided in Section 4, is how spatial patterns and spatial fluctuations can be used to infer underlying ecological processes and mechanisms from the indicators of the spatial fluctuation, spatial structure factor, spatial correlation functions, and LS growth behavior. Importantly, different ecological processes may lead a similar spatial patterns in some systems even for strikingly different ecosystems. Here, our proposed integrated indicators have potential for distinguishing the underpinning patterned behaviors.

#### **2. Theoretical Models**

#### *2.1. Allen–Cahn Model*

The original Allen–Cahn model is a reaction–diffusion equation describing the phase separation process in a multi-component alloy system [30], including from disordered to ordered transitions. Recently, McNally et al. [18] showed that Allen–Cahn model can well describe the spatial pattern formation of microbial community, where microbes use T6SS to kill neighbors to separate the initially mixed flora, forming cloned plaques that grow over time. Theoretically, the model is given as follows assuming the species density as variable *<sup>u</sup>* at the time *<sup>t</sup>* within a spatial domain of <sup>Ω</sup> (<sup>Ω</sup> ⊂ *<sup>R</sup>*2), i.e.,

$$\begin{aligned} \mu\_t &= -f(u) + \varepsilon \Delta u, u \in \Omega, \\ u \cdot \nabla u &= 0, u \in \partial \Omega, \end{aligned} \tag{1}$$

where *ut* is the partial derivative of *u* over time, *ε* is the diffusion coefficient, Δ is the Laplacian operator, *<sup>n</sup>* is the normal vector of the outer boundary, and <sup>∇</sup> is the gradient operator. Here, we take *f*(*u*) = *u*(*u* − *α*)(*u* − *β*), and 0 < *α* < *β*, then the zero points of *f*(*u*) is *u*<sup>0</sup> = 0, *α*, and *β*. Here, *u*<sup>0</sup> are the equilibrium points of the model in Equation (1).

It can be easily proved that 0 and *β* are stable equilibria of the model in Equation (1), whereas *α* is an unstable equilibrium. Over time, the solutions of the model in Equation (1) would be around 0 or *β*. In the spatial model in Equation (1), the solutions of phase separation will occur if *ε* > 0.

#### *2.2. Cahn–Hilliard Model*

The Cahn–Hilliard model describes the spontaneous separation of a mixed fluid into two pure phases [31], resulting in the formation of spatial patterns. This phase separation principle was used to elucidate the mechanism of pattern formation on small scales in a mussel bed ecosystem [13]. This self-organizing process results from the relationship between movement speed and local mussel density, which differs from the existing Turing principle. The general form of the Cahn–Hilliard model is given by

$$\begin{aligned} \mu\_t &= \varepsilon\_B \Delta (f(u) - \varepsilon\_A \Delta u), u \in \Omega, \\ u \cdot \nabla u &= 0, u \in \partial \Omega, \end{aligned} \tag{2}$$

where *ut* is the partial derivative of *u* over time, *ε<sup>A</sup>* and *ε<sup>B</sup>* are the local and nonlocal diffusion coefficients, <sup>Δ</sup> is the Laplacian operator, *<sup>n</sup>* is the normal vector of the outer boundary, and <sup>∇</sup> is the gradient operator. We take *f*(*u*) = *u*(*u* − *α*)(*u* − *β*) and 0 < *a* < *β*, then the zero points of *f*(*u*) are *u*<sup>0</sup> = 0, *α*, and *β*.

Mathematically, the stable equilibria of the Cahn–Hilliard model in Equation (2) are determined by the minimum values of the potential function *F*(*u*) with *F* (*u*) = *f*(*u*). Now, the potential function is expressed as *F*(*u*) = <sup>1</sup> <sup>4</sup>*u*<sup>4</sup> <sup>−</sup> *<sup>α</sup>*+*<sup>β</sup>* <sup>3</sup> *<sup>u</sup>*<sup>3</sup> <sup>+</sup> *αβ* <sup>2</sup> *<sup>u</sup>*<sup>2</sup> + *<sup>C</sup>* with constant *<sup>C</sup>*. *<sup>α</sup>* is the local maximum of *<sup>F</sup>*(*u*), 0 and *β* are the global minimum of *F*(*u*). Therefore, 0 and *β* are stable equilibria, and *α* is an unstable equilibrium. Note that Cahn–Hilliard model describes a mass-conservation ecosystem without mortality and birth processes.

#### *2.3. Cahn–Hilliard with Population Demographics Model*

It is natural to ask that what will happen when the population demographics are involved in the Cahn–Hilliard model. In general, we address the Cahn–Hilliard model with a population demographics (hereafter, called CHPD), which describes the growth and aggregation of a single species [16]. A similar model was used to describe patterns formation result from density-dependent motility rather than traditional chemotaxis behavior [20]. The general model is given by

$$\begin{aligned} u\_t &= \varepsilon\_B \Lambda (f(u) - \varepsilon\_A \Delta u) + ru(1 - \frac{u}{K}), u \in \Omega, \\ u \cdot \nabla u &= 0, u \in \partial \Omega, \end{aligned} \tag{3}$$

where *ut* is the partial derivative of *u* over time, *ε<sup>A</sup>* and *ε<sup>B</sup>* are the local and nonlocal diffusion coefficients, and Δ is the Laplacian operator. The last term is the logistic growth law, *r* (*r* > 0), is the net growth rate population, *K* is the environmental carrying capacity (*K* > 0), *n* is the normal vector of the outer boundary, and ∇ is the gradient operator. Now, the equilibria of the model in Equation (3) are 0 and *K*.

Linearizing the nonspatial model in Equation (3) using a first-order Taylor expansion at *u*<sup>0</sup> = 0 and *K*, respectively, yields

$$
\mu\_l = \lg(\mu\_0) + \lg'(\mu\_0)(\mu - \mu\_0),
\tag{4}
$$

where *g* (*u*) = *<sup>r</sup>* <sup>−</sup> <sup>2</sup>*ur <sup>K</sup>* . One sees that *u*<sup>0</sup> = 0 is an unstable equilibrium since *g* (0) = *r* > 0, whereas *g* (*K*) = −*r* < 0.

However, the spatial model in Equation (3) is linearized around *u*(*x*) = *u*<sup>0</sup> in Fourier space, i.e., *<sup>u</sup>*(*x*) = *<sup>u</sup>*<sup>0</sup> <sup>+</sup> <sup>∑</sup>*<sup>k</sup> <sup>δ</sup>u<sup>k</sup> <sup>e</sup>ik*·*x*+Λ*k<sup>t</sup>* yields

$$\begin{aligned} \delta u\_k &= \Lambda\_k \delta\_{k'}\\ \Lambda\_k &= -r - \varepsilon\_B f'(u\_0) \mathbf{k}^2 - \varepsilon\_A \varepsilon\_B \mathbf{k}^4, \end{aligned} \tag{5}$$

where *<sup>i</sup>* <sup>=</sup> √−1, *<sup>x</sup>* is position, and *<sup>k</sup>* is wave vector. If <sup>Λ</sup>*<sup>k</sup>* <sup>≤</sup> 0 is satisfied for all *<sup>k</sup>*, the system is stable. If Λ*k* > 0 is established, the system is unstable. Since *r* > 0, when *k* = 0, Λ*k* < 0. Hence, instability occurs if

$$\begin{aligned} \varepsilon\_A f'(u\_0) &< 0, \\ \varepsilon\_B f'(u\_0) f'(u\_0) &> 4\varepsilon\_A r. \end{aligned} \tag{6}$$

From this instability condition in Equation (6), one can see that the instability of the system is determined by the population growth rate (*r*) and the diffusion coefficients (*ε<sup>A</sup>* and *εB*). If the logistic growth term is ignored, then 0 and *β* are the stable equilibrium points and *α* is the unstable equilibrium point. The pattern of the model in Equation (3) will be phase separated at 0 and *β*. If both the diffusion term and the logistic term are considered, the phase separation is suppressed by the logistic growth term. Over time, the values of the patterned solution will be around 0 and *K*, which ascribes the Turing patterns.

#### **3. Materials and Methods**

#### *3.1. Numerical Simulations*

We use the finite difference methods [32] to solve the equation numerically in two dimensions. The Laplacian terms of the model in Equations (1)–(3) Care approximated by the following formula

$$\Delta u(\mathbf{x}, y) = \frac{u(\mathbf{x} + \Delta \mathbf{x}, y) + u(\mathbf{x} - \Delta \mathbf{x}, y) - 2u(\mathbf{x}, y)}{(\Delta \mathbf{x})^2} + \frac{u(\mathbf{x}, y + \Delta y) + u(\mathbf{x}, y - \Delta y) - 2u(\mathbf{x}, y)}{(\Delta y)^2}, \tag{7}$$

with step length Δ*x* = *Lx <sup>m</sup>* <sup>=</sup> 1.0 and <sup>Δ</sup>*<sup>y</sup>* <sup>=</sup> *Ly <sup>m</sup>* = 1.0. Here, *x* and *y* indicate the spatial coordinate. *Lx* and *Ly* are the physical length of the simulation in *x* and *y* directions and *m* is the grid number. The numerical simulations were implemented in two-dimensional space with grid numbers from *Lx* = *Ly* = 2048 to *Lx* = *Ly* = 20, 000, and periodic boundary conditions. For AC and CH models, the initial values are random numbers with a mean of 0.5. For CHPD model, the initial values are random numbers with a mean of 0.33.

The iteration over time uses the first-order Euler method, i.e.,

$$u(\mathbf{x}, y, t + \Delta t) = u(\mathbf{x}, y, t) + \Delta t f(\mathbf{u}),\tag{8}$$

with Δ*t* = 0.0025 that satisfied the stability condition [33].

All numerical simulations were implemented in Python 3.7 (Python Software Foundation, CWI, Amsterdam, Netherlands) and MATLAB 2019a (The MathWorks Inc., Natick, MA, USA). PyOpencl (https://documen.tician.de/pyopencl/) was used to accelerate the simulations; the code is available from Github (https://github.com/Kang-Zhang).

#### *3.2. Spatial Correlation Functions*

Spatial correlation function is widely applied to characterized spatial structure and its correlations biophysics [34]. It is usually expressed as

$$\langle \varrho(\mathbf{r}\_{l},\mathbf{r}\_{j})\rangle = \langle \left(\varrho(\mathbf{r}\_{l}) - \langle \varrho(\mathbf{r}\_{\uparrow})\rangle\right) \times \left(\varrho(\mathbf{r}\_{\uparrow}) - \langle \varrho(\mathbf{r}\_{\uparrow})\rangle\right)\rangle = \langle \varrho(\mathbf{r}\_{l})\varrho(\mathbf{r}\_{\uparrow})\rangle - \langle \varrho(\mathbf{r}\_{l})\rangle\langle\varrho(\mathbf{r}\_{\uparrow})\rangle,\tag{9}$$

where *r<sup>i</sup>* and *r<sup>j</sup>* are two different locations; *ϕ*(*ri*), *ϕ*(*ri*) are numerical solutions at *r<sup>i</sup>* and *rj*, respectively; and is the mathematical expectation. Thus, *<sup>g</sup>*(*ri*,*rj*) depends only on the displacement *<sup>r</sup>* between the two locations. That is,

$$\mathcal{g}(\mathbf{r}\_{i\prime},\mathbf{r}\_{j}) = \mathcal{g}(\mathbf{r}\_{i\prime}\mathbf{r}\_{i} + \mathbf{r}) = \mathcal{g}(\mathbf{r}).\tag{10}$$

Thus, the average of the correlation functions of all two points with a displacement of *r* is obtained as follows

$$\log(r) = \langle \sum\_{|r\_i - r\_j| \to r} [\langle \!\langle \boldsymbol{\varphi}(\boldsymbol{r}\_i) \!\!\boldsymbol{\varphi}(\boldsymbol{r}\_j) \!\rangle - \langle \!\boldsymbol{\varphi}(\boldsymbol{r}\_i) \!\rangle \langle \!\boldsymbol{\varphi}(\boldsymbol{r}\_j) \!\rangle] \rangle. \tag{11}$$

The spatial correlation function is a measure of the correlation of a variable between two spatial positions. If the value of the variable fluctuates in one direction at the same time, a positive value is taken; otherwise, a negative value is taken. If their fluctuations are completely unrelated, then the value is zero.

#### *3.3. Structure Factors*

Structure factor is an important indicator for analyzing long-range spatial correlations and scattering that is a function of scattering vector (*k*) and time (*t*). The normalized structural factor [35] is written as

$$s(\mathbf{k}, t) = \frac{\langle \frac{1}{N} | \sum\_{r} \exp^{-i\mathbf{k} \cdot r} [\varrho(r, t) - \langle \varrho \rangle] |^2 \rangle}{\langle \varrho^2(t) \rangle - \langle \varrho \rangle^2},\tag{12}$$

where *<sup>ϕ</sup>*(*r*, *<sup>t</sup>*) is the numerical solution at time *<sup>t</sup>* and position *<sup>r</sup>* and *<sup>ϕ</sup>* is the mean of all *<sup>ϕ</sup>*(*r*, *<sup>t</sup>*). Here, the structural factors *s*(*k*, *t*) are obtained by the Fourier transform on the solutions, and then the square of the modulus is taken and averaged in all directions.

#### *3.4. Density Fluctuation*

Density fluctuation is a measure to directly quantify observations of local crowd density, the rules that predict mass behaviors under new circumstances [36,37]. In a spherical window of radius *R*, the local variance *σ*<sup>2</sup> *<sup>ϕ</sup>*(*R*) associated with field fluctuations can be expressed in terms of the autocovariance function or the spectral function [29]. It is given by

$$
\sigma\_{\varphi}^{2}(R) = \frac{1}{v\_{1}(R)} \int\_{R^{d}} \psi(r) a\_{2}(r;R) dr = \frac{1}{v\_{1}(R)(2\pi)^{d}} \int\_{R^{d}} \tilde{\psi}(k) \tilde{a}\_{2}(k;R) dk,\tag{13}
$$

where *<sup>v</sup>*1(*R*) is the volume of a d-dimensional sphere of radius *<sup>R</sup>*; *<sup>r</sup>* <sup>=</sup> *<sup>r</sup>*<sup>2</sup> <sup>−</sup> *<sup>r</sup>*1. *<sup>ψ</sup>*(*r*) describes the autocorrelation functions, i.e., (*ϕ*(*r*1) <sup>−</sup>*ϕ*(*r*1))(*ϕ*(*r*2) <sup>−</sup>*ϕ*(*r*2)); and *<sup>ψ</sup>*(*k*) depicts the Fourier transform of *ψ*(*r*). *α*2(*r*; *R*) is the scaled intersection volume [37] of two spherical windows with radius *<sup>r</sup>* and *<sup>R</sup>* and *<sup>α</sup>*2(*k*; *<sup>R</sup>*) is its Fourier transform. One can call that *<sup>ϕ</sup>* is a hyperuniform state if its autocorrelation function subjects to the condition

$$\int\_{\mathbb{R}^d} \psi(\mathbf{r}) d\mathbf{r} = 0. \tag{14}$$

#### **4. Results and Discussion**

#### *4.1. Spatial Patterns*

Allen–Cahn model reveals a continual coarsening behavior, where cluster size become bigger and bigger and converge to one of its stable states, as shown in Figure 1. The system only displays a characteristic spatial scale at certain stages. A scale-free behavior appears at early stages. On the contrary, Cahn–Hilliard model displays a clearly characteristic spatial scale (Figure 2), and this scale follows a power law relationship with time. They are indistinguishable spatial patterns from Turing patterns (cf. Figures 2 and 3). It is notable that CH model describes the mass-conservation ecological processes that remarkably differ from the Allen–Cahn model despite both models reveal a coarsening behavior. The CHPD model displays classical Turing patterns with a stationary spatial scale when we consider a species birth–mortality process (Figure 3). In sum, it is difficult to infer the underlying ecological processes and potential mechanism from the merely observed spatial patterns.

**Figure 1.** The patterned dynamics of Allen–Cahn model in Equation (1) solved with the simple Euler algorithm starting a random initial condition: (**a**–**d**) four typical patterns at *t* = 0, *t* = 500, *t* = 5000, and *t* = 50, 000. Numerical simulation was implemented on the discrete 2048 × 2048 lattices with Δ*x* = Δ*y* = 1.0 and Δ*t* = 0.0025. It can be observed that the AC model has a relatively scattered wavelength on the spatial scales. Parameter values for the simulation: *α* = 0.5, *β* = 1.0, and *ε* = 0.5.

**Figure 2.** The patterned dynamics of Cahn–Hilliard model in Equation (2) solved with the simple Euler algorithm starting a random initial condition: (**a**–**d**) four typical patterns at *t* = 0, *t* = 5000, *t* = 50,000, and *t* = 500,000. Numerical simulation was implemented on the discrete 2048 × 2048 lattices with Δ*x* = Δ*y* = 1.0 and Δ*t* = 0.0025. It can be observed that the CH model has a more concentrated wavelength on the spatial scale than the AC model, which seems to be more regular on the time scales. Parameter values for the simulation: *α* = 0.5, *β* = 1.0, *ε<sup>A</sup>* = 0.5, and *ε<sup>B</sup>* = 1.0.

**Figure 3.** The patterned dynamics of Cahn–Hilliard with population demographics model solved with the simple Euler algorithm starting a random initial condition: (**a**–**d**) four typical patterns at *t* = 0, *t* = 5000, *t* = 50,000, and *t* = 500,000 . Numerical simulation was implemented on the discrete 2048 × 2048 lattices with Δ*x* = Δ*y* = 1.0 and Δ*t* = 0.0025. It can be observed that the pattern of the Cahn–Hilliard with population demographics model is dense at the initial moment, and the dots gradually become spatial regular patterns. The spatial patterns approximate a stable state after the critical time scales. At this time, the coarsening processes is suppressed by the population mortality and birth processes. Parameter values for the simulation: *α* = 0.5, *β* = 1.0, *ε<sup>A</sup>* = 0.5, *ε<sup>B</sup>* = 1.0, *r* = 0.0001, and *K* = 0.3.

#### *4.2. Spatial Correlation Functions*

Figure 4 illustrates the spatial correlation of the three models from *t* = 1 to *t* = 10,000. The spatial correlation functions reveal differently dampened oscillation behaviors at long time scales, where the location of minimal values increases with time (Figure 4). It is consistent with the coarsening processes of the spatial patterns. However, this behavior is not obvious in AC model, which means that the pattern has no characteristic scale. In sum, the spatial correlation function is not a significant indicator to distinguish three models.

#### *4.3. Structure Factors*

To investigate the characteristic wavelength of the patterns, we plot the structural factor *s*(*k*, *t*) versus wave numbers *k* at different times in Figure 5. A larger value of *s*(*k*, *t*) means a stronger dominant effect of its wave numbers.

For AC model, there is no dominant wave numbers of *s*(*k*), implying absence of spatial characteristic length on this system. It displays a power law decline with exponent of −3.12, predicting scale-free patterns (Figure 5d). However, *s*(*k*) has obvious dominant peak in the CH and CHPD models, indicating the regular spatial patterns. Furthermore, with time, the peak shifts to a smaller wave number, indicating that the spatial characteristic length is gradually increasing. For CHPD model, there is no significant change in the structure factor after a critical time scale (here 50,000), indicating that the spatial patterns do not change over time. This saturated behavior indicates the logistic growth term is playing a leading role (the phase separation is suppressed). In sum, the structure factors can serve as an excellent indicator to distinguish AC and CH/CHPD model.

**Figure 4.** Spatial correlation functions for the self-organized patterns at different times: (**a**) Allen–Cahn model; (**b**) Cahn–Hilliard model; and (**c**) Cahn–Hilliard with population demographics model. Different colors indicate the different times.

**Figure 5.** (**a**–**c**) The structure factors *s*(*k*) of the spatial patterns at different times. (**d**–**f**) The scaled structure factors *s*(*k*)*k*<sup>2</sup> <sup>1</sup> versus *k*/*k*<sup>1</sup> at different times. Different colors indicate the different times.

Moreover, all systems display a universal scaling behavior at high wave numbers (i.e., small wavelength). They collapse to master curves when the structure factors were normalized by *k*/*k*<sup>1</sup> for different times [38], as shown in Figure 5d–f. Unfortunately, to date we do not know the ecological significance of those slope values.

#### *4.4. Density Fluctuation*

The exponent of density fluctuation can be used to quantitatively measure the spatial distribution of the pattern [36,37]. As mentioned in the Method Section, we plot the variable-field fluctuations as a function of the window size *L* in Figure 6a. As expected, the exponent of density fluctuations predicts remarkable differences between AC model and CH/CHPD model despite it being impossible to directly distinguish from spatial patterns (cf. Figures 1–3). The exponent of the AC model is around −1.91, implying a random distribution of patterns at large-spatial scales. Interestingly, CH and CHPD models have the exponents of −2.90 and −2.60, respectively. Generally, this exponential value of the scaling law is used to describe the spatial-ordered features. It is called "hyperuniformity" when the exponent of the density fluctuation is less than −2.0, where the long-wavelength density fluctuations are suppressed [39]. In sum, the exponent of density fluctuation is also useful to quantify the spatial patterns between AC model and CH/CHPD models.

**Figure 6.** Spatiotemporal hallmarks of three different dynamics models (Allen–Cahn model -, Cahn–Hilliard model , and Cahn–Hilliard with population demographics model ). (**a**) The variable-field fluctuations as a function of the window size length *L* at *t* = 500, 000 associated with 20, 000 × 20, 000 systems. (**b**) The scaling behavior of the spatial wavelength *Rs*(*t*) versus increased times.

#### *4.5. Growth Law*

To measure how the spatial structure of the patterns evolves over time, we calculate the characteristic length, *Rs*(*t*), with location of the maximum of *s*(*k*, *t*) and its first moment at time *t*, i.e.,

$$k\_1(t) = \frac{\sum\_k k \, s(k, t)}{\sum\_k s(k, t)}.\tag{15}$$

Then, the spatial characteristic length (also called wavelength) change over time can be obtained as

$$R\_s(t) = \frac{2\pi}{k\_1(t)}.\tag{16}$$

The scaling behaviors of the spatial wavelength *Rs* (*t*) are shown in Figure 6b. The AC and CH models reveal power law relationship with exponents around 0.37 and 0.32, respectively. It should be noted that this scaling property is called the Lifshitz–Slyozov–Wagner law [40]. However, CHPD model has phase separation patterns at an early stage with scaling of 0.19, but it saturates at long-time scales. These patterns are equivalent to the Turing principle. In sum, this saturated behavior of growth rate can well distinguish CH model and CHPD model, which are impossible to be separated by structure factor, density fluctuation, and spatial correlation function (see Table 1 for summary).


**Table 1.** Summary of the different indicators on three mechanism models.

#### **5. Conclusions**

Although PDE models have been widely used to describe the spatial patterns in many ecological and biological systems, few studies focus on inferring the linkage between these visual patterned characters and their potential mechanisms [41–43]. To do this, one of the effective ways is to build the quantitative indicators of these macrocosm patterns, which can serve as the hallmarks of the underlying ecological processes. Thus, it is a key issue to explain the ecological mechanisms by comparing the temporal and spatial behavior of the spatial patterns [44–46].

Here, we use structure factors, density fluctuation, growth law, and saturation status (Table 1) to quantify the spatiotemporal dynamics of AC model, CH model, and CHPD model. There is a remarkable peak of structure factors in the CH model and CHPD model, but absent in AC model. The differences of spatial structures between the AC model and CH/CHPD model can be well distinguished by the exponent of density fluctuations. Both CH and CHPD models reveal a hyperuniform behavior with scaling of −2.90 and −2.60, respectively, but AC model displays a random distribution with scaling of −1.91. The differences in time scales among the patterns of the three models can be well distinguished by LS behaviors. Both AC and CH model display the LS behaviors with scaling of 0.37 and 0.32, respectively, but CHPD model with scaling of 0.19 at short-time scales and saturation at long-time scales.

We showed that the underlying mechanisms can well be identified by the spatial information, i.e., the structure factors, the density fluctuation, the growth law, and the saturation; however, seeking real ecological cases is still a challenge (see Table 2 for summary). To the best of our knowledge, the first indicator (structure factors) can be applied to calculate the spatial wavelength, with only spatial data. The second indicator (density fluctuation) can distinguish the random distribution and the hyperuniform state at large spatial scales. It also reveals the intensity of long wave suppression. The third indicator (growth law) can be obtained from the spatiotemporal observed sequences. It describes the change and the stability of spatial structure with time. There are still many other potential mechanisms to generate spatially self-organized patterns in the ecosystems [47,48] such as animal aggregation due to taxis and density-dependence [49], but they are beyond the scope of the discussion here. It is interesting for further comparison in the future research.

To sum, we find that the density-fluctuation exponent and LS behaviors can quantitatively infer the three underlying mechanisms from the spatiotemporal patterns. We need to integrate multiple indicators to distinguish the underlying ecological processes and mechanisms, such as Turing principle and phase separation principle.


**Table 2.** A summary of the three mechanisms on spatial self-organization in ecosystems.

**Author Contributions:** K.Z. performed research; W.-S.H. conceived the numerical analysis; Q.-X.L. constructed research; and K.Z., W.-S.H., and Q.-X.L. wrote the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Natural Science Foundation of China (41676084) and the National Key R&D Program of China (2017YFC0506001).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**


*Mathematics* **2020**, *8*, 112

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
