Postural instability via a loss of intermittent control in elderly and patients with Parkinson’s disease: A model-based and data-driven approach

Symptoms associated with postural instability are common in patients with Parkinson’s disease, and they become more prominent with disease progression, which causes falls and diminished quality of life.⁴ Contrary to popular belief, neuromechanical mechanisms of human upright posture during quiet stance are still subject to debate. Consequently, pathophysiology that causes postural instability in Parkinson’s disease is still poorly understood.

The intermittent time-delayed feedback control, referred to as the intermittent control in this study, is one of the novel hypothetical neural strategies for stabilizing upright posture during quiet stance.^5–7 Intermittency and discontinuity in the action of nonlinear postural controllers have been well established in the literature.^5–18 Although recent studies are beginning to associate intermittency with a healthy control strategy and loss of intermittency with deterioration in postural function,^19,20 it is still controversial whether such time-discontinuous, nonlinear controllers are more physiologically plausible compared to the traditional time-continuous linear controllers, i.e., the stiffness and impedance controllers as models of healthy strategy for postural stabilization.

Among many of the traditional models,^21–23 the majority have simply adopted a linear proportional (P) and derivative (D) feedback controller with a large signal transmission time delay and additive colored noise, i.e., the time-continuous delayed PD feedback control model.^24–28 On the other hand, the intermittent control model that we adopt in this study is a hybrid dynamical system that switches between two unstable subsystems in a state-dependent manner, driven by additive white noise.^6,7 Unlike other models of postural control with seemingly similar intermittent characteristics,^8–11,15 our intermittent control model exploits two types of instability. One is generated by the off-subsystem (open-loop control system), which exhibits unstable dynamics of purely mechanical inverted pendulum-like human body in the absence of the active feedback control, due to the fact that intrinsic (i.e., passive) ankle stiffness is not sufficient for stabilizing the upright posture,^29,30 whereby the upright posture is characterized by a saddle-type unstable equilibrium with stable and unstable manifolds in the phase space. The other is generated by the on-subsystem (closed-loop control system), which is exactly the same as the traditional continuous PD feedback control model, if it were adopted persistently without switching but with $P$ and $D$ gain values that exhibit delay-induced unstable oscillations. Note that typical $P$ and $D$ gains of the on-subsystem used for the intermittent control model are much smaller than those used for the continuous model, using which the on-subsystem is unstable for the former, and it is stable for the latter. In the intermittent control model, switching between those two unstable dynamics in an appropriate state-dependent manner makes the overall dynamics stable, where the switching function is implemented in such a way that the feedback controller is switched off when the state vector of the inverted pendulum is near the stable manifold of the saddle, and it is switched on otherwise. To avoid misunderstanding of our model, it is worthwhile to note that the switching between two unstable dynamics in the intermittent control model provides a stabilization mechanism qualitatively different from the well-known model of “open-loop and closed-loop control of posture”^8,31 that exploits another type of intermittency. In the case of the “open-loop and closed-loop control of posture”^8,31 and also in many other models with intermittency,^9–11,15 stability of the upright posture is ensured by the stable dynamics of closed-loop control mechanisms, in which the rate of convergence for every closed-loop interval is (statistically) greater than the rate of divergence for the corresponding open-loop interval. Contrastingly, in our intermittent control model,⁶ both closed- and open-loop dynamics are unstable, but the rate of convergence along the stable manifold of the open-loop dynamics until termination of the open-loop control can be (statistically) greater than the rate of divergence for the corresponding delay-induced oscillation in the closed-loop dynamics, leading to overall stability of the hybrid system. The important role played by the stable manifold in the intermittent control model is conceptually similar to that in the OGY chaos control that also exploits a stable manifold of a saddle point to make chaotic dynamics periodic.³² The difference between them is a way to exploit the stable manifold. A chaotic meandering behavior and a delay-induced unstable oscillation are exploited by the system to make the state point close to the stable manifold in the OGY and the intermittent control models, respectively.

The first objective of this study is to demonstrate that the intermittent control model is more physiologically plausible compared to the continuous PD feedback control model, as a healthy strategy for postural stabilization. This will be achieved using a Bayesian parameter inference to identify model parameters using which the model can statistically and optimally reproduce a diverse set of spontaneous postural sway data. Although a similar attempt has been made previously based on the continuous PD feedback control model,²⁷ the range of optimal parameter exploration of the current study is much wider compared to the previous study. This is because the continuous control model can be considered as part of the intermittent control model by the switching function of the intermittency parameterized by some parameters that determine state-dependent switching between the off- and the on-subsystems. That is, a region in the state space (referred to as the off-region) assigned for selecting the off-subsystem can be varied, and if we employ such parameter values that make the off-region null, the intermittent control model can never be switched to the off-subsystem, i.e., the intermittent control model becomes equivalent to the continuous control model. Therefore, our parameter exploration for the intermittent control model can examine whether behaviors of the intermittent control model are optimized by the parameter values that make the off-region null (or small) for representing the continuous control model or those to ensure a large area of the off-region for the intermittent control model.

The second objective is our main contribution, in which we examine a novel hypothesis that explores association between the loss of intermittency in postural control and deterioration of postural function.^19,20 Specifically, we will apply the Bayesian parameter inference on the intermittent control model using published data of spontaneous sway from healthy elderly subjects and patient with Parkinson’s disease and compare them to the data obtained from healthy young subjects. Previously, Maurer and colleagues have performed similar explorations based on the continuous PD feedback model.^27,33 They found that increases in stiffness and damping, i.e., $P$ and $D$ gains of the continuous controller and a large increase in the noise level with aging and postural symptoms in patients with Parkinson’s disease could account for a variety of sway measures reported in the literature for elderly and patients.^34,35 For example, increase in the sway amplitude due to therapeutic intervention³⁶ can be explained by the PD control model with smaller $P$ and $D$ gains compared to the corresponding parameters for patients with severe postural instability. This is contrary to the past view that regards large sway amplitude as a sign of decreased postural stability.³⁷ Moderately large sway amplitude, as a sign of preferable or improved postural function, might be counter-intuitive. However, it has been reported occasionally.^34,38,39 From a theoretical viewpoint, we are interested in alternative mechanisms (model-based interpretations) that can be obtained from the intermittent control model, instead of the continuous control model, for such severity-dependent variations in the measures of postural sway. Namely, we are interested in whether the loss of intermittency from the intermittent control model can be a major cause of postural instability due to aging and progression of Parkinson’s disease.

