ArticlesStochastic two-delay differential model of delayed visual feedback effects on postural dynamics
Abstract
We report on experiments and modelling involving the ‘visuo-postural control loop’ in the upright stance. We experimentally manipulated an artificial delay to the visual feedback during standing, presented at delays ranging from 0 to 1 s in increments of 250 ms. Using stochastic delay differential equations, we explicitly modelled the centre-of-pressure (COP) and centre-of-mass (COM) dynamics with two independent delay terms for vision and proprioception. A novel ‘drifting fixed point’ hypothesis was used to describe the fluctuations of the COM with the COP being modelled as a faster, corrective process of the COM. The model was in good agreement with the data in terms of probability density functions, power spectral densities, short- and long-term correlations (Hurst exponents) as well the critical time between the two ranges.
1. Introduction
Delays are an important feature of any biological control system. The feed- forward propagation of neural activity from sensors to the central nervous system, as well as from the central nervous system to muscles, can involve delays of hundreds of milliseconds. Delays complicate all control tasks involving sensors and effectors in any feedback configuration. This is because the controller, even if distributed, typically has a delayed response to errors detected between the desired and current states of the system. The study of such control systems with one or many feedback loops, and their robustness to sources of random fluctuations, is receiving ever more attention (Stepan 2009).
The postural control system is a canonical example of a neurological control system. Understanding its function is particularly important in the clinical realm, for example in the context of ageing. In fact, the control of posture degrades as we get older, and an understanding of the components of this control system and how they change with age is of paramount importance. In addition to delays in visual processing, there are delays of the order of hundreds of milliseconds in the purely proprioceptive control of posture based on mediated feedback from the muscles that apply torque to the ankle and other joints of the lower limb. Visual cues are important for stabilizing what is essentially an inverted (multisegmented) ‘pendulum’, and this visual processing also involves large delays.
It is essential to be able to perturb the healthy postural control system to gain a better understanding of its function, and of its dysfunction under pathological conditions. While it is not readily possible to extend or shorten the proprioceptive delay, it is possible to extend the visual delay, as has been done in the two recent studies of posture control (Rougier 2004; van den Heuvel 2009). This same strategy has also been used to study finger position control (Glass et al. 1988). It is known that symmetric systems with an unstable equilibrium, such as the inverted pendulum, can exhibit many different behaviours as parameters are changed. This is because their dynamics are organized around a co-dimension two bifurcation known as a Takens–Bogdanov bifurcation (in its symmetric version; see Redmond et al. (2002)), which arises when a pitchfork bifurcation meets a Hopf bifurcation. That same work showed how chaos may develop as nonlinearities are changed, adding a potential deterministic source of variability in the associated control problems. The addition of a second delay, as we do here to account for visual feedback, brings us into an unexplored dynamical territory.
In this paper, we first report new experimental measurements of posture control in healthy humans in which the visual delay is artificially increased by as much as 1 s. The variable that was measured and extensively analysed was the centre-of-pressure (COP) recorded using a force plate. We computed the probability densities of the fluctuations in anterio-posterior (AP) deviations from an upright position, as well as their spectral properties and Hurst exponents across delay conditions. This goes beyond the works of van den Heuvel et al. (2009), which concentrated on how the low- and high-frequency components of the sway behaved as a function of delay. Additionally, we summarize a novel dynamical model of postural control mediated by proprioceptive and visual feedback, which is based on the inverted pendulum. Thus, the model incorporates two delays associated with the corresponding control loops: proprioception and vision. The model also includes noise, which appears to be essential to simulate the ‘random walk’ executed by the COP. More noise is added to the slowly and randomly wandering fluctuation of the equilibrium position of the model to replicate the ongoing corrections to the position as the body shifts position during postural tasks. By making these fluctuations state-dependent, a better qualitative agreement between the model and the data is found. Postural fluctuations are known to exhibit different short- and long-term correlations (Collins & De Luca 1994). Here, our model was thus used to explore the stability of the system as a function of its delays, especially in terms of the cut-off between regimes of short- and long-range correlations.
Numerous models of postural dynamics have been proposed over the past two decades. Chow & Collins (1995) used an elastic-pinned polymer, modelled with a stochastic partial differential equation. Their model accounted for the characteristic mean-squared displacement and correlation functions typically seen in postural sway; however, it lacked an explicit control mechanism, which is our primary interest here. Kiemel et al. (2002) used high-order linear and optimal control theory models to estimate postural states. Higher order models require many parameters and thus are very complex. These models make use of Kalman filters and, therefore, states are evaluated instantaneously, lacking any minimal processing delay. Moreover, Frank et al. (2000) adopted the so-called vector-integration-to-end-point model to describe the COP dynamics, which was originally proposed by Bullock & Grossberg (1988) to explain the emergence of typical properties of reaching such as a speed-accuracy trade-off and a bell-shaped velocity profile. They identified the COP as a to-be-controlled (or to-be-minimized) difference between the actual position and the target position (in quiet stance the origin) and implemented the distinction between short- and long-term correlations by letting the drift coefficient of its dynamics depend on its correlation function (which led to a superdiffusive short-term regime).
A simpler approach to modelling postural control was based on the dynamics of the inverted pendulum. The validity of the inverted pendulum as a model for the upright quiet stance has been studied by Gage et al. (2004) and Winter et al. (1998), showing the relationship between COP and the centre-of-mass (COM). Eurich & Milton (1996) and Milton et al. (2009) have extended this approach to include noise, i.e. they treated the system as a stochastic delay differential equation. In both cases, the control was delayed and discontinuous since it was dependent on a constant threshold of ‘displacement’. This produces coexisting oscillations around each threshold, which were not evident in the data, but the more realistic global fluctuations due to noise-induced switches had Hurst exponents that agreed with the experiment. In the same context, a simple continuous control model was proposed by Yao et al. (2001). More recently, Stepan (2009) has worked on delay effects in the human sensory system while balancing. Also, recent developments show stabilizing effects of noise (Cabrera 2005; Milton et al. 2008). The analysis of stochastic delay differential equations in the context of posture control may benefit from techniques such as small delay approximations (Guillouzic et al. 1999), or two-state approximations if and when bistability is present (Tsimring & Pikovsky 2001).
To summarize, we focus on a two-delay stochastic differential equation to understand how generic features of the postural sway are related to several system parameters. Apart from providing a proper testing ground for the potential impact of vision on posture control, this approach allows for probing more deeply the boundary between the short-term, persistent (non-corrective) fluctuations and the more long-term, anti-persistent (corrective) fluctuations, for example determining how it depends on delays in the control loops and on the noise sources. Section 2 describes the experimental methods (data acquisition and analysis). Section 3 presents the empirical findings as a function of delay conditions: time series, standard deviations, power spectra, probability density and Hurst exponent estimates. Finally, the model is developed and analysed in §4 and we conclude in §5.
2. Methods
(a) Experiment
Eight men and six women aged 25.6±3.0 yr with masses and heights of 72±11 kg and 173±9 cm, respectively, participated in the experiment approved by the University of Ottawa Ethics Review Board. They signed an informed consent prior to participation. Participants had no cognitive or physical impairments. A participant’s task was to stand on a force plate (Advanced Mechanical Technology, Inc.) with their feet at a comfortable separation as straight as possible with their hands at their side in an upright stance. For visual conditions they looked at a 48 cm (19 inch) display, which showed a black background and a red dot at the centre representing the fixed target. A smaller white dot represented their COP and was displayed in real-time or with additional visual delay under experimental control. The experiment consisted of six conditions: eyes closed, and eyes open with delayed visual feedback of 0, 250, 500, 750 and 1000 ms. There were five trials for each condition, with a total of 30 trials per subject. Each trial lasted 120 s.
All footwear other than that with elevated heels or ankle constrictions was allowed. Participants stood on the same area on the force plate for all trials. Markers were placed on the force plate as a reminder of foot positions from trial to trial. Initial foot positions were determined before the first trial by finding the area on the force plate where participants stood upright and maximized the overlap of their COP with the target with minimal effort. Previous studies by van den Heuvel et al. (2009) and Rougier (2004) used instead the mean position of the COP’s first second of the trial as the target. The problem with that approach is that 1 s may not be enough time for initial COP transients to settle and consequently a less ‘balanced’ target might be used. That is, our technique aimed to find a better centred subject-specific target.
(b) Data analysis
Force plate data were sampled at 1000 Hz for AP and medio-lateral (ML) components. We will consider only the AP data in this article. These data were down-sampled offline to 100 Hz for analysis in Matlab (The Mathworks, Natick, MA, USA). Afterwards, 4 s at the beginning and 1 s at the end were removed to discard transients for a total of 115 s. Probability density functions (PDFs) were computed on the first three consecutive 30 s blocks, each with its mean removed. For every one of the six conditions, the PDFs were averaged over the three blocks for each subject, then across trials and finally across all participants. Power spectral densities (PSDs) were computed similarly via the Welch periodogram function with a non-overlapping Hamming window. Standard deviations were computed with the same partitioning. Mean-squared displacements were computed in 10 s blocks and averaged over 10 blocks per trial. They were averaged over all trials per participant and then over participants. A similar treatment of mean-squared displacement was given by Collins & De Luca (1993). The first quantity estimated from the mean-squared displacement was the critical time τc, defined as the first time that the Hurst exponent crossed H=0.5. We note that this also defines the boundary between the aforementioned two regimes of persistent and anti-persistent behaviour. A linear regression was then applied to the log–log of the mean-squared displacement K(τ) versus lag τ on [0.1,τc] and [τc,10] s, where the slopes of the fit were divided by 2 to give Hs and Hl, respectively. The slopes were divided by 2 for fractional Brownian motion (fBm) 
(Mandelbrot & Van Ness 1968; Newell et al. 1997).
3. Results
For the sake of legibility, results for the data and the model are presented simultaneously in the figures below, even though the model will only be presented in §4. First, let us refer to figure 1, which shows the COP for (i) no vision, (ii) 0 ms delay, and (iii) 750 ms delay, where all the time series were mean-centred. The condition with no vision was noticeably less stationary than the conditions with visual feedback, which suggested that visual feedback acted as an additional stabilizing and, probably, a corrective mechanism. By and large, figure 2 confirms this observation by displaying the mean (across participants) standard deviation of the COP time series for all experimental conditions; error bars represent the standard deviation over all participants and trials. It is interesting to note that, within error, the standard deviation varied little across visual delays.
Figure 1. COP time series for (a) no vision, (b) 0 ms delay and (c) 750 ms delay. Grey line, data; black line, model.
Figure 2. Standard deviation of COP times series for all conditions. Black line, data; grey line, model.
To disentangle effects in different time scales, figure 3 presents the mean PSD, Srr(f), averaged over participants for the same conditions as in figure 1. There was more power at lower frequencies for the no vision condition, since there appeared to be a slow, drift-like process (see also figure 1). A feature seen across all conditions was the ‘elbow’ in the middle of the curve, which divided the spectrum into two 1/fα-type regimes. This is characteristic of a multi-fractal one (Mandelbrot & van Ness 1968). If we proceed by assuming that the no vision condition has more drift, as suggested by the time series and the PSD, then its PDF, Pr, should have a relatively larger variance. Figure 4 clearly demonstrates this fact, whereas the two PDFs for visual feedback conditions were approximately the same regardless of delays (only two are shown here). We then computed the mean-squared displacement K(τ), and, from it, the short- and long-term correlations quantified via the corresponding Hurst exponents, Hs and Hl, as well as the critical time, τc. We note that for ordinary Brownian motion (oBm) H=0.5 holds, since 
. Figure 5 shows a schematic representation of the calculation of these quantities. Not only did we see multi-fractal behaviour in the PSDs, but also in the mean-squared displacement. Therefore, it can be detected from the time and frequency domains. Figure 6 shows the average values of Hs, Hl and τc across all participants. In fact, the Hurst exponents agreed with previous studies; i.e. Hs is between 0.5 and 1 and the range of Hl is from 0 to 0.5 for all conditions. Collins & De Luca (1995) have shown a significant difference between eyes open and eyes closed conditions for Hl. Similarly, figure 6 reveals a clear difference in the long-term Hurst exponent for the no vision and visual feedback conditions. The critical times, however, were the quantities of most interest since they define the boundaries of distinct control regimes (Collins & De Luca 1993). Across all conditions, even without vision, the critical times are approximately 0.5 s. Therefore, it appears that the artificial visual delay was not a major contributor to the value of τc.
Figure 3. PSDs for (a) no vision, (b) 0 ms delay and (c) 750 ms delay. Grey line, data; black line, model.
Figure 4. PDFs for (a) no vision, (b) 0 ms delay and (c) 750 ms delay. Grey line, data; black line, model.
Figure 5. Mean-squared displacement for the no vision condition. Dotted line, data; bold line, model.
Figure 6. (a) Short- Hs and (b) long-range Hurst exponents Hl and (c) critical times. Black line, data; grey line, model.
4. Model
(a) Inverted pendulum dynamics
In the models of Yao et al. (2001) and Eurich & Milton (1996) for the postural sway along the AP direction, the body is likened to an overdamped, inverted pendulum, which is considered as the ‘plant’ (see also below and Werness & Anderson (1984)). From figure 7, the inverted pendulum is a rigid mass-less rod connected to a point of mass m, and L is the distance from the origin to the COM. From this, one can obtain the moment of inertia I=mL2. There also exists a viscous damping force at the origin (not shown) proportional to the angular velocity. In combination, this leads to an equation for the rotational equilibrium in the yz-plane, where θ is the angle in that plane measured from the vertical direction
, which yields
. Because of θ=y/L, equation (4.2) reduces to 
. Henceforth, we will replace y with r and ε=mgL/β and write
Figure 7. Inverted pendulum free-body diagram where x represents ML and y, AP.
(b) Noise
Postural perturbations for quiet standing are traditionally modelled with Gaussian white noise denoted as ξ(t), which is a wide sense stationary process. This means that the first and second moments remain constant over time. Mathematically, its mean is 〈ξ(t)〉=0 and the autocorrelation function is 
, whereas the PSD is Sξ(f)=1; the latter follows via the Wiener–Khintchine relation. To avoid confusion with other works, it should be noted that, in this paper, all Gaussian white noise processes have a variance equal to 1.
(c) Feedback
The previously described inverted pendulum dynamics is intrinsically unstable at its fixed point. However, when observed experimentally, quiet standing appears fairly stable, and suggests that there must be a stabilizing mechanism, which exists to counteract this phenomenon (Eurich & Milton 1996). In the framework of continuous stochastic delay equations (see also Yao et al. 2001), a self-stabilizing model should be of the form
), where all parameters are positive (for the remainder of the paper as well) and κ>ε as follows:
(Mackey & Nechaeva 1995). In fact, since such oscillations are not manifested in the time series or PSDs of the experimental data, it is likely that this condition is satisfied. Equation (4.6), however, cannot reproduce both eyes closed and eyes open conditions since the eyes closed condition is non-stationary owing to a slow drift process as seen from the time series and PSDs in figures 2 and 3a. To explain this drift process, we turn to Zatsiorsky & Duarte (2000), who make reference to a migrating reference point (rambling), which we shall refer to dynamically as a ‘drifting fixed point’. This drifting fixed point may be modelled as an independent Wiener process, which is also a position, by adding it to equation (4.6),
We further wished to include visual feedback to explain both eyes closed and eyes open data in a single model. Comparing figures 2–4, we noticed that the drift process seemed negligible for the eyes open case when compared with that for the eyes closed one. Hence, we proposed, as earlier, to model the visual feedback similarly to the proprioceptive feedback. The visual feedback, however, must be corrective, i.e. directed towards the origin. Including this as part of equation (4.9) to minimize drift results in a doubly delayed Langevin equation. This inclusion reduces power in the PSDs at low drift frequencies, and the variance of the PDFs decreases, as seen experimentally. Equation (4.9) then becomes
and 
the final form of the model is
, r*=0, but if, temporarily, 
, r*=s/(κ−ε). Hence, the fixed point of r drifts as a function of the slow fluctuations of s. When γ=0 (no visual feedback), s is just the Wiener process and, thus, non-stationary; otherwise, vision acts as a corrective mechanism that keeps the drifting fixed point stationary. As one would expect, with relatively small visual delays, compared with the large delays, the system is more corrective as seen by the lower long-range Hurst exponents (figure 6a,b), yet the short-term Hurst exponents do not vary significantly across all visual delay conditions as well as without vision.(d) Comparison of the model with the data
Figure 1 shows how the model replicates the experimental data in terms of the fast corrective r-dynamics and the s-dynamics for the no visual and delayed visual feedback cases with non-stationarity and stationarity, respectively. The standard deviations of the realizations of r in figure 2 show a similar pattern to the data as well. All PSDs of the model in figure 3 display the same two regimes similar to the data, and although the quantitative agreement is not strong, it is acceptable for this coarse level of modelling. Figure 4 shows that the model has the same PDF as the data for the no vision case and for the delayed visual feedback. In particular, the variance is similar except for a higher peak at the origin.
Further, as figures 5 and 6 illustrate, the addition of the delayed visual feedback term to the s-dynamics produces lower values of the anti-persistent Hurst exponents (Hl) with respect to the no vision condition as seen in the experimental data. Additionally, in figure 6 the short- and long-term Hurst exponents are all in the same range for the model and the data. However, this is not true for the critical times, where the model has larger values than the data. Nevertheless, the values of τc for the model are less than 1 s, which agrees with previous studies. Moreover, the critical times are flat across all visual delays, similar to the data.
(e) Range of parameters
When simulating the model, we used subject-specific parameters for ε, κ and τprop. We attempted to model each participant in the experiment individually: ε=mgL/β, where m denotes the participant’s mass, L represents the distance from the force plate to the COM, which is approximated by dividing the participant’s height by 2, and the damping coefficient β could be given by Winter et al. (1998). To fulfil κ>ε, we fixed κ=1.5ε and τprop=0.6τcrit. As stated earlier, τvis=200+{0,250,500,750,1000} ms and γ=0.2 s−1 for all participants. Finally, the noise coefficients were ar=1.8×10−3 and as=1.2×10−3 m s−1/2. The choice of these parameters was made simply to fit the model to the data in terms of all results shown earlier. The stability of the two-delay model was, however, far from being known in any detail. Therefore, we performed a numerical stability analysis on r(t) in terms of τprop and τvis, where solutions with a standard deviation greater than 1 m are unbounded. Otherwise, we attempted to detect a fixed point or a marginally stable oscillation (since the system is linear). To detect a marginally stable oscillation we computed the mean of the PSDs in the range of ωcrit/2π±0.1 Hz. If this value was greater than the mean of the PSDs in the range [0,ωcrit/2π−0.1] Hz, and greater than some threshold value of 10−4 m2 Hz−1, it was labelled as a marginally stable oscillation; otherwise it was a fixed point solution.
Figure 8 shows the resulting numerically determined stability diagram. We chose typical parameter values for ε=2 and κ=1.5ε=3 Hz, and the noise coefficients were ar=as=0 m s−1/2 since we were investigating the deterministic case. Solutions were computed for 1210 s and the first half of each solution is discarded to remove transients. Initial history functions were constant and lasted for 1 s, where r0=0.001 and s0=0 m. We then decided to investigate the stability of the two-delay subspace within the experimentally realistic range of τprop=[0,0.5] and τvis=[0,1.2] s. It is clear that the model was stable in our experimental regime, since τprop did not exceed approximately 0.35 s and the maximum value of τvis was 1.2 s. One may also notice that the marginally stable oscillatory regime acts as a small boundary region separating the stationary and unbounded solution regimes.
Figure 8. For all time series τvis=0.5 s and (a) τprop=0.2 s, (b) τprop=0.365 s and (c) τprop=0.450 s. (d) Stability diagrams: white, fixed point solution; black, marginally stable solution; grey, unbounded solution.
(f) Further analysis
(i) Dependence of the critical time on delays
Let us recall figure 6c, where we have the critical time τc in relation to the experimental visual delay. This figure shows that τc is almost constant; however, we used the model to gain a better understanding of τc in terms of τprop and τvis. We averaged the values of τc over 25 realizations. Figure 9 reveals that the critical time is predominantly a function of the proprioceptive delay, τprop.
Figure 9. Critical time τc as a function of proprioceptive delay τprop and visual delay τvis. White on this figure indicates that τc does not exist between 0.1 and 10 s. The same parameters are used here as in figure 8, but with noise.
(ii) Another interpretation of the relationship between COM and COP
We next compared the model’s outcome with the empirical time series. Shown in figure 10a is a single realization of s(t). Contrasting this to a typical realization of r(t) in figure 1a, one may notice that s(t) does not oscillate as fast as r(t). If comparing the respective PSDs in figures 10b and 3c, Srr(f) appears greater than Sss(f) between approximately 0.1 and 1 Hz. Hence, it seems that s(t) is a low-pass filtered version of r(t).
Figure 10. (a) A sample realization of s(t) for experimental visual delay of 750 ms and (b) its corresponding PSD.
If we assume that r(t) represents the COP and s(t) the COM, then this phenomenon has already been reported by Winter et al. (1998), Caron et al. (1997) and Zatsiorsky & Duarte (2000) from a biomechanical perspective. Yet, our model is defined in the time domain by equations (4.11) and (4.12). Thus, because of the time-delayed term in equation (4.12), one may interpret s(t) as an exponentially decaying memory process of r(t).
(iii) Alternative models
During the process of model identification, we have tried numerous models, including
5. Concluding remarks
We have presented new experimental data on posture control in healthy subjects in the presence of visually delayed feedback. The data reveal that the visual feedback reduces the standard deviation of the postural sway, as expected from Collins & De Luca (1995), but it is not sensitive to the visual delay in the range investigated. This contrasts with the low-frequency behaviour, seen in van den Heuvel et al. (2009). We have also shown how the critical time separating the persistent and anti-persistent behaviours is well defined in the data, and, furthermore, that it is not very sensitive to the visual delay. We have proposed a novel modelling approach to such two-feedback posture control using a system of stochastic delay differential equations with two delays and noise, as well as a drifting fixed point meant to represent the slower fluctuation of the COM. The model qualitatively reproduces the main characteristics of our empirical findings with physiologically motivated parameter choices. The model further shows how, generically, the combination of the delayed feedback and noise produces the proper diffusive short- and long-term fluctuations. Finally, it enables us to show how the critical time is mainly a function of proprioceptive delay.
Acknowledgements
We thank James Jun and Maarten van den Heuvel for the very informative discussions which aided us with the experimental apparatus.
Footnotes
One contribution of 13 to a Theme Issue ‘Delayed complex systems’.
