Seasonal footprints on ecological time series and jumps in dynamic states of protein configurations from a non-linear forecasting method characterization

Abstract

We have analyzed phenology data and jumps in protein configurations with the non-linear forecasting method proposed by May and Sugihara [1]. Full plots of prediction quality as a function of dimensionality and forecasting time give fast and valuable information about Complex Systems dynamics.

The notion of prediction is at the heart of the scientific method and of some theories about the mind (the predictive brain) [2]. In this article we explore decades of ecological data [3] and data from a model for sampling protein configurations using a dynamic indicator $\rho(E,T_{p})$ introduced some time ago by May and Sugihara to distinguishing chaos from noise in time series [1] : $\rho^{2}$ is a number in the unit interval which quantify the quality of the predictions made by the non-linear forecasting method introduced by the authors. We have obtained valuable insights about such data when considering full plots of the quantity $\rho^{2}(E,T_{p})$. Here $E$ is a dimensionality parameter: in the forecasting method is the number of components of vectors constructed with the data. We are trying to make a prediction $T_{p}$ time steps into the future and the evaluation of that prediction is made with the quantity $\rho$. Characterizing the dynamics of a system is a basic goal in today’s big data dominated scientific practice. In this work we have applied the method proposed by May and Sugihara to analyze data spanning decades of population dynamics of several species from a Venezuelan forest [3] and to analyze some aspects of protein dynamics [9].

Figure 1: $\rho^{2}(T_{p})$ for $E=3$, $E=7$ and $E=12$ (see text). For $E=3$ we get seasonal scales from the signal. For $E=12$ we get chaotic dynamics from the signal for the case of flowers (top), and an el niño coupling for fruits (bottom). $T_{p}$ is given in months. $E=7$ is the first case outside the small dimensionality region in parameter space.

Given some data series $x_{n}$ the nonlinear forecasting model considers vectors of $E$ components of the form $(x_{i},x_{i-1},x_{i-2},\ldots,x_{i-(E-1)})$ and constructs a prediction for the next $x_{i+T_{p}}$ value in the series; the prediction’s quality is quantified with the correlation $\rho$ between observed and predicted values. Thus, for given data $x_{n}$, the quantity $\rho(E,T_{p})$ is the basic output of the method [1].

In figure 1 we show $\rho^{2}(T_{p})$ for Hieronyma moritziana Pax & K. Hoffm species for three values of the dimensionality parameter $E$, $E=3$, $E=7$ and $E=12$, for signals of flowers and fruits populations as a function of time [3]. It can be appreciated that the seasonal time scale (one year) manifests itself at $E=3$, where we obtain a periodic pattern. For $E=12$ there is a decay in prediction quality as $T_{p}$ increases, which is a signature of chaotic dynamics [1].

Full plots of $\rho(E,T_{p})$ can be very informative about the dynamic state of a system. In figure 2 we show $\rho(E,T_{p})$ for Hieronyma moritziana species, and some features emerge: only for $E<12$ the seasonal scale appears, and there are triangular patterns that repeat themselves each $12$ months. For flowers and $E>12$ we obtain a decay in $\rho$ with increasing $T_{p}$, while for fruits and $E>12$ we can observe a recurrence of period around $60$ months, which suggest an “el niño” coupling [3]. Other species have no structure in parameter space (like figure 5) or some distorted triangular patterns, see figure 3.

Figure 2: Dynamics as characterized by $\rho^{2}(E,T_{p})$ for Hieronyma moritziana species: flowers (left) and fruits (right).

Figure 3: Three cases for population dynamics with very different predictability patterns in parameter space. Up: flowers, Bottom: fruits. See also [3].

