Predictability of the onset of spiking and bursting in complex chemical reactions

For three complex chemical reactions displaying intricate dynamics, we assess the effectiveness of a recently proposed quantitative method to forecast bursting and large spikes, i.e. extreme events. Specifically, we consider predicting extreme events in (i) a copper dissolution model where Bassett and Hudson experimentally observed homoclinic (Shilnikov) chaos, (ii) a model derived from the mass action law of chemical kinetics, and (iii) an autocatalator model. For these systems, we describe how the alignment of Lyapunov vectors can be used to predict the imminence of large-amplitude events and the onset of complex dynamics in chaotic time-series of observables.


Introduction
2][3] Such extreme events are usually quite disruptive and manifest themselves as large-amplitude fluctuations in time-series of observables, either measured or computed.A plethora of methods to anticipate the onset of catastrophic events exist.Recently, we proposed a quantitative criterion to predict larger-amplitude events in chaotic systems. 4,5he criterion is obtained by studying the alignment of Lyapunov vectors (LVs), sometimes referred to as ''covariant'' Lyapunov vectors.The norms of the LVs, called Lyapunov exponents (LEs), measure the temporal expansion (or contraction) of the distance between two initially very close trajectories.Positive LEs are the signature of chaotic trajectories.We observed empirically 4,5 that a tangency between LVs occurs every time that an observable under consideration tends towards its maximum.This tangency is an alignment of Lyapunov vectors.In some cases the value of the maximum is inversely proportional to the degree of the vector alignment.
In the present paper, as a further application of the proposed criterion, 4,5 we investigate the reliability of the alignment between Lyapunov vectors as a predictive tool in the context of complex chemical phenomena.Specifically, we use Lyapunov vectors to predict the onset of bursting and spiking in three familiar reactions.The first reaction was considered in pioneering work in 1988 by Bassett and Hudson 6 to reproduce aspects of the experimental electrodissolution of a copper rotating disk in a H 2 SO 4 /NaCl solution.The second reaction is a mass-action model used to study the interplay of three interacting chemical species. 7,8The third reaction is a popular autocatalator model used to study complex oscillations in isothermal chemical systems. 9In all three examples, we find the alignment between Lyapunov vectors to be an effective predictor of extreme events.In addition, we point out some pitfalls of the method and discuss how to apply it effectively.
Lyapunov vectors provide a detailed step-by-step record of what exactly happens with the angles between stable and unstable manifolds during the whole evolution of a given dynamical process.1][12][13][14][15] In an Appendix, we summarize the steps involved in their computation.As discussed below, it was observed heuristically 4,5 that the onset of large peaks in a given physical variable corresponds to the alignment of LVs along the flow direction.Therefore, since Lyapunov vector alignment precedes large peaks, it can be used to predict the latter.It has been shown 4 that predictability using LVs is closely related to alignment conditions, obtained by defining thresholds of how close LVs become tangent to the flow direction.When the alignment conditions are satisfied at a given time, called the prediction time, an extreme spike is expected to occur.The appropriate definition of the alignment conditions depends sensitively on the specific dynamics of the system.
For systems with long quiescent intervals between large spikes, the alignment conditions are very useful to anticipate unexpected large spikes.However, when numerous smaller spikes occur between large spikes, the choice of the appropriate alignment conditions can be tricky and some preparatory trial-error checks may be needed.In general, while small and large spikes can be predicted by the alignment of LVs, large spikes correspond to the prediction times where LVs are strongly parallel to the flow direction.Thus, the amplitude of the predicted peaks is related to how close the LV direction aligns along the flow direction. 4part from the chemical phenomena described by the models, complex chemical reactions are governed by sets of differential equations which are significantly different from the equations so far considered and no knowledge concerning the interplay of their Lyapunov vectors exists.Thus, it is of interest to verify the predictive power of Lyapunov vectors in foretelling the imminence of large-amplitude events in chemical time-series.Here, we check predictions for three typical systems governed by nonlinear differential equations involving different combinations of polynomial terms.We computed LVs using the methodology explained in detail in ref. 4

