PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of ploscompComputational BiologyView this ArticleSubmit to PLoSGet E-mail AlertsContact UsPublic Library of Science (PLoS)
 
PLoS Comput Biol. 2010 December; 6(12): e1001033.
Published online 2010 December 16. doi:  10.1371/journal.pcbi.1001033
PMCID: PMC3002994

Rhythmic Dynamics and Synchronization via Dimensionality Reduction: Application to Human Gait

Jie Zhang, 1 , 2 , * Kai Zhang, 3 Jianfeng Feng, 1 , 4 and Michael Small 2
Jürgen Kurths, Editor

Abstract

Reliable characterization of locomotor dynamics of human walking is vital to understanding the neuromuscular control of human locomotion and disease diagnosis. However, the inherent oscillation and ubiquity of noise in such non-strictly periodic signals pose great challenges to current methodologies. To this end, we exploit the state-of-the-art technology in pattern recognition and, specifically, dimensionality reduction techniques, and propose to reconstruct and characterize the dynamics accurately on the cycle scale of the signal. This is achieved by deriving a low-dimensional representation of the cycles through global optimization, which effectively preserves the topology of the cycles that are embedded in a high-dimensional Euclidian space. Our approach demonstrates a clear advantage in capturing the intrinsic dynamics and probing the subtle synchronization patterns from uni/bivariate oscillatory signals over traditional methods. Application to human gait data for healthy subjects and diabetics reveals a significant difference in the dynamics of ankle movements and ankle-knee coordination, but not in knee movements. These results indicate that the impaired sensory feedback from the feet due to diabetes does not influence the knee movement in general, and that normal human walking is not critically dependent on the feedback from the peripheral nervous system.

Author Summary

Complex physiological rhythms arise from a large variety of biological systems that include natural pacemakers as well as feedback mechanisms, from the heartbeat to the rhythmic movement of human walking. Accurately extracting and characterizing the fluctuations underlying the biological rhythms is a fundamental problem which holds the key to understanding the mechanisms that govern the dynamics of biological systems. Usually such signals demonstrate certain oscillatory patterns, with each period displaying irregular fluctuation, or nontrivial dynamics, over time. This renders traditional spectral methods and nonlinear techniques less effective. We propose a novel approach to highlight the intrinsic fluctuations masked by the periodic component and noise through advanced dimension-reduction techniques, and apply it to human gait data from healthy subjects and diabetics. We find that this approach is capable of extracting the intrinsic dynamics and identifying the subtle synchronization pattern between knee and ankle. We find that although the two groups of individuals demonstrate remarkable differences in the dynamics of ankle movement and ankle-knee synchronization, the knee movement of both groups show similar dynamics. These results suggest that sensory feedback from a peripheral nerve system (like the feet) does not play an important role in regulating the motor control of human walking.

Introduction

Complex physiological rhythms and synchronization processes are ubiquitous in biological systems and are fundamental to life [1]. The human heartbeat [2], [3], walking [4], vocal cords vibration [5], blood pressure and respiration [6], white blood-cell count and tremor in patients [7], [8], epidemic dynamics [9] all demonstrate a stable, possibly nonlinear, oscillatory pattern along with highly irregular fluctuations from period to period. Such signals are variously known as semiperiodic, approximately periodic or pseudoperiodic time series. The fluctuation overlying the oscillatory pattern, or specifically, the cycle-to-cycle variability, arises from the combined effects from the changing environment, the nonlinear nature inherent to biological systems, and noise of various sources. It contains a wealth of information regarding the health or disease status of an individual subject. Usually, little or no a priori knowledge or models that govern the underlying system are available. Therefore accurately characterizing and quantifying such biological rhythms through data-driven approaches contributes significantly to our understanding of complex biological control systems [10] and have important applications in disease diagnosis.

Traditionally, rhythmic signals are fruitfully analyzed by linear methods like the Fourier transform and power spectrum analysis. However, physiological signals as outputs of complex biological systems are typically nonlinear and non-stationary, and can not be properly characterized by linear methods. A number of new techniques based on nonlinear dynamical system theory [11], [12] have also been intensively applied, like correlation dimension [13] and Lyapunov exponents [14]. Although the chaotic measures may provide new insights into the nonlinear nature of the system, they are severely hampered by the cyclic trend and noise in the system [15], [16]. Recent attempts include producing pseudoperiodic surrogate data [17], or performing a transformation from time domain to dual complex network domain [18][24]. Generally, there still lacks systematic and robust approaches to handle such oscillatory, possibly nonlinear time series. The Fourier analysis decomposes the signal into harmonics that span over the entire time-line, thus information about how each period, or cycle, changes over time has been averaged out. Similarly, nonlinear measures are also based on averaged properties of phase space attractor reconstructed from the data. This lack of discrimination among the individual cycles calls for more advanced signal processing techniques. Inspired by the recent advances in the field of dimensionality reduction [25], [26], we propose a novel and robust approach that can effectively capture the dynamics of cycle-to-cycle variation, which is especially suitable for analyzing approximately periodic, or semiperiodic data like human gait, ECG, respiration, and vowel data.

