Abstract
We are interested in understanding the mechanisms behind and the character of the sway motion of healthy human subjects during quiet standing. We assume that a human body can be modelled as a singlelink inverted pendulum, and the balance is achieved using linear feedback control. Using these assumptions, we derive a switched model which we then investigate. Stable periodic motions (limit cycles) about an upright position are found. The existence of these limit cycles is studied as a function of system parameters. The exploration of the parameter space leads to the detection of multistability and homoclinic bifurcations.
1. Introduction
In spite of recent advances in revealing the mechanisms responsible for balancing during quiet standing (see [1–11]) the physiological mechanisms are far from being understood. In recent years, there has been a growing interest among scientists to use mathematical modelling and numerical simulations to gain new insights into the problem of balancing. In the literature, it has been argued that forward body sway can be captured using models where the muscles and tendon–muscle complexes act as springs with certain stiffness and the neuromuscular system generates the corrective torque [12–15]. Models based on these assumptions have been tested experimentally [8,12,13,16]. In the studies of Maurer & Peterka [8] and Peterka [16], a singlelink inverted pendulum model with simple linear feedback was introduced. Different feedback laws were considered to obtain the best match to experimental data, and it was shown that a corrective torque proportional to the position and velocity signals combined with time delays can be used to account for the body sway observed experimentally.
Early work [17] suggested that intermittent control plays an important role in human motor control. Recent research [18,19] shows that intermittent control is a mechanism that can explain the results found in experimental tests. However, the authors also point out that there are alternative control models that can reproduce experimental results. Further research is needed to design experimental tests that can help to discriminate between different control models. An alternative model for human postural sway was considered in Eurich & Milton [20]. In this work, the authors use a reduced firstorder model with a timedelayed feedback. In place of intermittent control, Eurich & Milton [20] consider a switched control system depending on a threshold for the angular position. For the firstorder model, Eurich & Milton [20] show the existence of bistable limit cycles. In the study of Milton et al. [21,22], the switchlike balance controller was considered as a type of intermittent control. The idea of switched control was further explored in Milton et al. [23,24]. The authors considered a model with switched control, random perturbations that model noise, and time delay. It was shown that the interplay between noise and time delay may have a stabilizing property. The authors also argued that the control applied by humans is an adaptive switched control. The idea of switched control applied by humans during quiet standing serves as a starting point for the modelling approach adopted in this paper.
We approach the problem of human balancing by considering the subjects standing on both legs with eyes closed or open. We are interested in investigating the sway motion that occurs in the sagittal plane, i.e. we consider a forward–backward body sway.
The following assumptions are made:
— the body control is achieved using proprioception only (reception of internal signals such as posture and body sway through the length and elongation velocity of the muscles);
— we include time delays in the control system that represent the combined delays owing to sensory reception, neural transmission, neural processing, muscle activation and force development of the proprioception motorneural control system;
— we exclude the control effects of exteroception (sensing of external signals such as pressure) and vestibular system (it senses the angular velocity and the acceleration of the head). We believe that an inclusion of these effects should yield additional terms in our model such as, for instance, force feedback or a control term proportional to the measured acceleration;
— we assume a threshold value of the angle of the sway below which the motorneural control of the proprioception is not applied, similarly as in the work of Eurich & Milton [20]. This assumption is justified by the finite accuracy of sensing, as well as by the recent finding of impulsive like control muscle movements reported in Loram et al. [18];
— we assume that the motorneural control of the musculoskeletal system when sensing the error works like a PD (proportional–derivative) control system with the time delay in the position and velocity error signals. The letter ‘P’ refers to the corrective torque that is proportional to the error signal between the desired angle θ_{ref} (θ_{ref} = 0 in our case) and the measured angle multiplied by the proportionality factor, say K_{p}, that can be thought off as a stiffness factor. The letter ‘D’ refers to the corrective torque that is proportional to the velocity of the error signal between the desired velocity of the angle ( in our case) and the measured angular velocity with the proportionality factor, say K_{d}, that can be thought off as a damping factor.
Owing to the presence of a threshold value, our model is a hybrid (switched) system with time delay. The model we propose is similar to the one studied in Asai et al. [25]—the PD control is switched on or off depending on the value of state variables. The model studied in Asai et al. [25] is an extension of the studies conducted in Bottaro et al. [26,27] where the authors introduce intermittently switched on/off controller that allows for a bounded motion, and it is this type of control which Asai et al. [25] propose as the possible explanation for the sway motion during quiet standing. In the study of Asai et al. [25], it is also shown that the bounded dynamics is linked with a motion along the stable manifold of a saddle point. The control strategy used in Asai et al. [25] is based on the assumption that a bang–bang control does not allow sufficiently small bounded motion [27] to be produced, with the sway angle of about 1°, which could correspond to the sway motion during quiet standing. In our paper, we consider the simplest possible form of a dead zone, which is physiologically feasible, and we show in §4.4 that small scale stable periodic oscillations may be observed in bang–bang control systems in the presence of a dead zone.
We briefly summarize the results of our investigations. In order to maintain balance, the control of the proprioceptive system must use both: (i) the information on the elongation and (ii) the information on the velocity of the elongation of the muscles. This agrees with the fact that to stabilize an inverted pendulum (for small angles θ when the linear model of the pendulum is assumed) a linear controller must at least use the proportional and derivative signals of the error. Stabilization is achieved through the existence of stable oscillations about an offset angle ±θ_{0} where θ_{0} corresponds to the size of the dead zone. For a broad range of parameters, these oscillations are accompanied by stable larger scale oscillations. The difference in the amplitude of the angle of larger scale oscillations in relation to smaller scale oscillations is at least two orders of magnitude. In the context of upright standing, the small scale oscillations can be seen as a jitter about the nearly upright position; the ideal upright position cannot be achieved owing to our assumption that the motorneural control of the musculoskeletal system only reacts to values larger than the aforementioned offset angle θ_{0}, which is physiologically determined. By adding white noise to our model equations, we numerically show that close to a homoclinic bifurcation the system dynamics switches between two symmetric attractors. This scenario can explain noiseinduced switchings between two coexisting attractors reported in Eurich & Milton [20]. Finally, we should note that multistability as well as the presence of periodic oscillations have already been reported in the literature on switched delay differential equations, see [28–31]. However, the delay differential equations studied in these works contain discrete time delay in the position feedback only.
The rest of the paper is outlined as follows. In §2, a class of systems of interest is introduced. The system without the time delay in control function is first studied in §3. We then go on to explore the effects of time delay on system dynamics in §4. In the following §5, we introduce the modifications to the model that will be tested in the future. Finally, §6 concludes the paper.
2. The model equations
We simplify the biomechanics of the body by representing it as an inverted pendulum with the body sway occurring in the sagittal plane about the ankle joint axis. Gravity g acts on the centreofmass when the angle θ (measured in radians) between the vertical ankle joint axis and the body's position becomes nonzero; when there is no sway the body is vertical and θ = 0. The centreofmass m is located at height h above the ankle joint axis.
We assume that to control the upright position a corrective torque T is applied through a PD controller when some fixed, but nonzero, positive threshold θ_{0} is detected.
This leads to the following model equations: 2.1when there is no control applied to the system, and 2.2when the PD control is switched on; J is the moment of inertia of the body about the ankle joint axis. The delay terms are present in the applied torque generated by the PD controller. Namely, 2.3where K_{p} and K_{d} are negative constants, and Δ_{1}, Δ_{2} > 0 are time delays.
We make a number of observations regarding the physical interpretation of the model. For ϕ ≤ ϕ_{0}, the equations that govern system dynamics describe the dynamics of an inverted pendulum falling under the force of gravity. On the other hand, for ϕ > ϕ_{0}, we have . For Δ_{1} = Δ_{2} = 0 and sin(ϕ)≈ϕ, this equation is an equation of motion of a linear oscillator where the stiffness coefficient is given by (−K_{p}−mgh)>0 and the damping by −K_{d}>0. For Δ_{1} and Δ_{2} nonzero and positive, we obtain a secondorder linear delay differential equation whose dynamics is significantly more complex than that of a secondorder linear oscillator [32,33]. For sufficiently small Δ_{1} and Δ_{2}, and under appropriate continuity assumptions, the dynamics of the delayed system can be studied using a set of ordinary differential equations (ODEs) [34]. However, the presence of switchings makes this reduction impossible even for small values of time delays.
Making the approximation sinϕ ≈ ϕ, which is justifiable for small angles ϕ of the body sway, the model equations (2.1) and (2.2) become 2.4and 2.5
To reduce the number of parameters, we study system (2.4) and (2.5) in the nondimensional form; we introduce the nondimensional time and set , . Then and and we obtain 2.6and 2.7where , , with B = hK_{p}/gJ, , and . In what follows, we drop the bar symbol when referring to nondimensional time.
Setting and using firstorder representation, we arrive at a planar switched (Filippov) system of the form 2.8and 2.9where 2.10
3. System dynamics for τ_{1} = τ_{2} = 0, and τ_{1} ≠ 0, τ_{2} = 0
3.1. Switched system with no time delays
We start our investigations by discussing the dynamics of switched systems (2.8) and (2.9) for τ_{1} = τ_{2} = 0. We first make the following observations:

— let x = (θ, x) then, by the symmetry, we have that −F_{in}(−x) = F_{in}(x) and −F_{out}(−x) = F_{out}(x);

— the integral curves of F_{in} are symmetric with respect to the θ and axes.
Second observation implies that if an integral curve of F_{in} crosses the switching line Σ_{1} at some point a ≠ 0 then it either crosses Σ_{1} = {θ = θ_{0}}, at −a or Σ_{2} = {θ = −θ_{0}} at a. Similarly, if we consider a ∈ Σ_{2} then the integral curve either crosses Σ_{2} at −a or Σ_{1} at a.
The only equilibrium of equations (2.8) and (2.9) is the origin or a pseudoequilibrium. A pseudoequilibrium of equations (2.8) and (2.9) is the equilibrium of the full systems (2.8) and (2.9) that lies on the switching line Σ_{1} or Σ_{2}.
For A>0, the origin is an equilibrium of a saddle type characterized by the eigenvalues and with the associated eigenvectors and .
Lemma 3.1. Consider equations (2.8) and (2.9) with A > 0, and B and C negative , and such that A + B < 0. Then the only equilibrium of equations (2.8) and (2.9) is the origin. There are also two accumulation points (pseudoequilibria) which are located at (−θ_{0}, 0) and (θ_{0}, 0). The origin is an equilibrium of a saddle type and the pseudoequilibria are the only two global attractors of equations (2.8) and (2.9). The basins of attractions of these two pseudoequilibria are separated by the piecewisesmooth invariant manifolds of the saddle point of the full systems (2.8) and (2.9).
Proof. The first part of the lemma is trivial. The righthand sides (RHSs) of equations (2.8) and (2.9) consist of sets of linear equations. If the determinants of matrix L given by equation (2.10) and the matrix are nonzero, which is true for A>0 and A + B<0, then the only possible equilibria are located at the origin. It is vector field F_{in} that is defined in the neighbourhood of the origin. Then for A>0, the equilibrium is of a saddle type as the eigenvalues of L are and . The only other possible equilibria are pseudoequilibria which may only occur on the switching lines at points where the component of vector fields F_{in} and F_{out} is 0. It is easy to verify that this is the case only at (−θ_{0}, 0) and (θ_{0}, 0).
To prove the second part of the lemma, we first note that any point (θ, x) within the zone θ≤θ_{0} reaches some point on ± θ_{0} in finite time at some point, say P_{1}. Let us suppose that P_{1} lies on the positive part of Σ_{1}. From P_{1}, the flow follows the flow generated by F_{out} until some point P_{2} ∈ Σ_{1} is reached. The equations of motion along this segment of the flow are given by
Let and set and multiply both sides by . We then have
Thus, P_{2} must be closer to the origin than P_{1}. From P_{2}, the flow follows
Using the energy argument again, and owing to the fact that Σ_{1} and Σ_{2} are placed the same distance away from the origin, the energy at the next point of switching with Σ_{1} or Σ_{2} is the same as at P_{2}. There is a loss of energy only owing to the application of the flow F_{out} and after subsequent switches between F_{in} and F_{out} the trajectory will converge to either of the two pseudoequilibria since at these two points the system F_{out} reaches its lowest energy.
Finally, to determine which pseudoequilibrium is reached, we have to find the intersection points of the unstable and stable manifolds of the saddle point with Σ_{i} (i = 1,2). Concatenated segments joining subsequent intersection points form piecewisesmooth stable and unstable manifolds. By considering the evolution of a single trajectory, we can then show that the trajectory will be moving within a region bounded by the stable and unstable manifolds of the saddle point until it reaches either one of the two pseudoequilibria. Owing to system's symmetry if some point P_{in} converges to (θ, 0) then −P_{In} will converge to −(θ, 0).
Illustrative trajectory of switched systems (2.8) and (2.9) with A = 0.5, B =−0.6, C =−0.5, θ_{0} = 1 and the initial conditions (θ_{0}, x_{0} ) = (1,0.4) is depicted in figure 1.
If A + B < 0 and C is identically 0, it is easy to verify that our system exhibits an infinite number of periodic orbits which can be thought of as centres. Two representative periodic orbits for this case are depicted in figure 2.
If, on the other hand, B = 0 and C < 0 then the system trajectories will diverge to ± ∞.
We conclude this section with the physical interpretation of our results. When both the proportional and derivative control are applied instantaneously (no time delay owing to sensing, processing and actuation) then provided that A + B < 0 it is possible to achieve perfect stabilization. However, the angle θ at which the body is stabilized is offset from the vertical angle θ = 0.
3.2. Effects of the damping term and the position delay (τ_{1} ≠ 0 and τ_{2} = 0)
Our model system can be thought of as a generalization of the model studied by Eurich & Milton [20]. In this work, the authors model human postural sway using a secondorder ODE, which models inverted pendulum falling under the force of gravity, subject to timedelayed restoring force and noisy perturbation. Using an assumption that the postural sway is overdamped, for healthy subjects with eyes open, they reduce their model to a firstorder system. Our numerics suggest that this reduction removes stable pseudoequilibria present in the secondorder system. Thus, the reduced system does not have the possibility of producing multistable behaviour in the sense of having more than one pair of stable asymmetric orbits; by the system's symmetry stable asymmetric orbits come in pairs.
We observed the existence of stable pseudoequilibria in equations (2.8) and (2.9) with (2.8) having additionally a damping term . The switched system studied in this paper with the additional term in equation (2.8), and with τ_{2} = 0, can be thought of a system with damping and delayed position feedback and it is then the full model system studied in Eurich & Milton [20]. Numerical results with B varied and fixed A = 0.5, C = −10, time delay τ_{1} = 0.2 and τ_{2} = 0 are shown in figure 3. In figure 3a, for B = −0.49, the pseudoequilibrium located at (1,0) is unstable. We depict a representative diverging trajectory. A trajectory converging to a stable pseudoequilibrium, existing for B = −0.51, is depicted in figure 3b. The parameters are chosen so that the postural sway is overdamped. Reduction of this system to the first order will remove the possibility of having a stable pseudoequilibrium; the stable state may be of physiological importance.
4. Effect of time delay in the PD control
Time delay enters the model through the delay present in the PD controller. We assume that the time delay coming from the position and velocity signals are equal, that is, τ_{1} = τ_{2} = τ. This simplifying assumption is justified because in reality, we encounter distributed delays and there is no evidence that the delay coming from the position signal is significantly longer or shorter than the delay of the velocity signal.
To reduce the size of parameter space, for all our numerical experiments, we assume θ_{0} = 1. This assumption can be made without loss of generality as shown in appendix A.1.
4.1 Characteristic equation of the delay system
In order to understand the dynamics of our model, we have to understand the effects of time delay on the dynamics of the control system (2.9). We can rewrite equation (2.9) as a secondorder linear retarded differential difference equation. We have 4.1
Let , where is an arbitrary but nonzero constant. Then, equation (4.1) becomes 4.2
Equation (4.2) is the characteristic equation of the retarded differential difference equation (4.1). The eigenvalue solutions λ of this equation determine the character of the solutions of equation (4.1).
4.2. Dynamics for B = 0 and τ small, and for C = 0 and τ small
Let us consider system dynamics when C = 0 and τ is small. Then the characteristic equation (4.2) becomes 4.3
For τ = 0, we obtain a system that conserves energy, and the eigenvalues of the characteristic equation of the ODE that governs system dynamics are . For τ sufficiently small, we calculate the dominant eigenvalues of the characteristic equation (4.3) by expanding it in τ. We then have and after rearrangement we get 4.4
Approximating the RHS of equation (4.4) to zero the roots of the quadratic equation (4.4) are
By assumption B < 0, (A + B) < 0, and hence if τ > 0 and is sufficiently small there are at least two roots of equation (4.4) with positive real parts. This implies further that in the current case the trajectories of systems (2.8) and (2.9) will diverge to infinity for τ sufficiently small. A representative trajectory for this scenario is depicted in figure 4 with the parameters set to A = 0.5, B = −0.6, C = 0 and time delay τ = 0.1.
We can similarly derive the approximate characteristic equation when B = 0 and C < 0, and argue that trajectories will diverge to infinity for τ sufficiently small.
To conclude

— if the derivative control is not active (C = 0: the damping term not present) or if the proportional control is not active (B = 0: no additional stiffness), and the time delay is small enough, the switched systems (2.8) and (2.9) is not controllable, and all the trajectories diverge to ± ∞;

— numerics suggests that the above scenarios persist also for large values of τ.
4.3. Switched PD control: B and C nonzero
To create a stable limit cycle, it is necessary to move all the eigenvalues of the characteristic equation (4.2) to the left halfplane of the complex plane. Assume that both τ and C are small, that is both τ and C are 𝒪(ɛ). To leading order in τ and C the characteristic equation (4.2) is 4.5
The two roots of the quadratic equation (4.5) are 4.6
For C and τ sufficiently small these roots are complex conjugate, and it is the sign of the real part, that is the sign of C − Bτ, which determines whether all the roots of the characteristic equation (4.2) lie in the left halfplane of the complex plane.
Considering the previous numerical example if we set C =−0.07, we expect that a system trajectory will converge to an attractor since then C − Bτ < 0, and the dominant eigenvalues lie in the left halfplane of the complex plane. This is indeed the case as depicted in figure 5.
In figure 6, we depict a diverging trajectory of equations (2.8) and (2.9) for A = 0.5, B =−0.6, C =−0.05 and the time delay τ = 0.1. Hence the condition C − Bτ < 0 is violated and there exist eigenvalues of the characteristic equation (4.5) with a positive real part.
We should make a comparison here between our switched model and the delay differential equation that governs system dynamics outside of the dead zone. Assuming there is no dead zone (ideal sensing of the proprioceptive system) then the small scale stable oscillations present in the switched system correspond to the stable equilibrium states of the delay system. If, on the other hand, there are no attractors present in the switched system then the equilibrium of the delay system must be unstable, and the asymptotic dynamics of the switched and delay systems can be considered equivalent (the trajectories diverge to ± ∞).
4.4. Multiple stable oscillations
Having established the existence of stable periodic oscillations in our model, we wish to examine the region of parameter space where these oscillations persist, and also, what are the parameter values that are critical for the onset of different asymptotic dynamics. Initially, we set the parameters to these values for which we found stable limit cycles in §4.3; we set τ = τ_{1} = τ_{2} = 0.1, A = 0.5, B =−0.6, and we vary C. In figure 7a–h, we are depicting representative limit cycles and corresponding time series for C ∈ [−10, −0.1]. We note that there is little effect of the variation of C on the amplitude and the period of the oscillations for C ∈ [−10, −0.2] (figure 7c–h). When we vary C between −0.2 and −0.1, we observe that the value of at the switching between F_{in} and F_{out} increases from around 0.02 to 0.03 and the period from 0.5 to 0.7 (cf. figure 7c,d with figure 7a,b). To sum up, the effect of the variation of C on the existing periodic orbits is more pronounced for small values of C and there is little effect of the variation of C on the periodic orbits for −10 < C < −0.2.
To further investigate the system numerically, we obtained oneparameter orbit diagrams on which we depict the variation of a periodic point on the attractor versus parameter −C for distinct values of the time delay τ (figure 8). We found the coexistence of two families of attractors. Namely, we observe stable oscillations, e.g. figure 7a–h that coexist with stable oscillations of different period and amplitude. Note two curves obtained for τ = 0.25 one in a close proximity of x_{1} = 1 and the other further away from x_{1} = 1.
We found the coexistence of these two families of attractors also for τ = 0.5 and τ = 0.75 (not depicted in the figure). We conjecture that these two families of attractors are born in the limit as τ → 0. Their presence is due to the switched nature of the system. Representative examples of two coexisting attractors, each corresponding to one family of periodic orbits, found for τ = 0.25 are depicted in figure 9. The amplitude θ, say Amp(θ), of the limit cycle in figure 9c is Amp(θ) ≈ 0.02. This is about 50 times larger than Amp(θ) of the limit cycle in figure 9a. The period of the limit cycle in figure 9c is about 1.2 which is approximately six times longer than the period of the limit cycle in figure 9a; cf. the time series in figure 9d,b.
4.5 Homoclinic bifurcations and bistability
Another important feature of our system is the birth of stable symmetric orbits with long period through socalled homoclinic bifurcation. This phase space transition takes place under increasing values of τ when the switching between F_{in} and F_{out} occurs at points on Σ_{1} where the eigenvectors corresponding to the saddle point cross Σ_{1}.
Consider oscillations depicted in figure 10. Note that within the dead zone the orbit from figure 10a existing for A = 0.5, B =−0.6, C =−1 and τ = 1.1135 is very close to the stable and unstable manifolds of the saddle point (note the dashed lines in the figure superimposed on the orbit). For τ = 1.114 this orbit no longer exists, instead a stable symmetric orbit born in the homoclinic bifurcation is present (figure 10c). A homoclinic orbit existing at the bifurcation is formed from the stable and unstable manifolds of the saddle point. Within the dead zone, the shape of the orbit is given by the straight lines lying along the eigendirections of the saddle point and outside of it by the arches joining these straight lines. In other words, it is the stable and unstable manifolds that form a homoclinic orbit, and any trajectory starting its evolution on the orbit stays on it and reaches an unstable equilibrium of the saddle type in infinite time. This in turn implies that the periods of the orbits ‘before’ and ‘after’ the bifurcation are long—see time series in figure 10b,d, respectively. The straight lines emanating from the origin refer to the stable and unstable eigendirections of the saddle point. We believe that it is the homoclinic bifurcations that allow for the noiseinduced switchings between two coexisting attractors reported in Eurich & Milton [20]. We note that sufficiently close to the bifurcation there exist two stable asymmetric attractors as depicted in figure 11a. If we add white noise to the system then a typical trajectory will evolve around these two attractors as depicted in figure 11b. In the latter case, the switched systems (2.8) and (2.9) becomes a stochastic switched system with time delay, and its evolution within the dead zone is governed by 4.7and outside of the dead zone by 4.8where σ is the intensity of the applied white noise ζ (t). The parameter values are set to A = 0.5, B =−0.6, C = −1, σ = 0.1 and τ = 1.112. Both figures were generated using the firstorder Euler method with step size Δt = 0.001. The numerical integration scheme used to produce figure 11b is described in the appendix.
5. Modifications to the model
The investigations described here serve as an initial step to gain insight into the dynamics of upright balance. There is a number of modifications that could be considered to gain more insight into the problem of balancing. In the first instance, we could alter the model with the dead zone and time delay to a model with time delay and delay in the switching function. This modification reflects our assumption on the presence of time delay in the neural transmission and muscle activation when switching the PD control. In this paper, we consider the simplest possible model and that is why the dead zone is fixed. This can also be seen as modelling the time delay in the neural transmission only. In some papers, it is also argued that there is a time delay in the acceleration terms (see [35]). It could be also interesting to compare different control strategies and stable states that these can produce. Finally, we would like to develop the model to include the mechanism(s) of muscle tiring. This we believe to be a crucial element in the loss of balance in elderly people.
There are also other important issues. Namely, the dynamics of a multiple link inverted pendulum model could be more realistic. These modifications will be pursued in collaboration with experimentalists who work on human balance.
At this point, we would like to comment on the correspondence between the dimensional and nondimensional quantities of our model system. For example, if we take physiologically feasible values, similar to those used in Asai et al. [25], and set m = 60 kg, h = 1 m, g = 9.81 m s^{−2}, J = 60 kg m^{2}, time delay Δ_{1} = Δ_{2} = 0.1 s, the width of the dead zone θ = 0.02 rad, and the control coefficients K_{p} = 720 Nm rad^{−1}, K_{D} = 60 Nms rad^{−1} our nondimensional quantities are equal to A = 0.0001, B = 0.00012, C = 0.0031 and the time delay τ = 31. Our numerical simulations show that for these parameter values the system exhibits stable pseudoequilibria. We should note that these values are in a different regime from the one we explored numerically. The reason for our numerical exploration is that we want to discover all, or as much as possible, of qualitatively different dynamics that can be found in our model system. Having found multistability, small scale periodic oscillations, and a homoclinic bifurcation, all of which can be potentially important for human balance, it now remains to validate which of these scenarios can, indeed, be observed in experiments. This implies further that it may be necessary to use additional tools from dynamical systems analysis, for instance a numerical continuation of bifurcation curves, to obtain an exhaustive picture of the dynamics of the model system for parameter values which are physiologically feasible. However, this last point must be informed by experiments.
6. Conclusions
In this paper, we introduce a model of human balance during quiet standing following the idea that a human body can be modelled by a singlelink inverted pendulum, and balance is achieved using linear feedback control with time delay in the proportional and derivative error signals. We assume a threshold value of the angle of the sway below which the motorneural control system cannot detect any sway motion. We obtain a planar switched (hybrid) model. We find that to achieve stabilization, which is seen as ‘small’ oscillations about an upright equilibrium, it is necessary that both the proportional and derivative signals of the control system are used. These stable oscillations seem to represent closer to observation stable state for upright standing than the equilibrium points [23]. Therefore, we study the effects of parameter variations on their existence.
Our parameter study leads to the detection of a multiple number of stable oscillatory states existing for the same parameter values and for a wide range of the control parameters corresponding to the derivative term of the PD controller. We also find a homoclinic bifurcation that gives birth to a stable symmetric orbit with a long period. In particular, we show, using a numerical experiment, that close to a homoclinic bifurcation white noise introduced additively may result in the system switching between the two regions where symmetric stable solutions exist in the deterministic switched system leading to an apparent bistability; in other words, the switched system with added noise evolves for some time in the neighbourhood of each one of the two stable asymmetric limit cycles (present in the deterministic system) by switching between their regions of existence. This scenario can explain switchings between a pair of stable asymmetric attractors observed in the firstorder model in Eurich & Milton [20], which in turn was used to explain different scaling patterns that could be detected in human postural sway data.
Obviously, an important step in the investigations of human balance would be verification of the stable oscillatory states found in our model system. In particular, it would be interesting to see if, indeed, we can observe small scale oscillations close to an equilibrium against the larger scale stable oscillations and how are these linked with the physiology of human subjects.
Acknowledgements
Research partially funded by Engineering and Physical Sciences Research Council grant EP/E050441/1 (CICADA: Centre for Interdisciplinary Computational and Dynamical Analysis), the University of Manchester and Manchester Metropolitan University.
Appendix A
A.1. Model scaling
Consider the system for k > 0 given by A 1which is our standard control model with dead zone θ ≤ k. We aim to show that this is equivalent (after scaling) to a model with dead zone θ ≤ 1 of the form A 2
Thus having a smaller dead zone (smaller k) corresponds, in systems rescaled so that the dead zone is constant, to having a larger relative dissipation but a longer delay time; i.e. the effect is not obvious (larger relative dissipation implies greater stability, longer delay implies less stability). Note that the time t in the second equation is not the same as the time t in the first equation.
To see this, we will start from equation (A 2) and derive (A 1). First, let ψ(t) = kθ(t), so θ ≤1 becomes ψ ≤k and dψ/dt = k(dθ/dt) and (A 2) becomes A 3
Now introduce a new time variable u = kt, and let ϕ(u) = ψ(t), and y(u) = x(t). Note that A 4and similarly for y(u). Moreover and so using equation (A 3) A 5with a simpler equation if ϕ ≤k. Finally, using the definitions of ϕ and y, and the shift rule (A 4) gives (with primes denoting differentiation with respect to u) A 6which is equation (A 1) after the identifications ϕ → θ, y → x and u → t.
A.2. Euler's method for the switched system with the time delay and white noise
Switched systems (4.7) and (4.8) is a stochastic system with switching and time delay. The angle θ and the angular velocity are now random variables. The presence of the switching function implies that depending on the value of the random variable θ system equation is governed by a stochastic differential equation, or a stochastic delay differential equation. We switch between these two systems when the random variable θ is greater or smaller than the threshold value θ_{0}. It has been shown in the study of Kloeden & Platen [36] that stochastic delay differential equations A 7where τ is the time delay, ζ(t) is Gaussian white noise with intensity σ, can be approximated by A 8for h sufficiently small; τ is the time delay, h is the step size, k = τ/h and W_{n} is the standard Wiener process. Similarly, a stochastic differential equation A 9can be approximated by a discrete system A 10
Above numerical scheme was used to generate the trajectory of the stochastic switched systems (4.7) and (4.8). Switching occurs at the integration step n when θ_{n} − θ_{0} changes sign. The standard Wiener process is approximated numerically at each step t_{n} by a function that generates pseudorandom numbers with expected value μ = E[X] = 0 and standard deviation , where X is a random variable.
 Received April 7, 2011.
 Accepted June 1, 2011.
 This journal is © 2011 The Royal Society