Electrodissolution of copper
In 1988, Bassett and Hudson 6 performed a series of experiments on the electrodissolution of a copper rotating disk in a H 2 SO 4 /NaCl solution.They reported evidence of homoclinic (Shilnikov) chaos in their experimental time-series.Among this evidence, they presented a comparison of their time-series with those obtained from a mathematical model considered previously by Argoul et al., 16 in a discussion concerning homoclinic chaos in the Belousov-Zhabotinski reaction.The model is defined by the following three-dimensional flow Here, (x, y, z) represent chemical concentrations and ( : x, : y, z ˙) are their corresponding time derivatives.For m o 1.3 the origin is a stable steady state and for m = 1.3 a subcritical Hopf bifurcation occurs where the stable state becomes unstable and Shilnikov chaos appears.
Fig. 1 shows the chaotic attractor for m = 1.4 and initial conditions (À0.1, À0.1, 1.0).The attractor is obtained by solving eqn (1)-(3) numerically with the standard fourth-order Runge-Kutta algorithm and a fixed time-step h = 0.005.The first 10 7 steps were discarded as transient, with LVs computed subsequently for approximately 9.6 Â 10 6 steps.The Lyapunov exponents l i (i = 1, 2, 3) are l 1 = 0.234, l 2 = À0.0013, and l 3 = À2.739,where the subindex i = 1 refers to the unstable direction, i = 2 to the flow direction and i = 3 to the stable direction.Accordingly, y 12 , represented in colors in Fig. 1, is the angle between the unstable manifold and the flow direction.We observe that along the attractor the angle y 12 can change substantially.
In order to get a more quantitative insight concerning the ability of the alignment of LVs to anticipate spikes, Table 1 summarizes some statistical results.For this end, we consider the chaotic trajectory from Fig. 1 which has 2192 spikes over the integrated times, and say that there is a spike whenever x o À9.0 and an extreme spike when x o À15.0.To explain the quantities in Table 1, let us start by discussing the third line, which is the alignment condition displayed in Fig. 2. From the sample of 2192 spikes, N pp = 1693 of them were predicted, which is around 76.3% of the total number of spikes.However, N fa = 686 (41.0%) from the 1693 predicted spikes are false alarms, meaning that 1007 spikes are not extreme.False alarms are those alarms which predicted spikes with small amplitudes, i.e. predicted non-extreme spikes.From the total number 2192 À 1693 = 499 of non-predicted spikes, N ext = 181 of them were extreme spikes.Thus, for this alignment condition the method did not predict 181 extreme spikes.Such a large value of non-predicted extreme spikes has to be avoided.All predicted extreme spikes have an associated prediction time interval Dt p , which can be large or short.The prediction time interval is defined as the absolute value between the prediction time and the time for which the extreme spike occurs.While large prediction time intervals are most welcome to allow preparation to meet the extreme events, short prediction time intervals can bring difficulties for such preparations.In other   words, for very short prediction time intervals there may not be enough time to prepare for the impingent extreme events.
To quantify this, the last row displays the mean prediction time interval Dt p ¼ 15:2.
To better understand the relevance of y 12 , Fig. 2 shows a time-window displaying the variable x (black curve) from the chaotic attractor from Fig. 1.The angle y 12 is plotted in the red color.Apart from some small amplitude oscillations, y 12 tends to 0 or p before spikes appear in variable x.This means that the peaks in x can be predicted in time if we define thresholds for y 12 , namely y The overall message from Table 1 is that for stronger alignment conditions the number of non-predicted spikes increases, specially the extreme spikes, N ext .If the alignment conditions are relaxed, the number of predicted peaks increases but the number of false alarms also increases.Comparing the stronger alignment conditions with the more relaxed condition, N pp increases by around 65.2% but N fa increases by only 5.2%.In fact, more relaxed alignment conditions lead to longer prediction time intervals.The selection of a suitable threshold depends very much on what is needed for practical purposes.For example, if we want to avoid all extreme spikes, we should chose an alignment condition consistent with N ext -0.In this case the appropriate alignment condition would be y (min) = 2.44.On the other hand, if the large peaks provoke no major damage, one may accept some such peaks together with the condition to prevent a repetition of false alarms.This means, a compromise must be found between allowing some large peaks while at the same time not allowing the number of false alarms to dominate.To find such a compromise, one may consult Table 1 and select the suitable alignment condition providing an adequate balance.In such case, an adequate alignment condition would be y (min) 12 = 0.20 and y (max) 12 = 2.94.