In Sec. II, we introduce the intermittent control model, and dynamics of the model are quantified by several summary measures that are used for the Bayesian parameter inference. Moreover, stability analysis of the continuous control model is performed, which is helpful for understanding the intermittent control model. Section III describes methodologies of the Bayesian parameter inference. Results of the parameter inference are shown in Sec. IV, while we discuss our results in Sec. V.

Upright posture during quiet standing is often modeled by a single inverted pendulum in the sagittal plane.⁴⁰ Distal end of the pendulum is fixed in the plane by a pin joint as the ankle joint. The linearized motion equation for small angles and angular velocities can be formulated as

where $φ$ is the tilt angle, $I$ is the moment of inertia of the upright body around the ankle, $g$ is the gravity acceleration, $m$ is the body mass, $h$ is the distance from the ankle to the center of mass (CoM) of the pendulum, and $T$ is the total ankle torque. $mghφ$ represents the linearized gravitational toppling torque. For small values of $φ$ and $φ˙$⁠, $T$ may be modeled as

where the first two terms on the right-hand side represent passive feedback torques generated by the intrinsic mechanical impedance of the ankle joint (⁠$K$ and $B$ are the passive stiffness and viscosity, respectively) with no time delay. The third term $Tact$ represents the active neural feedback torque determined by the central nervous system as a function of delay-affected tilt angle $φΔ≡φ(t−Δ)$ and angular velocity $ωΔ≡ω(t−Δ)≡φ˙(t−Δ)$⁠, with $Δ$ being the feedback delay time due to the neural signal transmission. The last term $Tnoise$ is an endogenous noise modeled by the standard Gaussian white noise $ξ$ with standard deviation $σ$⁠. Using $ω≡φ˙$⁠, the state space representation of Eq. (1) can be rewritten as

where $k≡(mgh−K)/I$ and $b≡B/I$⁠. According to the experimental evaluations of the passive ankle visco-elasticity,^29,30 we set as $K=0.8mghN\,m/rad$ and $B=4.0N\,m\,s/rad$ and assume that they are constant throughout the study. Thus, $mgh−K>0$⁠, and the upright state $(φ,ω)=(0,0)$ is an unstable equilibrium of saddle type,⁶ with stable and unstable manifolds when no active control is provided, i.e., when $Tact=0$⁠.

The traditional continuous control and the intermittent control models employ the same PD controller but in a strongly different manner. The continuous control model uses the PD controller in a conventional manner,²⁷ i.e., continuously in time, which is why we call it the continuous controller. It generates the active feedback torque in Eq. (2), denoted by $Tactcont$⁠, formulated as

where $P$ and $D$ are the gains of proportional and derivative controllers, respectively. Stability of the continuous control model for a given delay $Δ$ is determined by the characteristic equation of Eq. (1),

when $K−mgh<0$⁠. Stability boundary in the $P$–$D$ parameter plane can be obtained using Stepan’s formula⁴¹ as a function of $Ω$ for $s=jΩ$ with $j2=−1$⁠,

For $Ω=0$⁠, Eq. (6) becomes

meaning that $P=mgh−K$ is part of the stability boundary. Depicting a curve on $P$–$D$ plane by varying $Ω$⁠, together with $P=mgh−K$⁠, we have a stability region of the continuous control model in the $P$–$D$ plane as Fig. 1 for several delays $Δ$⁠, showing that the stability region diminishes as $Δ$ increases.

In the intermittent control model, the PD controller is switched on and off intermittently, depending on the delayed state of the pendulum⁶ (Fig. 2). The corresponding active feedback torque in Eq. (2), denoted by $Tactint$⁠, is formulated as

where $Son$ and $Soff$ define two regions in the $φ$–$ω$ plane, by which the state-dependent activation and inactivation of the active torque are determined. In Fig. 2(b), $Soff$ is represented by gray regions and $Son$ by white regions. The parameters $P$ and $D$ are gains of the PD controller when it is activated. The PD controller is activated when the delayed state $(φΔ,ωΔ)$ is located in $Son$⁠, while it is inactivated when $(φΔ,ωΔ)$ is in $Soff$⁠. $Son$ and $Soff$ are separated by the boundaries defined by $φΔ=0$⁠, $ωΔ=lφΔ$⁠, and $φΔ2+ωΔ2=r2$⁠. The parameter $r$ determines the small circular sensory dead zone around the upright position. The parameter $l$ determines the slope of the on-off boundary line $ωΔ=lφΔ$⁠. In this study, we define the parameter $ρ$ that represents the ratio of $Son/(Son+Soff)$⁠, with neglecting the small sensory dead zone for simplicity. That is, by using the slope $l$⁠, $ρ≡0.5−arctan⁡(l)/π$⁠. $ρ=1$⁠, if the PD controller is continuously turned on (i.e., if S$off$ is null), which represents the continuous control model. $ρ<1.0$⁠, if the PD controller is switched off for a non-zero area of $Soff$⁠, which represents the intermittent control model. $ρ=0$ represents the inverted pendulum with no active controllers all the time. In this way, the continuous control model is part of the intermittent control model parameterized by the parameter $ρ$⁠. Our choice of the complicated geometry of $Soff$ has been encouraged by our preliminary study on reinforcement learning for stabilizing inverted pendulum based on the information of the time-delayed state point.⁴² That is, the energetically optimal choice of $Soff$ might not be necessarily composed of simple thresholds for the tilt angle of the inverted pendulum as in other controller models^8–10,15 with intermittent on-off switching mechanisms.

