Non-Gaussian center-of-pressure velocity distribution during quiet stance

In daily life, we perform several tasks that demand a constant maintenance and control of postural balance (e.g., walking, dancing, carrying shopping bags, and so on). From a biological point of view, the quiet upright stance requires of the body the control of the musculoskeletal system. Such control is possible when the central nervous system processes visual, vestibular and somatosensory information, and uses it to select and control the appropriate postural response [1–5].

The study of postural sway by using center-of-pressure (COP)-based measures has been useful to characterize the dynamics of the postural control system associated with maintaining balance during quiet standing. The dynamics of COP position, which is frequently related to a random walk, has been investigated using a variety of techniques [6–16]. In particular, some studies have focused COP velocities instead of trajectories. It has been argued that COP velocity is most accurate than position or acceleration to modulate human quiet stance [17–19]. Several characteristics of COP velocity, including temporal correlations, have been explored [19–25]. However, investigations related to velocity distributions are scarce. For instance, there are statistical significance tests indicating that COP velocity distributions are non-Gaussian [26,27], but more details about the shape of these distributions have not been studied. The present study is mainly directed to investigate COP velocity distributions.

In this work, sixteen healthy subjects participated in the experiment (eight females and eight males between twenty and twenty-eight years old). After the ethics committee of the Universidade Estadual de Maringá approved the study, all participants gave their informed consent prior to their participation. All participants received oral and written instructions about the experimental protocol. They signed a commitment consenting before performing the tasks of the experiment.

The volunteers stood barefoot and were asked to maintain quiet upright posture on the force platform at a distance of about 1.5 meters from a fixed point in the horizontal plane with their legs oriented in parallel at hip width apart. Each subject was submitted to this task for 60 s, and the same was repeated 10 times in an alternating fashion, i.e., participants took turns among themselves in order to alleviate signs of fatigue. Thus, a total of N = 160 time series for this experiment was obtained. Data were collected using a force platform (EMG SYSTEM do Brasil) to record the ground reaction forces with an acquisition rate of 100 Hz.

The collected COP trajectory time series were filtered by a low-pass filter with a cut-off frequency of 20 Hz. This procedure eliminates high-frequency fluctuations. We then considered a sampling frequency of 50 Hz. The velocity series were obtained by differentiating the position series. The velocity time series were visually examined and it was decided that the first 1.5 second of each series would be removed to avoid anomalies related to the destabilization of the subject at the beginning of the experiment.

A representative COP trajectory is displayed in fig. 1(A) for a subject. Although COP presents a two-dimensional behavior, we chose to emphasize the investigation on its velocity in the anteroposterior direction (forward and backward body movement) as each individual has a higher propensity to imbalance. By using Y_t and $Y_{t+1}$ to represent two consecutive terms of the COP position series in the anteroposterior direction, one has $\Delta Y_{t}=Y_{t+1} -Y_{t}$ as an element of the series of successive increments. The increment $\Delta Y_{t}$ will be referred to as COP velocity ($\Delta Y_{t}$ is proportional to the velocity). Figure 1(B) shows a typical position series and its corresponding velocity series for the anteroposterior direction. Moreover, standardized increments (standardized velocity) are written as $\Delta y_{t}=(\Delta Y_{t}-\mu)/\sigma$, where μ is the mean value and σ is the standard deviation of $\Delta Y_{t}$. Points above six standard deviations were ignored in all series because it was understood that these data include involuntary stimulation of the subjects.

Zoom In Zoom Out Reset image size

The asymmetrical character and "peakedness" for each COP velocity series were studied. It was also adopted another way to quantify the probability distribution of our total data set which is given by the moments of the distribution. Furthermore, in the direction of suggesting candidates to model COP velocity distributions, two families of non-Gaussians with fat tails related to generalized Fokker-Planck equations were considered.

To obtain the following results, the procedure taken was to group the ten velocity series of each subject in a single data set. The standard deviation, σ, of $\Delta Y_{t}$ for each subject depends a lot of the subject characteristics (fig. 2(A)), a result consistent with the literature [28]. The asymmetry measure of all COP velocity distributions does not present considerable deviations when compared to a symmetrical distribution as the normal one, i.e., each skewness is approximately equal to 0 (fig. 2(B)). The skewness measure employed here is $\langle(\Delta y_{t})^{3}\rangle$, where 〈⋯〉 denotes the average over the standardized COP velocities of a subject. Nevertheless, all measures of the "peakedness" of the probability distribution indicate a non-Gaussian behavior since the values of the kurtosis are significantly larger than 3 (fig. 2(C)). Here, the kurtosis measure is given by $\langle(\Delta y_{t})^{4}\rangle$.

