Abstract
The possibility of computing a cardiac age on the basis of spectral analysis of healthy individual tachograms was confirmed and facilitated by the use of a nonlinear technique: recurrence quantification analysis. The age of 112 subjects was predicted by this technique in terms of a progressive increase in the deterministic character of the heartbeat. This result confirms the “randomwalk” character of the heartbeat as predicted by the terminal dynamics paradigm, thus allowing for a simple and comprehensive model of the effect of aging on cardiac dynamics: as age progresses, heart rate dynamics become increasingly predictable (constrained) on a beattobeat basis. This implies a basically stochastic nature of heart rate dynamics, probably reflecting the continuous adjustments to an unpredictable internal environment.
 recurrence quantification analysis
 nonlinear dynamics
 heart rate variability
 terminal dynamics
in a recent work (5) we derived a “cardiac age” estimator based on the spectral analysis of the tachograms of healthy individuals in resting and tilt phases. This estimator was derived by means of a multiple linear regression model using the chronological age of individuals as dependent variable and 4 principal components extracted from 24 variables stemming from the Fourier analysis of the tachograms (5). The cardiac age estimator showed a remarkable correlation with the actual age of individuals (r = 0.74).
In another work (7) we demonstrated the possibility of a straightforward representation of cardiac dynamics in terms of a firstorder Markov model. According to this model, heartbeat dynamics could be considered a random walk, where the system at each beat was presented with three alternatives:1) remain in the same state (i.e., having a beat of a length very similar to the previous one),2) shift to the higher class of beat duration, or 3) shift to the lower class of beat duration.
Previously (7), we demonstrated that this model was effective for describing the nonoscillatory component of heart rate variability (HRV); in this work we wanted to prove the usefulness of the model to provide a simple description of the effect of senescence on HRV.
According to this model, we hypothesized that senescence, reducing the system plasticity, makes the dynamics selectalternative 1 overalternative 2 or3, with a consequent increase in predictability (i.e., deterministic, ruleobeying character) of the tachograms.
To assess this model, we analyzed the tachograms with a nonlinear analysis tool [recurrence quantification analysis (RQA)^{1}] that, in our opinion, given its independence of stationarity and its intrinsic discrete character, is particularly suited for the analysis of HRV (10, 15) and allows for a direct quantification of the deterministic character of the studied series. RQA allowed us to build a cardiac age indicator consistent with that derived from Fourier analysis, requiring only one component, which allowed for a straightforward interpretation of the effect of senescence on HRV in terms of an increase in the deterministic, ruleobeying character of beating. The terminal dynamics paradigm (7, 14, 17) offers a physical basis for the observed results.
METHODS
Data collection protocol.
To rely on a normal healthy population, the following exclusion criteria were adopted (5): history or actual evidence of cardiovascular, respiratory, renal, liver, or gastrointestinal disease; diastolic blood pressure >90 mmHg or systolic blood pressure >150 mmHg; body mass index >26 kg/m^{2}; smoking; diabetes; plasma cholesterol >5.7 mmol/l; arrhythmias or conduction abnormalities; echocardiographic evidence of wall motion abnormalities or valvular disease; and ultrasound evidence of significant carotid stenosis. Of the 455 subjects initially enrolled for this study, only 141 completed the study protocol; data from 112 of 141 subjects who completed the initial study (5) were analyzed in the present study (29 subjects were lost for the present study because of accidental loss of relative data resulting from a computer crash). Both populations had very similar statistical characteristics.
Heart rate recordings in all subjects took place at 8:30 AM in a quiet and comfortable environment according to the protocol described previously (5). The rest and tilt tachograms had an average length of 512 beats.
Electrocardiographic signals were digitized, stored on a hard disk, and sampled at a rate of 500 Hz, with 12 precision bits. Each QRS complex (lead II) and each trigger point were also subsequently checked by an expert cardiologist, who eliminated tracings with >1% premature beats or artifacts due to nonvoluntary motions.
The statistical analyses were carried out by means of the SAS package (6.12 release) for the IBM RISC 6000 (Unix system).
RQA.
RQA was first introduced in physics by Eckmann et al. in 1987 (6) as a purely graphic technique. Five years later, Zbilut and Webber (16) enhanced the technique by defining five nonlinear quantitative descriptors of the recurrence plot that were found to be diagnostically useful in the quantitative assessment of time series structure in fields ranging from molecular dynamics to physiology (8, 10, 12, 13).
This technique has been demonstrated to be particularly useful in quantifying transient behavior far from equilibrium (12); this feature is especially important in dealing with complex systems that can hardly be considered at equilibrium such as heart rate.
The RQA is based on the computation of a distance matrix between the rows (epoch) of the embedding matrix of the tachogram at unitary lag (given the discrete character of the heartbeat). This distance matrix is computed making use of Euclidean metrics; after this first step, the distance matrix, having as rows and columns the subsequent epochs of the series of length equal to the chosen embedding dimension (in this case set to 15), is colored, with darkening of the pixels located at specific (i,j) coordinates, which correspond to a distance betweeni andj rows (epochs) lower than a predetermined radius (for details see Refs. 8, 13, and 16).
The features of the distance function make the plot symmetrical (Di,j) =D( j,i) and with a darkened main diagonal corresponding to the identity line (Di,j = 0,j =i). The darkened points individuate the recurrences (recurrent points) of the dynamic process, and the plot can be considered a global picture of the autocorrelation structure of the system at hand (6, 16).
In other words, a recurrence plot visualizes the distance matrix between the epochs (rows) of the embedding matrix, which, in turn, represent the autocorrelation in the signal at all the possible time scales. In fact, it is important to note that the distance is computed for all the possible pairs of epochs: the elements near the principal diagonal of the plot corresponding to shortrange correlations (the diagonal marks the identity in time) and the longrange correlations corresponding to points distant from the main diagonal.
Besides the global impression given by the graphic appearance of the plot (see Fig. 1 for the RQA plot of 2 typical tachograms), the indexes developed by Trulla et al. (12) and Webber and Zbilut (13, 16) allow for a quantitative description of the recurrence structure of the plot. The RQA descriptors (4 of 5) used in this work are as follows:
1) The percentage of recurrence (REC) quantifies the percentage of the plot occupied by recurrent points. It corresponds to the proportion of recurrent pairs over all the possible pairs of epochs or, equivalently, the proportion of pairwise distances below the chosen radius among all the computed distances.
2) The percentage of determinism (DET) is the percentage of recurrent points that appear in sequence, forming diagonal line structures in the distance matrix. A line is defined a priori as each sequence that is equal to or longer than a predetermined threshold length. In this case, given the strong tendency of subsequent beats to be similar (7), the threshold was set to 5. (It is important not to confound this last threshold with the radius: although the radius defines a distance below which two epochs are considered recurrent, the DET threshold defines the minimum number of consecutive recurrent points that can be considered a line and then scored as deterministic.) DET corresponds to the number of patches of recurrent behavior in the studied series, i.e., to portions of the state space in which the system resides for a longer duration than expected by chance alone [quasiattractors (14, 17)]. This is a very crucial point: a recurrence can in principle be observed by chance whenever the system explores two nearby points of its state space. On the contrary, the observation of recurrent points consecutive in time (and then forming lines parallel to the main diagonal) is an important signature of deterministic structuring (6, 12, 13, 16). The formal superposition between determinism and Lyapunov exponents (6, 12) is proof of this point.
3) The entropy (ENT) is defined in terms of the ShannonWeaver formula for information entropy (11, 13) computed over the distribution of length of the lines of recurrent points and measures the richness of deterministic structuring of the series.
4) The maximal line (MAXLINE) is simply the length (in terms of consecutive points) of the longest recurrent line (sequence of consecutive recurrent beats) in the plot. MAXLINE accurately predicted (r = 0.93) the value of the maximum Lyapunov exponent in a logistic map from a regular to a chaotic regimen (12).
All these variables were obtained for resting and tilted phases. The two states are indicated by the suffixes 1 and 2 for rest and tilt phases, respectively.
Data analysis.
The RQA variables were supplemented by the computation of mean and standard deviation of the tachograms and by the differences between rest and tilt phases for all the computed indexes (the differences are shown by the prefix “dif” followed by the name of the corresponding variable).
As a first step in the analysis, the Pearson correlation coefficients (r) between the tachogram variables and the chronological age of patients (AGE) were computed. All the computed variables, except the dif variables, given their low correlation with age, were then submitted to a principal component analysis (PCA; see Refs. 3 and 5 for details of this widely used statistical technique). Here it is sufficient to stress that this technique allows us to build general synthetic indexes collecting all the relevant information present in the original variables. These indexes are by construction independent between them, allowing separation of the different aspects of the studied data set. The extracted components were then correlated with age.
Only the first, most informative component (PC1) was used to build a simple bivariate regression model to predict the actual age of individuals, with AGE as dependent variable and PC1 as regressor (5).
The estimate of chronological age given by this model was called cardiac age (ESTAGE).
To compare the results of this work with the previous study (5), the data set was classified by means of a cluster analysis procedure (see Ref. 5 for explanation of the kmeans clustering technique) into homogeneous classes of age. AGE was compared with ESTAGE at the cluster level to appreciate the general trend of cardiac regulation aging over the years.
RESULTS
Simple correlations of the measured parameters with age.
The entire population consisted of 112 healthy individuals 20–85 yr of age [52.42 ± 17.95 (SD) yr]. Table1 shows the Pearson correlation coefficients of the studied parameters with the AGE variable. All the RQA parameters and SD values exhibit high and significant correlations with age, whereas the heart rate (MEAN) is uncorrelated with age. The difference parameters are generally uncorrelated with age, except for weak correlations for difMEAN and difSD. Given this result, the dif variables were excluded from subsequent analyses.
The comparatively high correlations of DET and ENT variables with AGE are worth noting (Table 1); in the previous fast Fourier transform analysis, no single variable had a correlation coefficient with age >0.50 (5).
Given the high degree of correlation among the observed variables, the data set was analyzed by means of PCA to obtain a wellconditioned model for estimating chronological age and to obtain useful clues for the dynamic bases of the effect of aging on heartbeat regulation.
PCA.
The results of PCA are reported in Table 2in terms of factor loadings (correlation coefficients between experimental variables and components), percentage of variability explained by each single component, and correlation coefficients of each component with AGE as an external variable.
A fourcomponent solution explained 85% of total variability and represents an exhaustive summary of the data set structure. The first component (PC1) is by far the most relevant, explaining 52% of total variability. As shown in Table 2, PC1 is negatively correlated with SD and positively correlated with all the RQA parameters, all of them measuring the degree of structuring of the tachograms. The variables most heavily loaded on PC1 (and consequently most relevant for its interpretation) are the degree of determinism (DET) of the series and the richness of deterministic structure (ENT) in resting and tilted phases.
The heart rate is independent of PC1; in fact, MEAN1 and MEAN2 are loaded on PC2, which is orthogonal (independent) to PC1 and explains ∼14% of total variability. This is an important result: PCA allowed us to separate the mean value of the studied series (heart rate) from its shape parameters.
PC3 and PC4 account for peculiarities of rest and tilt phases, respectively. These peculiarities are expressed by the portion of information conveyed by REC1 and REC2 independent of the information carried by PC1. The only variable significantly loaded on PC3 is REC2, whereas REC1 is the only relevant parameter for the interpretation of PC4 (Table 2).
When examining a factorloading profile, we must consider that each component is orthogonal to the others and that the components are extracted by the algorithm in decreasing order of relative percentage of explained variance (3). This means that higherorder components explain the portion of information not accounted for by the lowerorder components. A useful analogy is thinking of the subsequent components as explaining the “residuals” of the first components.
Generally speaking, REC is linked to the amount of recurrent behavior of the studied dynamics, and, in the case of heartbeat, the amount of recurrent behavior is due to two wellseparated factors:1) the probability for a single beat to be similar to the previous one and2) the presence of cyclic behavior, i.e., the sympathetic and parasympathetic frequencies of HRV. Althoughtype 1 effects are explained by PC1 (as is evident by the DET, ENT, and MAXLINE high loadings),type 2 effects are explained by PC3 and PC4.
In summary, given the general role of “determinism quantifier” played by PC1, we can think of PC3 and PC4 in terms of “degree of cyclic behavior” of the tachograms.
Overall, PC1 can be considered the relative ruleobeying character of the heartbeat, PC2 the heart rate, and PC3 and PC4 the relative cyclic character of tilt and rest phases, respectively.
Only PC1 was significantly correlated with age (r = 0.74,P < 0.0001); thus the estimation of a cardiac age (5) was based on this single composite variable.
Setup of a cardiac age estimator.
The high correlation between AGE and PC1 allowed us to build a cardiac age index in terms of the best linear model (in a leastsquares sense) able to predict the actual chronological age of the individuals on the basis of the relative PC1 score.
The application of a leastsquares regression model to our data set generated the subsequent linear expression for cardiac age (ESTAGE)
However, an important result, from a practical and a theoretical point of view, is the ability of RQA to obtain the same results of Fourier analysis by only one regressor variable (instead of four). In fact, it is worth noting that PC1 was computed independently of its relation with age, which was introduced “a posteriori” as an external variable to be predicted. PC1 points to the direction of maximal variability of the data set relative to the tachogram description, so it is the best (most informative) cumulative index of the single tachograms. Its relation with age is then an “emergent” property that gives PC1 the character of general explanation of a senescence effect on cardiac dynamics. The advantage of using PC1 to estimate a cardiac age, instead of the original variables, is thus linked to statistical issues (using only one regressor instead of many) and to mechanistic issues, given that the principal components point to conceptually autonomous aspects of the investigated data set (3).
Given the observed factorloading profile, PC1 can be considered a global “inverse complexity score” measuring the degree of predictability of the tachograms.
This fact allowed us to greatly simplify the computation of cardiac age for the single subjects and to give a global explanation to the empirical results of the previous work in terms of increase in determinism as the basic effect of age (Fig. 1).
Clustering of ages.
The application of the kmeans algorithm to the age distribution of the population under study produced a fiveclass solution (Table 3) almost identical to that described previously, explaining ∼94% of the total variability of AGE.
The cardiac age (Table 3), as observed previously (5), decouples from chronological age at the level ofcluster D, reaching a plateau for the older population (cluster E), such that, an average chronological age of ∼81 yr scores a cardiac age of 64.
DISCUSSION
Besides the confirmation of the possibility to compute a cardiac age estimator for clinical and screening purposes and the existence of a saturation effect in cardiac age probably due to selection biases (i.e., the implicit selection linked to the fulfillment of the inclusion criteria at older ages), the main added value of this work is the elucidation of the dynamic character of the effect of aging on heart rate regulation and, more in general, in the proof of a “randomwalk” hypothesis for heartbeat dynamics. The emerging picture of heart rate regulation at the beattobeat level is that of a stochastic process. In this frame the agedependent increase of determinism has to be interpreted as an increase in the correlation time of the heartbeat (i.e., residence in the same quasistate) with a consequent increase in the general predictability of the signal. This picture comes from the consistency of the observed results with the Markov formalization introduced previously (7). It is important to stress that a Markov model is only a phenomenological description: what appears stochastic to us could be the result of complex adaptations to the internal state; moreover, any natural process has inertia, so two nearby states are, on average, more similar than two that are far apart. In these terms, the effect of senescence on cardiac dynamics can simply be expressed in terms of decreased “sensitivity” of the system to the changes of internal state.
In this work, age was by far the mainorder parameter of the total variability, and, in contrast to the previous work, its effect on HRV was not dispersed over different components but was concentrated on a unique factor. This factor (PC1) measures the influence of senescence on HRV and appears as a global determinism factor (Table 2).
With increasing age, heartbeat dynamics become increasingly deterministic and predictable. Recently, Giuliani et al. (7) proposed a “quantumlike” regulation of heartbeat with the heart oscillating between a few discrete states corresponding to different modal frequencies. The oscillation followed quantumlike dynamics with only two possible choices at each time step (beat):1) remaining in the last visited state and 2) shifting to an adjacent state, i.e., a simple firstorder Markov model. This kind of behavior, based on the theoretical framework of terminal dynamics [an innovative physical paradigm developed by Zak et al. (14), which models the dynamics of complex systems as short sequences of deterministic paths interspersed with purely stochastic singularities], was demonstrated for the first eigenvector of heartbeat sequences of rats and humans (7, 14, 17). The first eigenvector of the tachogram can be considered a filtered version of the original signal retaining the major dynamic features of the raw series (7). Given this general pattern of functioning, an increase in determinism corresponds to an increase in the relative probability for each single beat to remain in the same state (length class) as the previous one (increased probability for choice 1).
From an RQA point of view, this corresponds to an increase in the DET parameter (and, consequently, ENT and MAXLINE variables).
Terminal dynamics give us the physical rationale of how it is possible to exert such control. Randomness in such a paradigm is generated by the dynamics themselves and not from some error function. Furthermore, complex nonlinear interactions can be obtained at the singular points, whereas between these singular points the dynamics are constrained in a deterministic way. The advantage of this paradigm is that the organism is allowed to adjust the dynamics according to environmental influences at the singular points, whereas the deterministic character of a classical chaos paradigm would limit the adaptability of the system, given that, by definition, the deterministic character of the model makes the system strictly dependent on initial conditions. It may be argued that these objections can be countered by the existence of “control parameters,” which can retune the system but, given the tremendous amounts of noise in biological systems and the extreme sensitivity to initial conditions of chaotic dynamics, the energy expended to run adaptive controllers would be considerable (7, 14, 17).
In this work we demonstrated that aging acts on cardiac dynamics constraining the heart random walk on the last visited state, thus allowing us to generalize the Markovlike function to the heartbeat as a whole and not only to its first eigenvector.
This fact, in turn, implies that the optimal functioning of HRV regulation (younger ages) corresponds to a more stochastic character of HRV, paralleling a maximum adaptation power to an unpredictably changing environment. In older ages this adaptability is drastically reduced with a consequent increase in the deterministic character of the dynamics.
Our results parallel the findings of Mestivier and colleagues (10), who demonstrated an increase in heartbeat determinism in terms of MAXLINE in diabetic neuropathy.
More generally, our findings suggest that heart rate regulation is intrinsically stochastic and that the view of a basically deterministic (attractorlike) system with superimposed noise is basically incorrect. Besides the analysis of data and theoretical considerations (14, 17), our view is supported by the growing body of evidence (2, 4, 9) on the role of stochastic resonance effects in physiological systems.
From a practical point of view, our results may allow clinicians to compute in a very straightforward way (see ) the cardiac age of patients, which could be useful for general screening purposes. Screening situations in which the computed cardiac age could be useful include the epidemiologic and biomonitoring studies carried out to discover potentially neurotoxic effects of environmental hazards (1). In these cases, knowing the expected effect of age could greatly improve the power of the study. Obviously, the need to extend the method to the study of pathological situations remains crucial.
Acknowledgments
We thank Charles L. Webber and Joseph P. Zbilut for helpful discussions.
Footnotes

Address for reprint requests: A. Giuliani, Istituto Superiore di Sanità, TCE Lab, Viale Regina Elena 299, 00161 Rome, Italy.

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “advertisement” in accordance with 18 U.S.C. §1734 solely to indicate this fact.

↵1 RQA software, developed by C. L. Webber and J. P. Zbilut, can be downloaded from the following URL:http://homepages.luc.edu/∼cwebber/.
 Copyright © 1998 the American Physiological Society
Appendix
Computation of ESTAGE (cardiac age by means of RQA).
The subsequent steps indicate how to compute the cardiac age score. Obviously, all the registration parameters must be set as indicated inmethods. Otherwise it is best to generate another score following the procedure we outlined inmethods. The RQA software (RQD program, version 4.0) necessitates definition of the measuring parameters: in this case, lag = 1, emb = 15, norm = Euclid, dist = absolute, radius = 5, line = 5.
Given these premises, the first principal component is computed as
Having reckoned the PC1 value from the previous expression, the cardiac age results from