Human walking is a highly complex, rhythmic process which was found to exhibit long-range correlation and self-similarity, and has attracted sustained interest over the past decades [4], [27][33]. The fluctuations overlying the cyclic trend in human walking may reflect valuable information about the neuromuscular processes responsible for normal and pathological locomotor patterns. In particular, the stride interval (SI) (the duration of each gait cycle) has been intensively studied to quantify the physiological or pathological state associated with walking. For example, it has been found that the long range correlation properties are altered with aging and disease [10], [28]. The stride interval reflects the duration of each cycle, and a wealth of information contained in the waveform of the gait cycle is lost. Our approach utilizes the full waveform of each gait cycle, and is therefore expected to extract more information relevant to motor control of walking. We apply our method to the locomotion data collected from two groups of people — healthy subjects and neuropathic patients suffered from diabetes. We aim to find whether the extracted dynamical fluctuation of the knee and ankle movements as well as their synchronization pattern can vary between the healthy and diabetics group. Specifically, we want to find out whether the impaired sensory feedback due to diabetes can lead to different locomotion dynamics of the knee and ankle (and their synchronization pattern) compared with the healthy subjects.

Methods

Reconstructing Dynamics Underlying Cyclic Trend: Dimensionality Reduction

The general problem of dimension reduction has a long history. With advances in data collection, dimension reduction has reemerged as a prominent tool to unravel the high dimensional structure emerging in various disciplines. For example, it has been widely applied to gene and protein expression profiling for disease classification and prognostication [25], [26]. Generally, the large number of dimension reduction approaches can be categorized into linear methods, including the principle component analysis, multidimensional scaling [34], and nonlinear methods such as state-of-the-art Isomap [34], laplacian eigenmaps [35] and local linear embedding [34]. Usually biomedical data process nonlinear structures and that nonlinear dimensionality reduction methods might be more appropriate [26]. Here we use Laplacian Eigenmaps [35], which are based on spectral graph theory and projects the high-dimensional data into a low dimension so that two points nearby on the manifold are kept near to each other. We first illustrate with benchmark data from the chaotic Rössler system described by:

equation image
(1)

The time series from the An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e002.jpg-component (see Figure 1A), exhibits a strong periodic component along with irregular fluctuations, therefore it serves as an ideal example of an approximately periodic signal with nontrivial dynamics. Motivated by the fact that such data usually exhibit a highly redundant pattern in the form of repeated cycles, we can partition the time series into individual cycles An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e003.jpg (An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e004.jpg, which is the index of the cycle) at the peaks or troughs in the time series, see Figure 1A [36]. Each individual cycle can then be taken as a high dimensional vector An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e005.jpg (An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e006.jpg), whose dimension equals the number of points in that cycle. Our goal is to map these multiple, high-dimensional cycles to a set of new, low-dimensional (preferably one dimensional) representation, or embedding An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e007.jpg's, such that the proximity relations among An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e008.jpg's are maximally preserved in their low-dimensional counterparts An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e009.jpg's. In the case that each cycle An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e010.jpg is reduced to 1D (i.e., An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e011.jpg being a scalar), the derived An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e012.jpg constitute a new time series An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e013.jpg, which encodes the dynamics the original time series on the cycle scale.

Figure 1
Illustration of transforming a pseudo-periodic time series An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e014.jpg into a new series An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e015.jpg by reducing each cycle in An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e016.jpg to a point.

To achieve this, a weighted matrix An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e021.jpg is constructed, with each entry An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e022.jpg denoting the similarity between cycle An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e023.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e024.jpg, which can be chosen conveniently as the correlation coefficient An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e025.jpg (under the circumstance that An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e026.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e027.jpg differ in length, we shift the shorter vector along the longer one until An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e028.jpg maximizes). Then, the low-dimensional representation An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e029.jpg's can be cast as the solution of the following optimization problem, An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e030.jpg, which penalizes those mappings where nearby points An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e031.jpg's are relocated far apart in the space of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e032.jpg's. In case of univariate An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e033.jpg's, the objective can be written as An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e034.jpg, where An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e035.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e036.jpg is the graph Laplacian, and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e037.jpg is the diagonal degree matrix such that An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e038.jpg.