Zoom In Zoom Out Reset image size

In order to overcome the variability of the standard deviations, as in the cases of skewness and kurtosis analysis, and to try to identify more patterns related to COP velocities, the following analysis is based on the standardized COP velocity $\Delta y_{t}$. In this direction, consider the moments $\mu_{k}$ of the distribution of standardized COP velocity $\Delta y_{t}$, with $\mu_{k}\equiv\langle|\Delta y_{t} |^{k}\rangle$, where all the series were used as a single normalized data set. Figure 2(D) shows that the moments of the velocity distribution deviate from the Gaussian case as k increases. This result extends the kurtosis analysis for the global data because $\mu_{4}$ is the kurtosis.

In addition to the above analysis, tests of statistical significance were performed for data of each subject. The null hypothesis that the data is distributed according to the Gaussian distribution is rejected at the 5 percent level based on the Anderson-Darling and Cramér-von Mises tests. This analysis is consistent with previous studies where it was also indicated that COP velocity distributions are significantly different from Gaussians [26,27].

Next, the COP velocity distribution is explicitly calculated, i.e., $P(\Delta y_{t})$ is obtained by using all N = 160 series as a single normalized data set. From fig. 3, one verifies that $P(\Delta y_{t})$ remembers the normal distribution for $|\Delta y_{t}|\lesssim3$. In contrast, for $|\Delta y_{t}|\gtrsim3$, a deviation from a Gaussian is identified. A direct practical implication of this finding is that large velocities occur with larger probability than expected for a Gaussian process.

Zoom In Zoom Out Reset image size

Since COP movements can be viewed as random walks, candidates for COP velocity distributions could be related to Fokker-Planck equations. In the usual case of one-dimensional free-particle motion, the Fokker-Planck equation for the velocity v is given by [29]

$$\begin{equation} \frac{\partial P}{\partial t}=\frac{\partial}{\partial v}\left(D\frac{\partial P} {\partial v}\right)+\gamma\frac{\partial}{\partial v}(vP), \end{equation} \tag{ 1 }$$

where the noise strength D is constant, the friction coefficient γ is also constant and $P(v,t)$ is the probability velocity distribution at time t. In the stationary case $(\partial P/\partial t=0)$, with $P(v)\rightarrow0$ for $|v| \rightarrow\infty$, this equation reduces to

$$\begin{equation} D\frac{\text{d}P}{\text{d}v}+\gamma vP=0. \end{equation} \tag{ 2 }$$

The solution of this last equation is the Maxwell velocity distribution, i.e., the Gaussian $P(v)=N'\exp(-\gamma v^{2}/2D)$, with $N'$ being the normalization constant. As the COP velocity distribution exhibits non-Gaussian behavior, candidates for COP velocity distributions (generalized Maxwell velocity distributions) could be investigated from stationary solutions of generalized Fokker-Planck equations. Two cases of generalized Fokker-Planck equations that lead to a significant improvement in the description of the COP velocity distribution, when compared with a Gaussian, are focused here.

As a first generalization, consider $D=D_{q}P^{1-q}$ with D_q constant. This leads to a non-linear Fokker-Planck equation [30,31] which recovers the usual one when q = 1. The corresponding solution of eq. (2), when replacing v by $\Delta y_{t}$, is

$$\begin{equation} P_{q}(\Delta y_{t})=N_{q}\left[1-(1-q)\,b\,\Delta y_{t}^{2}\right]^{\frac{1}{1-q}}, \end{equation} \tag{ 3 }$$