We refer the system with the PD controller for the intermittent control model to as on-subsystem, whereas that without the PD controller to as off-subsystem. In the intermittent control model, the system switches between the on- and off-subsystems. Overall dynamics of the intermittent control model heavily depend on dynamics of the on-subsystem. Since the on-subsystem, if it were persistently activated, is identical with the continuous control model, dynamics (including stability) of the on-subsystem can be characterized by the stability map shown in Fig. 1. Typically, the intermittent control model operates with small $P$ and $D$ gains,⁶ for which the $(P,D)$ of the on-subsystem is located outside (below) the stability region as exemplified by the $(P,D)$ point indicated by the diamond-shaped symbol in Fig. 1. In this case, the system switches between two types of unstable dynamics as illustrated in Fig. 2(c). One is the saddle type instability of the off-subsystem. The other is delay-induced unstable oscillation of the on-subsystem due to the delay $Δ$ in active feedback control. The rationale of the intermittent control strategy is that the unstable dynamics of the on-subsystem can be exploited, in the sense that the actively driven spiraling out of the state point may provide an opportunity to transverse the stable manifold of the off-subsystem at some time, and inactivation of the active feedback controller at that time would trigger another transient behavior approaching the saddle point along or near the stable manifold of the off-subsystem. Thus, the intermittent controller can achieve bounded stability with limit cycle oscillation of the upright posture.⁵ 

Figure 3 exemplifies dynamics of the intermittent control model for $ρ=0.62$ (operated as the intermittent model) and for $ρ=1.0$ (operated as the continuous control model), where the CoM positions $xCoM$ in the anterior–posterior direction are calculated as $xCoM=hφ$ to be compared with experimental sway data. See  Appendix A for numerical integrations of the stochastic differential equations in Eqs. (3) and (4). For both cases, time-delay is set as $Δ=0.2$ s. As mentioned above, the $(P,D)$ parameters for the former case are set as $P=147=0.25mghN\,m/rad$ and $D=10N\,m\,s/rad$⁠, which results in unstable dynamics of the on-subsystem. This can be confirmed by the $(P,D)$ point located outside the stability region, as indicated by the diamond-shaped symbol in Fig. 1. Note that the noise intensity is very small at $σ=0.2$⁠, but the sway dynamics exhibit slow fluctuation with a relatively large amplitude. These parameters are determined ad hoc to mimic postural sway of healthy people under the intermittent control hypothesis, but they will be validated quantitatively in this study. For the latter case, $ρ=1.0$ (the continuous control model), where the gain values are set as $P=471(=0.8mgh)N\,m/rad$ and $D=270N\,m\,s/rad$⁠. In this case, the $(P,D)$ point is located inside the stability region, as indicated by the open circle in Fig. 1. These $P$–$D$ values have been used by the previous study²⁷ to mimic optimally postural sway of healthy people under the continuous control hypothesis. Note that the sway amplitude is small, in comparison with the former intermittent case because of the small noise intensity (⁠$σ=0.2$⁠). This value is about one tenth of the noise intensity used in the previous study.²⁷ That is, the continuous control model requires large noise intensity to show sway amplitudes similar to human postural sway. Moreover, it requires low-pass filtered colored noise for mimicking postural sway of healthy people.

In Fig. 3, sample trajectories on the $xCoM$–$vCoM$ plane are also depicted. Moreover, other characterizations of the sway time series are shown for each of the intermittent and the continuous control models. They include histograms of (1) $|xCoM|$⁠, (2) $|vCoM|$ representing the CoM velocity, (3) $|aCoM|$ representing the CoM acceleration, and power spectral densities (PSD) of (4) $xCoM$⁠, (5) $vCoM$⁠, and (6) $aCoM$⁠. These six histograms are used as the summary measures for the approximate Bayesian computation (ABC) as described below. Each histogram can also be characterized by its median value or by the slope of linear regression line for PSD of $xCoM$ at the low-frequency regime (⁠$−β$⁠), and peak frequencies in the PSDs for CoM velocity and acceleration. In particular, as in Fig. 3(a), postural sway of the intermittent control model exhibits a power-law-like behavior,^6,7,10,43 characterized by the scaling exponent $β$⁠, which is known to be close to 3/2 for healthy young people. $β$ values for the continuous control model are usually small [close to zero in Fig. 3(b)], implying that sway patterns of the continuous control model in the low-frequency regime are similar to white noise.

In this sequel, we summarize how the above-mentioned six summary measures are computed, for both experimental and simulated time series. There are some preprocessing procedures specifically for experimental and simulated time series. For experimental data, we need to obtain CoM time series from the center of pressure (CoP) data, which were measured by a force platform. See  Appendix B for this transformation. Then, a linear trend of the CoM time series was removed. We denote the linear-detrended CoM signal by $xCoM$ or simply $x[n]$ for the discretized time $n$⁠. Velocity $v[n]$ (or $vCoM$⁠) and acceleration $a[n]$ (or $aCoM$⁠) were approximated as

where $fs$ is the sampling frequency (⁠$fs=100Hz$⁠). The absolute values $|x[n]|$⁠, $|v[n]|$⁠, and $|a[n]|$ were represented as the histograms with 15 bins for each. The bin-width was adjusted for each experimental time series data, where we determined the maximum value of the histogram by taking a smaller value of either the maximum absolute value in the data or three times of standard deviation of the data, and divided it by 15. PSD for each of $x[n]$⁠, $v[n]$⁠, and $a[n]$ was calculated via 75% overlapped Fast Fourier Transform with 40 s rectangular time window for each time series. Then, powers at 10 frequencies equally spaced in logarithmic scale were selected up to 1.5 Hz. Each of all histograms and PSDs was normalized so that the sum of values of the histogram (15 values) or that of the PSD (10 values) becomes unity. In this way, we obtained the 75-dimensional summary measure vector, composed of all elements of the six summary measures.

Computation of the summary measures for simulated CoM data (⁠$xCoM=hφ$⁠) was performed in the same way as that for experimental data. For a given experimental data, we first computed its 75-dimensional summary measure vector, which determined the bin-width of each histogram. Those bin-width values were used in common for the histograms to characterize simulated data generated during exploring/inferring the parameter values, specifically for these experimental data. We employed a feat specifically to obtain the histogram of acceleration only for simulated time series. This was because CoM accelerations computed by Eq. (10) for simulated time series were excessively affected by the additive noise, due to the Euler approximation of the stochastic differential equation. That is, as in Eq. (A4) in  Appendix B, the discretized $ωn$ (⁠$v[n]=hωn$⁠) is influenced directly by the addition of white noise $Wn$⁠. Thus, computation of accelerations from the differences between successive $ωn$ estimates too large accelerations, which caused large differences between the histogram of acceleration for simulated and experimental time series. To avoid this problem, the angular acceleration, referred to as $α~n$⁠, was approximated using the vector field of the model as

without adding white noise $Wn$⁠, where $d$ is the discretized delay. In this case, $(φn,ωn)$ and $(φn−d,ωn−d)$ are affected by white noise, but they are smoothed by the integration for obtaining a solution. A sequence of $α~n$ values was used for computing the histogram of CoM acceleration for simulated data. However, we used $a[n]$ in Eq. (10) for computing PSD of CoM acceleration for simulated data.

We used 70 s long CoP time series data during quiet standing with eyes-open, measured from 155 PD patients (81 males and 74 females), age-matched 21 healthy elderly people (9 males and 12 females), and 21 healthy young adults (21 males). The UPDRS (Unified Parkinson’s Disease Rating Scale) Part III was used for quantifying the severity of motor symptom. Here, we used the sum of scores for items from 27 to 30 that are related to posture and gait functions. The zero point means no symptom, and the 16 points are for the most severe symptom. See Table I for the mean score of the patients. Most data have been acquired and analyzed in our previous studies. More specifically, data were taken from four studies.^39,44–46 See those literatures for full details of CoP measurements. Some data from PD patients were newly obtained by the same protocol with the previous studies at Osaka Toneyama Medical Center, for which all subjects gave informed consent approved by the ethical committees of Osaka Toneyama Medical Center. Moreover, some PD patients participated in the measurement several times at different days separated more than a few months. We regarded those data measured from identical individuals with PD at different days as the ones from different patients because they might show different sway characteristics depending on the severity of the patient at the measurement day. Therefore, the total number of PD patients was considered 272. Except for the healthy elderly, each subject performed more than two trials of quiet standing at each session, with sufficient resting time between trials. In this study, we randomly picked two trials data from all trials for each participant and used them for computing the summary measures. Then, averaged measures were used for the parameter inference for each subject. Since healthy elderly subjects performed only one trial, the summary measures for each subject were computed from the single trial. For all data, the sampling frequency of the CoP acquisition was $fs=100Hz$⁠.

We use a technique of parameter inference using the approximate Bayesian computation (ABC) based on Sequential Monte Carlo sampling (SMC)⁴⁷ with high performance computing resources for estimating multiple parameters of a highly nonlinear stochastic dynamical system model. More specifically, the Bayesian inference provides a framework for identifying values of parameters $θ$ of the model, which would lead to a model-generated data that is sufficiently similar to the observed data $y0$ by means of the likelihood function $L(θ)$⁠. That is, according to Bayes’ theorem,

where $p(θ)$ is the prior joint probability distribution of the parameters and $p(θ|y0)$ is the posterior joint probability distribution given (conditional) the observed data. However, it is almost impossible to compute the likelihood function due to the complexity of the system or the unobservable random quantities in the model.

To deal with this issue, approximate Bayesian computation (ABC) has been developed and widely used,⁴⁸ in which parameter values are inferred using summary measures, denoted by $Φobs$ for the observed data, without computing the likelihood function. In the ABC algorithm, a candidate parameter vector $θ$ is randomly generated according to initial beliefs of the parameter distributions $p(θ)$⁠, and dynamics of the model with the parameter vector $y(θ)$ is numerically simulated. Then, summary measures for the model-simulated data, denoted by $Φsim$⁠, are also computed. If the “distance” between $Φobs$ and $Φsim$⁠, which measures the discrepancy between the two data sets in the sense of the summary measures, is smaller than a threshold value $ϵ$⁠, the parameter vector $θ$ is accepted and registered as an element of posterior distribution of the parameters. Thus, ABC algorithm samples from the posterior distribution of the parameters by finding values that yield simulated data sufficiently resembling the observed data, and the sampling is continued until the number of registered element becomes $N$⁠.

The approximate Bayesian computation based on Sequential Monte Carlo sampling (ABC-SMC) is one of the extended algorithm of ABC, which iterates ABC algorithm sequentially with changing the threshold value $ϵ$ smaller gradually. In this study, we assimilated the intermittent control model into postural sway data using the ABC-SMC, as in Tietavainen et al.¹⁶ That is, we performed ABC-SMC to infer parameters of the intermittent control model. We considered six parameters $(P,D,ρ,Δ,σ,r)$ of the intermittent control model, and their values were inferred by the SMC-ABC. Other parameters were fixed in this study. Body mass $m$ of each subject was determined by the weight of the subject. However, $h$ was fixed as 1 m. The moment of inertia $I$ was calculated as $mh2$⁠. The passive elastic and viscosity coefficients were set as $K=0.8mgh$ and $B=4.0$ as mentioned before. The prior distribution of each parameter was assumed to be the uniform distribution within a given range as follows. $P∈[0,500]$⁠, $D∈[0,500]$⁠, $ρ∈[0.3,1.0]$⁠, $Δ∈[0,0.5]$⁠, $σ∈[0,3.0]$⁠, $r∈[0,0.01]$⁠.

For simulations, the initial state was set as $φ(0)=η$⁠, $ω(0)=0$⁠, and $φ(τ)=ω(τ)=0$ for $−Δ≤τ<0$⁠, where $η$ was a uniformly distributed random number in $[−0.01,0.01]$ rad. The stochastic differential equations of the model were numerically integrated using the Euler–Maruyama method with time step $Δt=0.001s$⁠. Each simulated data were resampled with the sampling frequency of 100 Hz. The first 20 s was discarded, and the remaining 70 s long simulation data were used for computing the summary measures. As described in Sec. II, the summary measure for each of experimental sway data $Φobs$ and simulated data $Φsim$ was 75-dimensional vector.

In this study, the distance between $Φobs$ and $Φsim$ was calculated using Jensen–Shannon divergence⁴⁹ $DJS$⁠, which is based on Kullback–Leibler divergence $DKL$⁠, defined as

where

and

where the index $i$ of summation indicates the number of elements consisting of the summary measure vector. ABC-SMC was performed with the initial threshold $ϵ=1$ and $N=500$⁠, and continued until the ABC-SMC session reached the threshold $ϵ=0.01$⁠, or the computation time was exhausted.

After the first ABC-SMC session was terminated, we examined whether the obtained posterior distribution of each parameter, i.e., each of six marginal distributions, exhibited multiple peaks (multimodality) or not. If none of six marginal distributions exhibited multimodality, the parameter inference and data assimilation were completed. However, if some of them exhibited multimodality with multiple peaks with sufficient valleys between them, we tried to improve the parameter inference, where ABC-SMC was performed again, but now starting from a non-uniformly distributed prior distribution. Namely, the posterior marginal distributions of the first ABC-SMC were divided into two sets of parameter values, according to the valleys of the posterior marginal distributions, and two prior joint distributions, each of which was composed of one of the two sets, were constructed. Then, the second ABC-SMC session was performed using each of the two prior joint distributions as in the first ABC-SMC session. By comparing the resultant two posterior distributions, we took one of them that showed smaller variance with single mode as the posterior distribution to complete the inference. However, if both of two posterior distributions still exhibited multimodality, we considered that the postural sway data for that subject could not be fitted by the intermittent control model, and the data were discarded from the subsequent analysis.

In this study, we used three High Performance Computers [one HPC5000 10Core (3.1 GHz) $×$ 2CPU, RAM 64GB and two HPC5000 12Core (3.1 GHz) $×$ 2CPU, RAM 64GB, HPC SYSTEMS]. We used 64 cores simultaneously. Eight cores were assigned for the parameter inference of sway data from a single subject, and thus the parameter inference using ABC-SMC was performed for eight subjects in parallel, which was completed by about 12 h, if all inferences were successfully completed only by the first session. The computation time could be much longer if second sessions were required.

For each experimental data with $Φobs$⁠, we obtained six posterior marginal distributions. We then characterized each of them by the mode value. Thus, each experimental data with $Φobs$ were characterized by a six-dimensional parameter point $(P,D,ρ,Δ,σ,r)$⁠. After completing the data assimilation for all data, we obtained a number of parameter points distributed in the six-dimensional parameter space. Each point characterizes individual postural sway based on the model.

Because we recognized that distributions of the parameter values of $σ$ and $r$ in the six-dimensional parameter space were relatively homogeneous, we performed the hierarchical cluster analysis to classify the parameter points in the four-dimensional space of $(P,D,ρ,Δ)$⁠. To this end, each of four parameters was standardized so that the mean and variance were zero and unity, respectively. Then, according to the centroid clustering method, hierarchical clusters were constructed using the $L2$ norm. A statistically meaningful number of groups (clusters), referred to as $n∗$⁠, was determined by the upper-tail method.⁵⁰ In this method, for a given number of clusters $n$⁠, distances between all possible combinations of centroids of the $n$ clusters denoted by $di,j$ for $i,j∈{1,…,n}$ were computed, by which the minimum distance $an$ was obtained as $an=mini,j(di,j)$⁠. Then, $n∗$ was determined as the largest number that satisfies $an∗<a¯+kσ$ for a positive integer $k$⁠, where $a¯=∑n=2Nan/(N−1)$ and $σ=∑n=2N(an−a¯)2/(N−2)$ with $N$ being the number of data points. We used $k=3$ and $k=4$⁠, which are said to be appropriate for $N$ of several hundreds. After clustering, all parameter points were plotted on the $P$–$D$ plane, separately for each cluster, with stability region of the on-subsystem for the inferred delay $Δ$⁠, by which we examined whether the parameter points were located inside or outside the stability region for the on-subsystem.

The parameters could be inferred successfully for sway data from most subjects (271 out of 314 subjects). Details of discarded data of 43 subjects are summarized in Table II. Figure 4 summarizes the hierarchical cluster analysis for the 271 inferred parameter points of $(p,D,ρ,Δ,σ,r)$ in the four-dimensional space of $(p,D,ρ,Δ)$⁠, where $p=P/mgh$⁠. Four to five clusters (groups) were determined as statistically meaningful by the upper-tail method, for which we discarded tiny clusters composed of only nine data in total as indicated by “ungrouped data” in Fig. 4. According to the result of the upper-tail method, we decided to consider five groups, which were clustered into two major groups. One major group (data from 150 out of 271 subjects) was composed of the groups A and B, and the other major group (data from 112 out of 271 subjects) was composed of the groups C, D, and E.

For the major group with groups A and B, the inferred $ρ$ values were distributed between 0.4 and 0.6, meaning that the model operates with the intermittent strategy (Fig. 5). As shown in Fig. 4, sway data from almost all healthy young participants (16 out of 18 subjects) were classified into this major group. In representative experimental sway and model-simulated sway with its inferred parameter values for each of the groups A and B, CoM sway exhibited slow variations with relatively large amplitudes. The $xCoM$–$vCoM$ planes for the model with the inferred parameters show that the PD controller is switched off when the state point is located near the second or the fourth quadrant of the plane. The mean value of posture and gait related UPDRS for the patients classified into group A was 4.6 out of 16 points, which was much smaller than that for the other patients. Specifically, mode and median values of the posterior marginal distributions for group A were $(p,D,ρ,Δ,σ)=(0.29,6,0.58,0.35,0.12,0.0006)$ and (0.32, 15, 0.58, 0.36, 0.13, 0.001), respectively. Those for group B were $(p,D,ρ,Δ,σ)=(0.49,72,0.42,0.33,0.14,0.0008)$ and (0.49, 84, 0.43, 0.34, 0.17, 0.0011), respectively. The inferred $(P,D)$ points of the on-subsystem for groups A and B were mostly located outside (below) the stability region in the $P$–$D$ parameter plane for delay $Δ$ about 0.2–0.4 s, as in Figs. 5(a) and 5(b), respectively. Particularly, $D$ values were very close to zero for group A [Fig. 5(a)]. This means that the model with those inferred parameter values operates as the intermittent control model that switches between two unstable dynamics of the off- and on-subsystems. The large amplitude characteristics of group A can be confirmed by broadly distributed shape of the posterior marginal distribution of $xCoM$ with the median value of 3.4 mm as shown in Fig. 5(a). Values of the scaling exponent $β$ for the PSD of $xCoM$ were close to 3/2 (1.71 for group A and 1.21 for group B).

The other major group was composed of groups C, D, and E. As shown in Fig. 4, $ρ$ values for groups D and E were close to 1.0. Thus, the model with the inferred parameter values for groups D and E operates with the continuous strategy. However, $ρ$ values for the group C were close to 0.5, indicating that group C operates with the intermittent strategy. As shown in Fig. 4, two healthy young subjects, other than 16 subjects in groups A and B, were classified in group C, but none of healthy young subjects were classified in groups D and E. Mean values of posture and gait related UPDRS for the groups D and E were 6.6 and 8.4 out of 16 points, which were larger than the mean value for patients in groups A, B, and C (Fig. 4). In representative experimental sway and model-simulated sway with its inferred parameter values for each of groups D and E, CoM sway amplitudes were relatively small. This can be confirmed by the less-distributed posterior marginal distribution of $xCoM$ with the median value of 2.1 and 2.3 mm for groups D and E, respectively, as shown in Figs. 6(b) and 6(c). Experimental and simulated sway data for group C were similar to those in groups A and B. The $xCoM$–$vCoM$ planes for the model with the inferred $ρ$ values for group C show that the off-region S$off$ was also similar to that of group A. Specifically, mode and median values of group C were $(p,D,ρ,Δ,σ)=(0.22,135,0.61,0.49,0.17,0.0011)$ and (0.24, 135, 0.64, 0.40, 0.18, 0.0014), respectively. Those for group D were $(p,D,ρ,Δ,σ)=(0.26,171,0.96,0.30,0.15,0.0006)$ and (0.26, 177, 0.94, 0.33, 0.24, 0.0011), respectively. Those for group E were $(p,D,ρ,Δ,σ)=(0.46,177,0.98,0.005,0.53,0.0008)$ and (0.44, 222, 0.96, 0.015, 0.66, 0.0011), respectively. The $(P,D)$ points of the on-subsystem for groups C, D, and E were located inside the stability region of the on-subsystem as shown in Figs. 6(a)–6(c), which is achieved by large values of the $D$ gain for groups C, D, and E. This means that the model for groups C, D, and E switches between stable dynamics of the on-subsystem and unstable dynamics of the off-subsystem. That is, the stabilizer of the model for these groups is the convergent dynamics of the on-subsystem. Values of the scaling exponent $β$ for the PSD of $xCoM$ were 1.08 for group D and 0.91 for group E, which were smaller than those for groups A, B, and C. Note that the inferred values of feedback delay $Δ$ for group E were almost zero. Moreover, the inferred noise intensities of group E were about four times larger than those of the other groups.

The intermittent delay feedback control model was assimilated into postural sway data from healthy young and elderly subjects as well as patients with Parkinson’s disease. A joint probability distribution of the parameter point $(P,D,ρ,Δ,σ,r)$ of the model was inferred for the sway data of each subject using the ABC-SMC method. Particularly, we focused on the parameters $ρ$ representing the ratio of on-region/(on-region+off-region) in the phase plane: $ρ∼1.0$ implies the model stabilized by the traditional continuous control strategy, while $ρ∼0.5$ implies the intermittent control model. Hierarchical cluster analysis based on the mode values of the inferred marginal distributions for four parameters $(P,D,ρ,Δ)$ uncovered two major groups within the parameter space. Based on the discussion which we will present next, we conclude that the one group represents stabilization by the intermittent control strategy (i.e., the intermittent control group) and the other by the continuous control strategy (i.e., the continuous control group).

For the intermittent control group (150 out of 271 subjects), $ρ$ values were close to 0.5 with small $D$ gains. More specifically, $(P,D)$ points for the on-subsystem were mostly located below the stability region in the $P$–$D$ parameter plane (Fig. 5), implying that people belong to this major group stabilize upright posture by the convergent dynamics along the stable manifold of the off-subsystem.^5–7 Sway data of almost all healthy young participants (16 out of 18 subjects) were better fitted by the intermittent control model rather than the continuous control model. Intermittent control for group A can be considered high efficiency based on the following characterizations. The $ρ$ values for group A were distributed between 0.5 and 0.6, which means that the PD controller is switched off when the state point is located in the second or the fourth quadrant of the phase plane, as shown in Fig. 4 (top). The $(P,D)$ points of the on-subsystem for group A were located completely below the stability region [Fig. 5(a)]. Particularly, $D$ values were very close to zero. Moreover, $P$ values were also small (around $0.3mgh$⁠). Thus, the intermittent control model with those parameter values is highly energetically efficient,^43,51 because energy consumption is zero (null) when the system is operated with the off-subsystem, and it is also close to zero even when the system is operated with the on-subsystem, because the $P$ and $D$ gains (and thus the active feedback torque) of the on-subsystem are very small. Intermittent control for group B might also be energetically efficiency but to a lesser degree compared to group A. This is because $ρ$ values for group B were distributed around 0.4, which means wider off-regions than group A. As discussed previously,⁶ the wide off-region for the intermittent control model requires large $P$ gains as well as large $D$ gains for stability, although the $(P,D)$ points for group B were located mostly below the stability region of the on-subsystem [Fig. 5(b)]. Posture and gait related UPDRS scores for the patients belonging to group A were smallest among those belonging to the other groups, followed by the scores for group B. This means that postural balancing capability of patients in groups A and B was less affected, while they are still capable of utilizing the intermittent control strategy.

Moreover, 112 out of 271 subjects were classified into the continuous control group (the groups belonging to C, D, and E), including two healthy young subjects in group C. $ρ$ values for groups D and E were close to 1.0, indicating that the system never utilizes the off-subsystem. $ρ$ values for group C were close to 0.5, implying intermittent switching between the off- and the on-subsystems. However, the $(P,D)$ points used for the on-subsystem of group C were located inside the stability region in the $P$–$D$ parameter plane [Fig. 6(a)]. Thus, we regarded the model with such stable on-subsystem as a type of continuous control model, although the model switches between the off- and the on-subsystems. This is because the stability of the upright posture in this case is achieved by stable dynamics of the on-subsystem, as in the typical continuous control model, not by the convergent dynamics of the off-subsystem. In other words, upright standing of people belonging to this major group is stabilized by the stable dynamics of the on-subsystem. People belonging to this major group are mostly subpopulations of elderly and Parkinson’s patients. Specifically, posture and gait related UPDRS scores for the patients belonging to the continuous control group were notably larger compared to those of the patients belonging to the intermittent control group.

Group D represents a typical continuous control strategy. The parameter points for group D are close to the one estimated by many studies that utilized the traditional continuous control model.²⁵ The inferred $P$ and $D$ gains, particularly $D$ gains, are notably larger than those for the intermittent control group, suggesting that continuous control strategy is energetically inefficient. Moreover, large $P$ and $D$ gains imply inflexible ankle joint (inflexible upright posture), as recognized by some previous studies.^34,38,39 Scaling exponent at the low-frequency regime of the PSD of the CoM sway was smaller than $β∼3/2$⁠, and relatively close to $0$ as discussed earlier.^6,52 That is, the PSD of CoM without scaling behavior at the low-frequency regime (closer to white noise) indicates a loss of intermittency in the feedback control. Posture and gait related UPDRS scores for patients belonging to the typical continuous control group were second largest, implying that these subjects exhibited severe postural symptoms.

Group E, in the continuous group, can be considered the non-reactive control group because the inferred values of the time delay $Δ$ were almost zero. This means that the feedback controller is not affected by time-lags for the signal transmission. This situation represents the traditional stiffness control, advocating that the stabilization of upright posture during quiet stance is achieved by stiffness of the ankle muscles.²² Stiffness control assumes that the central nervous system (CNS) does not actively change the level of motoneuron activation that determines stiffness of ankle muscles and the consequent stiffness of ankle joint during quiet stance even in the presence of small postural fluctuation (postural sway), but it is predetermined by the CNS in a feedforward manner. Thus, the feedback controller for the stiffness control is mathematically equivalent to the passive stiffness that does not involve any time delay. Because the inferred $P$ and $D$ gains for group E, particularly $D$ gain values, were very large, we conclude that the upright stance of people in this group is rigidly stabilized in a non-reactive manner using large values of the passive elasticity and viscosity at the ankle joints. Together with the fact that posture and gait related UPDRS scores for the patients belonging to the non-reactive control group were the largest, individuals in this group exhibited the most severe postural symptoms, and their ankle joints might also be too inflexible to induce postural fluctuations.^38,39,46

The results of this study, which are based on the intermittent control model, have similarities and discrepancies compared to the model-based characterization by Maurer and colleagues that utilizes the continuous control model.^27,33 A similarity is an overall qualitative tendency of large feedback gains to characterize postural instability with severe symptoms in patient with Parkinson’s disease. This is simply because the intermittent-control-model-based parameter exploration in this study inferred the continuous control model for sway data from patients with severe postural symptoms. Regardless, postural inflexibility in patient with Parkinson’s disease might be well characterized by continuous control strategy with larger feedback gains. Discrepancies become evident in quantitative characterizations of sway data from healthy people. The intermittent-control-model-based parameter exploration inferred an intermittent control model with small values of $P$ gain and very small values of $D$ gain for the on-subsystem. It also inferred a small noise intensity. On the other hand, continuous-control-model-based parameter exploration might infer relatively large values of $D$ gain and large noise intensity. Large values of the $D$ gain are inevitable for the continuous control model to avoid delay-induced instability, which imposes a large noise intensity to counteract strong damping force due to large $D$ gain for generating relatively large postural sway observed in healthy people. Note that the intermittent control model can generate oscillatory slow sway with relatively large amplitude even if the model is driven by white noise of very small intensity by the nature of the control strategy.⁵² 

For some groups obtained in this study, the posterior distributions showed relatively large deviations (i.e., less sharp peaks). For example, $D$ gain values in groups B and E and the delay parameter $Δ$ in group C were broadly distributed. Such uncertainty in the parameter estimation might be related to the sensitivity of the inverted pendulum balance control model. In this regard, it is important to recognize a critical difference in the sensitivity between the continuous and intermittent control models. The continuous control model is sensitive to the $D$ gain,⁵³ because a large delay $Δ$ makes the stability region small (the delay-induced instability), which enforces a tuned large value of the $D$ gain for maintaining stability. Contrastingly, the intermittent control model is quite insensitive to the $D$ gain in the wide area of small $D$ gain regime. This is because stability of the intermittent control model is very robust against the variation of $D$ value, as shown in our previous study,⁶ which is a major cause of the uncertainty in the $D$ value. Moreover, in the intermittent control model, a large value of the delay $Δ$ should be accompanied by a large area of the off-region (⁠$ρ$ value smaller than 1.0, typically around 0.5–0.6) and also a small value of the $D$ gain for making the switching system stable because the delay-induced instability of the on-subsystem should be compensated by the use of the off-subsystem. This way, values of the $D$ and $Δ$ are mutually dependent oppositely in the continuous and the intermittent control models.

ABC-SMC performed in this study for postural sway data is based on the work by Tietavainen et al.,¹⁶ with a critical improvement in the metric for quantifying distances between sets of summary measures. That is, the use of the Jensen–Shannon divergence (similar to Kullback–Leibler divergence) improved the accuracy of parameter inference, particularly for differentiating between continuous and intermittent control models. This improvement was needed since we considered a wider range of parameter exploration compared to Tietavainen et al.¹⁶ across the boundary separating the intermittent and continuous control models.

There are other methodologies useful for parameter estimations of nonlinear dynamical systems. The intermittent control model used in this study has been applied for parameter estimation by a Kalman filter based method¹⁷ and by Adaptive Gaussian process approximation.⁵⁴ In particular, McKee and Neale¹⁷ applied their method to postural sway data from healthy people and reported similar parameter values as our current study. Both methodologies reported relatively large values of time delay $Δ∼0.3$–$0.4s$⁠, which should be compared with a typical delay of about 0.2 s, as estimated previously.²⁵ Although this issue can be argued further, we will address it in the future.

Buza et al.⁵⁵ and Milton et al.⁵⁶ studied motor control during balancing tasks, in which they also tried to trace a path in the parameter space that is associated with a simple feedback controller and/or a prediction-based feedback controller, as human subjects become familiar with and expert of the tasks. In the case with moderate instability in the task, a linear delay feedback controller may be employed, in which the gain parameters may converge to an optimal point in the parameter space, leading to high linear stability with energetically efficiency.⁵⁵ In the case with large instability, such as during short-stick balancing, the gain parameters tend to converge to the edge of stability for either with continuous controllers⁵⁶ or with time-discontinuous intermittent controllers.^57,58 The latter situation for expert subjects may be analogous to the high-efficiency intermittent control strategy implemented with the $(P,D)$ points near (but outside) the edge of linear stability region for the on-subsystem [Fig. 5(a)].

Finally, some of the recent studies that used the intermittent postural control model suggest that intermittency in the feedback control is harmful and that it can lead to postural instability.^59,60 However, the results of our current study are apparently contradictory to those previous studies. Such discrepancy may be due to inadequate (non-systematic) parameter exploration, by which one might overlook the beneficial aspects of the intermittent control strategy.

Overall, the results of our current study may be interpreted by postulating that increased postural instability of most Parkinson’s patients and some elderly persons might be characterized as a dynamical disease.^61–63 Indeed, it has been demonstrated that a sequence of period-doubling or period-doubling-like bifurcations occurs as a parameter that determines the (size of) off-regions in the phase plane, corresponding to $ρ$ in this study, changes.⁵² This means that dynamics of the intermittent control model and those of the continuous control model are qualitatively different and that transitions between them are accomplished through the sequence of bifurcations.

In this study, inference of postural control strategy was performed using sway data acquired from a number of people acquired at a single or a few different instances of time. Such a “transverse study” that could provide hard evidence of a transition between different control strategies (with an increasing $ρ$-parameter) is almost impossible to carry out in practice. However, one might consider observing sway data from the same subject over very long periods of time with an increasing physical or mental load, expecting to observe transitions and an increase in the value of $ρ$⁠.

This study was supported by the following grants from the Ministry of Education, Culture, Sports, Science and Technology (MEXT)/Japan Society for the Promotion of Science (JSPS) KAKENHI: No. 20H05470 (T.N.), No. 19H04181 (T.N.), No. 17K13016 (Y.S.), No. 20K11989 (Y.S.), and No. 20K19412 (M.M.), and by the Uehara Memorial Foundation (T.N.).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The state space representation of the model equations (3) and (4) is rewritten as

The Euler approximation of this simultaneous stochastic delay differential equations with a fixed time step $Δt=tn+1−tn(n=0,1,…)$ is then described as

where $d=Δ/Δt$⁠. $Wn$ is independent Gaussian random variables with mean $E[Wn]=0$⁠, variance $E[Wn2]=1$⁠, and autocorrelation $E[WnWm]=δnm$⁠.

The transformation from CoP to CoM was performed as follows. First, anterior–posterior component of the CoP signal was low-pass filtered using the fourth-order zero-phase-lag Butterworth filter with the cutoff frequency of 10 Hz. Then, CoM positions were estimated from the filtered CoP using the filter proposed by Morasso et al.⁶⁴ defined as

where $Ω$ is the angular frequency and $he$ is the parameter of “effective distance.” In this study, we set $he=1$ as a typical value, because it has been shown that the transform from CoP to CoM is less dependent on the value of $he$⁠.