Visualization of shortterm heart period variability with network tools as a method for quantifying autonomic drive
Acknowledgements
The authors DM, ZRS, BG, DW and AK acknowledge the financial support of the National Science Centre, Poland, UMO: 2012/06/M/ST2/00480.
This work was partially realized under the SKILLS project of the Foundation for Polish Science (172/UD/SKILLS/2012) and was cofinanced by the European Union through the European Social Fund.
Abstract
Introduction and purpose
Signals from heart transplant recipients can be considered to be a natural source of information for a better understanding of the impact of the autonomic nervous system on the complexity of heart rate variability. Beattobeat heart rate variability can be represented as a network of increments between subsequent intervals, which makes possible the visualization of shortterm heart period fluctuations.
Methods
A network is constructed of vertices representing increments between subsequent intervals, and edges which connect adjacent increments. Two modes of visualization of such a network are proposed. The method described is applied to nocturnal Holter signals recorded from healthy young people and from cardiac transplant recipients. Additionally, the analysis is performed on surrogate data: shuffled RRintervals (to display shortrange dependence), and shuffled phases of the Fourier Transform of RRintervals (to filter out linear dependences).
Results
Important nonlinear properties of autonomic nocturnal regulation in shortterm variability in healthy young persons are associated with increments: accelerations and decelerations of a size greater than about 35 ms. They reveal that large accelerations are more likely antipersistent, while large decelerations are more likely persistent.
Changes in increments in a heart deprived of autonomic supervision are much lower than in a healthy individual, and appear to be maintained around a homeostatic state, but there are indications that this dynamics is nonlinear.
Conclusion
The method is fruitful in the evaluation of the vagal activity  the quantity and quality of the vagal tone  during the nocturnal rest of healthy young people. The method also successfully extracts nonlinear effects related to intrinsic mechanisms of the heart regulation.
Since the vagal part of autonomic regulation is considered to be responsible for nocturnal heart rhythm with large accelerations and decelerations, we hypothesize that vagal activity is a crucial source of the complexity in shortterm heart rate variability.
1 Introduction
Network methods have been successfully used to capture and represent properties of multilevel complex manmade systems [13] and living organisms [3]. The use of network representations in the characterization of time series complexity is a relatively new but quickly developing branch of time series analysis ([10], [12]). The most direct method is to map a time series into a graph in which vertices represent signal values, while edges link values that are consecutive in a signal. The correspondence between time series formed by consecutive cardiac interbeat intervals, so called intervals, and such networks was studied by Campanharo et al. [7]. The topology in these networks appeared as a clique, i.e., each state is reachable from each other in a single step. Understanding interval dynamics arising from a network with such a structure is not straightforward. It appears, in general, that information provided by a network graph strongly depends on the nature of sequences and our knowledge about the underlying dynamics ([12], [13]). Therefore the use of network methods, e.g., visualization or/and structure decomposition, is effective only if they are used in conjunction with other sources of learning. We show how tools developed within the scope of complex networks can be fruitfully applied to the qualification and quantification shortterm heart period dynamics.
Fluctuations in intervals are known to have a scaleinvariant structure which demonstrates fractal ([17], [31], [22]) and multifractal ([15]) properties. These fluctuations appear as a result of many component interactions acting over a wide range of time and space scale. Competing stimuli from the autonomic nervous system are assumed to be the reason for the fractal organization observed in intervals [28]. By observing subsequent changes in intervals, we obtain beattobeat information about the resulting force of these interactions. We have found that important dynamical aspects about the autonomic competitive regulation can be described by changes in intervals, namely by increments ([20], [21]).
Signal increments can be decomposed into their magnitude (absolute value) and their direction (sign). Magnitude and sign analysis has been used to investigate the scaling properties of RRintervals [2]. It has been found that magnitude series are longrange correlated, while sign series are anticorrelated. Furthermore, it has also been shown that during sleep, the strength of these correlations varies depending on the stage of the sleep: rapideyemovement (REM) or other (nonREM) sleep stages [16]. It appears that both the strongest anticorrelations in the sign signals, and largest exponents for longrange correlations for the magnitude signals are in REM sleep. Furthermore, it is assumed that during REM sleep, the nonlinear properties of the heartbeat dynamics are more pronounced.
During sleep, the heart rate is mostly regulated by the autonomic nervous system and is less influenced by physical or mental activity. Moreover, during the night, vagal (parasympathetic) predominance is present, which makes this period a useful state to observe autonomic activity. The classical tools applied to measuring heart rate variability during physiological sleep have shown that the REM stage is characterized by a likely sympathetic predominance associated with a vagal withdrawal, while the opposite trend is observed during nonREM sleep (see [29] for a review). Previous studies have also shown that alternations in nocturnal heart rate variability have clinical importance, e.g., may explain why sudden cardiac death in many cases occurs during the night [30].
Heart transplantation surgery destroys the nerve connections between the organism and the graft — the donor heart is completely denervated, the vagal ganglia at the sinus node are cut off from medulla signals. The regulation is driven by the intrinsic heart mechanisms. The sympathetic control works indirectly at the level of circulating adrenergic hormones in the blood. The lack of vagal activity has the effect, for example, that heart transplant recipients have a resting heart rate higher than the average in healthy people, and their heart rate variability is significantly reduced [5]. The exception is a small respiratory sinus arrhythmia ([25], [11]), which is assumed to be an effect of the intracardiac reflex ([1], [32]) or mechanical stretch of the sinus node. Cardiac reinnervation has been demonstrated in longterm heart transplant recipients ([6], [24], [9]), but it seems that it is limited to the sympathetic nerves. Therefore, the comparison of the nocturnal heart rate variability in healthy young individuals and heart transplant patients gives a unique opportunity to show the impact of autonomic (especially vagal) activity on heart rate regulation.
Following [2], to discover in which way properties of networks constructed from increments demonstrate nonlinear or/and linear dependences among consecutive intervals, we investigate properties of artificially modified interval data [26]. In the following, we argue that network methods are successful in detecting nonlinear properties in the dynamics of autonomic nocturnal regulation in shortterm variability. Two modes of visualization of networks constructed from increments are proposed. The first is based on the handling of a state space. The state space of increments can be modified by a bin size used to code a signal, and by the role of a given vertex as the representation of events occurring in a signal. The second mode relies on the matrix representation of the network on the twodimensional plane. This approach is similar to the accepted method, known as the Poincaré Plot representation of time series for evaluation of heart rate variability. The methods introduced will be applied to nocturnal Holter signals recorded from healthy young people and from cardiac transplant recipients. Thus we obtain a way to filter out the intrinsic heart variability from the autonomic drive and then to quantify complexity in the shortterm interval variability related to nocturnal rest. Changes in increments in a heart deprived of autonomic supervision provide us with insight into beattobeat dependences in forces governing the intrinsic heart dynamics.
2 Method
2.1 Groups and signals studied
Twentyfourhour Holter ECG recordings during a normal sleepwake rhythm were analyzed in two study groups. The first group, the Young, consisted of healthy young volunteers (18 female, 18 male, age 1932). The second group, the HTX, comprised heart transplant patients (surgery at ages 28 to 65). Data from the HTX group was constructed of 20 recordings obtained from 10 patients without any signs of heart graft rejection and who had undergone surgery more than 12 months previously. The Holter recordings were first analyzed using Del Mar Reynolds Impresario software and screened for premature, supraventricular and ventricular beats, missed beats and pauses. Finally, the signals were thoroughly manually corrected and annotated.
Since the analysis concentrates on hours of sleep, the intervals were analyzed from 24:00 to 04:00 in the case of the group and from 22:00 to 05:00 in the case of signals from the group.
As the method of signal preprocessing may impact the results, only long, goodquality fragments of ECG were analyzed. Namely, the parts which contained more than five hundred normaltonormal intervals were extracted and then joined together in a sequence of 10000 intervals. What is worth noting is that all the signals constructed were built from less than seven consistent parts. The number of 10000 points was chosen to ensure proper statistical relevance.
2.2 Signal preprocessing
Our Holter equipment provides data with a 128 Hz sampling frequency. Therefore, the intervals have a limited resolution: 1s/128=7.8125 ms, which can be approximated as ms. This value, denoted as , is accepted as the signal resolution. To decrease the number of different values appearing in a sequence of intervals, we use a binning procedure based on multiples of . Namely, we always set the bin size to as for . The bin quantization described has the effect that intervals take values which are multiples of the bin size . As a consequence, increments are also multiples of , namely, .
The two types of artificially modified cardiac signals were constructed for their further use in statistical tests:

shuffled signals, which were obtained by random shuffling of intervals;

surrogate signals, which were calculated by randomization of phases in the Fourier transform of intervals.
The analysis of shuffled signals tests the presence of dependencies in the signals studied, while the analysis of surrogate signals provides information about whether these dependences are linear or not [26]. Signals of both types were prepared with the help of the TISEAN software [14]. For each cardiac signal, we prepared ten shuffled and ten surrogate signals.
A network was constructed separately for each signal analyzed. Then the mean network for each of the groups of signals studied was established by collecting networks corresponding to the same class of subjects. The confidence interval (CI) for each element of the mean network was also estimated. Calculations were performed with the special software prepared by us.
2.3 Transition network for increments
Let be a time sequence of intervals binned in a . Let be a time sequence of increments, i.e., . Discrete values of the set serve as states in the state space of the transition network indexed by the bin value .
Let denote the number of different states in the network state space, and let us arrange them as follows. If the smallest state is and the greatest state is , then the vertices of a network are labeled consequently as:
(1) 
A directed edge from a vertex to a vertex is established if and represent a pair of consecutive events in a time sequence . Namely, there is a moment in time , for which . If a given pair of increments occurs many times in , the weight of this edge increases accordingly to represent counts of occurrences.
Note that the weight of the edge measures the size of a set consisting of the following events:
(2)  
This means that:

if , hence both increments are negative or both are positive, we observe a run of accelerations or decelerations, accordingly;

if , we observe an alternation between an acceleration and a deceleration or vice versa.
This completes the construction of the transition network from a given time series. The resulting network is directed and weighted. The sums of weights of edges adjacent to a given vertex (total number of and total number of edges) provide the basic network characteristics (called  degree and  degree, respectively), which quantify the role of the vertex in a network. But a network constructed from time series is specific in that each edge from a given vertex is accompanied by an edge to this vertex (with the exception of vertices representing the first and last events in consistent parts of a signal), which implies that the and degrees of each vertex can be considered to be equal to each other. This degree, if normalized by the length of time series, is directly related to the probability that an event represented by occurs in a signal.
The modular structures, also called the community structure, in networks have been shown to be relevant to the understanding of the structure and dynamics of the system studied [13]. However, this problem has been found to be difficult and has not yet been satisfactorily solved ([18], [12]). Here we propose to investigate modularity in the transition network by the socalled core graph ([27], [18]). The core graph is constructed from a given network by the removal of all the vertices with a probability less than . Then all the edges which connected these deleted vertices with the other parts of a network are removed. The sum of normalized weights in the resulting subgraph is called the volume of the core graph. A decay in this volume with an increasing value is known as the network disintegration [20].
2.4 Transition network graph
Visualization of a transition network is challenging because usually a transition network consists of many vertices which are densely, often completely, interconnected. The plot of such a network may be barely readable. Therefore the graph organization requires a special effort.
There are parameters in the method which have to be thoroughly tuned:

— the bin size which is used in preprocessing RRintervals and which determines the number of states in the state space;

— the probability of neglected events, which also allows a reduction in the number of states.
Since states in the state space are ordered according to the values of their labels, see Eq. 1, we plot them in a circle arranged clockwise according to increasing value of the vertex label from to . Moreover, if we call two vertices neighboring when the magnitude of difference between their labels is equal to , then we can code transitions between neighboring vertices by colors.
Here we use the following color code:

violet to mark neighboring vertices, i.e., loops describing events of two adjacent accelerations or decelerations of the same value; the case of the loop denotes the situation when three consecutive intervals have the same value;

green to mark neighboring vertices, i.e., transitions to the nearest neighbors in the state space; they denote the smallest possible observable changes in subsequent accelerations and/or decelerations within a given binning;

blue to mark transitions for neighboring vertices;

red to mark transitions between neighboring vertices;

yellow to mark the transitions linking neighboring vertices;

black to mark transitions of a size larger than which, for example, in the case of ms means changes of at least of 40 ms.
Moreover, we use also the width of an edge to visualize the weight of a given transition.
In the following, we use the popular software PAJEK [4] to plot graphs of transition networks.
2.5 Matrix representation of a transition network
Adjacency matrices and transition matrices are standard representations of any network [12]. For a transition network with vertices, the adjacency matrix is a matrix. The number of the outgoing edges from vertex to vertex is counted and designated as . If there is no edge between these vertices, then . Hence
In the following, we normalize counts by the total number of events. As a results, means the probability of a given transition. Referring to a signal with increments stands for the probability that value occurs after in a signal. The transition matrix is obtained by dividing elements of each row of the matrix by the total weight of vertex . Thus,
Therefore describes a Markov walk on a network where a walker being in vertex moves to with a probability .
It appears that the contour plots of adjacency and transition matrices provide a readable visualization of transition networks obtained from increments even in the case when the bin is equal to the resolution of signals. Manipulations in the range of the axes allow one to pass through the whole range of values obtained. However, when departing from , less probable events with larger standard errors are estimated. The large variations between neighboring points lead to an unclear picture if the plots are constructed from signals binned in a small bin size. Therefore, in the following, we limit our interest to the ranges of increments which contain the most probable events.
Each point of the contour plots can be resolved into the three interval patterns as described by formula (2). Moreover, these events can be translated into codes of shortterm variability proposed by Porta et al. [23]:  0 variation,  1 variation,  2 likely variations and  2 unlikely variations. The relation between three interval patterns and their description by 0, 1 or 2 variations is shown in Fig. 1.
3 Results
3.1 Graphs of transition networks for the Young
In the presentation of our results, we first refer to some graphs which demonstrate the possibilities of the method of visualization introduced. In particular, these graphs clarify the influence of the parameters and on the graph shape.
In Fig. 2, there are six graphs which represent networks obtained from signals of the Young group (A)(D) and their surrogates (E), (F). The left column networks were prepared with signals binned in ms, while the right column graphs come from signals binned in ms. The first and third row graphs show vertices which correspond to deceleration/acceleration events appearing in a signal with probability . The second row graphs contain vertices representing more rare events because only vertices with are plotted.
All the graphs are complete, which means that each vertex is connected to all others. The color of the edges corresponds to size of the change in a way defined in Sec. 2.4. The weight of an edge is represented by their width. The global weights of the most important transitions are described additionally by giving their probability value in per cents and CI.
We see that the transition is the most frequent in all the graphs, namely the transition from nochange to nochange occurs the most often. However the probability of observing such an event is different depending on the binning applied to the signals. The bin size works like a magnifying glass, making it possible to perceive a single event in greater detail. For example, vertex in the network constructed from signals binned at 64 ms (B) represents the whole core graph (A) obtained at .
Fig. 2 also gives graphs obtained from networks constructed from surrogates of cardiac signals of the Young group. They look similar to graphs obtained from original signals, which might indicate that dynamical relationships among the most important events are of a linear type. We do not present graphs obtained for the shuffled signals of the Young group because they are barely readable.
Moreover, when signals are binned in , the core graph of consists of only a few vertices (A), while the number of vertices grows sharply if we decrease the probability to (C). Hence, graph (A) rather superficially describes the dynamics of the system, since a slight change in any parameter strongly influences the graph shape. This is different in the case of signals binned in ms; compare (B) to (D). Therefore the volume of a given core graph gives us an important message about how dynamical forces presented by a graph are meaningful for the overall dynamics of the system studied. Here, the volumes of the cores presented in Fig. 2 are (A): , (B): , (C): , (D): , (E): , (F): (mean ).
3.2 Graphs of transition networks for Htx
The signals with intervals obtained from patients after are plain in the sense that consecutive increments do not differ much. Therefore, it becomes possible to present readable core graphs even when is low, e.g., , and the bin size is equal to the signal resolution ms. Moreover, graphs obtained from modified HTX signals, their surrogates and when shuffled, are also clear. All these graphs are shown in Fig. 3.
Figure 3 demonstrates the following:

How different the graph obtained from HTX_shuffled is from all the other graphs. This indicates that adjacent intervals resulting from intrinsic mechanisms of variability are strongly correlated;

There is similarity between graphs representing signals HTX and HTX_surrogates. However we should note that there is a noticeable difference in the probability of transitions between nochange vertex and vertices . This observation may indicate that dependences between intervals cannot be approximated by some linear relations;