where N_q is a normalization constant and b is a function of q, D_q and γ. Note that the Gaussian distribution corresponds to $P_{q}(\Delta y_{t})$ in the limit $q\rightarrow1$ and $P_{q}(\Delta y_{t})$ has fat tails because $P_{q}(\Delta y_{t})\propto|\Delta y_{t}|^{2/(1-q)}$ for large $|\Delta y_{t}|$ and q > 1. In addition, $P_{q}(\Delta y_{t})\approx N_{q}(1-b\,\Delta y_{t}^{2})$ for $\Delta y_{t}^{2}\ll1/b$. Therefore, eq. (3) has two aspects of COP velocity: a central part resembling a Gaussian and long tails when q > 1.

Equation (3) has been employed in many contexts, from physics and engineering to economics and biology [32,33]. For instance, it was applied with success to model the non-Gaussian velocity distribution of Hydra viridissima cells in two kinds of cellular aggregates [34]. Because the connection of eq. (3) with the Tsallis generalized statistical mechanics, it is usually referred as q-Gaussian. This statistics has been a useful theoretical tool to discuss anomalous behaviors, in particular, when related to non-Gaussians [32,33]. For instance, the usual Maxwell distribution is replaced by eq. (3). In this scenario, many aspects related to the usual (Boltzmann-Gibbs) statistical mechanics have been generalized. An example is the use of the non-linear Fokker-Planck equations [30,31,35,36] replacing the usual ones.

As well as the non-linear Fokker-Planck equation was first introduced in connection with a non-linear anomalous diffusion equation (porous medium equation) [30], another generalized Fokker-Planck equation can be related to other anomalous diffusion. In this direction, the same functional form of Richardson anomalous diffusion equation [37] is employed here but replacing the distance r by the velocity v. This proposal corresponds to using $D=D_{s}|v|^{s-2}$ in eq. (1), with D_s and s being constants. In this case, the solution of eq. (2), a generalized Maxwell distribution, becomes

$$\begin{equation} P_{s}(\Delta y_{t})=N_{s} \exp(-a\, |\Delta y_{t}|^s), \end{equation} \tag{ 4 }$$

where v was replaced by $\Delta y_{t}$, N_s is a normalization constant and a is a function of s, D_s and γ. Equation (4) reduces to a Gaussian for s = 2 and becomes stretched for s < 2. For this reason, eq. (4) is referred to here as stretched Gaussian. Analogously to the q-Gaussian, the stretched Gaussian is common in many scenarios, such as diffusion in fractals [38], two-dimensional turbulence [39], predator-prey encounters [40], atom deposition into a porous substrate [41], distribution of sound amplitude in music [42], and ergodicity breaking in heterogeneous diffusion processes [43,44].

Since the standard deviation of $\Delta y_{t}$ is equal to one, standardized q-Gaussian and stretched Gaussian distributions must be employed. As consequence, one has $b=1/(5-3q)$, with $q<5/3$, in eq. (3) as well as $a= [s\,\Gamma(1+1/s)]^{-s/2}\,\Gamma(3/s)^{s/2}$ in eq. (4), with $\Gamma(x)$ being the Gamma function. By using the maximum likelihood method, the COP velocity distribution for our global data was approximated by a q-Gaussian and a stretched Gaussian. Figure 3 suggests that the use of q-Gaussians and stretched Gaussians improves substantially the adjustment of the COP velocities in comparison with a Gaussian.

The same procedure was repeated to analyze the velocity distribution for each subject, i.e., $P^{(i)}(\Delta y_{t})$, with $i=1,2,\dots,16$. In this case, as previously pointed out, the ten series for each volunteer were grouped into a series and each $P^{(i)}(\Delta y_{t})$ was approximated by standardized q-Gaussian and stretched Gaussian distributions. An illustration of an individual COP velocity distribution for a representative subject is presented in fig. 4(A). The results of the individual analysis are shown in fig. 4(B) ($1.165\leq q\leq1.309$ and $1.179\leq s\leq1.488$), indicating that the q parameter remains larger than 1 and the s parameter is smaller than 2. Moreover, these parameters are close to their global values obtained in the previous analysis, $q=1.244$ and $s=1.301$. Tests of statistical significance (Anderson-Darling and Cramér-von Mises tests) indicated that more than 70% (80%) of 160 anteroposterior COP velocity series are not rejected at 1% level for q-Gaussians (stretched Gaussians).

Zoom In Zoom Out Reset image size