For comparison and reference we show in figure 4 (left) $\rho^{2}(E,T_{p})$ plots for chaotic time series with different Lyapunov exponents [5]. It can be appreciated that the intersection points on the $E$ and $T_{p}$ axes are the same, and that the larger the Lyapunov exponent the smaller those intersection points. Similar triangular patterns were obtained for Hieronyma moritziana like species [3] but with a major difference (see figure 2): for several species we have obtained a kind of inversion and repetition of the triangular patterns obtained for chaotic series, which prompts the question: how can we obtain low prediction quality at short times and high prediction quality at longer times?. We have found that a random decay process with resinsertion [4] produces the same kind of triangular patterns obtained for Hieronyma moritziana like species, this result is shown in figure 4 (right). We have generated a random decay signal with the explicit functions of González et al: $X_{n}=\sin^{2}(\theta\pi z^{n})$, with parameter $z=p/q$, $p$ and $q$ being relative primes ($q>p$, $z<1$) [5]. A random decay with reinsertion sounds very ecosystemic indeed. We get poor prediction at short times because is a random process, and we obtain a high prediction quality at longer times because once the signal has decayed to near zero values the dynamic is very predictable. It can be observed in figure 2 that for Hieronyma moritziana and within one year we obtain a kind of superposition of two triangles, corresponding to two overlapping random decay processess: wet and dry periods of the year in the Venezuelan forest [3]. A simpler plot like the one shown in figure 4 (right) highlights the underlying proposed mechanism. Since we have random decay processess of size $6$, it can be said that the case $E=7$ (see figure 1) is the smaller dimensionality value $E$ outside the small dimensionality region in parameter space.

Figure 4: left: $\rho^{2}(E,T_{p})$ for chaotic dynamics with different Lyapunov exponents $\lambda$. From top to bottom: $z=1.6,\ 1.75,\ 2.0$, with $\lambda=\log z$ [5]. Similar results were obtained for the logarithmic map. right: $\rho^{2}(E,T_{p})$ for synthetic data emulating a random decay process with reinsertion. We have used $z=2/3$ in the explicit functions of references [5] for the random decay process. See text.

In figure 5 we verify, with the help of the GHWS model [6], that an ideal gas type of dynamics would produce very small values of $\rho^{2}$, for most $E$ and $T_{p}$ values. The signal used in this case was the activity in the system as a function of time [6, 7, 8].

Figure 5: $\rho^{2}(E,T_{p})$ for an ideal gas like dynamics. The processed signal was generated from the GHWS model [6] in a zone of parameter space that produces zero average correlations of activity (see [8], section $3.6$), which results in the ideal gas equation of state [7].

In figures 6 and 7 we show a $\rho^{2}(E,T_{p})$ plot for the atoms coordinates of a protein obtained with molecular dynamics simulations [9, 10, 11]. As an example, we have chosen a typical intrinsically disordered protein such as alpha synuclein associated with neurodegenerative Parkinson’s disease, which is known for its complex dynamic conformational behavior and tendency to form aggregates [12]. This protein is made up of three domains or regions known as the N-terminal, NAC, and C-terminal domains [13]. The signal to be processed in this case is the position of the atoms in one of the three regions of the alpha synuclein protein, for a given time. As time passes towards equilibrium the protein jumps between very different dynamic states: fully disarmed ($\rho\sim 0$ for most $E$ and $T_{p}$), rigid (high values of $\rho$ for most $E$ and $T_{p}$), a chaotic-like structure (triangular pattern), see figure 6, and several more complicated $\rho(E,T_{p})$ patterns. It is important to note that structural disorder has been highlighted recently for alfa synuclein wild type [14], and that we are associating predictability (large values of $\rho^{2}$) with rigidity through lost of effective degrees of freedom. A particular pattern emerges for an alpha synuclein dimer scenario: we obtain a particular dimension value $E_{0}$ for which no prediction can be made about the position of the atoms along the polymer chain, for any value of $T_{p}$, see figure 7 (left). Regarding $E_{0}$, we have interpreted this result as implying that there are structures of lenght $E_{0}$, involving $E_{0}$ number of atoms, that appear in the protein (for a given time).