Mass-action kinetics
In the next example, we revisit a chemical model introduced by Gaspard and Nicolis 7,8 to study homoclinic phenomena.The dimensionless governing equations are ez ˙= x À az 3 + bz 2 À cz.Accordingly, the angle y 12 is the angle between the unstable manifold and the flow direction.The trajectory from Fig. 3 and all subsequent simulations in this section were obtained by solving eqn ( 4)-( 6) numerically with the standard fourth-order Runge-Kutta algorithm with a fixed time-step h = 0.00005.The first 10 7 steps were discarded as transient, with LVs computed subsequently for approximately 10 7 steps.
The trajectory forming the attractor in Fig. 3 starts parallel to the xy plane with a chaotic spiraling-out motion, with the angles y 12 oscillating between p/4 (cyan) and p/2 (green), while the reinjection loops are shown in black and red, associated with y 12 approaching either 0 or p, respectively.To better understand the time evolution of y 12 , Fig. 4 presents the evolution of z (black lines) and the associated LV angle y 12 (red lines) for a regular trajectory in Fig. 4(a) and for a chaotic trajectory in Fig. 4(b).By 'regular' trajectory we mean a trajectory which has no positive Lyapunov exponent.In both cases, regular and chaotic, the LV angle y 12 oscillates fast in time but tends to approach 0 or p before the spikes in z are observed.This means that the peaks in z can be predicted in time when introducing suitable thresholds for y 12 , namely y (min) 12 and y (max) 12 .As exemplified in Fig. 4, the two horizontal blue lines indicate thresholds for these trajectories, namely y (min) 12 = 0.1 and y (max) 12 = 3.04 for the regular case and y (min) 12 = 0.5 and y (max) 12 = 2.64 for the chaotic case.In order to choose adequate thresholds for these trajectories, a statistical analysis of the spikes must be performed, as done in the next paragraph.When y 12 crosses these thresholds, spikes in z are to be expected.We note that y 12 tends to approach 0 or p not only before the large spikes, but also before the small spikes.The crucial distinction is that, preceding large spikes, the angle comes closer to 0 or p than before the small spikes.From a physical point of view, our results show that each time the unstable (chaotic case) and the stable (regular case) manifolds align along the direction of motion, a spike in the z concentration is expected.
Fig. 4 also shows that, far from spikes, the angle y 12 oscillates in time.Such oscillations occur when the eigenvalues of the Jacobian matrix are complex and no invariant manifolds exist.In fact, manifolds rotate and the angles oscillate in time.This oscillatory behavior was observed and explained for specific parameter values of the He ´non map. 5 Analogous angle oscillations were reported recently for interesting critical transition problems, 17 although no explanation concerning the origin of the oscillations was given there.
The results above suggest that for times called prediction times, as exemplified in Fig. 4 by black arrows, it is possible to anticipate the imminent occurrence of a spike.It also becomes clear that prediction time intervals, defined as the absolute values between the prediction times and the times for which the extreme spike occurs, depend on the thresholds chosen.
Similar to before, Table 2 summarizes some statistical results.For this we consider the chaotic trajectory which has 549 spikes during the integration interval.We say that there is a spike when z 4 0.9 and an extreme spike when z 4 3.0.Let us start to discuss the first line, which has a very strong alignment condition.From the total of 549 spikes, N pp = 526 of them were predicted.However, N fa = 33 from the predicted spikes are false alarms, meaning that the 33 spikes are not extreme (z o 3.0).From the 549 À 526 = 23 non-predicted spikes, N ext = 16 are extreme spikes.The last row displays the mean prediction time intervals Dt p .
The general behavior presented in Table 2 is similar to the one in Table 1.For strong alignment conditions the number of non-predicted spikes increases, especially the number N ext of extreme spikes.If the alignment conditions are relaxed, there is an increase in the number of predicted peaks as well as in the number of false alarms.However, the percentage of false alarms is much smaller here.Comparing the stronger alignment conditions with the more relaxed condition, one sees that N pp increases by 4.2% while N fa increases by only 1.0%.More relaxed alignment conditions require longer prediction time intervals, excluding almost all shorter (o1.0) prediction time intervals.
As discussed in the previous section, the appropriate alignment condition depends on what is sought in practical situations.If we again want to find a compromise between allowing some large peaks while at the same time not allowing the number of false alarms to dominate, consulting Table 2 the suitable alignment condition providing an adequate balance is y (min) 12 = 0.03 and y (max) 12 = 3.11.