A similar analysis was performed to the mediolateral direction (left-right movement as X_t in fig. 1(A)). One found $q=1.177$ and $s=1.443$ for the global velocity distribution. For the individual cases, $1.104\leq q\leq1.226$ and $1.296\leq s\leq1.619$ were obtained. The same tests of statistical significance (Anderson-Darling and Cramér-von Mises) pointed out that more than 50% (60%) of 160 velocity series are not rejected at 1% level for q-Gaussians (stretched Gaussians).

In addition to the analysis of COP velocity distributions, we present below results concerning COP position distributions. Following a similar procedure to the velocity series, the fitting parameters obtained for the global position distribution were q close to 1 when employing q-Gaussians and s close to 2 for stretched Gaussians for both directions, anteroposterior (AP) and mediolateral (ML) ones. Explicitly, these values were $q=1.01\ ({q=1.059})$ and $s=2.023\ (s=1.928)$ for AP (ML). Skewness in both directions was near to zero, $-0.02\ (-0.07)$, and kurtosis 3.02 (3.22) for AP (ML) direction. For each individual, the values of q, s, skewness and kurtosis remained close to the global values. For the tests of statistical significance (Anderson-Darling and Cramér-von Mises) less than 18% (22%) of 160 position series were not rejected at 1% level for q-Gaussians (stretched Gaussians) in both directions. These results show that the non-Gaussian behavior for the position series is subtle or almost imperceptible.

To conclude our discussion about COP velocity distribution, the main achievements are summarized and possible implications of our findings are raised. Consistently with the literature [28], fig. 2(A) indicated that standard deviations exhibit large variability from a subject to another and therefore cannot be used as a typical characteristic of all subjects. To investigate more deeply aspects unrelated to the standard deviation, a standardized COP velocity was employed here. In contrast to the investigation of standard deviations, the skewness analysis (fig. 2(B)) pointed out in the direction that COP velocity distributions are well approximated by symmetric ones. In addition to this common aspect for all subjects of our analysis, the kurtosis values indicated a robust fat-tail pattern for COP velocity distributions when compared with a Gaussian (fig. 2(C)) and reinforced by a moment study (fig. 2(D)). The non-Gaussian profile was also verified by tests to decide if the velocity data comes from a Gaussian distribution, consistently with previous studies [26,27]. In addition, for standardized velocities, very similar distributions were identified for all subjects. And, considering generalized Fokker-Planck equations for COP velocities, the non-Gaussians distributions with fat tails can be interpreted as generalized Maxwell velocity distributions. These distributions were well approximated by q-Gaussians and by stretched Gaussians (figs. 3 and 4).

The empirical results obtained here for COP velocity distributions, together with those about correlations [19–21,24], compose a useful guide to model COP dynamics. However, others studies are still required to obtain a broader scenario for COP velocity distributions. Future works could investigate the sensitivity of COP velocity distributions to changes in postural steadiness related to age, neurological and orthopedic diseases or other pathologies affecting the postural control system.

Concerning the models employed here, they represent a first view of the scenario of COP velocities. In the model related to q-Gaussians, the noise strength D is proportional to a negative power of the probability distribution, $P^{1-q}$ with q > 1, leading a non-linear stochastic process. For the case of stretched Gaussians, the noise strength is proportional to $|v|^{s-2}$ with s < 2, that is, a heterogeneous model. If one is interested in a model that accommodates q-Gaussians and stretched Gaussians in a unified way [45], $D\propto P^{1-q}|v|^{s-2}$ could be used. Another possibility would be to consider other models that have stretched Gaussians as solutions. This is the case of the fractional Fokker-Planck equation [46,47], which is equivalent to continuous time random walk models with scale-free waiting times. Such a model presents a degree of memory and could be useful in connection with some regulatory mechanism in the maintenance of the upright standing. In turn, this Fokker-Planck equation and heterogeneous models exhibit weak non-ergodic behavior [43,44,48]. Motivated by a possible connection between these models and COP dynamics, the time dependence of the mean square displacement, measured from single COP trajectories, could be investigated. This procedure would give direct experimental verification of an anomalous stochastic process and of a non-ergodic character. Furthermore, to employ this procedure is useful because probability distributions could not be the best way to identify anomalous stochastic processes. All these aspects illustrate possible perspectives for further analysis involving upright quiet stance data.

This research was supported by the Brazilian agencies Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES).