## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# Identifying the dynamics of complex spatio-temporal systems by spatial recurrence properties

Edited by H. Eugene Stanley, Boston University, Boston, MA, and approved March 25, 2010 (received for review September 22, 2009)

## Abstract

Complex spatio-temporal systems may exhibit irregular behaviors when driven far from equilibrium. Reaction-diffusion systems often lead to the formation of patterns and spatio-temporal chaos. When a limited number of observations is available, the reconstruction and identification of complex dynamical regimes become challenging problems. A method based on spatial recurrence properties is proposed to deal with this problem: generalized recurrence plots and generalized recurrence quantification analysis are exploited to show that detection of structural changes in spatially distributed systems can be performed by setting up appropriate diagrams accounting for different spatial recurrences. The method has been tested on two prototypical systems forming complex patterns: the complex Ginzburg–Landau equation and the Schnakenberg system. This work allowed us to identify changes in the stability of spiral wave solutions in the former system and to analyze the Turing bifurcations in the latter.

- complex Ginzburg-Landau equation
- spatial recurrence plots
- spatially distributed systems
- Turing patterns
- reaction-diffusion systems

Reaction-diffusion equations are widely used for modeling a large variety of observed natural and physico-chemical spatial phenomena, such as chemical reactions (1), biological systems and population dynamics (2), and pattern formation in physical systems (3, 4). In fact, when driven outside from equilibrium, spatially extended systems may exhibit irregular behavior both in space and time. This phenomenon, commonly known as spatio-temporal chaos (5, 6), consists in the spontaneous emergence of spatial patterns like Turing structures, traveling and spiral waves, and turbulence. Reaction and diffusion have been widely investigated both from the theoretical (7, 8, 9) and experimental viewpoints (10).

One important and still open issue arises in studying an unknown dynamical spatio-temporal system when only the observation of one or a part of its spatial variables (or a combination of them) is known and observation data is poorly available. The problem of space-state reconstruction and model identification of a spatio-temporal dynamical system has been investigated in lattice dynamical systems (11), whereas a method for spatial forecasting from single snapshots has been proposed in (12). In such cases one has to cope with the problem of understanding the dynamics of a system by using only a limited amount of data.

This paper attempts to provide a methodology for extracting information on the dynamics of a distributed dynamical system from a limited number of spatial observations and, in particular, for detecting structural changes in its dynamics.

## Discussion

### Spatial Recurrence Strategies.

Let us consider the images shown in Fig. 1. We may consider these images as snapshots of an unknown two-dimensional system, and we want to extract as much information as possible about the underlying system. Indeed, they represent the numerical solutions at equal times of a spatially distributed system (which will be described later) for different values of its characteristic parameters. One can notice that some images have similar patterns, whereas others are substantially different.

Of course, the solution of the inverse problem is a formidable task, which cannot be solved in general, except in very special cases. For this reason, we focus on a more limited scope: to devise a systematic methodology based on experimental data to detect changes in the system dynamics and in its qualitative behavior.

We start by considering that the formation of a pattern in a certain region depends on the spatio-temporal dynamics of the system. Hence, spatial state recurrences may give insight in the evolution of the system.

Indeed, the concept of recurrence is strictly related to that of dynamical systems, as originally stated by Poincaré. For time series, temporal recurrence of states has been widely investigated by means of recurrence plots and recurrence quantification analysis (13), whereas the problem of spatial recurrence in biomedical images was investigated by Marwan et al. (14), who developed the generalization for spatial systems, leading to generalized recurrence plots (GRP) and generalized recurrence quantification analysis (GRQA) (see *Materials and Methods* for details).

Recently, the authors of this paper proposed the use of GRP and GRQA for the analysis of spatially distributed systems by introducing the determinism-entropy (*D*-*E*) diagram, a tool accounting for the recurrence properties of images generated by a dynamical system (15). In this paper, we intend to go beyond this task and answer the following question: Is this framework capable of extracting all the useful information from a sequence of images for detecting structural changes in the dynamics of the generating system?

Following ref. 15, determinism is an appropriate measure accounting for the global structure of the image, whereas entropy represents a measure accounting for the patterns at small scales. In fact, the combination of both indicators in the determinism-entropy diagram provides a suitable description of the image at micro and macro scales. Table 1 reports the values of both indicators for the images of Fig. 1. One can notice that *D* ranges approximately from 10.6 to 27.5 and *E* from 1.8 to 3.2. We observe that the points corresponding to the pair (a),(b) are quite close in the *D*-*E* plane. The same holds for the pair (e),(f), whereas the pair (c),(d) shows values different from the other two pairs and in addition (c) is somewhat far from (d). Such differences may be related to structural changes (see ref. 16 for similar behavior in one dimensional systems). The hypothesis of abrupt changes in the recurrence indicators is consistent with the nature of the system used to generate the patterns; i.e., the Complex Ginzburg–Landau Equation (CGLE) [1]where , *x*, *y* are the spatial variables and (17, 18). This is one of the most important and studied nonlinear partial differential equations describing pattern formation in a spatially distributed system near a Hopf bifurcation. The interest in Ginzburg–Landau equations has been sustained because its solutions display a very rich spectrum of dynamical behaviors when the parameters are varied. This allows one to describe a large variety of physical systems, such as nonlinear waves, second order phase transitions, superconductivity and superfluidity, Bose–Einstein condensation and population dynamics (19). An important feature of the CGLE is the structure, nature, and role of defects, i.e. points in the space where |*A*| = 0. For two-dimensional systems, defects are points that appear and disappear in pairs, and, for small enough values of *α*, they are able to generate (stable or unstable) spiral waves (further details in *Materials and Methods*). The behavior of the solution in the parameter space (*α*,*β*) is very complex and still under investigation.

To give an answer to the problem posed at the beginning of this section, in the following we present a methodology based on recurrence properties of the images for a detailed analysis of the CGLE and reaction-diffusion systems.

### Detection of Structural Changes.

The behavior of the CGLE in the plane (*α*,*β*) has been extensively studied by Huber (17, 20, 21), Brito (22), and Chaté and Manneville (23, 24). On the basis of extensive simulations, they define a set of transition lines in the region *β* > 0. In particular, they identify the convective instability limit, the absolute instability limit, and the regions in which the phase turbulence and defect turbulence are observed.

In our investigation, a systematic analysis of the CGLE has been performed in a square domain of the parameter space *α* × *β* = [-1.5,1.5] × [-1.5,1.5]^{*}. The parameter plane has been sampled by varying *α* by steps of 0.25, and *β* of 0.1. Specifically, the real part of solutions *A*(*x*,*y*) of the CGLE has been considered for each pair (*α*_{i},*β*_{i}).

To guarantee the stationarity of the measures, we first look at the temporal evolution of both determinism and entropy starting from random initial conditions. Fig. 2 shows the dynamics of *D* for *α* = -1 and *β*∈{-1,0.1,1}. For each *β*, *D* reaches a different saturation value, after the transient. For *E* the same behavior has been observed. This analysis allows for the selection of an appropriate stopping time for the integration and guarantees the stationarity of both *D* and *E*.

A preliminary information on the variation of the GRQA measures with respect to *β* is obtained by introducing a sensitivity indicator defined by [2]where and are the averages of entropy and determinism, respectively. A small value of *K*(*β*) indicates that both *E* and *D* do not vary significantly with *β*. Fig. 3 reports *K*(*β*) for three different values of *α*: One can notice that, depending on *α* and *β*, the three curves show two clearly distinct behaviors characterized by low and high sensitivity.

High values of *K*(*β*) point out that small variations of *β* may produce strongly different patterns: This variability is related to the position and number of defects, allowing for the formation of structures with different spatial distribution. This is consistent with the fact that defects arise, collide, and annihilate in the unstable wave solution case. On the contrary, in the stable case, the number and location of defects is robust to small variations of the parameter *β*, giving small values of the sensitivity.

After this preliminary analysis, the values of *D* and *E* for each image are computed and reported in the *D*-*E* diagram. Fig. 4 shows the position of all the analyzed images in the diagram. Observe that the images are mainly organized in two regions, indicated by *A* (*D* ∼ [8,14], *E* ∼ [1.6,2]) and *B* (*D* ∼ [20,39], *E* ∼ [2.6,3.9]). The two regions are separated by a narrow zone depicted in gray (*G*).

It is worth noticing that points in region *A* are much more condensed with respect to region *B*. This is explained by the very low sensitivity of the GRQA measures, as mentioned before (Fig. 3). The role of the transition zone *G* becomes clear when looking at Movie S1, showing the dynamics of the points in the *D*-*E* diagram as functions of *α* and *β*. One can notice a systematic behavior: The points (i.e., the steady state solutions of the CGLE for a fixed pair (*α*,*β*)) are confined in one of the two regions (*A* or *B*), then move toward the other region by crossing the transition zone, and remain inside it. For example, consider the value *α* = 0. For increasing *β*, the points initially belong to *A* and remain in *A* until a critical value is reached; after this, the points cross *G* and jump in *B* for a new value . By increasing *β*, the scenario is observed again in the opposite direction: Starting from *B*, the points remain inside *B* until a critical value is reached, then they cross again *G* and jump in *A* after the value . The arrows in Fig. 4 indicate the directions along which the motion of the points for increasing *β* takes place.

The behavior shown in Movie S1 and the regions with high and low sensitivity with respect to *β* suggest that it is possible to detect structural changes in the dynamics of the CGLE by looking at the region jumps. For each value of *α*, following the evolution of the points in the *D*-*E* diagram, we collect the values for which the transitions are observed and we report them in Fig. 5, where the blue circles correspond to the values of *β* for which the jumps *G*–*A* and *A*–*G* are observed, whereas the red diamonds correspond to the values of *β* for which the jumps *G*–*B* and *B*–*G* take place.

One can notice that the blue circles lie with good agreement on the Eckhaus stability limit (the line *S*_{1} in paper (23). The region below the Eckhaus limit is mapped in the region *A* of the *D*-*E* diagram. Analogously, the unstable region of the parameter space is mapped to the region *B*, and the transition zone is mapped to the light gray zone *G* separating *A* and *B*. By looking at the values of *β* for which the jumps *B*–*G* and *G*–*B* happen, we identify with very good accuracy the transition region that was previously observed for some values of *α* and *β* by Huber et al. (21) and by Bohr et al. (25). They analyzed the onset of turbulence in the CGLE for finite domain size calling it transient turbulent behavior.

#### Reaction-diffusion systems

In this subsection the application of the proposed method to reaction-diffusion equations is discussed. Recently, the authors of the present paper have shown in (26) that determinism and entropy are suitable indicators for detecting different routes to Turing pattern formation in a Belousov–Zhabotinsky reaction performed in a water–oil microemulsion (27). Here we focus on one of the prototypical reaction-diffusion systems showing Turing bifurcations: the Schnakenberg model (28) (see *Materials and Methods* for mathematical details). Referring to Eq. **9** and according to ref. 29, we choose *γ* = 800, *k*_{2} = 1, whereas *k*_{1} varies in the range [0.1,0.6]. As expected, the Turing bifurcation can be easily detected by computing the recurrence indicators *D* and *E*. In fact, for diffusion coefficient smaller than the critical value *d*_{c}(*k*_{1},*k*_{2}), the solution converges to a homogeneous state yielding the maximum feasible values of *D* and *E*. Fig. 6 shows the evolution of both *D* and *E* for *k*_{1} = 0.1 and 9.8 < *d* < 15. Notice that near the critical value *d*_{c} ∼ 10 an abrupt change in both values of *D* and *E* indicates the occurrence of the bifurcation. Furthermore, for *d* > *d*_{c}, both indicators show some variations before reaching the saturation values *D*^{∗}(*k*_{1} = 0.1) = 20.64, *E*^{∗}(*k*_{1} = 0.1) = 1.868. This may be explained by the formation of transient patterns (Fig. 6, *Insets*) reported in ref. 29. Let us now collect the saturation values *D*^{∗}(*k*_{1}) and *E*^{∗}(*k*_{1}) for *k*_{1}∈[0.1,0.6]. Fig. 7 shows that *D*^{∗}(*k*_{1}) increases quadratically, whereas the same behavior is not observed for the entropy, whose regime value *E*^{∗}(*k*_{1}) ∼ 1.8 is independent of *k*_{1}. The behavior of both *D*^{∗} and *E*^{∗} can be explained by considering the Turing patterns: The increase of determinism may be related to the decrease of the spatial frequency of the patterns (Fig. 7, *Insets*). On the other side, entropy does not show important changes because the information does not vary significantly by reducing the spatial frequency.

## Conclusions

Reconsidering the experiment proposed at the beginning of the paper, we recognize that the pair of images (a),(b) falls in region *A* of Fig. 1, and the pair (e),(f) falls in region *B*, whereas the pair (c),(d) falls in region *G*. Then, with reference to Fig. 5, we can conclude that (a),(b) come from stable spiral wave solutions of the CGLE and (e),(f) come from unstable spiral wave solutions. Notice that, despite the different visual aspect, (a) and (b) refer to the same dynamical structure. This conclusion holds also for (e) and (f) but with different levels of determinism and entropy. Following the analysis performed in the paper, this indicates that the system underwent a transition by changing the parameter *β*. Moreover, the transition giving rise to images (c) and (d), takes place according to a path that is clearly visualized in Movie S1. The transition from (a),(b) to (e),(f), which has been clearly detected through the *D*-*E* diagram analysis, is not easily detectable by a visual inspection [e.g., notice that image (b) appears visually more different from (a) than from (d)]. It is worthwhile to remark that using different classical methods for measuring image complexity (like, for example, pixel-based entropy (30) fails to discriminate the different pattern structures in Fig. 1 (*A*)–(*F*) (15).

Summarizing, we have shown that the generalized recurrence quantification indicators are suitable for characterizing the spatial patterns of the complex Ginzburg–Landau equation. Specifically, following the evolution of the image positions in the *D*-*E* diagram, we were able to identify, for each value of *α*, a set of threshold values of *β* for which a region jump takes place. By collecting these values, it was possible to reconstruct, with good agreement, the Eckhaus limit of convective instability of spiral waves in the (*α*,*β*) plane. A transition zone of transient turbulence was reproduced extending the preceding results.

Regarding reaction-diffusion systems, we have considered the Turing patterns formed by the Schnakenberg model. As in the CGLE case, the quantification of spatial recurrences allowed for the detection of the different regimes observed in the pattern formation. In particular, the determinism was found to increase quadratically with the parameter *k*_{1} of Eq. **9**. Moreover, the transition from homogeneous to fully formed patterned states was easily detected.

The two examples addressed may be considered as prototypes for covering a wide range of phenomena. Furthermore, the technique proposed can be usefully applied either for detecting structural changes in unknown systems or for uncovering bifurcations in dynamical spatio-temporal systems, whose complexity prevents the application of classical bifurcation analysis (8, 31).

Indeed, there are several application areas where the proposed approach could help to understand the dynamical features of biological and physical phenomena from sporadic spatial data, such as spatial ecology (32), plankton turbulent patterns (33), physiology of tissues (2), and functional MRI (34). Further applications of the proposed method to these phenomena are the subject of ongoing work.

## Materials and Methods

### Complex Ginzburg-Landau Equation.

In the following only basic information is provided (for an exhaustive treatment of the CGLE the reader is referred to ref. 18). The CGLE reads [3]where and . The first term of the right hand side is related to the linear instability mechanism leading to oscillation. The second term accounts for diffusion and dispersion, whereas the cubic term insures, for *β* > 0, the saturation of the linear instability and is involved in the renormalization of the oscillation frequency. When *α* = 0 and *β* → ∞ one has the real Ginzburg–Landau Equation, which possesses a Lyapunov functional and exhibits relaxational dynamics (35). For *α* → ∞ and *β* = 0 one recovers the nonlinear Schrödinger equation. In two dimensions, the solutions of the CGLE are families of plane waves *A* = *a*_{k} exp *ı*(*kx* + *ω*_{k}*t*) with and *ω*_{k} = 1/*β* - (*α* + 1/*β*)*k*^{2}, where *k* is the wave number.

We have integrated the equation using a pseudospectral code (36), in a square domain of *L* = 512 points with periodic boundary conditions^{†}.

### Generalized Recurrence Plots.

The recurrence plot (RP) (37) is essentially a two-dimensional binary diagram indicating the recurrences that occur in an *m*-dimensional phase space within a fixed threshold *ϵ* at different times *i*, *j*. For time series the RP is easily expressed as a two-dimensional square matrix with ones and zeros representing the occurrences (ones) or not (zeros) of states and of the system: , where , *i*,*j* = 1,…,*N*, *N* is the number of the measured states , Θ(·) is the step function, and ||·|| is a norm. In the graphical representation, each nonzero entry of **R**_{i,j} is marked by a black dot in the position (*i*,*j*). Because any state is recurrent with itself, the RP matrix fulfills **R**_{i,i} = 1 and hence contains the diagonal line of identity.

An RP is characterized by typical patterns, whose structure is helpful for understanding the underlying dynamics of the system investigated. A homogeneous distribution of points is usually associated with stationary stochastic processes; e.g., Gaussian or uniform white noise. Periodic structures, like long diagonal lines parallel to the line of identity indicate periodic behaviors, whereas drifts in the structure of the recurrences are often due to a slow variation of some parameter of the system and white areas or bands indicate nonstationarity and abrupt changes in the dynamics. For an extensive discussion of RPs and recurrence quantification analysis (RQA) measures the reader may refer to ref. 13. Recurrence plots may be exploited for the analysis of systems showing complex patterns in time and space. In ref. 14 the GRP for a *d*-dimensional data-set as the 2*d*-dimensional RP is defined by [4]where is the *d*-dimensional coordinate vector and is the associated phase-space vector. This GRP accounts for recurrences between the *d*-dimensional state vectors. The line of identity is replaced by a linear manifold of dimension *d* for which , .

The discretized solution of the two-dimensional CGLE at a fixed time can be represented by images (one for *ℜ*{*A*} and one for *ℑ*{*A*}); in this case the state or *ℑ*{*A*(*x*_{i1},*y*_{i2})} and the GRP reads as follows: **R**_{i1,i2,j1,j2} = Θ(*ϵ* - ||*x*_{i1,i2} - *x*_{j1,j2}||), *i*_{k}, , where each black dot represents a spatial recurrence between two pixels, and every pixel is identified by its coordinates (*i*_{1},*i*_{2}), being *i*_{1} and *i*_{2} the row and the column index respectively. In this case, the recurrence plot is a four-dimensional RP and contains a two-dimensional identity plane, defined by setting *i*_{1} = *j*_{1} and *i*_{2} = *j*_{2}.

Because the GRP of an image is four-dimensional, its visual inspection is possible only by projections in three or two dimensions. Although this is possible, relevant information is hard to extract, and one must cope with the fact that GRPs lose their visual appeal. Despite this drawback, RQA can still be performed because the structures described before can be easily extracted, and in the following we describe how to generalize the structures formed by the recurrences. A line structure of length *l* is composed by *l* consecutive recurrent points, according to the following formulation: [5]On the basis of the above definition, we build the histogram *P*(*l*) of the line lengths and define the GRQA measures as in the one dimensional case. In particular, we focus on recurrence rate (*R*), determinism (*D*) and entropy (*E*), defined as [6][7][8]where *l*_{min} is the minimum length considered for the diagonal structures.

The recurrence rate *R* is a density measure of the RP, accounting for the fraction of recurrent points in the spatial domain with respect to the total number of possible recurrences. Determinism (*D*) is a measure of the predictability of the system; in fact, it is the fraction of recurrent points forming diagonal structures in the RP with respect to all the recurrences. The entropy (*E*) is a complexity measure of the distribution of the diagonal lines in the GRP because it refers to the Shannon entropy with respect to the probability to find a structure of exactly length *l*.

The computation of the measures based on the diagonal lines and their distribution provides valuable information about the structure of the RP and the underlying structure of the solution under investigation. In particular, the determinism is the measure of the global pattern recurrence in the image: A value of determinism > 60–70% indicates that the image has strong recurrent components. In this sense, the measure fits the need to describe globally the patterns showed by the image. On the other side, the entropy provides a measure of the complexity of the GRP with respect to the diagonal structures: A low entropy indicates a poor organization of the line structures and is related to the small scale patterns.

### Schnakenberg System.

The Schnakenberg system describes a simple chemical reaction with limit cycle behavior (28). It has been widely used for investigating Turing instabilities in biological and ecological systems (2). The dimensionless system of equations reads [9]where *u*(*x*,*y*,*t*), and *x*, *y* are the spatial variables; *γ* is proportional to the spatial domain size and represents the relative strength of the reaction terms, *k*_{1} and *k*_{2} depend on the reaction rates and *d* is the ratio of the diffusion coefficients of the two reactants. The critical diffusion coefficient *d*_{c} is a function of the parameters *k*_{1} and *k*_{2}. For further details, the reader is referred to ref. 2, page 78.

The system has been integrated in a square domain of *L* = 250 points using a finite element method by means of `Comsol` `Multiphysics` software.

## Footnotes

^{1}To whom correspondence should be addressed. E-mail: mocenni{at}dii.unisi.it.Author contributions: C.M., A.F., and A.V. designed research, performed research, contributed new reagents/analytic tools, analyzed data, and wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.0910414107/-/DCSupplemental.

↵

^{*}We have also considered negative values of*β*because, according to ref. 23, this region of the parameter space is still not well investigated.↵

^{†}The analysis was performed by considering increasing domains: For sizes*L*= 256,512,1024 the GRQA measures did not change significantly.

## References

- ↵
- Epstein IR,
- Pojman J

- ↵
- Murray J

- ↵
- Cross M,
- Greenside H

- ↵
- Rutherford A,
- Aronson DG,
- Swinney HL

- ↵
- ↵
- Cross M,
- Hohenberg PC

- ↵
- Grindrod P

- ↵
- Mei Z

- ↵
- Banks H,
- Kunish K

- ↵
- Vanag V,
- Epstein IR

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Huber G

- ↵
- ↵
- Sherratt J,
- Smith M,
- Rademacher J

- ↵
- Cladis PE,
- Palffy-Muhoray P

- Huber G

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Vanag VK,
- Epstein IR

- ↵
- ↵
- ↵
- Gonzalez R,
- Woods R,
- Eddins S

- ↵
- Hirsch MW,
- Smale S

- ↵
- ↵
- ↵
- Huettel SA,
- Song AW,
- McCarthy G

- ↵
- ↵
- ↵
- Eckmann JP,
- Kamphorst SO,
- Ruelle D

## Citation Manager Formats

### More Articles of This Classification

### Physical Sciences

### Physics

### Biological Sciences

### Related Content

- No related articles found.

### Cited by...

- No citing articles found.