MMultiplex Recurrence Networks
Deniz Eroglu , , ∗ Norbert Marwan , Martina Stebich , and J¨urgen Kurths , , Potsdam Institute for Climate Impact Research (PIK), 14473 Potsdam, Germany Department of Physics, Humboldt University, 12489 Berlin, Germany Senckenberg Research Station of Quaternary Palaeontology Weimar, Am Jakobskirchhof 4, 99423 Weimar, Germany Institute for Complex Systems and Mathematical Biology,University of Aberdeen, Aberdeen AB24 3UE, United Kingdom (Dated: March 9, 2020)(Preprint version of doi:10.1103/PhysRevE.97.012312)We have introduced a novel multiplex recurrence network (MRN) approach by combining recur-rence networks with the multiplex network approach in order to investigate multivariate time series.The potential use of this approach is demonstrated on coupled map lattices and a typical examplefrom palaeobotany research. In both examples, topological changes in the multiplex recurrencenetworks allow for the detection of regime changes in their dynamics. The method goes beyondclassical interpretation of pollen records by considering the vegetation as a whole and using theintrinsic similarity in the dynamics of the different regional vegetation elements. We find that thedifferent vegetation types behave more similar when one environmental factor acts as the dominantdriving force.
I. INTRODUCTION
In order to understand the dynamical behavior of sys-tems in a broad range of scientific fields such as physics,biology, medicine, climatology, economy etc., time seriesanalysis provides crucial techniques. Although investi-gation of time series can be done by various techniques,such as basic statistics, symbolization, power spectra orsimilarity analysis, phase space based methods have be-come an important role in dynamical systems’ analysis.Recurrence of a trajectory in its phase space is one ofthe most important fundamental features of dynamicalsystems. Related approaches have been used for severaldecades. The first recurrence based approach is known asthe Poincar´e recurrence theorem introduced in 1890 [1].The theorem indicates that almost all trajectories of dy-namical systems will turn back infinitesimally close totheir previous positions after a sufficiently long but finitetime [1]. Continuous dynamical systems can be definedby a set of ordinary differential equations. These equa-tion sets are called as flows and if the phase space of aflow has a bounded volume, then the Poincar´e recurrencetheorem is always valid [2]. Among the many differentapproaches of analyzing dynamical systems by their re-currence, the recurrence plot (RP) is a multifaceted andpowerful approach to study different aspects of dynami-cal systems. Introduced by Eckmann et al. in 1987, theRP is a matrix to show the times of recurrences of a tra-jectory in its phase space [3]. Afterwards, many statisti-cal quantification methods based on RP were developedfor characterizing dynamical properties, regime transi-tions, synchronization, etc. [4].Understanding underlying dynamics and detectingpossible regime changes in the evolution of dynamical ∗ [email protected] systems are important problems studied by time seriesanalysis. For instance, we assume a dynamical system (cid:126) ˙ x = f ( (cid:126)x, r ) (1)where (cid:126)x ∈ R m , f : R m → R m and a control parame-ter ( r ) that has not to be constant on time r = g ( t ).The purpose of the analysis is to detect possible dynam-ical regime changes in the time series caused by the timedependence of r . Transitions in the dynamics can be de-tected by different RP based measures, which in generalare powerful to study complex, real-world systems [5–7]. Examples of their successful application to real-worldsystems have been found in medicine [6, 8–12], Earthscience [7, 13–17], astrophysics [18, 19], electrochemistry[20], and others [21–23].In the last decade, transformation of a time series to acomplex network has become a very powerful approachto analyze complex dynamical systems. There are sev-eral ways to convert a time series to a network such assymbolic dynamics based techniques [24, 25], visibilitygraphs [26–28], cycle networks [29], or recurrence net-works (RN) [30–32]. In this work we consider recurrencebased approaches, since it is well known that recurrencesare thumbprint of characteristic properties of dynamicalsystems [22, 33–37]. Moreover, recently we have pre-sented a work on fundamental behavior of recurrenceplot measures which proves that the Huberman-Rudnickunique scaling is valid for the RP measures, while chang-ing the control parameter of a given dynamical system[38]. The adjacency matrix of a complex network rep-resents the structure of the system and thus determinesthe links between the nodes of a network. For unweightedand undirected networks, the adjacency matrix is binaryand symmetric, hence very similar to an RP. As a con-sequence, we know that there is a similar unique scalingbehavior for RN measures and RNs have been used toinvestigate real-world systems such as the climate [7] orthe cardio-respiratory system [39]. a r X i v : . [ phy s i c s . d a t a - a n ] M a r Naturally many real systems possess many degrees offreedom and such systems can be described by multivari-ate time series. Each component of these systems canbe considered as a time series and we can reconstruct aphase space with them. Meanwhile, when the number ofcomponents are huge, we need to have longer time seriesfor enough occurrence of recurrences in the phase spacein order to analyze system by the recurrences. How-ever, in several disciplines like astrophysics, earth sci-ences and economy, having long time series cannot beensured. Therefore the components are analyzed one byone or some dimension reduction is applied, but thesemight result in further information loss from the system.An increasing number of dimensions requires longertime series for traditional recurrence networks. Usu-ally, the period of a trajectory in the phase space is ex-tending with an increasing number of dimensions. Withother words, if the dimension of the phase space is in-creasing, we need longer and longer time series to haveenough number of recurrences to analyze the given sys-tem. Therefore, traditional RNs are not efficiently appli-cable for short multivariate time series to interpret thebehavior of the system, because of the scarcity of recur-rences. In order to gain the maximum yield from datasets which has a many degrees of freedom, in this workwe introduce a multiplex network-based approach to an-alyze multivariate time series. Network modeling is amain tool in characterizing a broad range of problems inphysics, biology, social interactions, Earth science, econ-omy etc. If there are interactions between componentsof the system they can be modeled by a coupled net-work structure. For example, the RN approach has beenused to identify the direction of inter-systems relation-ships between bivariate time series [40]. Now we repre-sent each component of the system as a separated RNand interpret it as a layer of a multiplex network [41–45] that we call multiplex recurrence network (MRN). Asimilar approach based on visibility graphs was recentlydiscussed by Lacasa et al . [28], resulting in multiplex visi-bility graphs [28]. Using the measures given by Lacasa etal. [28], we will show a comprehensive comparison of reg-ular RN and multiplex recurrence networks on coupledchaotic systems. By employing the recurrence proper-ties of the dynamical system for the multiplex networkconstruction, we focus on the dynamics, whereas the visi-bility graph is more restricted to the statistical propertiesof the system (e.g., Hurst exponent) [27, 32]. We demon-strate our new method’s efficiency by investigating highdimensional systems which are not possible to be ana-lyzed by the traditional RN approach directly.As a real-world application, we analyze a palaeoclimaterecord which is a multivariate data set of pollen taxa rep-resenting the variability of past vegetation in NE Chinaover the Holocene and found significant variations in thecongruence of the vegetation dynamics of the consideredtree species.
II. METHODS
Recurrence based techniques have been successfullyused for time series analysis of physical, biological, eco-nomical, climate systems. Multi-layer networks havebeen recently introduced as a powerful representation ofa specific network of networks. In this section, we dis-cuss RNs and how to use them in multi-layer networksin order to analyze multivariate time series.
A. Recurrence Networks
A time series { u i } Ni =1 can be reconstructed as a trajec-tory (Eq. 1) in its phase space with time delay embedding[46] x i = ( u i , u i + τ , ..., u i + τ ( M − ) , (2)where M is the embedding dimension which can be foundby a false nearest neighbours approach and τ is the em-bedding delay which can be computed by mutual infor-mation or auto-correlation [47].Two state vectors of the reconstructed time series areconsidered to be recurrent if the second vector falls intothe neighbourhood (an (cid:15) -radius sphere) of the first vec-tor. For the trajectory (cid:126)x i ( i = 1 , . . . , N, (cid:126)x i ∈ R M ), theadjacency matrix of RN, A , is defined as A i,j ( (cid:15) ) = Θ (cid:0) (cid:15) − (cid:107) (cid:126)x i − (cid:126)x j (cid:107) (cid:1) − δ i,j , i, j = 1 , . . . , N, (3)where N is the trajectory length, Θ( · ) is the Heavisidefunction, (cid:107) · (cid:107) is the Euclidian norm and δ i,j is the Kro-necker delta ( δ i,j = 1 if i = j , otherwise δ i,j = 0) [31, 48].A RN is constructed in the following way: Consider thetime points of a time series as nodes of a network; if thenodes are sufficiently close to each other, in other words,if the space vectors are neighbours, then a link betweenthem is drawn. A represents the network, where A i,j = 1if i and j are connected, otherwise A i,j = 0.In all recurrence based applications, the threshold (cid:15) isan arbitrarily selected small number. The selection of (cid:15) can affect the results easily. In order to have a reason-able analysis, some threshold selection techniques wereproposed [4, 48, 49]. However none of them is a cer-tain way to choose the threshold. To be consistent, inthis work we always used the way which depends on thestandard deviation of the time series [50]. B. Multiplex and Weighted Recurrence Networks
Multiplex structure.
In this work, an m -layer mul-tiplex network is constructed by RNs. For m -dimensionalmultivariate time series, we can create m different RNswhich have the same number of nodes and each node islabeled by its associated time. These networks will formthe different layers of a multilayer network. The lay-ers are connected each other with the same time labeled timexyz L x L y L z t - t - t t + t + t - t + t - t - t t + t + t - t + t - t - t t + t + t - t + L x L y L z I x,y I y,z I x,z multivariate time series recurrence network multiplex network weighted network of recurrence networks FIG. 1. Illustration of the procedure for generating network structures from a multivariate time series:multivariate time series → recurrence networks → multiplex recurrence network → weighted network of recurrence networks.Dashed lines between multiplex network’s layers connect all layers to all. nodes. This procedure requires that the time samplingis the same for all of the used time series. If a multi-layer network consists of m -layers which has the samenumber of nodes and the connections between layers areonly between a node and its counterpart in the other lay-ers, then we call such networks “multiplex”. Networks,transformed from multivariate time series, are compat-ible with the definition of multiplex networks, becauseeach node is uniquely assigned to a certain time point ofthe multivariate time series, i.e., we will find the equallytime-labelled nodes in all layers.We consider an m -dimensional multivariate time series { u ( t ) } Nt =1 , with u ( t ) = ( u ( t ) , u ( t ) , . . . , u m ( t )) ∈ R m forany value of t . Then, the RN of the κ th component of u ( t ) is created and located into the associated layer κ of the multiplex network. For an example with N = 3time series these steps are illustrated in Fig. 1. First, wehave three time series and construct a phase space foreach signal. The RN of the κ th component of the timeseries is calculated by Eqs. (2)–(3) and then placed in tolayer κ . We denote the adjacency matrix of the κ th layeras A [ κ ] = a [ κ ] ij and a κij = 1 if nodes i and j are connectedin layer κ , a κij = 0 otherwise. The giant adjacency matrixdescribing the entire multiplex network is denoted by A = A [1] I N . . . I N I N A [2] . . . ...... . . . . . . I N I N . . . I N A [ m ] (4)where I N is the identity matrix of size N .In order to measure the similarity between layer κ and γ of the MRN, we use the interlayer mutual information I κ,γ [28]: I κ,γ = (cid:88) k [ κ ] (cid:88) k [ γ ] P ( k [ κ ] , k [ γ ] ) log P ( k [ κ ] , k [ γ ] ) P ( k [ κ ] ) P ( k [ γ ] ) (5)where P ( k [ κ ] ) and P ( k [ γ ] ) are the degree distributionsof RNs at layer κ and γ respectively, and P ( k [ κ ] , k [ γ ] )is the joint probability of the existence of nodes whichhas k [ κ ] degree at layer- κ and k [ γ ] at layer- γ . The de-gree distribution P ( k ) is a probability distribution func-tion which holds a general structure information that howmany nodes have each degree. The mutual informationmeasures how much a system is similar to another. In-stead of computing mutual information between originaltime series, degree distributions are considered in Eq. 5.The difference to the commonly used direct mutual in-formation of time series is that I κ,γ does not comparethe probability of states (estimated from the time series)but the topological structure in the phase space basedon the recurrences (using the recurrence matrix). Thissimilarity measure I κ,γ quantifies the information flowbetween the multiplex networks and, thus, the charac-teristical behaviour of the system. The average of thequantity of I κ,γ over all possible pairs of layers of MRN,gives a scalar variable I = (cid:104) I κ,γ (cid:105) κ,γ which captures themean similarity of the degree distributions of the RNs,i.e., the order of coherence in the system.Another measure to quantify the coherence of the orig-inal multivariate system by the MRN is the average edgeoverlap : ω = (cid:80) i (cid:80) j>i (cid:80) κ a [ κ ] ij m (cid:80) i (cid:80) j>i (cid:16) − δ , (cid:80) κ a [ κ ] ij (cid:17) (6)where δ ij is the Kronecker delta symbol. This measurerepresents the average number of identical edges overall layers of the multiplex network [28]. Like the inter-layer mutual information Eq. (5), ω estimates the similar-ity and coherence with averaged existence of overlappedlinks from node- i to j between all layers κ and γ . Notethat ω can take values in the interval [1 /m, i and j occurs in only different layers ω = 1 /m ,i.e., a [ κ ] ij = 1 and a [ γ ] ij = 0 , ∀ γ (cid:54) = κ . If all links are identicalin all layers then ω = 1. Weighted structure.
Another approach is the pro-jection of all layers onto one weighted network represen-tation. Now we consider each single layer of an MRN asa node and weighted edges between nodes κ and γ aredetermined by the quantity I κ,γ , Eq. (5).In order to quantify high-dimensional systems, con-verting multilayer systems to weighted structures is com-putationally a very efficient approach. The adjacencymatrix of the multiplex network A is a giant matrix( N m × N m ), but in the weighted network case the sizeof the matrix is only N × N . Furthermore, the quan-tification of weighted networks is a well-developed analy-sis [51, 52]. Among many measures of weighted networks,we use the clustering coefficient C w and the average pathlength L w in order to detect the transitions between dif-ferent dynamical regimes. The weighted network cluster-ing coefficient is given by C w = 1 m (cid:88) κ k κ ( k κ − (cid:88) γ (cid:88) β ( I κ,γ I κ,β I γ,β ) (7)where k κ is the degree of node κ . The average shortestpath length of the weighted network is L w = (cid:88) κ,γ ∈ V d ( κ, γ ) m ( m − , (8)where V is the set of nodes (layers of multiplex network)in the weighted network and d ( κ, γ ) is the weighted short-est path length between nodes κ and γ .Both MRN measures, I and ω , represent the similari-ties in the linking structures of the RNs and have highervalues for more similar ones. For example, for periodicsystems, even if there is phase difference, the correspond-ing RNs would be rather similar, resulting in high valuesfor I and ω . When the systems are more chaotic and havefinally a more different recurrence structure, I and ω willdecrease. This is similar for C w , since in periodic casesthe number of triangle structures in networks is increas-ing. However, it is opposite for L w because the diameterof a denser network is in general smaller. III. COUPLED MAP LATTICES
As a first application, we consider multi-componentdynamical systems, namely coupled map lattices(CMLs), which are discrete-time models of diffusively coupled oscillators on a ring model of m sites. CMLs arewell-studied dynamical systems that model the behaviorof non-linear systems and exhibit a variety of phenom-ena [53]. x [ κ ] t +1 = (1 − ε ) f ( x [ κ ] t ) + ε (cid:0) f ( x [ κ − t ) + f ( x [ κ +1] t ) (cid:1) , (9)CMLs have been analyzed with various techniques toquantify their very rich dynamics. Now we compare mul-tiplex and regular RN approaches for CMLs of m = 5diffusely coupled chaotic logistic maps f ( x ) = 4 x (1 − x ),which show interesting dynamics in the range of the con-trol parameter ε ∈ [0 , .
4] studied here with an incre-ment of ∆ ε = 0 . N = 15 ,
000 for each value of ε . Inorder to discard transients, we delete the first 10,000 val-ues, resulting in time series consisting of 5,000 values thathave been used for all analyses of the CMLs in this paper.All simulations of CMLs are repeated and averaged over100 realizations, since the system strongly depends onthe initial conditions. In order to compare MRN and RNtechniques, we construct a five-dimensional phase spacefor regular RN and five separated single RN for each layerof MRN.For consistency, we use a recurrence threshold that isproportional to the standard deviation σ of the data [4].In this work, we use (cid:15) = 0 . σ . The logistic map is aone-dimensional dynamical system, therefore embeddingis not required for applying recurrence networks in thelayers of MRN.Although both techniques could detect the transitionsfrom fully developed turbulence (FDT), ε ∈ [0 . , . ε ∈ [0 . , . ε ∈ [0 . , . ε ∈ [0 . , .
4] (Figs. 2, 3). Multiplex network’smeasures ( I and ω ) can recognize every single transitionand, especially ω , can distinguish splitting of trajecto-ries into two attractors in STI at ε = 0 .
36 very clearly(Fig. 2). The MRNs are more sensitive to regime changesthan the RNs as well as the MRNs are applicable onlarge system sizes when RNs are not suitable to deal withthem.
Large systems.
Simultaneous analysis of interact-ing components of a system is very important for deepunderstanding of the underlying dynamics. The RN ap-proach is not an appropriate technique to analyze suchsystems, because while the number of degrees of free-dom is increasing, the size (volume) of its phase spaceis getting larger as well. Therefore, in order to analyzesuch systems, we need very long data sets proportionalto the size of the dimension of the system for establish-ing recurrences. Although the recurrence based analysesnicely deals with analyzing short time series, finding longdata is a common problem in time series analysis. How-
FDT FDTFDT FDTPPSPPS STISTI (a)(b)
FIG. 2. Results for m = 5 CMLs (Eq. 9), (a) average mutualinformation ( I ), and (b) average edge overlap ( w ) of MRN. ever, MRNs overcome this problem since they use eachcomponent of the system as a layer of the giant network.In this application, we use the same parameters of theRN and the isolated dynamics as in the small system sizeone. Fig. 4 presents the results of the MRN for m = 200diffusely coupled logistic maps, which possess one moredifferent regime than the small system, called Brownianmotion of defect (BMD), ε ∈ [0 . , . ε ∈ [0 . , .
12] and ε ∈ [0 . , . ε ∈ [0 . , .
16] and PPS for ε ∈ [0 . , . ε changes the dynamical regime of the maps. According tothe high chaos in the FDT, each RN has a rather differenttopology compared to the others, i.e., all the single layers FDT FDTPPS STIFDT FDTPPS STI (a)(b)
FIG. 3. Results for m = 5 diffusively coupled map lat-tices CMLs, (a) clustering coefficient C , and (b) average pathlength L of the (regular) recurrence network. of the MRN are quite different. The BMD is a less chaoticregime, thus, the single layers of the MRN are less dif-ferent than in the FDT regime (increasing the inter-layersimilarity). The PPS and STI regimes are closer to eachother as both regimes exhibit periodic dynamics but atdifferent scales. Therefore, distinguishing the transitionbetween PPS and STI is the hardest task. Neverthe-less, I and ω are able to detect the corresponding regimetransitions. However, ω performs better, because it is aspatio-temporal measure, i.e., it checks the existence ofedges at the same locations in the different layers wherethe calculation for I does not take the nodes’ identitiesinto account. IV. REAL-WORLD EXAMPLE –MULTIVARIATE PALAEOCLIMATE ANALYSIS
So far we have been testing the performance of thenew MRN method using prototypical models. As a real-world application with changing dynamics and with sev-eral variables, we treat a palaeoclimate problem. Theinvestigation of the linkage between the climate condi-tions and specific environmental responses, as well as ofchanges in these relationships represent an important sci-entific challenge in palaeoclimatology and ecology in or-der to improve the understanding of climate impacts andfeedback mechanisms. Information on past vegetation isuseful to understand the palaeoecological processes link-ing climate and biosphere changes. Vegetation dynamicsdepends strongly on the environmental conditions and (a)(b)
FDT FDTPPSSTIBMDFDT FDTPPSSTIBMD
FIG. 4. Results for m = 200 diffusively coupled map latticesCMLs (a) average mutual information ( I ), (b) average edgeoverlap ( w ) of multiplex network. changes in the vegetation dynamics can, therefore, indi-cate critical changes in the environment (anthropogenicor climate impact).Pollen assemblages collected from lake sediments arecommonly used proxies for palaeoecological and palaeo-climate investigations. Here we focus on the pre-cisely dated Late Pleistocene-Holocene pollen record (last16,500 years (16.5 kyr BP)) from the Sihailongwan Lake(42 ◦ (cid:48) N, 126 ◦ (cid:48) E), located in the Longgang vol-canic field, Jilin Province, NE China (Fig. 6) [56–58].The Sihailongwan Lake is situated near the north-ern edge of the Asian summer monsoon system. Orig-inally the Sihailongwan pollen record consists of 103 dif-ferent pollen taxa, whereby several of them contain re-dundant information. We, therefore, selected seven treepollen taxa (
Pinus koraiensis , Betula , Carpinus , Juglans , Quercus , Ulmus , Fraxinus and one herbaceous genus(
Artemisia ), representing typical regional vegetation ele-ments with different sensitivities on specific environmen-tal conditions during the past 16.5 kyr BP (Fig. 7).Since, the pollen record is irregularly sampled, we firstinterpolate the data into N = 3 ,
500 points leading tothe time resolution ∆ t = t i +1 − t i ≈ .
95 years ∀ i ∈ [0 , N − (a)(b) FDT FDTPPSSTIBMDFDT FDTPPSSTIBMD
FIG. 5. Results for m = 200 diffusively coupled map latticesCMLs (a) clustering coefficient of weighted network C w , (b)average path length of weighted network L w Lake Sihailongwan
Bay of Bengal P a c i f i c O c e a n FIG. 6. Location of the Lake Sihailongwan from where thepollen record for this study was sampled. tivariate pollen data set reveals distinct fluctuations of ω (Fig. 7). A bootstrapping statistical test is appliedto evaluate these variations with respect to the null-hypothesis that multivariate dynamics is constant overtime within 95% confidence [59].The considered measure ω quantifies the amount ofsimilar dynamics in the different pollen taxa. The higher ω , the more similar vary the different vegetation types inthe sense of recurrence properties. This means, it is notnecessary that the pollen abundances go into the same di-rections, but that, e.g., periodical variations are similar. Pinus k.
Betula
Carpinus
Juglans
Selected pollen taxa of Sihailongwan Lake (%) Average edge measure (ω)
400 20
Quercus
Ulmus
Fraxinus
Artemisia competing disturbance factors vegetation dynamics without related climate changevegetation dynamics triggered by warmingno significant impact of the climate event A ge ( k r y B . P . ) additional disturbances byanthropogenic activity successional re-expansion of treestransformation glacial-interglacial forest finished0461012821416 0.18 0.2 0.22 0.24begin Late-glacial begin Holocene8.2 eventbegin human influencebegin Younger Dryas vegetation dynamics likely triggeredby minor climate shift (stronger winter monsoon) FIG. 7. Pollen data set from Sihalongwan Lake (left) and derived average edge measure ω calculated in moving windows(right). The basevalue of the shading (0.199) marks the upper 95% confidence interval, i.e., exceeding this value correspondswith highly significant increase of ω . The grey shaded bands mark important environmental and climate transitions (see text). Low ω values indicate, in contrast, less similar behavior.This could happen, e.g., when the system is disturbedand some vegetation populations recover faster than oth-ers.The pollen assemblages reveal relatively open, cold-dry adapted vegetation at Sihailongwan at the end ofthe pleniglacial. Starting at 16.0 kyr BP, we find a de-crease of ω until 15.0 kyr BP. The apparent decrease ofthe similarity in the floristic response likely results fromrather instable environmental conditions or competingfactors influencing the vegetation dynamics. The insta-ble environment at Sihailongwan is strongly evidenced bythe sediment composition during this time interval. Inparticular, the occurrence of graded event layers with re-worked soil material implies prevailing permafrost condi-tions with reduced seepage and erosion events during thefinal stage of the last glacial period [56]. Moreover, theyoungest Heinrich event (H1; 16.0 and 14.6 kyr BP), mayhave resulted in stronger winter monsoon [60]. By con-trast, the ongoing change of CO from the glacial to theinterglacial level might have individually affected boththe physiology of the prevailing plants, and the vegeta-tion dynamics in an increasingly less moisture-stressedenvironment [61, 62].Subsequent to the H1 event, the environmental con-ditions during the Late-glacial interstadial changed towarmer and moister climate, obviously related with a general change in the regional vegetation to one of moresimilar vegetation dynamics. The short positive excur-sion of ω at 13.6 kyr BP is, however, clearly not relatedto a climate shift, but most likely results from a changein interspecific competition after reaching a critical pop-ulation density of Fraxinus [56].The Younger Dryas temperature decline about 12.7 kyrBP ago is one of the most prominent climate changesobserved in Northern Hemisphere palaeoclimate records.Noticeable changes in the Sihailongwan pollen recordprovide a clear indication of a Younger-Dryas-like coolingin NE China beginning at 12.7 kyr BP. Unlike the abrupttermination of the Younger Dryas in the North Atlanticevidence, the more gradual changes of pollen assemblagesindicate a successive substitution of the prevailing borealwoodland by species-rich broadleaf deciduous forests atthe Sihailongwan lake during the first millennium of theHolocene (11.7 and 10.7 kyr BP).Our MRN analysis reveals significant high values of ω in the period between 12.7 and 10.5 kyr BP, which are in-terrupted by an abrupt drop in ω between about 12.3 and11.7 kyr BP. The timing of this prominent change in veg-etation dynamics corresponds to the recently detected bi-partition of the Younger Dryas at Suigetsu Lake (Japan)[63], but neither the composition of pollen assemblagesnor the sediments of the Sihailongwan sequence reveal asubstantial environmental shift at 12.3 kyr BP. It is likely,that the ecosystem at Sihailongwan Lake may not be assensitive to changes in winter conditions as the Japanesesite [63], so that the climate shift may have influencedthe intrinsic vegetation dynamics at Sihailongwan, butcould not trigger substantial floristic and/or landcoverchanges.At approximately 10.7 kyr BP, the dynamic spreadof thermophilous Juglans and a more gradual increaseof
Quercus mark a shift in the vegetation compositionto oak- and walnut-rich mixed conifer hardwood forests.Subsequent to this transition, the MRN analysis dataindicate persistent vegetation dynamics, but the generalshift to lower ω values implies less similar behaviour ofthe pollen taxa. Interestingly, a shift in ω can be observedaround 8.2 kyr BP, which might be related to the mostprominent drop of the global temperatures during theHolocene. However, as its amplitude does not exceedthat of other Holocene changes in vegetation dynamics,the shift does not prove a significant impact of the 8.2event on the vegetation in NE China.While substantial human impact cannot be tracedby the pollen assemblages, monsoon changes, forestsuccession dynamics and flowering activity may there-fore primarily drive the observed dynamics in thepollen values/vegetation at Sihailongwan throughout theHolocene. Nevertheless, weak indication of farming ac-tivity in conjunction with fires and modest grazing mayexplain the minor shift to lower values from 1.8 yrs BP.Taken together, in our example we find that the pat-terns of vegetation dynamics are related to both intrin-sic vegetation dynamics and external disturbance factorslike climate changes, erosion or human impact. Obvi-ously, the different vegetation types behave more similarwhen one environmental factor acts as the dominant driv-ing force. This is particularly evident at the change tocooler and drier conditions at the beginning of YoungerDryas. Also the rapid mass expansion or the collapse ofone or more tree species can have a similar effect on thedynamics of vegetation, as detected in our example at13.6 kyrs BP and throughout the first millennium of theHolocene. On the other hand, less similar behaviour ofthe vegetation results from weaker or competing distur-bance factors. Our example demonstrates that the newdeveloped MRN technique goes beyond the classical in-terpretation of the pollen amplitude variation as a proxyof environmental conditions. A better understanding ofthis part of the climate-biosphere interaction is of crucialimportance as non-linear feedback mechanisms and tip-ping points cause high uncertainty and an unpredictablefuture for humankind [64, 65]. V. CONCLUSIONS
In this study we have introduced a novel multiplex re-currence network (MRN) approach which allows us toanalyze m -dimensional multivariate data simultaneouslyvia joining m recurrence networks (RNs) together. An-alyzing and understanding the characteristics of singleor low dimensional data with RNs is a known approach[7, 37]. Nevertheless, RN analysis of high dimensionalsystems is not trivial, because the increasing dimensionof the phase space would require longer time series. How-ever, studies of real-world problems are often linked withshort time series and long time series are often not avail-able. Our new approach considers each observation vari-able in its own phase space and combines them later on.Therefore, it is possible to analyze relatively short mul-tivariate time series with MRNs.Our extensive tests of the method have demonstratedthat the MRN approach enables us to detect abruptregime transitions accurately in large dimensional sys-tems. It can be used to understand the underlying dy-namical properties of prototypical models and real-worldapplications.By analysing a pollen data set for the last 16 kyr,we are able to identify periods when different vegeta-tion types behave more similar and periods of less simi-lar behaviour. The changes in the vegetation dynamicscoincide well with known climate transitions and the in-creasing human impact. This analysis goes beyond theclassical interpretation of the pollen amplitude variationas a proxy of environmental conditions. This new viewcan provide insights into plant communities and their dy-namics with respect to climate responses.Especially for some research fields such as palaeocli-mate (pollen taxa, planktonic samples etc), neuroscience(short time crises in EEG), economics (stock markets,currencies), where many dependent signals are collectedat the same time, our method provides a quantitative andobjective way to investigate their dynamics and detectshidden regime changes. ACKNOWLEGMENTS
This work has been supported by German-Israeli Foun-dation for Scientific Research and Development (GIF),GIF Grant No. I-1298-415.13/2015 and the EuropeanUnion’s Horizon 2020 Research and Innovation pro-gramme under the Marie Sk(cid:32)lodowska-Curie grant agree-ment No. 691037 (QUEST). [1] H. Poincar´e, Acta Mathematica , 1 (1890).[2] A. Katok and B. Hasselblatt, Introduction to the Mod-ern Theory of Dynamical Systems (Cambridge UniversityPress, Cambridge, 1995). [3] J.-P. Eckmann, S. Oliffson Kamphorst, and D. Ruelle,Europhysics Letters , 973 (1987).[4] N. Marwan, M. C. Romano, M. Thiel, and J. Kurths,Physics Reports , 237 (2007). [5] L. L. Trulla, A. Giuliani, J. P. Zbilut, and C. L. WebberJr., Physics Letters A , 255 (1996).[6] N. Marwan, N. Wessel, U. Meyerfeldt, A. Schirdewan,and J. Kurths, Physical Review E - Statistical, Nonlinear,and Soft Matter Physics , 1 (2002).[7] J. Donges, R. Donner, M. Trauth, N. Marwan,H. Schellnhuber, and J. Kurths, Proceedings of the Na-tional Academy of Sciences , 20422 (2011).[8] B. H. Jansen, International Journal of Bio-Medical Com-puting , 95 (1991).[9] D. Kaplan, M. Furman, S. Pincus, S. Ryan, L. Lipsitz,and A. Goldberger, Biophysical Journal , 945 (1991).[10] M. A. Riley, R. Balasubramaniam, and M. T. Turvey,Gait & Posture , 65 (1999).[11] Y. Neuman, N. Marwan, and D. Livshitz, Complexity , 28 (2009).[12] S. Carrubba, A. Minagar, A. L. Chesson Jr., C. Frilot II,and A. A.Marino, Neurological Research , 286 (2012).[13] G. R. Richards, Quaternary Science Reviews , 709(1994).[14] N. Marwan, M. H. Trauth, M. Vuille, and J. Kurths,Climate Dynamics , 317 (2003).[15] T. Matcharashvili, T. Chelidze, and J. Peinke, NonlinearDynamics , 399 (2008).[16] I. Ozken, D. Eroglu, T. Stemler, N. Marwan, G. B. Bagci,and J. Kurths, Physical Review E , 062911 (2015).[17] D. Eroglu, F. H. McRobie, I. Ozken, T. Stemler, K.-H. Wyrwoll, S. F. M. Breitenbach, N. Marwan, andJ. Kurths, Nature communications , 1 (2016).[18] N. Asghari, C. Broeg, L. Carone, R. Casas-Miranda,J. C. C. Palacio, I. Csillik, R. Dvorak, F. Freistet-ter, G. Hadjivantsides, H. Hussmann, A. Khramova,M. Khristoforova, I. Khromova, I. Kitiashivilli, S. Ko-zlowski, T. Laakso, T. Laczkowski, D. Lytvinenko,O. Miloni, R. Morishima, A. Moro-Martin, V. Paksyu-tov, A. Pal, V. Patidar, B. Pecnik, O. Peles, J. Pyo,T. Quinn, A. Rodriguez, M. C. Romano, E. Saikia,J. Stadel, M. Thiel, N. Todorovic, D. Veras, E. V. Neto,J. Vilagi, W. von Bloh, R. Zechner, and E. Zhuchkova,Astronomy & Astrophysics , 353 (2004).[19] N. V. Zolotova, D. I. Ponyavin, N. Marwan, andJ. Kurths, Astronomy & Astrophysics , 197 (2009).[20] D. Eroglu, T. K. D. M. Peron, N. Marwan, F. A. Ro-drigues, L. F. Costa, M. Sebek, and Z. Kiss, PhysicalReview E , 1 (2014).[21] N. Marwan, European Physical Journal – Special Topics , 3 (2008).[22] P. Grassberger and I. Procaccia, Physica D , 189 (1983).[23] P. Grassberger and I. Procaccia, Physica D , 34 (1984).[24] J. Zhang and M. Small, Physical Review Letters , 3(2006).[25] M. Small, 2013 IEEE International Symposium on Cir-cuits and Systems (ISCAS2013) , 2509 (2013).[26] L. Lacasa, B. Luque, F. Ballesteros, J. Luque, andJ. C. Nu˜no, Proceedings of the National Academy of Sci-ences of the United States of America , 4972 (2008),arXiv:0810.0920.[27] L. Lacasa, B. Luque, J. Luque, and J. C. Nu˜no, EPL(Europhysics Letters) , 30001 (2009).[28] L. Lacasa, V. Nicosia, and V. Latora, Scientific Reports , 15508 (2015).[29] J. Zhang, J. Sun, X. Luo, K. Zhang, T. Nakamura, andM. Small, Physica D , 2856 (2008). [30] X. Xu, J. Zhang, and M. Small, Proceedings of the Na-tional Academy of Sciences , 19601 (2008).[31] N. Marwan, J. F. Donges, Y. Zou, R. V. Donner, andJ. Kurths, Physics Letters A , 4246 (2009).[32] R. V. Donner, M. Small, J. F. Donges, N. Marwan,Y. Zou, R. Xiang, and J. Kurths, International Jour-nal of Bifurcation and Chaos , 1019 (2011).[33] P. Grassberger and I. Procaccia, Physical Review Letters , 346 (1983).[34] V. Afraimovich, Chaos , 12 (1997).[35] B. Saussol, S. Troubetzkoy, and S. Vaienti, Journal ofStatistical Physics , 623 (2002).[36] B. Saussol and J. Wu, Nonlinearity , 1991 (2003).[37] R. V. Donner, J. Heitzig, J. F. Donges, Y. Zou, N. Mar-wan, and J. Kurths, European Physical Journal B ,653 (2011).[38] O. Afsar, D. Eroglu, N. Marwan, and J. Kurths, EPL(Europhysics Letters) , 10005 (2015).[39] G. M. Ram´ırez ´Avila, A. Gapelyuk, N. Marwan,T. Walther, H. Stepan, J. Kurths, and N. Wessel,Philosophical Transactions of the Royal Society A ,20110623 (2013).[40] J. H. Feldhoff, R. V. Donner, J. F. Donges, N. Mar-wan, and J. Kurths, Physics Letters, Section A: Gen-eral, Atomic and Solid State Physics , 3504 (2012),arXiv:1301.0934.[41] G. Bianconi, Physical Review E , 062806 (2013).[42] V. Nicosia, G. Bianconi, V. Latora, and M. Barthelemy,Physical Review Letters , 058701 (2013),arXiv:1302.7126.[43] M. De Domenico, A. Sol´e-Ribalta, E. Cozzo, M. Kivel¨a,Y. Moreno, M. A. Porter, S. G´omez, and A. Arenas,Physical Review X , 041022 (2013).[44] M. Kivel¨a, A. Arenas, M. Barthelemy, J. P. Gleeson,Y. Moreno, and M. a. Porter, Journal of Complex Net-works , 203 (2014), arXiv:1309.7233.[45] V. Nicosia and V. Latora, Physical Review E , 032805(2015).[46] N. H. Packard, J. P. Crutchfield, J. D. Farmer, and R. S.Shaw, Physical Review Letters , 712 (1980).[47] H. Kantz and T. Schreiber, Nonlinear Time Series Anal-ysis (University Press, Cambridge, 1997).[48] D. Eroglu, N. Marwan, S. Prasad, and J. Kurths, Non-linear Processes in Geophysics , 1085 (2014).[49] C. Ahlstrom, P. Hult, and P. Ask (2006) pp. III688–III691.[50] S. Schinkel, O. Dimigen, and N. Marwan, EuropeanPhysical Journal: Special Topics , 45 (2008).[51] A. Barrat, M. Barth´elemy, R. Pastor-Satorras, andA. Vespignani, Proceedings of the National Academyof Sciences of the United States of America , 3747(2004), arXiv:0311416 [cond-mat].[52] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, andD. U. Hwang, Physics Reports , 175 (2006).[53] K. Kaneko, Chaos , 279 (1992).[54] K. Kaneko, Physica D: Nonlinear Phenomena (1989),10.1016/S0378-4371(00)00431-3.[55] K. Kaneko and I. Tsuda, Complex Systems: Chaos andBeyond (Springer-Verlag Berlin Heidelberg, 2001).[56] M. Stebich, J. Mingram, J. Han, and J. Liu, Global andPlanetary Change , 56 (2009).[57] M. Stebich, K. Rehfeld, F. Schl¨utz, P. E. Tarasov, J. Liu,and J. Mingram, Quaternary Science Reviews , 275 (2015).[58] M. Stebich, K. Rehfeld, F. Schl¨utz, P. E. Tarasov, J. Liu,and J. Mingram, PANGAEA (2015), 10.1594/PAN-GAEA.852704.[59] N. Marwan, S. Schinkel, and J. Kurths, EurophysicsLetters , 20007 (2013).[60] Y. J. Wang, Science , 2345 (2001).[61] S. A. Cowling and M. T. Sykes, Quaternary Research ,237 (1999).[62] E. Monnin, E. J. Steig, U. Siegenthaler, K. Kawamura,J. Schwander, B. Stauffer, T. F. Stocker, D. L. Morse,J.-M. Barnola, B. Bellier, D. Raynaud, and H. Fischer,Earth and Planetary Science Letters , 45 (2004).[63] G. Schlolaut, A. Brauer, T. Nakagawa, H. F. Lamb, J. J.Tyler, R. A. Staff, M. H. Marshall, C. Bronk Ramsey, C. L. Bryant, and P. E. Tarasov, Scientific Reports ,44983 (2017).[64] T. M. Lenton, H. Held, E. Kriegler, J. W. Hall, W. Lucht,S. Rahmstorf, and H. J. Schellnhuber, Proceedings ofthe National Academy of Sciences of the United Statesof America , 1786 (2008).[65] J. Rockstrom, W. Steffen, K. Noone, A. Persson, F. S.Chapin III, E. F. Lambin, T. M. Lenton, M. Scheffer,C. Folke, H. J. Schellnhuber, B. Nykvist, C. A. de Wit,T. Hughes, S. van der Leeuw, H. Rodhe, S. Sorlin, P. K.Snyder, R. Costanza, U. Svedin, M. Falkenmark, L. Karl-berg, R. W. Corell, V. J. Fabry, J. Hansen, B. Walker,D. Liverman, K. Richardson, P. Crutzen, and J. A. Fo-ley, Nature461