Figure 6: Protein dynamics, for data obtained with molecular dynamics simulations [9, 10, 11]. On the left a dynamic state, at a particular time, for which $\rho^{2}\sim 0$ for most $E$ and $T_{p}$: the protein is disarmed. At the center a dynamic state, at a particular time, for which $\rho^{2}$ is large for most $E$ and $T_{p}$: the protein is rigid. To the right we show a chaotic like structure obtained for a particular time (see figure 4 (left) and text).

Figure 7: left: $\rho^{2}(E,T_{p})$ from protein dynamics, a particular dynamic state for data obtained with molecular dynamics simulations: $\rho^{2}\sim 0$ for $E=E_{0}\sim 12$, for any $T_{p}$ considered. right: $\rho^{2}(E,T_{p})$ for $\alpha$-sinuclein protein with mutation A53T, after $100ns$. The input signal was the $y$ coordinate of the atoms in the N domain.

So far we have shown transient protein dynamics states. In figure 7 (right) we show $\rho^{2}(E,T_{p})$ for the N region of a mutated protein (A53T) after $100ns$ of simulation. The N-terminal region of the alpha synuclein protein was chosen given its important biological function of binding to the cell membrane [15, 14]. No such state were achieved for the same protein in the wild type case, and this particular dynamic state, a near equilibrium one, fits a chaotic structure scenario: good predictability for small $T_{p}$ and eventually unpredictability for large $T_{p}$, but with an approximately horizontal frontier between large and near zero values of $\rho^{2}$. These results will be further investigated elsewhere, but let’s consider the case of an horizontal frontier (in parameter space) between high and low values of $\rho^{2}$. If we look at figures 1 and 2, we can see that when $T_{p}$ is a multiple of $12$ months $\rho^{2}$ shows a weak dependence on dimensionality $E$, which in turn implies that variations of $\rho$ in parameter space (see below) points only in the $T_{p}$ direction. This suggest that an approximately horizontal frontier like the one shown in figure 7 (right) could be associated with an interaction with large parts of the protein (the equivalent of seasonal scales in this case). Under these considerations we can say that the mutated protein gains some (pathological) rigidity because of self-interaction. Besides, it appears that this pathology could be associated to a lack of cost to make predictions, since the dimensionality parameter $E$ can be directly related to memory capacity which certainly it is expected to have a “cost”.

Figure 8: For data obtained from the logistic map we show the vector field $\nabla\rho$ (left) and the modulus of the gradient of $\rho$ in parameter space $|\nabla\rho|$ (right). See text.

Complex adaptive systems. If we were interested in places of parameter space where adaptation could be done more efficiently, by changing the dynamic state of the system with minimal variation of parameters, then our ideal spots in parameter space would be places where the gradient (in parameter space) of $\rho$ is large. As an example, in figure 8 we show the modulus and direction of $\nabla\rho$ for the logistic map. We have estimated $\nabla\rho$ from finite differences of parameters values.

From the triangular patterns found for Hieronyma moritziana like species in the case of phenology (figure 2) and for chaotic dynamics (figure 4 (left)) it is clear that lines of slope $-1$ play a role in parameter space $(E,T_{p})$. The relation $T_{p}=E-1$ is significant in the method introduced by May and Sugihara: it limits the amount of information that in principle could be found in vectors of the form $(x_{i},x_{i-1},x_{i-2},\ldots,x_{i-(E-1)})$ when trying to make a prediction about $x_{i+T_{p}}$. The lines of slope $-1$ are perpendicular to lines of the form $T_{p}=E-1$, which suggests that the direction in parameter space given by $\nabla\rho$, as can be seen in figure 8, is very relevant in this method, which also highlights the natural framework it provides for considering adaptation in Complex Systems, as commented in the previous paragraph.

We have presented results for several types of system-dynamics based on the forecasting method introduced by May and Sugihara [1]. Our results highlight the powerfulness of full plots of $\rho(E,T_{p})$ obtained with the nonlinear forecasting method to analyze and characterize the dynamics of Complex Systems, for scales as large as forest Ecosystems and as small as Proteins.