Constructing ordinal partition transition networks from multivariate time series

A growing number of algorithms have been proposed to map a scalar time series into ordinal partition transition networks. However, most observable phenomena in the empirical sciences are of a multivariate nature. We construct ordinal partition transition networks for multivariate time series. This approach yields weighted directed networks representing the pattern transition properties of time series in velocity space, which hence provides dynamic insights of the underling system. Furthermore, we propose a measure of entropy to characterize ordinal partition transition dynamics, which is sensitive to capturing the possible local geometric changes of phase space trajectories. We demonstrate the applicability of pattern transition networks to capture phase coherence to non-coherence transitions, and to characterize paths to phase synchronizations. Therefore, we conclude that the ordinal partition transition network approach provides complementary insight to the traditional symbolic analysis of nonlinear multivariate time series.

The basic idea of ordinal partition network method can be traced back to identifying ordinal patterns of time series 31,32 . Considering a one-dimensional time series = … x t { ( )} t L 1, , comprising of L points from a dynamical system, the original phase space can be reconstructed by time delay embedding )] x with dimension D x 33,34 . The next step is to compute the rank order of τ x based on relative amplitudes, which is conveniently represented by a symbol π x (t). When sliding windows from t = 1 to N = L − (D x − 1)τ in the embedded space, a symbolic representation of the trajectory π x (t) is produced. One traditional approach, following the symbolic representations, is to compute permutation entropy based on the frequency plot of order pattern which yields very well established statistical measures in nonlinear time series analysis 31 . In the recent decades, ordinal symbolic representation of time series has found a number of interesting applications in science and engineering, for instance, biomedical recordings 35 , finance 36 , climate sciences 37 . Some recent progress has been comprehensively reviewed in ref. 35. However, the transition behavior between ordinal patterns remains largely unclear. The recent ordinal partition network representations capture the evolutionary behavior of the ordinal patterns 25,26 , which sheds novel insight on the standard ordinal symbolic analysis of time series.
For a given embedding dimension D x , there are a total of D x ! unique ordinal patterns that can possibly occur in a time series, neglecting equality. A stochastic process of stationary increments fulfils P(x t = x t+τ ) = 0 and therefore the probability to have ties x t = x t+τ is zero. For empirical time series, we can avoid ties by adding a tiny white noise with continuous distribution 38 . Therefore, the original phase space is decomposed into D x ! equivalent partitions 31 . It is intuitive that all D x ! patterns almost occur with equal frequencies in a time series generated by a stochastic process for N → ∞. However, a set of patterns may never occur in a time series produced by deterministic dynamics. Therefore, it is possible to quantify determinism in time series data by counting the forbidden patterns. However, complications arise in real time analysis. For instance, missing ordinal patterns might be related to finite time length during the period of observation and correlated stochastic processes, which require some revised methods for the detection of determinism in relatively short noisy data [39][40][41][42][43] . From the viewpoint of ordinal partition networks, both the frequencies of order patterns and the transitions between different patterns are inhomogeneous. Therefore, network properties thus obtained are sensitive to different system dynamics, which successfully characterize the difference between healthy and patients from EEG data 25,26 .
Most of the recent works have focused only on univariate time series {x(t)}. The embedding dimension D x and time delay τ are two important parameters for constructing ordinal partition networks, in particular having crucial impacts on the appearance of forbidden order patterns [27][28][29] . However, the generalization to multivariate time series remains largely untouched 5,44 . Most of the observable phenomena in the empirical sciences are of a multivariate nature. For instance, assets in stock markets are observed simultaneously and their joint development is analyzed to better understand tendencies. In climate science, multiple observations (temperature, pressure, precipitation, human activities etc, from different locations) are the basis of reliable predictions for the future climate conditions. We propose to construct ordinal partition transition networks from multivariate data.