The above constrained minimization is solved by the generalized eigenvalue problem An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e039.jpg, where An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e040.jpg's (An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e041.jpg) are eigenvalues sorted in an ascending order, and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e042.jpg's are the corresponding eigenvectors. The minimum eigenvalue An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e043.jpg is zero, corresponding to an eigenvector (An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e044.jpg) whose entries are all An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e045.jpg. Therefore it is degenerate and the optimal solution is actually provided by An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e046.jpg, the eigenvector of the second smallest eigenvalue [35]. As is shown in Figure 1B, the eigenvector An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e047.jpg provides a cycle-scale representation of the original time series by reducing each full cycle to a single point. We use a general notation An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e048.jpg (An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e049.jpg for “cycle”, and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e050.jpg indicates the An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e051.jpgth cycle) for this simplified representation of the original time series. It is worthwhile to note that other dimension reduction schemes can also be adopted. For example, we can compute the Euclidian distance among cycles [36], and multidimensional scaling can be readily applied in this case to reduce the cycles to scalars. The An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e052.jpg derived from multidimensional scaling and the Laplacian Eigenmap in fact yield quite similar results, see Figure 2.

Figure 2
Return maps of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e053.jpg series and Poincaré section points An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e054.jpg.

The applicability of dimension reduction techniques is generally justifiable, considering the low correlation dimension of most real world pseudoperiodic data. For this kind of data, the trajectories of nearby cycles in phase space usually have similar orientations. Such redundancy can be effectively removed through dimension reduction, leaving only the useful degree of freedom. Finally, it is worthwhile to mention that for long time series with large number of cycles, the nyström method can be adopted to solve large scale spectral clustering problem [37], [38].

The Advantage of Characterizing Dynamics on the Cycle Scale

A popular method in nonlinear time series analysis is to reduce a continues flow to a series of discrete points, called a Poincaré section. The Poincaré section is the intersection of flow data in the state space with a hyperplane transversal to the flow. Thus each cycle in the data is simplified into a single point on the Poincaré section, which preserves many properties of periodic or pseudoperiodic orbits. Now we compare An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e059.jpg obtained by dimension reduction and the Poincaré section points An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e060.jpg obtained by collecting the local minimum points in the data. As can be seen Figure 2, the return plot (i.e., plot of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e061.jpg versus its previous values An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e062.jpg) of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e063.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e064.jpg show similar quadratic form. Further calculation of the chaotic measure such as correlation dimension indicates that they have the same dynamical origin.

One problem with the An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e065.jpg series is that it is highly susceptible to noise that is inevitable in biological data. To see this, we plot the return map for the An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e066.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e067.jpg series obtained from the noisy An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e068.jpg data, see Figure 3. We find that although both return maps display a clear quadratic form intrinsic to the chaotic Rössler system, the return map of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e069.jpg is more vulnerable to noise as the points are more dispersed than that of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e070.jpg. We use the variance An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e071.jpg of the least-square fit to the quadratic function (An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e072.jpg) to quantify the influence of noise, with An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e073.jpg, almost An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e074.jpg times larger than that of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e075.jpg (An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e076.jpg). Now we explain why An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e077.jpg is more robust to noise. First we take the “cycle” as the basic unit rather than a discrete point in the time series. Obviously the latter is more vulnerable to noise. Second, the acquisition of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e078.jpg is based on an optimization that preserves the proximity relation among all cycles simultaneously, while An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e079.jpg is obtained by treating each cycle independent of one another. Our approach utilizes the richer information of pairwise cycle correlation, therefore it not only excavates the inherent dynamics obscured by the cyclic trend, but also offers an extra robustness to noise due to the global nature of this method.

Figure 3
Return plots of Poincaré section points An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e080.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e081.jpg.

Detecting Degree of Synchronization from Bivariate Oscillatory Data

Another interesting phenomena associated with rhythmic process is synchronization between self-sustained oscillators, which plays an important role in understanding coordination or cooperation in biological systems [39], [40]. Several different types of synchronization have been observed, such as complete synchronization [41][43] generalized synchronization [44], and phase synchronization [45].

The evaluation of degree of synchronization from the outputs of coupled systems is of particular interest to study the interactions in biological systems. For example, the consistency of mutual nearest neighbors and the peakness of the phase difference distribution are used to characterize the dynamical interdependence [46] and phase synchronization [45], respectively. Here we are interested in the case where two processes are phase synchronized, but the synchronization strength is hard to estimate due to noise and non-phase-coherence, or is too subtle to be differentiated due to the mask of strong phase synchronization. For example, the knee and ankle move perfectly in phase during human walking. Under such circumstances the phase synchronization index will take on high values for both healthy subjects and diabetics and it may not probe the subtle difference in the degree of synchronization masked by the strong phase synchronization and noise.