There is a similarity between graphs of HTX and those of the Young which are binned at an interval 8 times longer than the HTX series; see Fig. 2(D). This approximate similarity may give rise to the conjecture that the variability driven by autonomic regulation enhances eight times the magnitude of fluctuations of the intrinsic heart period variability.
3.3 Network disintegration
The decay of the volume of the core when is increasing can provide us information about the collective structures in the system dynamics. In Fig. 4, we show these decays for all signals studied binned in ms.
It appears that the network formed from signals of the HTX group is significantly more resistant to the vertex removal than all the other networks. This network decays the slowest. Moreover its core volume stays unchanged at for . This firmed core is built from transitions between the three vertices and . A similar network core also emerges from the network constructed from surrogate signals of the group. Its volume of is significantly lower than the volume of the network constructed from original cardiac signals (by MannWhitney test, ) for all in the interval described. This fact may indicate that the network organization resulting from the cardiac signals relies on important nonlinear dependences.
On the contrary a similar property does not hold in the case of the disintegration of the network constructed from the Young group signals. The volume of the core decays in the same way in both networks: the network made from cardiac signals and the network made from surrogate signals. The disintegration of both networks goes fast, for example at , the network volume is about in both cases.
Obviously, the networks obtained from shuffled signals decay more quickly when compared to the networks produced from cardiac series.
3.4 Adjacency and transition matrices for the Young
In Fig. 5, we show contour plots of both matrices: adjacency and transition obtained from signals recorded from the Young group. Together, the plots obtained from surrogates of these signals and shuffled intervals are presented. All the plots represent signals binned at and for changes smaller than 100 ms for and smaller than 70 ms in case of .
From the plots in Fig. 5, it is evident that shuffled intervals provide different matrix representations. The meaningful smaller probability of events corresponding to small accelerations or decelerations (i.e., probability of events around ) is a noticeable effect of the independence of intervals. Moreover, the symmetry in these matrices is different from symmetries which are present in other plots. These symmetric features can be explained by elementary counting of the sets built from the three interval events of the specific type. Since the state space of intervals is discrete and limited, the values of for any take one of the values from the set of different values:
(3) 
For simplicity, let us assume that all events of Eq. 3 are equally probable, . Then the number of possible monotonically growing threeelement sequences (, , ) with constructed from these values, namely patterns of the 2LV type, is . On the other hand, from each sequence of the 2LV type, one can construct two different threeelement sequences of the 2UV type. Hence the total probability of events of the 2UV type is twice as large as of the 2LV type.
In the case when events from the list described in Eq. 3 are not distributed uniformly, the calculations are more demanding but finally they lead to the conclusion that if , then events are more probable than or .
In addition, comparing matrices obtained from cardiac signals and their surrogates, Fig. 5 first row, we see:

there is a small deficiency of events close to in cardiac signals when compared to signals with surrogates;

there are three regions in in which cardiac signals dominate their surrogates. They can be described as three interval events of the type:

2LV but related only to large decelerations;

1V when after a small change, a large acceleration occurs;

2UV but only for a large deceleration after a large acceleration.
By large changes above, we mean increments greater than 30 ms.

3.5 Adjacency and transition matrices for Htx
In the absence of any influence of the autonomic nervous system, the network representation of increments consists of considerably fewer vertices than for a typical healthy person, as has been already shown in Fig. 3. Fig. 6 presents results aimed at widening our understanding of the nonlinear effects of the intrinsic mechanisms controlling the heart contractions. Note that the contour plots in Fig. 6 are in different scales from the plots representing signals of the Young group in Fig. 5.
We see in Fig. 6 similar symmetric features in all three plots with matrices, namely that alternating changes with are more dominant than monotonic changes. Following our discussion in the previous subsection about the imprints of randomness, this observation may imply stochastic independence of the underlying dynamics. However, the cardiac dynamics is more concentrated around transitions from a nochange event to the smallest increments possible, namely to ms, than if it resulted from random sources. Furthermore, a closer analysis of (compare Figs 3(A) and 3(B)) reveals that the system represented by the cardiac signals is less likely to stay in the nochange vertex, and changes of size ms between vertices representing transitions of and ms in size occur more often in the cardiac signals [] than in their surrogate series [].
Additionally, a comparison between the corresponding matrices provides important distinctions between cardiac signals and their surrogates in the system reaction after the larger accelerations, namely, if ms. It appears that when accelerating, the system is more resistant to a pendulumlike reaction. This has the effect that the interval is able to retain the shorter rhythm for the next contraction.
In Fig. 7, we show plots of differences between matrix graphs arising from cardiac signals of the Young group when the signals are binned in ms, and the HTX group. We see that a similarity is apparent between the dynamics underlying these two systems. The basic distinction relies on events of the 2LV type, where two subsequent accelerations or two subsequent decelerations are involved. Thus it can be hypothesized that independently of the scale, the purpose of the autonomic regulation seems to be to speed up or slow down the heart period. Since such persistence could be involved in some overall purpose like actual bodily needs, the next supposition can be formulated as follows. While the intrinsic heart control mechanisms are devoted to keeping the homeostasis, the control of the autonomic nervous system aims at satisfying physical demands.
4 Conclusions
Network structure methods are able to visualize, describe and differentiate heart rate dynamics in healthy young subjects and HTX patients. Using these methods, the general dynamical properties of heart rate can be defined as correlated or not correlated (employing the comparison of raw and shuffled signals) and linearly or nonlinearly correlated (comparing raw signals and signals with shuffled phases of Fourier Transform).
The essential feature of complex dependencies in nocturnal heart rhythm in our group of healthy young persons is related to large increments, both decelerations and accelerations. This feature manifests itself in that large accelerations are more likely antipersistent, while large decelerations are more likely persistent. This observation also seems to be an important indicator of healthy heart rate.
Moreover, since the vagal part of autonomic regulation is considered responsible for large increments, we may hypothesize that vagal activity is a crucial source of complexity in shortterm heart rate variability. In healthy young individuals, the change in vagal tone during sleep (e.g. change from high vagal activity to its withdrawal between nonREM and REM sleep stages) allows us to observe the specific patterns of heart rate dynamics. We interpret the nonlinear relationship observed between consecutive accelerations and decelerations in the case of bigger changes (accelerations and decelerations of more than 35 ms) as an effect of vagal activity.
Although in HTX patients’ heart rate regulation is mostly intrinsic with no autonomic control, the relationship between consecutive accelerations and decelerations is also observed, but in this case the scale of changes is much lower. increments vary as fluctuations around a homeostatic state. However, the organization of this homeostatic state in the case of raw signals shows that it involves dynamical forces more strongly than if the dynamics were driven by linear forces only. In posttransplant patients, the nonlinear dependencies are also characterized by the appearance of sequences made of bigger ( ms) accelerations followed by smaller decelerations (less than 10 ms). This means that an increase in heart rate is not so effective as in healthy individuals but is still possible. We hypothesize that this pattern of heart rate in HTX patients may be a result of gradual sympathetic reinnervation.
References
 [1] Armour, J.A., 2008. Potential clinical relevance of the ’little brain’ on the mammalian heart. Experimental Physiology 95:165–76.
 [2] Ashkenazy, Y., Ivanov, P.Ch., Havlin, S., and Peng, C.K., 2001. Magnitude and sign correlations in heartbeat functions. Phys. Rev. Lett. 86:1900–03.
 [3] Bashan, A., Bartsch, R.P., Kantelhardt, J.W., Havlin, S., and P.C. Ivanov. 2012. Network physiology reveals relations between network topology and physiological function. Nature Communications 3:702.
 [4] Batagelj, V., and A. Mrvar. 1998. Pajek: a program for large network analysis. Connections 21:47–57.
 [5] Bigger Jr, J. T., Steinman, R. S., Rolnitzky, L.M., Fleiss, J.L., Albrecht, P., and R.J. Cohen. 1996. Power law behavior of RRinterval variability in healthy middleaged persons, patients with recent acute myocardial infarction, and patients with heart transplants. Circulation 93:2142–51.
 [6] van De Borne, P., Neubauer. J., Rahnama, M., et al. 2001. Differential characteristics of neural circulatory control: early versus late after cardiac transplantation. Circulation 104(15):1809–13.
 [7] Campanharo, A.S.L.O., Sirer, M.I., Malmgren, R.D., Ramos, F.M., and L.A.N. Amaral. 2011. Duality between time series and networks. PLoS ONE 6(8):e233378.
 [8] Cohen, M. C., and J.A. Taylor. 2002. Shortterm cardiovascular oscillations in man: measuring and modelling the physiologies. J. Physiology 542(3):669–83.
 [9] Cornelissen, V. A., Vanhaecke, J., Aubert, A. E., and R. H. Fagard. 2012. Heart rate variability after heart transplantation: A 10year longitudinal followup study. J. Cardiology 59:22024.
 [10] Donner, R.V., Zou, Y., Donges, F.F., Marwan, N., and J. Kurths. 2010. Recurrence networks  a novel paradigm for nonlinear time series analysis. New J. Physics 12:140. http://iopscience.iop.org/13672630/12/3/033025.
 [11] Eckberg, D.L. 2003 The human respiratory gate. J. Physiol. 548:339–52.
 [12] Fortunato, S. 2010. Community detection in graphs. Physics Reports 486:75–174.
 [13] Havlin, S., Kenett, D.Y., BenJacob, E., et al. 2012. Challenges in network science: Applications to infrastructures, climate, social systems and economics. Eur.Phys.J. Special Topics 214:273–293.

[14]
Hegger, R., Kantz, H., and T. Schreiber. 1999.
Chaos 9:413–34. Software accessible from
w.mpipksdresden.mpg.de/ tisean/Tisean_3.0.1 .
 [15] Ivanov, P.C., Amaral, L.A., Goldberger, A.L., Havlin, S., Rosenblum, M.G., Struzik, Z.R., and H.E. Stanley. 1999. Multifractality in human heartbeat dynamics. Nature 399:461–465.
 [16] Kantelhardt, J.W., Ashkenazy, Y., Ivanov, P.Ch., et al. 2002. Characterization of sleep stages by correlations in the magnitude and sign of heartbeat increments. Phys. Rev. E 65:051908.
 [17] Kobayashi, M., and T. Musha. 1982. 1/f fluctuation of heartbeat period. IEEE Trans Biomed Eng 29:456–457.
 [18] Kumpula, J.M., Kivelä, M., Kaski, K., and J. Saramäki. 2008. A sequential algorithm for fast clique percolation. Phys. Rev. E 78:026109.
 [19] Makowiec, D., Struzik, Z.R., Graff, B. et al. 2013. Complexity of the heart rhythm after heart transplantation by entropy of transition network for RRincrements of time intervals between heartbeats. in Conf. Proc. IEEE Eng. Med. Biol. Sci. 2013:6127–30.
 [20] Makowiec, D., Struzik, Z.R., Graff, B., et al. 2013. Community structure in networks representation of increments in beattobeat time intervals of the heart in patients after heart transplantation. Acta Phys. Pol. B 44:1219–33.
 [21] Makowiec, D., Struzik, Z.R., Graff, B., ZarczynskaBuchowiecka, M., and J. Wdowczyk. 2014. Transition network entropy in characterization of complexity of the heart rhythm after heart transplantation. Acta Phys. Pol. B: to appear .
 [22] Peng, C.K., Havlin, S., Stanley, H.E., and A.L. Goldberger. 1995. Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series. Chaos 5:82–87.
 [23] Porta, A., Tobaldini, E., Guzzetti, S., Furlan, R., Montano, N., and T. GnecchiRuscone. 2007. Assessment of cardiac autonomic modulation during graded headup tilt by symbolic analysis of heart rate variability. Am J Physiol Heart Circ Physiol 293:H702–08.
 [24] Porta, A., Catai, A.M., Takahashi, A.C.M., et al. 2011. Causal relationship between heeart period and systolic arterial pressure during graded headup tilt. Am. J. Physiol. Regul. Integr. Comp. Physiol 300:R378–86.
 [25] Radaelli, A., Valle, F., Falcone, C. et al. 1996. Determinants of heart rate variability in heart transplanted subjects during physical exercise. European Heart J. 17:462–71.
 [26] Schreiber, T., and A. Schmitz. 2000. Surrogate time series. Physica D 142:346–82.
 [27] Seidman, S.B. 1983. Network structure and minimum degree. Social Networks 5:269–87.
 [28] Struzik, Z.R., Hayano, J., Sakata, S., Kwak, S., and Y. Yamamoto. 2004. 1/f scaling in heart rate requires antagonistic autonomic control. Phys Rev E Stat Nonlin Soft Matter Phys 70:050901.
 [29] Tobaldini, E., Nobili, L., Strada. S., Casali, K.R., Braghiroli, A., and N. Montano. 2013. Heart rate variability in normal and pathological sleep. Frontiers in Physiology 4:294.
 [30] Vanoli, E., Adamson, P.B., BaLin, M.P.H., Pinna, G.D., Lazarra, R., and W.C. Orr. 1995. Heart rate variability during specific sleep stages. A comparison of healthy subjects with patients after myocardial infarction. Circulation 91(7):1918–22.
 [31] Yamamoto, Y., and Hughson, R.L., 1991. Coarsegraining spectral analysis: new method for studying heart rate variability. J Appl Physiol 71:1143–1150.
 [32] Zarzoso, M., Ryseavaite, K., Milstein, M.L. et al. 2013. Nerves projecting from the intrinsic cardiac ganglia of the pulmonary veins modulate sinoatrial node pacemaker function. Cardiovascular Research 99:566–75.