Results
Ordinal pattern definitions. Given a scalar time series {x(t)} which is produced by a deterministic dynamical system, the order structure of the time series depends on the embedding dimension D x and time delay τ 33,34 . Let us start with D x = 2. Neglecting equality, we have two relations between x(t) and x(t + τ), namely, two symbol sequences representing order patterns π x : x For dynamical systems with continuous distributions of the values, we can neglect equality because the Lebesgue measure of ties is zero 31 . In addition, a large amount of numerical simulations suggest that the results do not change qualitatively and it does not matter whether we count 45 . In the practical application, we may easily test for < and ≤. Therefore, we follow the routine as suggested in ref. 45 and do not separately consider equalities. Time delay τ is chosen as 1 in this work. By this choice, the order pattern π x 1 captures the increasing trend, respectively, π x 0 corresponds to the decreasing trend of the time series. This definition is equivalent to considering the signs of the increments Δx(t) = x(t + 1) − x(t) by a first-order difference of the original series. In other words, the associated order patterns capture the variations of x(t) in its velocity space, showing dynamic rather than static information based on the displacement directly.
When generalizing the above idea to two dimensional time series (x(t), y(t)), we restrict our discussion on embedding dimension D x = D y = 2 for individual variable. Therefore, we have four different combinations of order patterns depending on the signs of increments (Δx(t), Δy(t)) ( Table 1): In the phase space of (x(t), y(t)), we have order pattern Π(t) ∈ (π 1 , π 2 , π 3 , π 4 ) which captures the increasing or decreasing behavior. The example of Fig. 1(a,c) shows the construction of ordinal patterns for two dimensional series (x(t), y(t)). In a full analogy, based on increments (Δx(t), Table 1. Order patterns in two dimensional time series (x(t), y(t)).
sional time series (x(t), y(t), z(t)) is enumerated in Table 2 and visualized in Fig. 1(b,d). Therefore, the dimension of order pattern Π(t) for an n-dimensional time series is D = 2 n since each component has either increasing or decreasing trend at time t.
Note that, in this work, we do not apply time delay embedding technique to obtain the multi-dimensional phase space from one univariate time series. Instead, given multi-variate time series, we consider the increments between two consecutive time points of each measurement in the space of multi measurements. In other words, our approach captures the dynamic properties of the multi-variate time series in its associated velocity space (difference space). Therefore, time delay τ in the order pattern definition (Eq. (1)) has rather a different interpretation with the time delay that is often used in embedding. Traditionally, one chooses appropriately an embedding dimension and time delay to reconstruct phase space from a given univariate time series. We can certainly generalize the discussion to the case of time delays larger than 1 (i.e., τ > 1) and embedding dimension D x > 2 for each variable (measurement), but we think that the physical meaning in terms of dynamics becomes ambiguous for multivariate time series.
Ordinal partition transition networks. Given a multi-variate time series, for instance, the two dimensional case of (x(t), y(t)), we denote the frequency of the i-th pattern π i as p(π i ) which is computed over the time interval = … t N 1, , . One important property is that ordinal patterns from a deterministic process have different frequencies p(π i ). Permutation entropy  O is then introduced to characterize the inhomogeneous appearance of ordinal patterns as following where the sum runs over all D = 2 n permutations. We use log 2 and hence the units of O  are bits. For a n-dimensional independent identical distributed stochastic process, one obtains the largest entropy  = n O since each of D = 2 n ordinal patterns is expected to have the same frequency.
We illustrate the above algorithm by using a toy model of three dimensional identical independent periodic time series in Fig. 2. Different combinations of periods of the 3D periodic series often yield different values of permutation entropy  O . In addition, the limited number of possible ordinal patterns (forbidden patterns) are widely observed reflecting the determinism of the series.
Most of the current studies focused on the computation of permutation entropy  O considering the frequencies of order patterns, which do not disclose the transition behavior between order patterns. Therefore,  O is static. For instance, the details of the dynamics remain unclear in Fig. 2 Fig. 3(a,b).

Figure 1.
Order pattern definitions for (a) two dimensional series (x(t), y(t)) and (c) its increment series (Δx(t), Δy(t)). (b) Three dimensional series (x(t), y(t), z(t)) and (d) the corresponding increment series (Δx(t), Δy(t), Δy(t)). Signs of the increment series and the ordinal patterns are respectively indicated in (c,d). Comparing to a 3-dimensional uncorrelated independent identical distributed random uniform noise, the ordinal partition network is a complete connected graph ( Fig. 3(c)). In addition, we indicate each directed link by its transition frequency w ij = p(π i → π j ), following the time iterations of the series. Finally, we come up with a weighted directed network characterized by a weighted adjacency matrix Here, based on W, the regularity of the order pattern transition properties is quantified by the Shannon entropy  T , which is where the sum runs over all possible 2 2n transitions. In a full analogy to O  , for a n-dimensional independent identical distributed stochastic process, one obtains the largest entropy = n 2 T  . It is rather straightforward to compute T  for time series produced by stochastic processes. However, we need to pay attention to the case of continuous systems, which yield a larger proportion of self-loops in the resulted networks as will be illustrated below. Here we take the chaotic Rössler system as an example which reads ω ω where a = 0.165, and ω = 1.0. The Eq. (4) are numerically integrated by the fourth-order Runge Kutta method with integration step h = 0.01. The first 10000 transient data points are discarded and time series consisting of N = 500000 data points are analyzed. Short segments of time series (x, y, z) are shown in Fig. 4(a). Due to the continuity property of the system, there are many plateaus reflecting the invariance of the order patterns during certain time intervals (Fig. 4(b)). These plateaus are reflected by self-loops in the resulting transition networks ( Fig. 4(c)). In most of the existing studies of complex networks, self-loops are avoided because of both  , where the transition route π 1 → π 5 → π 6 → π 8 → π 4 → π 3 → π 1 is observed. The values on links represent the corresponding transition frequencies of the ordinal patterns. Note that N = 500000 data points are used in obtaining (c,d).
the computational simplicity and theoretical concerns 46 . In the case of the Rössler system, there is about 99% self-loops and only about 1% non-self-loops (as indicated by the arrows in Fig. 4(c)).
The weighted adjacency matrix W can be split into diagonal and off-diagonal entries. Taking into account the numerical observations that the diagonal elements of W (self-loops) are much larger than the off-diagonal ones, Therefore, we obtain   ≈ T O for a continuous system when self-loops are considered. In this case, the transitions between different ordinal patterns are hard to be captured by  T .
In order to emphasize the importance of non-self transitions between ordinal patterns, we remove the self-loops as shown in Fig. 4(d) by setting the diagonal values to be 0 in W. This is typical of most research work on complex networks 46 . Furthermore, we remove self-loops before computing the weighted matrix W to keep the Note that self-loops should not be expected with large amounts in stochastic processes.
Ordinal pattern partitions of phase space. Ordinal pattern transition networks provide us with an alternative for phase space partitions, which utilizes nullclines of the systems. Here we show two examples covering discrete and continuous dynamical systems.
Example (1): the Hénon map x t y t x t y t x t 2 is chosen as an example for a chaotic two-dimensional map. The order pattern partitions of the attractor are shown Fig. 5(a), which is color coded by different order patterns. A segment of time series is shown in Fig. 5(b). The histogram of order patterns (Fig. 5(c)) discovers that π 4 are forbidden patterns of the system, which yields 1 76 T . According to the order pattern definitions in Table 1, the phase space partitions of the Hénon map are delineated by nullclines as follows: These two lines are shown in Fig. 5(a), where we find no points of the attractor lying in the region of order pattern π 4 no matter with the iteration steps. The disappearance of π 4 pattern suggests that there is no intersection between the π 4 partition and the attractor, except for the unstable fixed point of (0.63, 0.19) (intersection of L 1 and L 2 ). Example (2): the Rössler system, Eq. 4 is chosen as a continuous dynamical system. When the parameter a = 0.165, the attractor is shown in Fig. 6(a) with phase space points being further color-coded by ordinal patterns. The boundaries of each partition are determined by the corresponding nullclines, i.e., dx/dt = 0, dy/dt = 0 and dz/dt = 0. The ordinal pattern transition network is shown in Fig. 4(d). In this case, neither π 2 nor π 7 appears, which is explained as follows 47 . When transforming phase space into ordinal partition transition network, we associate to each state (x, y, z) with an order pattern such as (+, −, +) (as listed in Table 2). This 3-dimensional ordinal pattern describes which variables of (x, y, z) are increasing and which are decreasing at a given time.
Taking variable x as an example, both y and z are repressors to x because of the negative signs in the Jacobian (−ω,  Fig. 4(d). Panel (c) is the same as (a) with a = 0.26, where a significant number of π 2 patterns are highlighted, and (d) is the ordinal partition transition network (self-loops are excluded), where an alternative transition from π 4 → π 2 → π 1 has been observed. and −1). However, x is an activator for variable y because the element of the Jacobian matrix is ω being positive everywhere in phase space. The scheme of activation and repression of the Rössler system is shown in Fig. 6(b).
The condition that the trajectory can cross a nullcline meaning a change from increasing to decreasing (or vice versa) are equivalent to determining local maximum or minimum. Because of the continuity of the system, we have the following rules to have a maximum or minimum 47 : (i) a variable cannot have a maximum if all its repressors are decreasing and all its activators are increasing; (ii) a variable cannot have a minimum if all its repressors are increasing and all its activators are decreasing. These two rules yield all possible transitions between different order patterns as shown in Fig. 6(b). However, the transition network observed for a typical trajectory in phase space is determined by the given set of parameters a and ω. For the case a = 0.165 and ω = 1, we only find the transition route π 1 → π 5 → π 6 → π 8 → π 4 → π 3 → π 1 (shown in Fig. 4(d)), and meanwhile π 2 and π 7 are forbidden patterns. Increasing the value a to 0.26, the Rössler system presents screw-type chaotic oscillations with irregular kicks, which yield an alternative transition from π 4 → π 2 → π 1 as highlighted in Fig. 6(c). Therefore, we observe two transition routes of the patterns from π 4 to π 1 (Fig. 6(d)), while π 7 remains to be absent. In other words, the appearance of π 2 pattern suggests that the changes of the ordinal patterns are sensitive to the geometric changes of the attractor.

Identifying dynamical transitions.
We apply ordinal partition transition networks to identify dynamical transitions in two different cases: (i) phase coherence to non-coherence transition which is a weak chaos-chaos transition, (ii) paths to phase synchronization transitions. In both examples, we show frequency plots of ordinal patterns (without self-loops), complexity entropy measures of O  (with self-loops) and T  (without self-loops). In addition, we compare the case (i) to coherence index (CI), and case (ii) to mean rotation frequency of each oscillator Ω i .
Example (1) shows phase coherence to non-coherence transitions in the chaotic Rössler system (Eq. (4)), where the parameter a is systematically varied in the range [0.15, 0.25]. As it has been systematically shown in ref. 48, this parameter range comprises different kinds of dynamics, including periodic windows, phase coherent chaos (existence of a well-defined rotation center in phase space) as well as non-phase coherent chaotic oscillations (lack of a distinct center of rotation). The transition between phase coherence and non-phase coherence chaos occurs at a c ≈ 0.206. More specifically, for a < a c , the chaotic attractors are always phase coherent, whereas they are non-phase coherent for a > a c . We refer readers to ref. 48 for further discussion on various measures to detect this chaos-chaos transition as well as periodic windows, ranging from traditional measures of phase coherence factor, phase diffusion coefficient, recurrence quantification based discriminators, and recurrence network based measures. In this work, in order to avoid repetitions we only discuss the capabilities of ordinal pattern changes and entropies  O and  T in detecting the transition from phase coherent to non-coherent chaos by comparing to the measure of phase coherence index (see Methods). Figure 7 shows the bifurcation diagram when the parameter a is changed. First, the frequency of ordinal pattern π 2 is zero (f(π 2 ) = 0) when a < a c , and becomes positive when a > a c . In contrast, f(π 3 ) decreases when a > a c (Fig. 7(a)). Much smaller changes are observed for the other ordinal patterns π 1 , π 4 , π 5 , π 6 and π 8 . Pattern π 7 does not appear in the entire interval of a. Ordinal patterns π 2 and π 3 are sensitive to the geometric changes of the attractor, capturing the transition from phase coherent to non-coherent regime.
In addition,  O shows rather small changes while  T is a constant value when a < a c (Fig. 7(b)). The behavior of  O and T  has been confirmed by the coherence index (see Methods) as shown in Fig. 7(c). Meanwhile, we find some discrepancies for these measures when the control parameter a increases towards the transition point between phase coherence and non-coherence regimes. In particular, both O  and the coherence index increase slightly before the transition point a c , while T  increases sharp at a c , as highlighted by vertical dashed lines in Fig. 7. This is because of the homo-clinic point at the origin. As the control parameter a increases within the phase coherent regime, the attractor successively grows and finally extends to the vicinity of the origin shortly before the transition to the funnel regime, where a unique rotation center of trajectories in phase space is lost. The dynamics in the (x, y)-plane becomes very slow whenever a trajectory gets close to the homo-clinic point. As a consequence, there is a high density of sampled points on a trajectory in the neighborhood of the origin. By the same time, these re-injection to and ejection from the origin are rather irregular events, which introduce fluctuations in the computations of  O and the coherence index. In contrast, when computing  T , the local velocities change direction from increasing to decreasing (or vise versa) only after the transition to the non-phase coherent regime. Therefore, T  shows good sensitivity on the change of the local velocity space when the control parameter a passes the transition from phase coherent to non-coherent regime, which are shown in Fig. 7.
Note that all measures of  O , T  and coherence index show pronounced local maxima in periodic windows (for instance, at a = 0.227 and a = 0.245) 48 .
Example ( where k = 1, 2, 3 and κ is the coupling strength. We consider non-identical oscillators by choosing ω 1 = 0.98, ω 2 = 1.02, ω 3 = 1.06. The parameter a is chosen as 0.165 such that the subsystems are in the phase coherent regime ( Fig. 6(a)). The oscillator k = 2 is bidirectionally coupled to both k = 1 and k = 3, whereas there is no direct coupling between k = 1 and k = 3. The Eq. (8) are numerically integrated by the fourth-order Runge Kutta method with integration step h = 0.01. The first 10000 transient data points are discarded and time series consisting of 150000 data points are analyzed. We construct ordinal pattern transition networks from the x k components, namely, (x 1 , x 2 , x 3 ), following the same pattern definitions as shown in Table 2.
Our motivation of using Eq. (8) is to study the variations of ordinal patterns on the paths to phase synchronization, focusing on the evolutionary process of the transition networks with different regimes of synchronization, which is more complex than the case of a single Rössler system. The results are shown in Fig. 8. Furthermore, the results of Fig. 8 have been averaged over 50 random initial conditions when integrating Eq. (8).
In the regime of no synchrony (κ < κ c1 = 0.036), three oscillators evolve almost independently such that all ordinal patterns have the same frequencies of 0.125. There are rather small gradual changes only when κ approaches to κ c1 (Fig. 8(a)). The entropy value  T is more sensitive to these gradual changes showing a pronounced decreasing trend, while O  seems to be a constant ( Fig. 8(b)). The average rotation frequencies Ω k (see Methods) of each oscillator are shown in (Fig. 8(c)), which confirms no synchrony in this coupling regime.
In the regime that phase synchronization appears between oscillators k = 1 and k = 2, but not with k = 3 (κ κ κ ∈ = . . ), we observe monotonic increasing trends for order patterns π 1 , π 2 , π 7 , and π 8 ( Fig. 8(a)). In addition, we find relatively slower increasing trends for patterns of π 4 and π 5 . In contrast, some monotonic decreasing trends are found for π 3 and π 6 . The changes in the frequencies of order patterns are captured by both entropy values O  and T  , showing gradual decreasing trends ( Fig. 8(b)). The average rotation frequencies Ω k are shown in Fig. 8(c), where k = 1 and k = 2 are phase locked to the same rotation frequency but not with k = 3.
From the viewpoint of high dimensional systems of coupled oscillators, in the process from non-synchrony to phase synchronization we find that the transition networks have experienced rather random transitions between all possible ordinal patterns to a state of transitions between a limited number of ordinal patterns as shown in Fig. 9. In addition, we find π 3 and π 6 are forbidden patterns if all three oscillators are synchronized.

Conclusions
In this work, we have proposed to construct ordinal partition transition networks from multivariate time series, which help us to analyzing the interaction patterns between different components. The basic idea is to capture the directions of changes in the associated velocity space, which yields dynamic instead of static information in the original phase space. The resulting ordinal partition transition networks are weighted directed networks, which are fundamentally different to recurrence networks 19 and visibility graphs 10 . For time series from both discrete and continuous dynamical systems, we find that the frequencies of observed ordinal patterns are inhomogeneous which is quantified by the entropy  O . In addition, the transition frequencies between different ordinal patterns are inhomogeneous as well, which is characterized by the entropy T  . Note that no essential difference between  O and T  is expected for discrete systems, however for continuous systems,  T is a better way to characterize the ordinal partition transition networks because O  is more influenced by self-loops as shown in Fig. 4. The ordinal partition transition network utilizes nullclines to generate partitions, resulting in a Markov chain representation of the time series in phase space. The transition between two ordinal patterns is determined by the changes of signs of the increments of the variables. As we have demonstrated in the chaotic Rössler system, this definition is sensitive to capture the geometric changes in phase space, for instance, from phase coherence to non-coherence transition.
Note that our ordinal partition transition network generation algorithm is different to the recent work on constructing temporal networks to capture the memory effects 50 . It will be a future subject to characterize the memory effects by means of ordinal partition transition networks. In addition, we have focused on embedding dimension D x = 2 and delay τ = 1 for each variable which captures either increasing or decreasing trends of a time series in differenced space. One can certainly generalize the algorithm to higher values of D x and τ, however, it is computationally more demanding. For instance in a n-dimensional multivariate series is used. Additionally, the dimension of the transition matrix W is (3!) n × (3!) n . Meanwhile, the increase of dimension D x requires longer time series in order to estimate the transition frequencies of ordinal patterns more reliably. From the view point of the algorithm, no computational complexity is introduced by using a large time delay τ > 1, which, however, lacks a proper interpretation in terms of velocity of the variable. There is one open problem to estimate the ordinal partition transition matrix reliably from short time series, especially when noise plays a significant role. We have applied ordinal partition transition networks to investigate the paths to phase synchronization, showing that a dominant transition route emerges in the networks when the coupling strength is increased undergoing different synchronization transition regimes. As the degree of synchronization is strengthened, dynamics of the coupled systems are locked to the synchronization manifold which yield a dominant transition route in the resulting ordinal partition networks. Before the appearance of synchronization, another challenging task is to distinguish the indirect from direct coupling directions 49,51 , which is very common in climate data analysis, i.e., extracting network interaction patterns from multi-channel time series from distant places 52 . In the case of three coupled Rössler subsystems as we considered (Eq. (8)), the oscillator k = 2 is bidirectionally coupled to both k = 1 and k = 3, whereas there is no direct coupling between k = 1 and k = 3. It is possible to introduce ordinal recurrence plots 45 to tackle this problem.
Traditionally, the computation of permutation entropy based on ordinal symbolic representation of time series does not include pattern transition behavior following the trajectory in phase space. In contrast, ordinal partition transition network approaches take into account the time evolution information explicitly, which hence provide much complementary insights to the traditional symbolic analysis, showing high potentials for experimental time series analysis.

Methods: Phase coherence index
Here we summarize the major steps in computing phase coherence index as we have done in ref. 48. We restrict our attention in this work to the standard analytical signal approach. Here, a scalar signal x(t) is extended to the complex plane using the Hilbert transform . . denotes Cauchy's principal value of the integral, which yields the phase The above definition is straightforward for phase coherent dynamics. In the regime of non-phase coherent dynamics, an alternative phase definition has been proposed based on the local curvature properties of the analytical signal 53 , i.e., φ = .  t dy t dt dx t dt ( ) arctan ( )/ ( )/ ( 11) Since in the standard Hilbert transform-based definition, the phase variable φ(t) does not necessarily increase monotonously in time, we quantify this monotonicity in order to obtain a simple heuristic order parameter for phase coherence, which we will refer to as the coherence index T 0 with φ φ =  t d t dt ( ) ( )/ . Furthermore, the instantaneous frequency of a chaotic oscillator is then defined as the derivative of the phase variable with respect to time. Averaging this property over time yields the mean frequency π φ Ω = . d t dt 1 2 ( ) (13)