To solve this problem, we propose to quantify the degree of synchronization between two noisy, phase-synchronized oscillatory processes on the cycle scale through the reduced representation An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e090.jpg of the original data, thereby minimizing the influence of phase synchronization. We illustrate with the An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e091.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e092.jpg components of the noisy Rössler system, which are perfectly in phase and therefore serve as ideal benchmark data. The two time series are first segmented into cycles according to the local minimums of either An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e093.jpg or An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e094.jpg time series, then we apply Laplacian eigenmap to both and obtain reduced representations, i.e., An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e095.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e096.jpg for An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e097.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e098.jpg time series. Finally we calculate the linear correlation coefficient An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e099.jpg between An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e100.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e101.jpg as an indicator of the degree of synchronization between An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e102.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e103.jpg data. We find that the extracted An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e104.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e105.jpg series can successfully reveal the synchronization pattern in presence of noise, which is demonstrated by an increasing trend in the corresponding scatter plot (Figure 4B). The Poincaré section points, however, are less informative of the synchronization degree due to the presence of noise, as An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e106.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e107.jpg do not demonstrate clear correlation (Figure 4A).

Figure 4
Correlation between An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e108.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e109.jpg component of noisy Rössler system as is revealed by Poincaré section points An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e110.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e111.jpg.

Results

Data Description

Now we apply the method proposed in previous section to human gait data collected from two groups: the healthy controls (CO) and neuropathic patients (NP, with significant diabetic neuropathy), each with 10 subjects [4]. The kinematic data were collected from a portable data-logger equipped on the subjects during continuous overground walking for 10 minutes (sampled at 66.7 Hz). Three electrogoniometers were placed on the approximate joint centers of the hip, knee, and ankle joints of the right leg to measure their sagittal plane movements. Here we consider the signals measured from knee and ankle joints movement.

Characterizing Human Locomotion Dynamics

Human locomotion is a highly complex, rhythmic process that involves control from subcortical locomotor brain regions and feedback from various peripheral sensors. Typically, the human gait time series (see Figure 5 for knee and ankle movement) exhibit a stable frequency while irregular stride-to-stride fluctuation. For biological signals with a strong periodic component, vital information regarding the pathology and pathophysiology of the subject is hidden in the cycle-to-cycle variation. Accurately extracting this fluctuation and characterizing its dynamics are expected to provide important insights into the underlying neuromuscular control of walking and yield greater diagnostic information. For example, the stride interval (SI), defined as the time duration of each gait cycle, has been widely used to study the human gait [27], [28], [31][33]. It was shown that SI series displays long-range correlation intrinsic to the healthy locomotor system.

Figure 5
Time series (upper panel) and the corresponding phase space reconstructions (lower panel) of knee and ankle locomotion from a healthy subject.

The SI series contains the information of the duration of each gait cycle. Another source of information consists in the waveforms of the gait cycles, which is not reflected in SI series [29]. As can be seen in Figure 6, the stride interval from the gait data seems to contain insufficient information to reflect the intrinsic dynamics due to digitization. Therefore it is natural to expect that the An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e118.jpg series, which is obtained by comparing the waveforms of each cycle (therefore preserves the full dynamical pattern within each cycle) will contain more information relevant to the locomotion dynamics. Meanwhile, An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e119.jpg successfully removes the periodic trend that may obscure the underlying dynamics, and has the same advantage as SI series. In the following we will demonstrate that An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e120.jpg obtained by Laplacian Eigenmap can reveal the dynamical fluctuation underlying the cyclic trend more effectively than the SI series, so that we can distinguish clearly between the healthy and pathological groups and make inference about the neuromuscular control, especially on the role of sensory feedback from the feet in regulating dynamics of human walking.

Figure 6
Time series of the extracted An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e121.jpg and the stride interval (SI) from ankle movement data of a healthy subject.

First we check the ankle movement (see Figure 5B), and use the Laplacian eigenmap to extract the fluctuation An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e123.jpg on cycle scale for the two groups. To quantitatively characterize the time evolution of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e124.jpg, we furthermore compute its power spectrum density (PSD), see the top row in Figure 7. We find that most CO subjects demonstrate broad band spectrums (i.e., An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e125.jpg noise) that scale as An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e126.jpg, with An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e127.jpg (see Figure 8A), indicating the presence of long range correlation (i.e., the strides separated by a large time span are still statistically correlated). In comparison, the power spectrum of the diabetic patients are mostly flat resembling white noise processes (An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e128.jpg), which means that the strides at different times are mostly uncorrelated. The values of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e129.jpg in the two groups are statistically different (An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e130.jpg). This difference, however, has not been found with either the stride interval (SI) series (see the middle row in Figure 7) or the raw data (see the bottom row in Figure 7).

Figure 7
Power spectrum density (PSD) for a typical healthy subject (the left panel) and a diabetes patient (the right panel).
Figure 8
Mean and standard deviation of the derived statistics for control (CO) and neuropathic (NP) groups.

The absence of long range correlation in the ankle kinematics of the NP group suggests the alteration of the locomotor pattern in the lower limbs of neuropathy patients. This is due to the loss of peripheral sensation in the lower limbs, which arises from the gradual dying back of nerves from the fingers and toes typical of diabetes. Despite the deterioration in peripheral nervous system, the knee movement of the NP group is still found to demonstrate a stable long range correlation indistinguishable from the CO group (An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e138.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e139.jpg, which are statitcally identical with An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e140.jpg), see Figure 8B. These results suggest that the impaired peripheral feedback caused by the dying nerves in the feet does not influence the upper-limb dynamics, which leds us to another fundamental problem in human walking, i.e., what is the role of sensory feedback in adjusting the global locomotor dynamics? To understand this, we examine further the degree of synchronization between the knee and ankle movement using dimension reduction.

Assessing Synchronization between Knee and Ankle Movement

Human walking involves the coordination of two major joints, i.e., the knee and the ankle, whose movements during continuous walking are obviously in phase due to the physical connection between them. However, we find that correlation between knee and ankle movement for the two groups can hardly be distinguished by the phase index of the signal due to the presence of strong phase synchronization. Also, noise tends to destroy the local structure in phase space and thus hampers the dynamical dependence measures [47]. To circumvent these difficulties, we propose to compare the dynamics of the two time series by using their Laplacian eigenmap An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e141.jpg's. Note that each time series can be segmented by either its own local maximums, or those of its partner series (shown in Figure 5). Therefore we will segment each time series twice and compute the averaged correlation coefficients An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e142.jpg's between An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e143.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e144.jpg for these two segmentation schemes.

Figure 9A shows the typical synchronization pattern between ankle and knee movement for healthy subjects. The scatter plot between An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e145.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e146.jpg demonstrates a significant increasing trend, indicating that the knee and ankle movements are highly synchronized. The correlation coefficient An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e147.jpg between An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e148.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e149.jpg for the healthy group takes on a high value: An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e150.jpg, see Figure 8C. For diabetics, however, there is little correlation between An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e151.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e152.jpg, as is manifested in randomly distributed points in Figure 9B. The correlation coefficient in this case is also low: An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e153.jpg, and the values of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e154.jpg is statistically different for the two groups (An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e155.jpg). All the results of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e156.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e157.jpg values for the two groups of individuals are summarized in Supplementary Table S1. Again, the discrimination between CO and NP groups cannot be achieved by SI series, which always exhibits a strong correlation between the two joints (Figure 9, lower panel), corresponding to high degree of phase synchronization. Finally, we point out that a more comprehensive description of synchronization can be achieved by examining more Laplacian eigenvectors. In the current case the single eigenvector An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e158.jpg already encodes the primary variability and is thus sufficient for the discriminative task.

Figure 9
Synchronization pattern between knee and ankle locomotion revealed by An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e159.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e160.jpg(upper panel) and stride interval series An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e161.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e162.jpg (lower panel).

The lack of significant synchronization between ankle and knee movements observed in diabetic patients suggests the “incoordination” between the ankle and knee movements. This may arise from the gradual deterioration of the nerves in foot and toes of the diabetics, which is unable to produce sufficient neural feedback for the lower limb to be coordinated with the upper limb. In section “Characterizing Human Locomotion Dynamics” we have found that both the ankle and knee movements for healthy subjects demonstrate long range correlation, while for patients, only the knee movement show long range correlation. This finding is consistent with the result obtained here, i.e., the ankle and knee movement are more synchronized for healthy people than for the diabetics. Finally it should be noted that our study is limited by the relatively small sample size (each group has 10 subjects). Therefore significance tests are performed to verify the differences observed between the two groups of individuals. The current conclusion will be further validated on a larger data base available in the future.

Discussion

Central Nerve Control over Peripheral Nerve System

A fundamental question concerning human walking is the origin of the long range correlation (or An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e166.jpg noise) found in human gait data, the mechanism of which are not exactly clear [28]. Generally, the locomotor system incorporates inputs from both the central nervous system such as the motor cortex and basal ganglia, and peripheral inputs and sensory feedbacks. Both these two kinds of inputs are suggested to be possible reasons for the presence of the long range correlation in normal human walking [4].

In section “Characterizing Human Locomotion Dynamics”, we found that although the locomotion dynamics of the ankle shows significant difference between the normal persons and the patients in terms of long range correlation, their knee movements demonstrate similar scaling properties. These results support the belief that the impaired peripheral feedback from the sensors in the feet of diabetics influences only the lower limb locomotion while not that of the knees. We therefore conclude that human walking is not critically dependent on the feedback from peripheral feedback of the lower-limb, and that the central nervous system is playing a major role in regulating locomotor dynamics. In fact it has been found that pathology in central nervous system, such as Huntington's disease, can result in a loss of long range correlation in the gait dynamics [28]. For diabetics, although the peripheral sensory feedback is weakened, their central nervous system is not damaged and still plays an important role in adjusting the locomotion dynamics. This is why their knee locomotion still demonstrate An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e167.jpg dynamics. Finally, it was pointed out that diabetics may still retain proximal somatosensory inputs, and visual or vestibular feedback information [33]. Further study need to be done to clarify the role of these factors in regulating the dynamics of human walking.

Possible Application to Other Areas

Our approach may be of great relevance, and is expected to provide more accurate and robust characterization and diagnostics to the complex oscillatory data observed in general biological and engineering fields. Reconstructing the dynamics on the cycle scale also brings new vitality to a number of other approaches which are otherwise not suitable for analyzing rhythmic data directly, such as detrended fluctuation analysis [48], recurrence plot [20][22], [49], entropy measures [50][52], surrogate data method [17], causality analysis[53][56], and so on. The inherent periodicity in time series can cover up the intrinsic dynamical fluctuation, thus extracting the dynamics on a cycle scale is crucial for subsequent analysis. For example, human vowel data and sunspot number variation are typical oscillatory time series in biomedical and astrophysical fields. We can easily extract the corresponding An external file that holds a picture, illustration, etc.
Object name is pcbi.1001033.e168.jpg series from such signals and feed them into DFA algorithm. By this we are able to tell if the original data has long range correlation on the cycle scale more reliably than other approaches.

Another striking example of physiological rhythms and their interaction is the complex, human cardiovascular system (CVS) [57], [58]. For example, the human heart is driving the blood circulation, in which one heartbeat corresponds to one cycle of blood pressure variation. Meanwhile, the heart rate is modulated by respiration through the so called Respiratory Sinus Arrhythmia (RSA). Analyzing the interaction among these biological rhythms from simultaneously measured ECG, blood pressure, and respiration force is crucial for understanding the cardiorespiratory control and disease diagnostics. Due to noise and non-phase-coherence of ECG, blood pressure and respiration signal, traditional measures may not be able to capture the subtle changes in the degree of synchronization, which could possibly be probed more accurately by our approach.

Finally, it is worthwhile to note that our approach can also be applied to time-course microarray data [59], as well as Functional Magnetic Resonance Imaging (FMRI) time series [61][63]. For example, genes related to cell-cycle control are always expressed periodically and subject to fluctuation, thus reliably characterizing the correlation or synchronization pattern among such genes using microarray data is expected to provide more insights into the corresponding regulatory network responsible for cell-cycle transcription and regulation. Similarly, FMRI experiments, which are unusually designed with periodic stimuli, always lead to periodic BOLD responses demonstrating large fluctuations [64][66]. Our approach is therefore expected to evaluate more reliably the functional connectivity among different regions in the brain for a deeper understanding of cerebral function.

Supporting Information

Table S1

β and ρ values for Control (CO) and Neuropathic (NP) groups.

(0.03 MB PDF)

Acknowledgments

We thank Dr. J. B. Dingwell for providing the human gait data.

Footnotes

The authors have declared that no competing interests exist.

JZ is supported by Hong Kong Polytechnic University(G-YX0N) and Fudan University. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

References

1. Glass L. Synchronization and rhythmic processes in physiology. Nature. 2001;410:277–284. [PubMed]
2. Bezerianos A, Bountis T, Papaioannou G, Polydoropoulos P. Nonlinear time series analysis of electrocardiograms. Chaos. 1995;5:95. [PubMed]
3. Babloyantz A, Destexhe A. Is the normal heart a periodic oscillator? Biol Cybern. 1988;58:203–211. [PubMed]
4. Dingwell J, Cusumano J. Nonlinear time series analysis of normal and pathological human walking. Chaos. 2000;10:848. [PubMed]
5. Little M, McSharry P, Moroz I, Roberts S. Testing the assumptions of linear prediction analysis in normal vowels. J Acoust Soc Am. 2006;119:549. [PubMed]
6. Small M, Judd K, Lowe M, Stick S. Is breathing in infants chaotic? Dimension estimates for respiratory patterns during quiet sleep. J Appl Physiol. 1999;86:359. [PubMed]
7. Haurie C, Dale D, Mackey M. Cyclical neutropenia and other periodic hematological disorders: a review of mechanisms and mathematical models. Blood. 1998;92:2629. [PubMed]
8. Edwards R, Beuter A. Using time domain characteristics to discriminate physiologic and parkinsonian tremors. J Clin Neurophysiol. 2000;17:87. [PubMed]
9. Olsen L, Schaffer W. Chaos versus noisy periodicity: alternative hypotheses for childhood epidemics. Science. 1990;249:499. [PubMed]
10. Goldberger A, Amaral L, Hausdorff J, Ivanov P, Peng C, et al. Fractal dynamics in physiology: alterations with disease and aging. P Natl Acad Sci Usa. 2002;99:2466. [PubMed]
11. Kantz H, Schreiber T. Cambridge University Press; 2004. Nonlinear time series analysis.369
12. Small M. World Scientific Pub Co Inc.; 2005. Applied nonlinear time series analysis: applications in physics, physiology and finance.245
13. Grassberger P, Procaccia I. Characterization of strange attractors. Phys Rev Lett. 1983;50:346–349.
14. Wolf A, Swift J, Swinney H, Vastano J. Determining Lyapunov exponents from a time series. Physica D. 1985;16:285–317.
15. Cazelles B. How predictable is chaos? Nature. 1992;355:25–26. [PubMed]
16. Abarbanel H, Brown R, Sidorowich J, Tsimring L. The analysis of observed chaotic data in physical systems. Rev Mod Phys. 1993;65:1331–1392.
17. Small M, Yu D, Harrison R. Surrogate test for pseudoperiodic time series data. Phys Rev Lett. 2001;87:188101.
18. Zhang J, Small M. Complex network from pseudoperiodic time series: Topology versus dynamics. Phys Rev Lett. 2006;96:238701. [PubMed]
19. Zhang J, Sun J, Luo X, Zhang K, Nakamura T, et al. Characterizing pseudoperiodic time series through the complex network approach. Physica D. 2008;237:2856–2865.
20. Donner R, Zou Y, Donges J, Marwan N, Kurths J. Recurrence networks¡aa novel paradigm for nonlinear time series analysis. New J Phys. 2010;12:033025.
21. Donner R, Zou Y, Donges J, Marwan N, Kurths J. Ambiguities in recurrence-based complex network representations of time series. Phys Rev E. 2010;81:015101. [PubMed]
22. Marwan N, Donges J, Zou Y, Donner R, Kurths J. Complex network approach for recurrence analysis of time series. Phys Lett A. 2009;373:4246–4254.
23. Yang Y, Yang H. Complex network-based time series analysis. Physica A. 2008;387:1381–1386.
24. Lacasa L, Luque B, Ballesteros F, Luque J, Nuño J. From time series to complex networks: The visibility graph. P Natl Acad Sci Usa. 2008;105:4972. [PubMed]
25. Tarca A, Carey V, Chen X, Romero R, Draghici S. Machine learning and its applications to biology. PLoS Comput Biol. 2007;3:e116. [PMC free article] [PubMed]
26. Lee G, Rodriguez C, Madabhushi A. Investigating the efficacy of nonlinear dimensionality reduction schemes in classifying gene and protein expression studies. IEEE ACM T Comput Bi. 2008:368–384. [PMC free article] [PubMed]
27. Hausdorff J, Peng C, Ladin Z, Wei J, Goldberger A. Is walking a random walk? Evidence for long-range correlations in stride interval of human gait. J Appl Physiol. 1995;78:349. [PubMed]
28. Hausdorff J, Mitchell S, Firtion R, Peng C, Cudkowicz M, et al. Altered fractal dynamics of gait: reduced stride-interval correlations with aging and Huntington's disease. J Appl Physiol. 1997;82:262. [PubMed]
29. Bartsch R, Plotnik M, Kantelhardt J, Havlin S, Giladi N, et al. Fluctuation and synchronization of gait intervals and gait force profiles distinguish stages of Parkinson's disease. Physica A. 2007;383:455–465. [PMC free article] [PubMed]
30. Dingwell J, John J, Cusumano J, Diedrichsen J. Do Humans Optimally Exploit Redundancy to Control Step Variability in Walking. PLoS Comput Biol. 2010;6:e1000856. [PMC free article] [PubMed]
31. Dingwell J, Kang H, Marin L. The effects of sensory loss and walking speed on the orbital dynamic stability of human walking. J Biomech. 2007;40:1723–1730. [PubMed]
32. Kang H, Dingwell J. Separating the effects of age and walking speed on gait variability. Gait Posture. 2008;27:572–577. [PubMed]
33. Gates D, Dingwell J. Peripheral neuropathy does not alter the fractal dynamics of stride intervals of gait. J Appl Physiol. 2007;102:965. [PMC free article] [PubMed]
34. Roweis S, Saul L. Nonlinear dimensionality reduction by locally linear embedding. Science. 2000;290:2323. [PubMed]
35. Belkin M, Niyogi P. Laplacian eigenmaps for dimensionality reduction and data representation. Neural Comput. 2003;15:1373–1396.
36. Zhang J, Luo X, Small M. Detecting chaos in pseudoperiodic time series without embedding. Phys Rev E. 2006;73:016216. [PubMed]
37. Zhang K, Kwok J. Density-weighted Nyström method for computing large kernel eigensystems. Neural Comput. 2009;21:121–146. [PubMed]
38. Zhang K, Kwok J. Clustered Nyström Method for Large Scale Manifold Learning and Dimension Reduction. IEEE T Neural Networ. 2010;21:1576–1587. [PubMed]
39. Pikovsky A, Rosenblum M, Kurths J. Cambridge University Press; 2001. Synchronization A Universal Concept in Nonlinear Sciences.432
40. Arenas A, Daz-Guilera A, Kurths J, Moreno Y, Zhou C. Synchronization in complex networks. Phys Rep. 2008;469:93–153.
41. Pecora L, Carroll T. Synchronization in chaotic systems. Phys Rev Lett. 1990;64:821–824. [PubMed]
42. Lin W, He Y. Complete synchronization of the noise-perturbed Chuas circuits. Chaos. 2005;15:023705. [PubMed]
43. Lin W, Ma H. Synchronization Between Adaptively Coupled Systems With Discrete and Distributed Time-Delays. IEEE T Automat Contr. 2010;55:819–830.
44. Kocarev L, Parlitz U. Generalized synchronization, predictability, and equivalence of unidirectionally coupled dynamical systems. Phys Rev Lett. 1996;76:1816–1819. [PubMed]
45. Rosenblum M, Pikovsky A, Kurths J. Phase synchronization of chaotic oscillators. Phys Rev Lett. 1996;76:1804–1807. [PubMed]
46. Schiff S, So P, Chang T, Burke R, Sauer T. Detecting dynamical interdependence and generalized synchrony through mutual prediction in a neural ensemble. Phys Rev E. 1996;54:6708–6724. [PubMed]
47. Quiroga R, Arnhold J, Grassberger P. Learning driver-response relationships from synchronization patterns. Phys Rev E. 2000;61:5142–5148. [PubMed]
48. Peng C, Havlin S, Stanley H, Goldberger A. Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series. Chaos. 1995;5:82. [PubMed]
49. Marwan N, Carmen Romano M, Thiel M, Kurths J. Recurrence plots for the analysis of complex systems. Phys Rep. 2007;438:237–329.
50. Amigó J, Zambrano S, Sanjuán M. Permutation complexity of spatiotemporal dynamics. Europhys Lett. 2010;90:10007.
51. Amigó J, Zambrano S, Sanjuán M. True and false forbidden patterns in deterministic and random dynamics. Europhys Lett. 2007;79:50001.
52. Amigó J, Zambrano S, Sanjuán M. Combinatorial detection of determinism in noisy time series. Europhys Lett. 2008;83:60005.
53. Brovelli A, Ding M, Ledberg A, Chen Y, Nakamura R, et al. Beta oscillations in a large-scale sensorimotor cortical network: directional influences revealed by Granger causality. P Natl Acad Sci Usa. 2004;101:9849. [PubMed]
54. Kamiński M, Ding M, Truccolo W, Bressler S. Evaluating causal relations in neural systems: Granger causality, directed transfer function and statistical assessment of significance. Biol Cybern. 2001;85:145–157. [PubMed]
55. Guo S, Wu J, Ding M, Feng J. Uncovering interactions in the frequency domain. Plos Comput Biol. 2008;4:e1000087. [PMC free article] [PubMed]
56. Ge T, Kendrick K, Feng J. A Novel Extended Granger Causal Model Approach Demonstrates Brain Hemispheric Differences during Face Recognition Learning. Plos Comput Biol. 2009;5:172–181. [PMC free article] [PubMed]
57. Chon K, Mullen T, Cohen R. A dual-input nonlinear system analysis of autonomic modulation of heart rate. IEEE T Bio-Med Eng. 1996;43:530–544. [PubMed]
58. Luchinsky D, Millonas M, Smelyanskiy V, Pershakova A, Stefanovska A, et al. Nonlinear statistical modeling and model discovery for cardiorespiratory data. Phys Rev E. 2005;72:21905. [PMC free article] [PubMed]
59. Rustici G, Mata J, Kivinen K, Lió P, Penkett C, et al. Periodic gene expression program of the fission yeast cell cycle. Nat Genet. 2004;36:809–817. [PubMed]
60. Wichert S, Fokianos K, Strimmer K. Identifying periodically expressed transcripts in microarray time series data. Bioinformatics. 2004;20:5. [PubMed]
61. Friston K, Holmes A, Poline J, Grasby P, Williams S, et al. Analysis of fMRI time-series revisited. Neuroimage. 1995;2:45–53. [PubMed]
62. Worsley K, Friston K. Analysis of fMRI time-series revisited¡aagain. Neuroimage. 1995;2:173–181. [PubMed]
63. Friston K, Harrison L, Penny W. Dynamic causal modelling. Neuroimage. 2003;19:1273–1302. [PubMed]
64. Bullmore E, Brammer M, Williams S, Rabe-Hesketh S, Janot N, et al. Statistical methods of estimation and inference for functional MR image analysis. Magn Reson Med. 1996;35:261–277. [PubMed]
65. Marchini J, Ripley B. A new statistical approach to detecting significant activation in functional MRI. Neuroimage. 2000;12:366–380. [PubMed]
66. Guy C, ffytche D, Brovelli A, Chumillas J. fMRI and EEG responses to periodic visual stimulation. Neuroimage. 1999;10:125–148. [PubMed]

Articles from PLoS Computational Biology are provided here courtesy of Public Library of Science