The autocatalator model
Next, we consider the use of LVs to predict outliers in the complex oscillations generated by the isothermal chemical Table 2 Statistics for the chaotic trajectory with 549 peaks (z 4 0.9) for distinct values of thresholds y (min) 12 and y (max)

12
. The number of extreme spikes (z 4 3.0) for this trajectory is 509.The quantities are: N pp , the number of predicted spikes; N fa , the number of false alarms, i.e. predicted spikes with smaller amplitudes (z o 3.0); N ext , the number of non-predicted extreme spikes (z 4 3.0) and Dt p , the average prediction time interval where a, b, g are dimensionless concentrations and m, k, s, d are control parameters of the reaction.Fig. 5 displays the attractor for the typical parameters k = 2.5, d = 1, m = 0.2983 and s = 0.013 using the initial condition (a, b, g) = (0.5, 2.5599, 0.5).As before, colors represent the values of the angle y 12 for each point on the attractor.The attractor is formed by a chaotic trajectory with Lyapunov exponents l 1 B 0.755, l 2 B 0.002 and l 3 B À11.117.The attractor in Fig. 5 and in subsequent simulations was obtained by solving eqn ( 7)-( 9) numerically with the standard fourthorder Runge-Kutta algorithm and a fixed time-step h = 0.05.The first 10 7 steps were discarded as transient, with LVs computed during the next 6 Â 10 6 steps.For better visualization, only 12% of these points were plotted in Fig. 5.The angle y 12 changes along the attractor.The rightmost green ridge clearly visible in Fig. 5 is related to angles close to p/2 and corresponds to strong hyperbolic motion.In sharp contrast to the reaction considered in Sections 2 and 3, here the variable of interest, namely the concentration b, never relaxes to random quiescent intervals.It always increases towards a peak (large or not) and decreases very fast to a minimum.After each minimum, b promptly starts to increase towards the next maximum, while y 12 and y 13 tend to approach Fig. 5 Attractor for the autocatalator, eqn ( 7)-( 9), with colors representing the angle y 23 .0 or p, usually right after the minimum.This is the alignment of both stable and unstable manifolds along the direction of motion and precedes the appearance of a peak.This is in agreement with our recent proposal 4 that the alignment of LVs along the flow direction precedes peaks in a chaotic dynamical system.Therefore all peaks in the concentration b can be easily predicted by the angles y 12 and y 13 .However, since the concentration b always goes, after the minima, towards the peaks, the alignment of LVs to predict peaks is not necessary in this case.In fact, the dynamics is apparently strongly nonhyperbolic for all the duration of the integration.For all three models considered here, the alignment of LVs is able to predict all peaks during the time evolution, large or not.Therefore, while the comparatively large number of false alarms found in Section 2 (B46.7%) is compatible with this detection power, false alarms found in Section 3 (B7.3%)are acceptable.However, this seemingly negative aspect of the method can be essential to predict spikes in systems where the quiescent motion does not give any indication that a spike will appear.For example, in the quiescent intervals shown in Fig. 2 and 4, the plotted quantities x and z, respectively, present some smaller spikes with increasing amplitudes which indicate that an extreme spike might appear.All these smaller spikes are detected by the alignment of LVs and lead to false alarms.Sometimes, the time evolution of the observed quantities x, z and b in the examples above provides clues that a large spike will appear.But the difficult situation is when such clues are not found in the time evolution.In this circumstance, LV alignment is quite helpful to predict large spikes.Comparing results for the two systems in Sections 2 and 3, one sees that when the alignment conditions are relaxed, the number of predicted peaks increases much more than the number of false alarms.In addition, the number of unpredicted extreme spikes tends to zero very fast.All in all, for the reactions considered, the LV alignment is able to foretell the imminence of spikes.

Conclusions and outlook
This work assessed the alignment of LVs as a quantitative tool to predict large-amplitude spikes in three chemical models governed by distinct polynomial dynamics.For the models in Section 2 and 3, the time evolution of the concentrations presents random intervals of quiescent motion.In these cases, large events are unpredictable by analyzing just the time evolution of the concentrations.Here, LVs are very helpful to predict spikes.In the model presented in Section 4 the concentrations oscillate continuously between minima and maxima.There is no relaxation and the arrival of spikes can be recognized by the time evolution of the concentrations.The three examples considered show that the alignment of Lyapunov vectors along the flow direction provides a straightforward and reliable quantity to predict small and large events for typical chemical reactions.
The LV alignment as a spike prediction method is expected to be very general since it is able to detect all spikes.It is particularly useful when random intervals of quiescent motion are present in the time evolution of the concentrations.Due to the lack of analytical results, the relationship between the LV alignment and extreme events in dynamical systems must be tested for each particular case.An open challenge is to devise a way to apply LVs (i) to observational data, namely, for data which lack an underlying mathematical framework, and (ii) to realistic systems where the presence of noise is inevitable.

Conflicts of interest
There are no conflicts to declare.

Appendix -numerical procedure
The most unstable direction of a chaotic trajectory can be determined from the direction of the eigenvector related to the largest LE, obtained as an eigenvalue of the Jacobian matrix of the system.The procedure for obtaining it is well known and widely used in the nonlinear and complex systems community and an explicit numerical code for the determination of Lyapunov exponents is given by Wolf et al. 18 It involves a Gramm-Schmidt (GS) orthonormalization procedure which, unfortunately, destroys the information about the direction of all eigenvectors not related to the largest Lyapunov exponent.To compensate for this loss, a method proposed by Ginelli et al., 14 based on forward and backward time integrations, allows one to recover the direction of all unstable and stable invariant manifolds.This procedure is somewhat time consuming but very precise.In order to obtain stable and unstable directions in phase space we need to perform (i) forward and (ii) backward time evolutions: (i) The procedure summarized here is well known and uses ordinary differential equation (ODE) integrators and the aforementioned numerical code, 18 based on the GS orthonormalization procedure.For the forward time evolution we integrate the equations of motion for variables ( : x, : y, z ˙) (eqn (1)-(3), for example), together with their corresponding linearized equations with the Jacobian matrix Eqn (10) governs the time evolution in the tangent space.Choosing Dx(0) = a (1) 0 , Dy(0) = a (2) 0 and Dz(0) = a (3) 0 , as an initially orthonormal basis G 0 = (a (1) 0 ,a (2) 0 ,a (3) 0 ), the evolution in tangent space up to a time t = nh (n = 0, 1, 2, 3,. ..) is then obtained from where J n = J nÀ1 ,J n . ..J 0 is the product of the Jacobian matrix of the map evaluated for every orbital point at times t = nh.The subscript represents the number of iterations with time step h.It is implicit here that at the initial iteration (n = 0) we have already computed the orbit for a sufficiently long time (forward transient time) so that the orthonormal basis of the tangent space converged already to the asymptotic Gram-Schmidt vectors. 11In the numerical code, 18 there are three sets of linearized equations, one for G 0 = (1, 0, 0), one for G 0 = (0, 1, 0) and the last one for G 0 = (0, 0, 1).To avoid divergences, the matrix G ˜n is renormalized for each step.As usual, this can be done by the QR decomposition The matrix R n is upper-triangular and contains the information obtained in the GS orthonormalization procedure of G ˜n. Here, G n = (a (1) n ,a (2) n ,a (3) n ) are the GS vectors before the orthonormalization at iteration n and (g (1) n ,g (2) n ,g (3) n ) are the GS vectors after the orthonormalization.Since g (1) n , g (2) n and g (3) n are orthonormal by construction, they can only provide information about the local rates of expansion (contraction) of the vectors, stored in the diagonal elements of R n , from which the standard Lyapunov vectors can be obtained.The correct direction of the eigenvectors of the Jacobian matrix were corrupted by the orthonormalization procedure.They can be recovered by the backward integration, as shown below.
(ii) The essential property is that for the backward time evolution the unstable direction is the stable direction from the forward time evolution.Thus, in the backward evolution the Lyapunov vector related to the largest Lyapunov exponent points along the stable direction from the forward evolution.In this way we define the Lyapunov vectors as v (1) n = c (1,1) n g (1) n , ( 15) n g (1) n + c (2,2) n g (2) n , ( 16) 3) n g (1) n + c (2,3) n g (2) n + c (3,3) n g (3) n . ( Since g (1) n , g (2) n and g (3) n were computed in the forward time evolution, the coefficients c ( j,i) n can be determined from the dynamics in tangent space for the backward direction.In the matrix form we have C jÀ1 = R j À1 C j , where R j À1 is the inverse of the R j matrix.The time j starts to count after the backward transient, which is the time in the backward evolution sufficiently long to converge the tangent initial conditions close to the LVs.The LVs have normalized length so that the columns of C j must be normalized to 1.The initial condition for C n , before starting the backward evolution, can be a generic nonsingular upper triangular matrix, which is the GS basis.Finally, the angles between LVs are calculated from y ij = cos À1 (v (i) n Áv ( j) n ).
Since in this work we use a fourth-order Runge-Kutta integrator, the error in the variables is around h 5 , which is negligible for the time steps h that we use.As convenient rules of thumb, one has to chose time steps small enough compared to the relevant oscillations of the system.For example, the peaks in Fig. 4 occur for small intervals of time, around Dt B 0.4.Thus, to compute many points along such a short peak, we used h = 5 Â 10 À5 so to have around 8000 for each peak (not all plotted here).Once the above compromises are satisfied, the predictions are independent of the time step.On the other hand, the number of steps is more related to how much computational time is available or, in our case, how many peaks are needed to have a reliable statistical analysis.The precision of Lyapunov vectors is more dependent on the transient times than on other parameters.
It is worth mentioning that in general there is no need to realize the backward time evolution.It is enough to determine the unstable direction from the forward evolution and compare it to the flow evolution direction.Below, we present in schematic form the main steps used to determine Lyapunov vectors.

Fig. 1
Fig. 1 Representative attractor for the model of dissolution of copper, eqn (1)-(3), with colors representing the values of the angle y 12 between the unstable invariant manifold and the flow direction.

Table 1 12 .
Statistics for the chaotic trajectory with 2192 peaks (x o À9.0) for distinct values of thresholds y (min) 12and y (max) The quantities are: N pp , the number of predicted spikes; N fa , the number of false alarms, i.e. predicted spikes with smaller amplitudes (x 4 À15.0);N ext , the number of nonpredicted extreme spikes (x o À15.0) and Dt p , the average prediction time interval y(min)

12 .
The two horizontal blue lines in Fig.2are illustrative thresholds for this trajectory, namely y (min) 12 = 0.14 and y (max) 12 = 3.0.Instants of time where y 12 crosses these thresholds (i.e., y 12 o 0.14 or y 12 4 3.0) are called prediction times and a large spike in x is expected to occur.Three such examples are indicated with black arrows in Fig.2, the arrows chosen subjectively, for illustration purposes.For the first two arrows on the left we have y 12 4 3.0 and for the arrow on the right y 12 o 0.14, and they precede the large spikes with x o À15.The thresholds define the alignment condition for the LV along the flow direction and allows us to predict the spikes.The effectiveness to predict spikes may change for different alignment conditions, as discussed below in more detail.

( 6 )
As before, of interest is the temporal evolution of the concentrations (x, y, z), which depend on the control parameters a, b, e, f, g, s, c, a, b.As usual, we kept fixed the constants a = 0.5, b = 3, e = 0.01, f = 0.5, g = 0.6, s = 0.3 and c = 4.8 and change the pair (a, b).Fig. 3 displays the chaotic attractor for a = 0.7825 and b = 0.55, starting from the initial condition (x, y, z) = (1, 1, 1).Colors represent the values of the angle y 12 for each point on the attractor.The attractor consists of a chaotic trajectory with Lyapunov exponents l 1 B 0.13, l 2 B 0.0005 and l 3 B À544.

Fig. 2 A
Fig. 2 A time window displaying the evolution of x (black curve), of y 12 (red), and two straight lines (blue) indicating the values of y (min) 12 = 0.14 and y (max) 12 = 3.0 for the alignment condition.

Fig. 4
Fig. 4 Time window showing the variable z (black lines) and y 12 (red lines) for part of (a) a regular trajectory with parameters a = 0.7825 and b = 0.39213 and (b) the chaotic trajectory from Fig. 1.Two horizontal blue lines indicate thresholds to predict spikes and the two black arrows show two prediction times where y 12 crosses the thresholds.

Fig. 6
presents a time-window of the evolution of b (black lines) and associated LVs for the chaotic trajectory from Fig. 5. Fig. 6(a) presents the angle y 12 between the unstable manifold and flow direction, while Fig. 6(b) displays the angle y 23 between the stable manifold and flow direction.Besides the very fast variations of these angles, they assume values close to 0 and p for almost all integrated times.Next, we consider a regular trajectory obtained for k = 2.5, d = 1, m = 0.29776 and s = 0.013 and initial conditions (a, b, g) = (0.5, 0.25599, 0.5).The Lyapunov exponents for this case are l 1 B 0.003, l 2 B À0.914 and l 3 B À16.16.Therefore the angles y 12 and y 13 are the angles between stable manifolds along the flow direction.The time evolution of the concentration b is shown in Fig. 7 (black lines) with colors representing the angle y 12 obtained along the trajectory.Results are analogous to the chaotic case discussed above, also for y 23 .Since for the autocatalator the variables y 12 and y 23 show essentially the same behavior, we use y 23 to show that the prediction method also works when using angles between other manifolds.

Fig. 6
Fig. 6 Time evolution of the concentration b (Â3) (black lines) for a chaotic trajectory and (a) the angle y 12 (red line) and (b) y 23 (green line).

Fig. 7
Fig. 7 Time evolution of the concentration b (black lines) for a regular trajectory and (a) the angle y 12 (red line) and (b) y 13 (green line).