Generation of the Human Biped Stance by a Neural Controller Able to Compensate Neurological Time Delay

Abstract

The development of a physiologically plausible computational model of a neural controller that can realize a human-like biped stance is important for a large number of potential applications, such as assisting device development and designing robotic control systems. In this paper, we develop a computational model of a neural controller that can maintain a musculoskeletal model in a standing position, while incorporating a 120-ms neurological time delay. Unlike previous studies that have used an inverted pendulum model, a musculoskeletal model with seven joints and 70 muscular-tendon actuators is adopted to represent the human anatomy. Our proposed neural controller is composed of both feed-forward and feedback controls. The feed-forward control corresponds to the constant activation input necessary for the musculoskeletal model to maintain a standing posture. This compensates for gravity and regulates stiffness. The developed neural controller model can replicate two salient features of the human biped stance: (1) physiologically plausible muscle activations for quiet standing; and (2) selection of a low active stiffness for low energy consumption.

About the Authors:

Ping Jiang

* E-mail: [email protected]

Affiliation: Department of Precision Engineering, School of Engineering, The University of Tokyo, Tokyo, Japan

Ryosuke Chiba

Affiliation: Research Center for Brain Function and Medical Engineering, Asahikawa Medical University, Asahikawa, Japan

Kaoru Takakusaki

Affiliation: Research Center for Brain Function and Medical Engineering, Asahikawa Medical University, Asahikawa, Japan

Jun Ota

Affiliation: Research into Artifacts, Center for Engineering (RACE), The University of Tokyo, Kashiwa, Japan

Introduction

Stance postural control (SPC), which allows individuals to maintain an upright stance, is one of the most important and basic requirements for a comfortable life [1]. The physiologically inspired neural controller (NC) model for SPC has potential applications in a variety of fields. For example, the NC model could be used to assist device development and design robotic control systems. To realize such applications, the model must successfully replicate the functionality of the human neural controller (e.g., the ability to maintain a posture via muscle coordination, while compensating for the neurological time delay). The realized model must also replicate the salient features (e.g., physiologically plausible muscle activations) of the human biped stance in computer simulations.

System identification is one approach to developing an NC model based on experimentally measured data [2–13]. Applications of this method aim to develop an NC model that can simulate human movement that agrees with experimental data. To achieve this, the variables in the NC model are typically identified by minimizing the discrepancies between the simulation results and the corresponding experimental data. Using system identification, researchers have studied the influence on SPC of joint and muscle stiffness [3–6], sensory information [7–10, 12, 13], feedback (fb) gains [14–17], and the muscle torque generation process [11]. Several of the studies cited above [2–19] have proposed potential NC models, such as the stiffness control [3], sensory fb control [17], and disturbance estimation and compensation (DEC) [20] models. An optimal postural control model was proposed by Qu et al. [21] to identify the change of balance control mechanism due to aging. This model was further adopted to investigate the role of passive and active toque for maintaining postural balance [22]. A sliding model [23] developed by Zhang et al., enables the passive and active joint stiffness to be estimated. He also proposed a novel identification approach to estimate passive and active joint stiffness without any additional perturbations [24]. However, the efficacy of these models is solely determined by their ability to reproduce experimental data, and it is difficult to confirm whether they truly reflect the mechanism behind human postural control. In addition, it is unclear whether these models can be employed to control a more complex human body model (e.g., a musculoskeletal model with many muscles).

Forward modeling is another effective approach for developing NC models. Different from system identification, forward modeling does not use any experimental data as input. Instead, the variables in the NC model are typically designed by optimizing an assumed performance criterion (e.g., the joint angle displacement) [25]. It is possible to validate an NC model developed using forward modeling by examining whether the simulated results reflect the salient features of human movement, and comparing those results with experimental data. To date, three types of NC model that generate the biped stance have been proposed. The first is based on conventional feed-forward (ff) control in conjunction with fb control [26–30], whereby the neural controller estimates the state of the human body internal model [31] (the joint angle and angular velocity) using an optimal estimator [32]. Torque is added directly using a proportional-differential (PD) controller that utilizes the estimated joint kinematics information [26]. For this model, van der Kooij et al. [26] have reported that the controller can compensate for a delay of approximately 80 ms, and that it can simulate the salient features of sensation integration observed in experiments. The second model type is based on fb control only [18, 33–35]. For an inverted pendulum model and assuming fb control only, Masani et al. [18] reported that a PD controller can compensate for a delay of approximately 185 ms, provided the gain is sufficiently high. The designed model can reproduce the position and velocity of the center of mass (CoM), along with the joint torque data. The third model type employs intermittent control [36–42], whereby the PD fb controller is activated intermittently based on the status of the joint angle.

To generate a human-like biped stance that reflects the appropriate salient features, an NC model should: (1) successfully maintain a musculoskeletal model with muscles in a standing posture; and (2) compensate for the neurological time delay. Many previous studies [4, 5, 7–19, 26, 28–30, 33–38, 40–44] have utilized simplified musculoskeletal models with fewer joints and torque actuators in place of muscle actuators. They have also neglected the human musculoskeletal anatomy (e.g., the CoM of each body segment, muscle attachment points, and skeletal inertial properties). These previous studies have employed one [4, 5, 7–9, 11, 13–16, 18, 30, 36–38], two [10, 12, 17, 19, 29, 40–42, 44], or three degrees of freedom (DoF) [26, 28, 43] inverted pendulum models with no muscles, a one-DoF inverted pendulum model with one [33] or three muscles [34], or a three-DoF inverted pendulum model with nine muscles [35]. However, Günther et al. [43] and Hsu et al. [45] have suggested that all the joints should contribute to the biped stance. Further, Horlings et al. [46] identified the importance of muscle strength to postural stability. Therefore, a musculoskeletal model that more accurately represents the human anatomy than the above inverted pendulum models should be developed. For example, Clark et al. [47] have used a reflex controller to hold a musculoskeletal model with 92 muscles in a standing position; however, neurological time delays similar to those of humans was not considered, and simulated muscle activations were not evaluated. Further, a DEC model [10, 20, 48] has been successfully implemented in a robot featuring eight muscles, incorporating a realistic time delay. However, a robot with more muscles and more physiologically accurate anatomy is needed to approximate the human anatomy. To date, no NC model that can compensate for the neurological time delay in the case of a more human-like, complex musculoskeletal model in a standing posture has been developed.

In this study, we employ a forward modeling method to develop a computational model of a neural controller that can simulate a human-like biped stance and yield physiologically plausible features for this stance.

Specifically, we aim to develop an NC model (red block in Fig 1) capable of: (1) realizing human SPC (Fig 1); and (2) reflecting the features that generate muscle activations within physiologically plausible ranges for the human biped stance. In realizing this postural control, we derive a musculoskeletal model (green block in Fig 1) with 70 muscles and seven joints in a standing position under the influence of a 120-ms neurological time delay (flesh-colored blocks in Fig 1).

[Figure omitted. See PDF.]

Fig 1. Human stance postural control.

Red block: neural controller to actuate muscles during quiet standing. Green block: musculoskeletal model. Flesh-colored blocks: 120-ms neurological time delay. The delay is comprised of a 40-ms feedback delay, 40-ms transmission delay, and (at most) 40-ms activation dynamics delay. The activation dynamics delay corresponds to a 20-ms activation time delay during the muscle activation process and a 40-ms deactivation time delay during the muscle deactivation process.

https://doi.org/10.1371/journal.pone.0163212.g001

Methods

We developed an NC model (red area in human SPC (Fig 2)) to simulate the biped stance, considering a 120-ms neurological time delay (flesh-colored blocks in Fig 2). A musculoskeletal model was created to approximate the human musculoskeletal system anatomy (green block in Fig 2). The NC model consists of ff and fb controls, and was designed to actuate the muscles in order to hold the musculoskeletal model in an upright standing posture. The ff control output, u_ff, is a set of predetermined constant muscle activations that are used to maintain an objective posture. Details are given in the “Feed-forward control” section. Feedback control is modeled as a PD control incorporating fb information on the muscle length and lengthening velocity.

[Figure omitted. See PDF.]

Fig 2. Proposed NC model for human stance postural control.

u, u_ff, and u_fb are the total, ff, and fb controls, respectively; τ_fb, τ_trans, and τ_act are the feedback, transmission, and activation dynamics delay, respectively. L^MT and are the length and lengthening velocity of the muscular-tendon actuator, respectively. L^M and are the length and lengthening velocity of the muscle fiber, respectively. a indicates the muscle activation, k_p and k_d are the fb gains, and and are the muscular-tendon actuator length and lengthening velocity for the objective posture, respectively. Each of the symbols listed above represents a 70-dimensional vector.

https://doi.org/10.1371/journal.pone.0163212.g002

The large number of variables (53) in the NC model poses a challenge. The variable design process was divided into two stages. The ff control variables were designed in the first stage, after which various u_ff candidates were collected. Among them, several u_ff were selected and used to design the fb control variables via an optimization process.

Note that the selected u_ff must, in conjunction with fb control, simulate muscle activations that reflect physiologically plausible features of the biped stance in order to be considered valid. Therefore, we compared the simulated muscle activations with experimental data to verify which u_ff could simulate muscle activations that were within the physiologically plausible range.

Musculoskeletal model

We created a musculoskeletal model (Fig 3a) in OpenSim 3.3 (SimTK.org) [49], which is an open-source biomechanical simulator. The model is based on the Gait2392 Model [50] and ToyDropLanding Model [51], both of which are provided by the OpenSim 3.3 software, and is designed to simulate the biped stance. The created model is comprised of eight body segments, 70 muscular-tendon actuators, and seven joints, with seven DoFs, as shown in Fig 3a. In this study, we primarily focused on the motion of the lower extremities; therefore, the upper extremities (head, trunk, and arms) are modeled as a single body segment connected to the pelvis by the lumbar joint. The lower extremity section is composed of the femur, shank, and foot, which are connected by the hip, knee, and ankle joints, respectively. As we primarily studied the stance motion in the sagittal plane, all joints are modeled as 1-DoF rotational joints for flexion and extension. The 70 muscular-tendon actuators in the ToyDropLanding Model [51], which are assumed to be important for quiet standing, are employed. These muscles are modeled as 70 Millard 2012 Equilibrium muscular-tendon actuators [52]. Each actuator incorporates a compliant tendon because compliant tendons such as the Achille’s tendon plays an important role in stance postural control [53]. To apply an accurate human musculoskeletal anatomy to the model, the kinematic (e.g., the body segment lengths, joint positions, and muscle attachment points) and dynamic (e.g., the body segment inertial properties, muscle isometric forces, and optimal muscle lengths) information of the Gait2392 Model is used. Note that the Gait2392 Model has rather accurate human anatomy data and has been widely employed in gait analysis [54–57].

[Figure omitted. See PDF.]

Fig 3. Musculoskeletal model.

(a) Musculoskeletal model. HAT: head, arms, and torso. The OpenSim musculoskeletal model consists of seven DoFs and 70 muscular-tendon actuators. q₁–q₇ are the seven joint coordinates, each of which has one DoF for joint flexion and extension. (q₁–q₇: Lumbar, right hip, right knee, right ankle, left hip, left knee, and left ankle flexion-extension, respectively.) (b) CoM position and projection at base of support (BoS). The green-and-black ball is the CoM position, and the green ball is the vertical CoM projection at the BoS. (c) Foot-ground contact. The foot-ground contact is represented by four compliant spheres (heel, toe, and two metatarsals) on each foot, with Coulomb friction.

https://doi.org/10.1371/journal.pone.0163212.g003

The foot-ground contact is modeled as contact between four compliant balls (heel, toe, and two metatarsals in Fig 3c) with a plane (Ground in Fig 3a). The Hunt–Crossley Contact Model [58] defined in the ToyDropLanding Model is employed to calculate the contact force. Specifically, the contact force (ground reaction force) is calculated when the feet begin to penetrate the ground, according to(1) (2)where R_x and R_y are the frictional force and the vertical ground reaction force, respectively, and μ(v) is the frictional coefficient function related to the relative velocity v between the ground and foot-ground contact point. Note that μ(v) is defined in terms of μ_s, μ_d, μ_v, and v_t, which are the static, dynamic, and viscous friction coefficients and the speed at which the static friction reaches its peak value, respectively. The equation for μ(v) was defined by Sherman et al. [59]. Further, E is the contact stiffness, h is the penetration height (the foot soft-tissue deformation displacement), and b is the dissipation coefficient. More details on implementing the Hunt–Crossley model in OpenSim are given by Sherman et al. [59].

In this study, μ_s, μ_d, and μ_v are set to 0.9, 0.9, and 0.6, respectively; v_t is set to 0.1 (m/s), b is set to 0.5 (m/s)⁻¹, and E is set to 10⁸ (Nm^−1.5). The contact-model variable configuration is determined to avoid unrealistic contact and maintain a reasonable computation speed [59, 60].

Neurological time delay

In realizing the biped stance, we considered a maximum neurological time delay of 120 ms, which includes a 40-ms feedback delay τ_fb, 40-ms transmission delay τ_trans, and maximum 40-ms activation dynamics delay τ_act. Here, τ_trans is the delay affecting the neural controller’s transmission of the neural controls to the muscle-motion neurons, whereas τ_fb is the delay affecting the sensory receptors’ receipt of the sensory information. Both delays are modeled using a pure time delay of 40 ms, in accordance with Masani et al. [18]. The third delay, τ_act, affects the muscle activation dynamics, which is the process through which the muscles generate force following receipt of the control signals. τ_act is 20 or 40 ms, depending on which motion process influences the neurons. The muscle activation dynamics can be divided into activation and deactivation processes, as shown in the breakdown of the activation dynamics delay given in Fig 1. In the activation process, the muscle-motion neurons fire to a higher activation level, as commanded by the neural controller. In contrast, for the deactivation process, the muscle-motion neurons deactivate to a lower activation level. In accordance with Eq (4), τ_act is set to be equivalent to a 20-ms activation time delay t_act during the activation process and a 40-ms deactivation time delay t_dact during the deactivation process. Further, t_act and t_dact are set to 20 and 40 ms, respectively, in accordance with Zajac [61], Winters [62], and Jacobs [63]. The muscle activation dynamics are modeled as a first-order differential equation in OpenSim, with(3) (4)where u_i is the total output from the NC model to the ith muscle and a_i represents the muscle activation of the ith muscle (i = 1, 2, 3, …, 70). δ is the delay coefficient.

In this paper, the maximum neurological time delay of 120 ms is referred to as the “120-ms neurological time.”

Neural controller

Both the ff and fb controls are assumed to be necessary for human SPC. It is widely acknowledged that fb sensory information plays a vital role in postural control [20]; therefore, sensory fb control is indispensable during quiet standing. In addition to fb control, previous physiological studies have indicated the possible existence of ff control [64, 65]. For example, Fitzpatrick et al. [64] suggested that sensory fb control is important, but not sufficient, to stabilize posture in the case of perturbations. Further, Gatev et al. [65] reported that muscle contraction occurs prior to stabilization of the CoM position, implying that ff control may be employed to predict the future CoM position and achieve postural stabilization.

In accordance with the findings of these previous physiological studies [20, 64, 65], our proposed neural controller incorporates ff and fb controls, with outputs labeled u_ff and u_fb, respectively. Thus, the total output of the NC model, u, which is a combination of u_ff and u_fb, excites the muscle motion neurons to drive the skeletal system, such that(5)where u_ff and u_fb are the ff and fb output to 70 muscles. Thus, u, u_ff, and u_fb are 70-dimensional vectors.

Feed-forward control.

Unlike the ff control proposed in previous studies [26–30], in this study, u_ff is the set of activations needed to keep the musculoskeletal model standing in the defined objective posture, i.e., the upright standing posture shown in Fig 3a.

We determine this upright standing posture by tuning q₁, q₂, q₄, q₅, q₇ (Fig 3a) using the Covariance Matrix Adaptation Evolution Strategy (CMA-ES) optimizer available in OpenSim 3.3 (SimTK.org) [49]. Details of the optimizer are described in the “Variable design” section. Note that q₃ and q₆ are fixed to 1°, because we assume that the knees are mildly flexed in the objective posture. The initial solutions for q₁, q₂, q₄, q₅, and q₇ are set to q₁ = −10, q₂ = −5, q₄ = 0, q₅ = −5, and q₇ = 0 for the CMA-ES optimizer. The optimizer generates posture candidates and evaluates their objective function J_pos for each iteration. The optimization terminates when the convergence tolerance reaches the default value defined in the optimizer.

The objective function J_pos that is minimized by the CMA-ES optimizer can be expressed as(6) (7)where CoM_0,x is the vertical projection of the body CoM at the base of support (BoS) (green ball in Fig 3b), CoP_0,x is the center of pressure (CoP) position under each candidate posture. BoS_x is the BoS center (yellow ball in Fig 3c).

J_pos is used to evaluate whether the posture satisfies the static stability conditions. The well-known condition for standing stability in static situations is that the vertical projection of the CoM should be within the BoS [66]. The CoP should also be within the BoS, to ensure that the ground reaction force can be transmitted from the ground to the feet [67]. Hence, both “(CoM_0,x − BoS_x)²” and “(CoP_0,x − BoS_x)²” are used to adjust the objective posture, so as to allow both the vertical projection of the CoM and the CoP to be close to the BoS center. This adjustment ensures that the objective posture satisfies the stability condition. Note that the CoP position is not fixed. Rather, the objective-posture CoP position is simply the initial CoP position, which changes with body movement during quiet standing.

There are many posture candidates that satisfy the static condition. To obtain a posture that can be maintained with minimal torque, the following term is also incorporated:(8)where Ω_n is the joint torque necessary for the nth joint to maintain this objective posture. We use w_static and w_Ω to represent the weights of the corresponding terms: w_static is set to 10000, so as to reject any postures that do not satisfy the static stability condition, and w_Ω is empirically set to 0.1. The specific values of the coordinates for the obtained objective posture are listed in Table 1.

[Figure omitted. See PDF.]

Table 1. Coordinate values for objective posture.

https://doi.org/10.1371/journal.pone.0163212.t001

Once the objective posture has been determined, a specific constant muscle activation, c_i, which is independent of the fb information, is sent directly to the ith muscle. Thus,(9)where c_i is the ith element of u_ff, and ensures that the ith muscle maintains balance in the objective posture (i = 1, 2, 3, …, 70).

The u_ff components compensate for gravity and increase the joint stiffness resulting from muscle contraction. Note that human joint stiffness contributes to balance during SPC. In this paper, two sources of joint stiffness are considered: the passive mechanical source (tendons), which is modeled as a spring-like tendon component in a muscular-tendon actuator [52], and the active source, called active stiffness, resulting from u_ff . u_ff is similar, but not identical, to muscle co-contraction, which stimulates the activation of antagonist muscles around a joint so as to fixate that joint. Muscle co-contraction increases the joint stiffness only. The net torque generated by muscle co-contraction is zero; therefore, it does not contribute to gravitational compensation. However, the net torque generated via u_ff compensates for the gravitational torque.

In a previous muscle co-contraction study [68], the active stiffness level was typically quantified based on the muscle activation level. Similarly, in this paper, the square norm of u_ff, ||u_ff|| (), is used to quantify the active stiffness level.

Feedback control.

Sensory information is critical for realizing the biped stance. Human SPC is realized via an fb mechanism that actuates muscles to generate appropriate corrective torques, so as to counter the destabilizing torque due to gravity [8]. The corrective torques are generated through the integration of multisensory inputs, including the visual, vestibular, and proprioceptive somatosensory inputs.

Proprioceptive somatosensory input is assumed to make the largest contribution to the SPC during quiet standing. Winter et al. reported that vestibular input may not make a significant contribution to the control of the upright stance [3], and Sousa et al. argued that, as in normal conditions, proprioceptive information has more relevance than other sensory sources. Hence, we employed the proprioceptive sensory input as fb information only, despite the existence of multisensory inputs [69].

The fb control is approximated by a PD controller (Eq (10)). This controller functions based on proprioceptive sensory information concerning the muscular-tendon length and lengthening velocity. Thus,(10)where is the maximum muscle lengthening velocity, , and are the delayed muscular-tendon length and lengthening velocity, respectively, is the muscular-tendon length reference, which is the length for the objective posture, and is the lengthening velocity reference. Further, k_p,i and k_d,i are the PD gains of the ith muscle (i = 1, 2, 3, …, 70). is set to 0 m/s in order to approach a stable stance. Note that this controller is different from stretch reflex control which adopts muscle fiber length as feedback information. Instead, muscular-tendon length, sum of muscle length and tendon length, is adopted as feedback information. The muscle fiber length moves in the opposite direction of human body. For example, when musculoskeletal model moves forward, soleus muscle fiber length contract to generate the bias of the tendon and torque to move the musculoskeletal model backward. This paradoxical muscle movement characteristics coincides with that reported in previous study [53].

Variables in neural controller.

In total, the neural controller has 210 variables, including a 70-dimensional u_ff (70 constant c values in Eq (9)) and 140 PD gains (70 proportional k_p and 70 derivative k_d gains, Eq (10)). In this paper, we mainly focus on the anterior-posterior body movement in the sagittal plane during quiet standing (joint displacements are symmetrical: q₂ = q₅, q₃ = q₆, and q₄ = q₇). Further, we assume that muscles can be grouped based on their function on the joints. Hence, the following assumption-based simplifications are employed:

1. The control laws for the ff and fb control on the muscles are taken to be symmetrical (e.g., the left and right soleus have the same c and the same k_p and k_d);

2. Muscles located around the same joint and having similar functions on the joint are assumed to have the same PD gains (e.g., the pectineus and psoas are taken to have the same PD gains, because they are both positioned around the hip and have hip flexion functionality).

Specifically, we divide the 70 muscular-tendon actuators into nine groups in accordance with the muscle functions and attachment point positions, obtaining the lumbar extensor, lumbar flexor, hip extensor, hip flexor, knee extensor, knee flexor, ankle extensor, ankle flexor, and biarticular muscle groups. The biarticular muscle group is introduced because the biarticular muscles generate torque on two joints and may have different effects on those joints. The k_p and k_d of each muscle are not unique, and members of the same muscle group have the same values for these fb gains. Thus, the gain of each muscle is dependent on its muscle group. For example, the pectineus and psoas are assigned to the hip extensor group, and both have fb gains of [k_{p_l_ex}, k_{d_l_ex}]. The following assignments are made:(11)where “group” indicates the muscle group.

The muscle gains are grouped for muscles in the same group, and each member of a given group has the same PD gains but different c. Note that the u_ff inputs are not grouped, because u_ff is added to the muscles directly without any fb information. It is difficult to balance the net torque generated by the extensors and flexors, as each muscle has unique muscle properties. Muscle elongation can actually differ among muscles of the same group because of their different moment arm. We assume that moment arm is approximately constant during quiet standing. In general, the muscle fiber length and muscle moment arm are positively correlated [70]. Hence, in Eq (10), we use to normalize muscular-tendon length feedback information to reduce the effect of moment arm. We also normalize the velocity information by to obtain dimensionless derivative gains. This normalization allows muscles in the same group to be controlled by the same gains. As a result of the normalization, [k_p,i, k_d,i] is dimensionless.

Thus, the number of essential variables for our proposed neural controller design has been reduced from 210 to 53, including a 35-dimensional u_ff and 18 PD gains (as listed in Eq (11)).

Variable design

As explained above, 53 variables are required to model the controller. This is very challenging as regards determining a suitable solution for many variables to keep a musculoskeletal model standing under the influence of the 120-ms neurological time delay. As u_ff is constant and independent of the neurological time delay, the various u_ff candidates (the 35 constant c values in Eq 9) are first determined (Fig 4, red area); τ_trans and τ_fb are both 0 ms for these variables. Details of the calculation are given in the “u_ff calculation” section. Among the various calculated u_ff candidates, several are selected based on the value of ||u_ff||. Subsequently, for each selected u_ff, the PD gains are optimized (Fig 4, blue area) so as to maintain the standing posture under the influence of the 120-ms neurological time delay.

[Figure omitted. See PDF.]

Fig 4. Variable design algorithm.

The red and blue areas indicate the u_ff calculation algorithm and the PD gain optimization algorithm, respectively.

https://doi.org/10.1371/journal.pone.0163212.g004

u_ff calculation.

This section describes the calculation of the many 35-dimensional u_ff candidates (the c values in Eq (9)) needed to hold the musculoskeletal model in an objective standing posture (Table 1).

To calculate u_ff, which is independent of the neurological time delay, the muscle activations during a stable biped stance in the objective posture must be obtained, with τ_trans and τ_fb both set to 0 ms. The obtained muscle activations define u_ff, and can be used as the u_ff input. A PD controller receives the same feedback information but different gain types in Eqs (10) and (11) is used to generate a biped stance in the objective standing posture. The desired muscle activations are then obtained by calculating the integrated muscle activations over the stable stance period (Eq (15)). These are used as the u_ff input.

Note that the PD controller above is only used to calculate u_ff, and is not related to the PD controller that copes with feed-forward control to compensate for the 120-ms neurological time delay used in the following “PD gain optimization” section.

To obtain physiologically plausible u_ff, the gains of the flexors and biarticular muscles are scaled with respect to the extensors. We allow all of the flexors and extensors to have identical gains, and scale the gains of the flexor and biarticular muscles by factors of 0.5 and 0.2, respectively, with respect to those of the extensors. The scaling of PD gains is based on the physiological knowledge that the extensors are the dominant mechanism during quiet standing, whereas the biarticular muscles contribute little to this posture [71]. Therefore, we assume that the extensors have larger gains than the flexors and biarticular muscles, and that the biarticular muscles have the lowest gains. Thus, the gains in Eq (11) can be described by P and D as follows:(12) (13) (14)where P and D are the values of the extensor proportional gain and derivative gain, respectively. Note that the scaling of gains is only conducted for the u_ff calculation, and P and D are the only variables in the u_ff calculation.

To calculate various u_ff candidates, as shown in Fig 4 (red area), we conduct an extensive search for the P and D values needed to maintain a standing posture for 60 s, when τ_trans and τ_fb are both 0 ms, using a forward dynamics simulation. Note that τ_act is not neglected, because the muscle dynamics must be incorporated to obtain physiologically plausible muscle functionality. The simulation end time was set to 60 s because this is assumed to be sufficiently long to allow the musculoskeletal model to achieve a stable stance state. The initial posture was set to the same as the objective posture (Table 1), and the initial muscle activations a_i(t)|_{t = 0} were set to zero.

We search for P and D in the range 0.0–2.0 at increments of 0.1. The muscle activations are only recorded if the musculoskeletal model is capable of standing (i.e., if the CoM was higher than 0.4 m). The integration of the muscle activations when the musculoskeletal model achieves a stable standing posture is adopted as u_ff. In other words, the c_i in Eq (9) are calculated according to(15)where t₁–t₂ is the period for which the musculoskeletal model maintains a stable posture (t₁ = 3 s and t₂ = 5 s).

||u_ff|| is then computed to quantify the active stiffness level.

PD gain optimization.

Among the various u_ff candidates obtained as described above, we select several based on ||u_ff||. For each selected u_ff, the PD gain variables are designed based on an optimization procedure (Fig 4, blue area) to cope with the selected u_ff and compensate for the 120-ms neurological time delay. CMA-ES is employed to optimize the 18 PD gain variables that act to hold the musculoskeletal model in a standing posture for 60 s, using a forward dynamics simulation that incorporates the selected u_ff along with the 40-ms τ_trans and τ_fb. The initial posture is updated to the stable posture that can be maintained with the selected u_ff, and the initial a values are updated in accordance with u_ff (a_i(t)|_{t = 0} = u_ff,i). CMA-ES is an evolution algorithm for solving nonlinear black-box optimization problems [72, 73], and has been successfully app lied by Dorn et al. [60] to optimize a complicated controller for gait generation. This algorithm does not calculate the gradient of the objective function but, rather, estimates the covariance matrix. The variables are the population size λ, initial standard deviation σ, and the initial solution and termination criteria. In this study, the optimizer was initialized by setting λ = 20 and σ = 0.005 for fast convergence. The initial solution was generated from a seed. In addition to the default termination criteria, a maximum iteration number of 750 was defined. The simulation was conducted so as to evaluate 20 candidate solutions generated by the CMA-ES in parallel in each iteration. The number of child threads created for parallel execution is equivalent to the number of computer cores.

CMA-ES is used to minimize the objective function J, where(16) (17)with T_simu = 5 s being the simulation end time and T_fail being the time at which the musculoskeletal model begins to collapse or the heel or toe leave the ground. This event is monitored by an event trigger that terminates the simulation when the height of the CoM is less than 0.4 m, or when the heel or toe contact force is 0 N. T_fail is the time at which this event occurred. Note that, if the failure event did not occur until T_simu, then T_fail = T_simu and J_fail = 0. The event trigger decreases the computation time of the optimization during the early iterations. J_stability is used to evaluate the stability of the biped stance and the deviation of the posture from the objective posture. This term is expressed as(18)where q_n(t) is the nth coordinate value (angle displacement of one joint DoF) and q_n(0) is the nth initial coordinate value.

The failure weight, w_fail, was set to 500 000 in order to reject any solutions in which the musculoskeletal model failed to stand. The stable weight, w_stability, was set to 50 so as to rapidly discover solutions that hold the musculoskeletal model in a stable standing position and in a posture close to the objective. The values of 500 000 and 50 were determined based on a study reported by Dorn et al. [60]. In addition, if or is smaller than a minimum value (1.0e⁻⁶), they will be set to 0.

Evaluation of simulated results

To evaluate our simulated results, we compared them with experimental data, including the CoM anterior-posterior (AP) displacement range, CoP AP displacement range, joint correlations and, the muscle activation range for a human biped stance [71, 74].

To investigate which selected u_ff could function with the fb control to generate physiologically plausible muscle activations, we compared the simulated muscle activations against the experimental muscle activation range [71] to determine whether the simulated activations were in the range of the experimental data. Further, the deviation of the simulated muscle activations from the range data was calculated. Experimental muscle activation range and mode value data were reported by Panzer et al. [71], who studied 24 normal subjects, including 12 young subjects (age: 21–57 years; mean age: 38.4 years) and 12 elderly subjects (age: 63–77 years; mean age: 68.1 years). In that study, the subjects stood in an erect posture on a platform with a fre ely chosen foot position. Electromyographic (EMG) data normalized by the maximal voluntary contraction (MVC) of eight muscles were collected.

Further, to investigate the physiologically plausibility of postural sway patterns, we evaluated the following two aspects: 1. CoM AP and CoP displacement trajectory, and 2. multi-joint coordination. We plotted the CoM AP and CoP displacement trajectory. The CoM and CoP AP range are used as indicators to evaluate the physiological plausibility. The simulated CoM and CoP AP displacement range were compared against experimental data reported by Warnica et al. [74], who studied 16 young adults (whose age, height, and body mass (mean (SD)) were 22.6(1.4) years, 173.3(11.1) cm, and 70.7(12.9) kg, respectively). In that study, CoM and CoP displacement data for each adult were measured during quiet standing, and then the CoM and CoP AP displacement range were calculated and analyzed. We calculated the hip-ankle, hip-knee, and knee-ankle angle correlation coefficient, and compared against experimental data [43] to evaluate the multiple-joint coordination.

In addition, to evaluate whether the generated biped stance motion was stable, the CoM AP displacement versus CoM AP velocity was plotted.

Results

Variable design

We conducted an extensive search for u_ff candidates, and obtained a total of 402 u_ff that successfully held the musculoskeletal model in a standing posture. Among them, we selected nine u_ff based on ||u_ff|| (||u_ff|| = 0.06, 0.89, 2.07, 3.44, 4.00, 5.04, 6.00, 7.00, and 8.02), and conducted further PD gain optimization with a neurological time delay of 120 ms. The PD gains were optimized for each u_ff using the CMA-ES algorithm. As a result (Table 2), the gains for all u_ff (except ||u_ff|| = 0.06) were found to successfully maintain the musculoskeletal model in a standing posture (S1 Video).

[Figure omitted. See PDF.]

Table 2. Variable design results.

https://doi.org/10.1371/journal.pone.0163212.t002

Evaluation of simulated results

To investigate which u_ff could function with the PD control to generate physiologically plausible muscle activations, we compared the simulated muscle activations (red dots in Fig 5) against the experimental muscle activation range (gray shaded boxes in Fig 5). The deviations of the simulated activations from the range data are listed in Table 3. The muscle activations when ||u_ff|| = 2.07 were generally within the range of the experimental data. Only the result for rectus abdominus 2 was slightly higher (0.025) than the higher limit of the experimental data range. As for the other ||u_ff||, two muscles fell outside the experimental range for ||u_ff|| = 0.89, three muscles fell outside the experimental range for ||u_ff|| = 4.00, 6.00, 7.00, and 8.02, and four muscles fell outside the experimental range for ||u_ff|| = 3.44, and 5.04.

[Figure omitted. See PDF.]

Fig 5. Simulated muscle activations.

The red dots are simulated activations, the black dots are mode values, and the gray shaded boxes indicate the experimental activation data ranges. The experimental mode values and activation data ranges were reported in [71].

https://doi.org/10.1371/journal.pone.0163212.g005

[Figure omitted. See PDF.]

Table 3. Deviations from value ranges.

https://doi.org/10.1371/journal.pone.0163212.t003

To investigate the physiological plausibility of postural sway patterns during a 60-s simulation, the CoM AP, CoP AP, and CoM height displacements were plotted, as shown in Figs 6 and 7. The CoM oscillates around a stable state value, and the CoP oscillates around the CoM trajectory. From the CoM height displacement, it was determined that all the selected u_ff allowed the musculoskeletal model to maintain a standing position via a periodical height sway. In addition, the CoM AP and CoP AP displacement ranges (maximum displacement minus minimum displacement over 30–60 s, which is the period in which the musculoskeletal model achieved a stable-stance state) were obtained; see Table 4. We compared the ranges with those reported by Warnica et al. [74], who reported the mean±SD of CoP AP and CoM AP to be 20.48±6.97 mm and 17.36±5 mm, respectively. We confirmed that the simulated CoM AP had a smaller range for all ||u_ff||, whereas the CoP AP range for all ||u_ff|| except ||u_ff|| = 6.00 had a smaller range than those of a human.

[Figure omitted. See PDF.]

Fig 6. CoP AP, CoM AP, and CoM height displacement for ||u_ff|| = 0.89, 2.07, 3.44, and 4.00.

The green, red, and blue solid lines are the CoP AP, CoM AP, and CoM height displacements, respectively. The green, red, and blue dotted lines are the CoP position, CoM horizontal position, and CoM height position values for the objective posture, respectively. The positive direction of the “CoP and CoM AP displacement (m)” axis represents the anterior direction.

https://doi.org/10.1371/journal.pone.0163212.g006

[Figure omitted. See PDF.]

Fig 7. CoP AP, CoM AP, and CoM height displacement for ||u_ff|| = 5.04, 6.00, 7.00, and 8.02.

The green, red, and blue solid lines are the CoP AP, CoM AP, and CoM height displacements, respectively. The green, red, and blue dotted lines are the CoP position, CoM horizontal position, and CoM height position value for the objective posture, respectively. The positive direction of the “CoP and CoM AP displacement (m)” axis represents the anterior direction.

https://doi.org/10.1371/journal.pone.0163212.g007

[Figure omitted. See PDF.]

Table 4. CoP and CoM AP displacement ranges.

https://doi.org/10.1371/journal.pone.0163212.t004

Note that the CoM AP, CoP AP, and CoM height values for the objective posture (the dotted line in Fig 6) differed for each ||u_ff||. This is because the objective posture was updated after each u_ff calculation. That is, the same objective posture was employed in each u_ff calculation; however, the final stable posture (which was used as the objective posture during the PD gain optimization) was different, because of the muscle force-generation capability. The muscles worked to realize the objective posture for the musculoskeletal model; however, the maximum isometric force rendered the generation of sufficient force difficult and, as a result, the muscles could only maintain a reasonably similar posture to the objective.

In addition, joint correlation coefficients were calculated as shown in Table 5. when ||u_ff|| < 7.00, hip-ankle angle are negatively correlated, whereas the knee-ankle angle are positively correlated. This result fits the experimental data that hip and ankle angle are negatively correlated (-0.91±0.054); Both in-phase (0.88±0.000) and anti-phase (-0.87±0.054) of correlations between ankle and knee were observed [43]. However, hip-knee correlation are anti-phase when ||u_ff|| = 0.89, 3.44, and 6.00, which differs from the experimental data that hip and knee exhibit a positive correlation (0.89±0.053) [43]. When ||u_ff|| = 7.00 and 8.02, the hip-ankle angle has positive correlation, whereas both the hip-knee and knee-ankle angle are negatively correlated. In this case, musculoskeletal model has high overall joint stiffness and sways like an inverted pendulum. The rotations of the hip, knee and ankle all contribute to move CoM in the same direction, i.e., the hip and ankle flex, and the knee extend to move CoM backward.

[Figure omitted. See PDF.]

Table 5. Joint angular correlations.

https://doi.org/10.1371/journal.pone.0163212.t005

To evaluate whether the generated biped stance motion was stable, CoM AP displacement versus CoM AP velocity was plotted, as shown in Fig 8. As shown in the figure, the attractor of each ||u_ff|| is a limit cycle, which indicates that all of the ||u_ff|| capable of keeping the musculoskeletal model standing stably. The CoM oscillation may result from the non-linear dynamics of the system (e.g. muscular-tendon dynamics).

[Figure omitted. See PDF.]

Fig 8. CoM AP displacement vs. CoM AP velocity.

https://doi.org/10.1371/journal.pone.0163212.g008

Discussion

Our first goal was to develop an NC model capable of compensating for a 120-ms neurological time delay, so as to allow a musculoskeletal model to stand. The postural control model proposed in this study is the first to have been reported to successfully maintain a musculoskeletal model in a standing posture for the case of a physiologically plausible human anatomy (multiple joints, 70 muscles, and human-like skeletal inertial and muscle dynamics properties) and a 120-ms neurological time delay. In previous studies, researchers have generally employed an inverted pendulum model to investigate the delay compensation mechanism; however, such a model cannot be used to investigate the muscle activation contributions. Hence, we replaced this simplified human model with a human-like musculoskeletal model. Our simulation results indicate that the u_ff, in conjunction with the fb control, can compensate for the neurological time delay. In this study, the fb control employed proprioceptive sensory information only, because this was assumed to have more relevance than other sensory data such as visual and vestibular inputs during quiet standing [69]. However, Peterka has reported that multisensory inputs are important for humans to perform a task, and the contributions of vestibular sensory inputs increase with an increase in the level of external disturbance [7]. The present simulation results indicate that, for an unperturbed stance, a certain ||u_ff|| functioning with proprioceptive sensory feedback can compensate for the loss of visual and vestibular sensory input. For the case of a perturbed stance, however, we believe that other sensory input should be incorporated in order to counter the external disturbance. Further, an estimation and anticipation mechanism [20, 30], which would estimate the disturbance acting on the human body and “the internal model” (usually the human body orientation and posture), may also be necessary to counter such a disturbance.

Our second goal was to develop an NC model capable of generating a human-like stance and exhibiting the salient features of the human biped stance. One feature reflected by our NC model is physiologically plausible simulated muscle activations. In particular, the muscle activations for ||u_ff|| = 2.07 are considered to be physiologically plausible activations for the human biped stance. That is, the majority of the muscle activations for ||u_ff|| = 2.07 were within the physiologically plausible range, based on data obtained via an experimental study, although the activation result for rectus abdominus 2 was slightly higher than the experimental data range. Note that rectus abdominus 2 has higher activations because the functions of all the muscles surrounding the lumbar that contribute to lumbar flexion are condensed into only two lumbar flexors (rectus abdominus 1 and rectus abdominus 2) in the musculoskeletal model. Hence, rectus abdominus 1 and rectus abdo minus 2 require slightly elevated activations to generate the forces that are actually produced by all of the lumbar flexors. The second feature reflected by the NC model is that humans employ a strategy involving a low muscle active stiffness during quiet standing, so as to achieve low energy consumption. The present simulation results indicate that human beings may be capable of standing using various active stiffness levels. ||u_ff|| = 2.07 is a physiologically plausible active stiffness level, and is low compared with the highest level (||u_ff|| = 8.02) among the selected u_ff. Therefore, for a normal person standing with a relaxed posture, a low active stiffness level may be sufficient to compensate for the neurological time delay and maintain a standing posture. This coincides with the well-known physiological result that humans select a low active stiffness level during quiet standing to reduce energy consumption.

The two features discussed above coincide with current physiological knowledge on the human biped stance; therefore, our proposed NC model and variable design framework successfully generated a physiologically plausible human-like biped stance and may be used in various potential applications, such as to assist with device development and design robotic control systems.

Furthermore, the NC model can also partly represent the characteristics of sway patterns. The postural sway patterns observed when musculoskeletal model stands under physiologically plausible muscle activations (||u_ff|| = 2.07) partly resemble the features of that of human beings. The hip and ankle angle are negatively correlated; Both hip-knee and knee-ankle angle have a positive correlation (Table 5). This result coincides with the correlation coefficient obtained from experimental data [43]. However, the simulated CoM AP and CoP AP ranges are smaller than those indicated by the experimental data. This difference is likely to be because the sensory noise level, which is assumed to be one of the variants accounting for the postural sway [75], is not incorporated in the fb control. The sensory noise may induce larger CoM AP sway range, as well as higher anti-phase coupling of the hip and the ankle to maintain the balance. In addition, the neglect of factors such as age and heart rate may have caused the smaller postural sway. However, it would be difficult to develop an NC model that incorporates all the factors that influence postural sway. In this study, we primarily focused on whether or not the employed NC model could generate muscle activations for a stable biped stance. The NC model utilized in this study provides a foundation for the development of a more complex NC model, which could reflect more physiologically plausible features such as postural sway.

Limitations and Future Work

One limitation of our study is that some rational simplifications were made to the PD controller. More gains should be included to achieve a more natural and stable biped stance simulation, which would enable more physiologically plausible muscle activations.

In addition, our current postural control model can only be used to simulate quiet standing; how it will deal with disturbances should also be investigated. Before that, however, a more sophisticated model that incorporates more physiologcally plausible components or mechanism should be created. Firstly, the fb controller should incorporate more fb information, such as visual and vestibular sensory inputs. Secondly, Mergner has noted that sensory inputs may have a threshold and be affected by noise [20], and a sensory integration and disturbance anticipatory mechanisms may exist. Such a mechanism should therefore be employed in our postural control model to investigate the ff input response to this mechanism as a means of maintaining balance. In addition, other important indicators for stance postural control such as the sway path, sway density, and power spectral density should be investigated in the future to evaluate the model.

Further, whether the model is overfitting or not should be validated. The NC model coordinates 35 muscles (left and right muscles are controlled symmetrically) to maintain the musculoskeletal model in a stance posture. However, only 8 of 35 simulated muscle activations were compared to experimental data. More muscle activation data should be measured in the experiment, and then be used to validate the NC model. Moreover, the complexity of the NC model and musculoskeletal model necessary to study bipedal stance postural control should also be investigated.

A final limitation is that the comparison between the simulated results and experimental data was inaccurate, owing to the difference in the weights, heights, and stance posture of the musculoskeletal model and the experimental subjects. The comparison was also affected by the EMG data normalization method employed by Panzer et al. [71]. In future, experiments on human quiet standing should be conducted as part of the controller design project, with the experimental setup matching the simulation conditions.

Conclusion

An NC model was developed to generate a human-like biped stance. Rather than an inverted pendulum, a musculoskeletal model was used to approximate the human anatomy. The NC model utilized in the postural control model consisted of ff and fb controls. Further, a variable design framework was developed for the NC model so as to maintain the musculoskeletal model in a standing position under the influence of a 120-ms neurological time delay. The NC model generated physiologically plausible muscle activations for the biped stance. The NC model also reflected a salient physiological feature, i.e., that humans select a low active stiffness level during standing so as to achieve low energy consumption.

Supporting Information

[Figure omitted. See PDF.]

S1 Video. Biped stance simulation video.

https://doi.org/10.1371/journal.pone.0163212.s001

(MP4)

Acknowledgments

The authors are grateful to the OpenSim community members for their helpful comments on code design and memory-leak fixes. We would also like to thank Editage (www.editage.jp) for English language editing.

Author Contributions

1. Conceptualization: PJ RC KT JO.

2. Data curation: PJ.

3. Formal analysis: PJ.

4. Funding acquisition: RC KT JO.

5. Investigation: PJ.

6. Methodology: PJ.

7. Project administration: JO.

8. Resources: PJ.

9. Software: PJ.

10. Supervision: JO.

11. Validation: PJ.

12. Visualization: PJ.

13. Writing – original draft: PJ.

14. Writing – review & editing: PJ RC KT JO.

Citation: Jiang P, Chiba R, Takakusaki K, Ota J (2016) Generation of the Human Biped Stance by a Neural Controller Able to Compensate Neurological Time Delay. PLoS ONE 11(9): e0163212. https://doi.org/10.1371/journal.pone.0163212

References

1. Chiba R, Takakusaki K, Ota J, Yozu A, Haga N. Human upright posture control models based on multisensory inputs; in fast and slow dynamics. Neuroscience research. 2016;104:96–104. pmid:26746115

2. van der Kooij H, van Asseldonk E, van der Helm FCT. Comparison of different methods to identify and quantify balance control. Journal of Neuroscience Methods. 2005;145(1–2):175–203. pmid:15922036

3. Winter DA, Patla AE, Prince F, Ishac M, Gielo-Perczak K. Stiffness control of balance in quiet standing. Journal of Neurophysiology. 1998;80(3):1211–1221. pmid:9744933

4. Morasso PG, Schieppati M. Can muscle stiffness alone stabilize upright standing? Journal of Neurophysiology. 1999;82(3):1622–1626. pmid:10482776

5. Morasso PG, Sanguineti V. Ankle muscle stiffness alone cannot stabilize balance during quiet standing. Journal of Neurophysiology. 2002;88(4):2157–2162. pmid:12364538

6. Loram ID, Lakie M. Direct measurement of human ankle stiffness during quiet standing: the intrinsic mechanical stiffness is insufficient for stability. The Journal of Physiology. 2002;545(3):1041–1053. pmid:12482906

7. Peterka RJ. Postural control model interpretation of stabilogram diffusion analysis. Biological Cybernetics. 2000;82(4):335–343. pmid:10804065

8. Peterka RJ. Sensorimotor integration in human postural control. Journal of Neurophysiology. 2002;88(3):1097–1118. pmid:12205132

9. Peterka RJ, Loughlin PJ. Dynamic regulation of sensorimotor integration in human postural control. Journal of Neurophysiology. 2004;91(1):410–423. pmid:13679407

10. Hettich G, Assländer L, Gollhofer A, Mergner T. Human hip–ankle coordination emerging from multisensory feedback control. Human Movement Science. 2014;37:123–146. pmid:25150802

11. Masani K, Vette AH, Kawashima N, Popovic MR. Neuromusculoskeletal torque-generation process has a large destabilizing effect on the control mechanism of quiet standing. Journal of Neurophysiology. 2008;100(3):1465–1475. pmid:18596181

12. Hettich G, Fennell L, Mergner T. Double inverted pendulum model of reactive human stance control. In: Multibody Dynamics Conference; 2011. Available from: http://www.posturob.uniklinik-freiburg.de.

13. van der Kooij H, De Vlugt E. Postural responses evoked by platform pertubations are dominated by continuous feedback. Journal of Neurophysiology. 2007;98(2):730–743. pmid:17460106

14. Masani K, Popovic MR, Nakazawa K, Kouzaki M, Nozaki D. Importance of body sway velocity information in controlling ankle extensor activities during quiet stance. Journal of Neurophysiology. 2003;90(6):3774–3782. pmid:12944529

15. Welch TDJ, Ting LH. A feedback model reproduces muscle activity during human postural responses to support-surface translations. Journal of Neurophysiology. 2008;99(2):1032–1038. pmid:18094102

16. Welch TDJ, Ting LH. A feedback model explains the differential scaling of human postural responses to perturbation acceleration and velocity. Journal of Neurophysiology. 2009;101(6):3294–3309. pmid:19357335

17. Hidenori K, Jiang Y. A PID model of human balance keeping. Control Systems, IEEE. 2006 Dec;26(6):18–23.

18. Masani K, Vette AH, Popovic MR. Controlling balance during quiet standing: Proportional and derivative controller generates preceding motor command to body sway position observed in experiments. Gait & Posture. 2006;23(2):164–172.

19. Boonstra TA, Schouten AC, Van der Kooij H. Identification of the contribution of the ankle and hip joints to multi-segmental balance control. Journal of Neuroengineering and Rehabilitation. 2013;10(1):23. pmid:23433148

20. Mergner T. A neurological view on reactive human stance control. Annual Reviews in Control. 2010;34(2):177–198.

21. Qu X, Nussbaum MA, Madigan ML. A balance control model of quiet upright stance based on an optimal control strategy. Journal of biomechanics. 2007;40(16):3590–3597. pmid:17628566

22. Qu X, Nussbaum MA. Evaluation of the roles of passive and active control of balance using a balance control model. Journal of biomechanics. 2009;42(12):1850–1855. pmid:19640541

23. Zhang H, Nussbaum MA, Agnew MJ. Development of a sliding mode control model for quiet upright stance. Medical engineering & physics. 2016;38(2):204–208.

24. Zhang H, Nussbaum MA, Agnew MJ. A new method to assess passive and active ankle stiffness during quiet upright stance. Journal of Electromyography and Kinesiology. 2015;25(6):937–943. pmid:26547842

25. Ackermann M, van den Bogert AJ. Predictive simulation of gait at low gravity reveals skipping as the preferred locomotion strategy. Journal of Biomechanics. 2012;45(7):1293–1298. pmid:22365845

26. van der Kooij H, Jacobs R, Koopman B, Grootenboer H. A multisensory integration model of human stance control. Biological Cybernetics. 1999;80(5):299–308. pmid:10365423

27. Allender L, Lockett MJ III, Martin MJ, Savage-Knepshield P. Designing Soldier Systems: Current Issues in Human Factors. Ashgate Publishing, Ltd.; 2013.

28. Kuo AD. An optimal control model for analyzing human postural balance. Biomedical Engineering, IEEE Transactions on. 1995;42(1):87–101.

29. Kuo AD. An optimal state estimation model of sensory integration in human postural balance. Journal of Neural Engineering. 2005;2(3):S235. pmid:16135887

30. van der Kooij H, Jacobs R, Koopman B, van der Helm F. An adaptive model of sensory integration in a dynamic environment applied to human stance control. Biological Cybernetics. 2001;84(2):103–115. pmid:11205347

31. Morasso PG, Baratto L, Capra R, Spada G. Internal models in the control of posture. Neural Networks. 1999;12(7):1173–1180. pmid:12662652

32. Borah J, Young LR, Curry RE. Optimal estimator model for human spatial orientationa. Annals of the New York Academy of Sciences. 1988;545(1):51–73. pmid:3071213

33. Elias LA, Forner-Cordero A. Postural control model of human stance with force feedback and passive ankle torque. In: 3rd National Meeting of Biomechanics Engineering (ENEBI) 2011; 2011. p. 153–154.

34. Elias LA, Watanabe RN, Kohn AF. Spinal mechanisms may provide a combination of intermittent and continuous control of human posture: Predictions from a biologically based neuromusculoskeletal model. PLoS Computational Biology. 2014 11;10(11):e1003944. pmid:25393548

35. Jo S, Massaquoi SG. A model of cerebellum stabilized and scheduled hybrid long-loop control of upright balance. Biological Cybernetics. 2004;91(3):188–202. pmid:15372241

36. Bottaro A, Casadio M, Morasso PG, Sanguineti V. Body sway during quiet standing: is it the residual chattering of an intermittent stabilization process? Human Movement Science. 2005;24(4):588–615. pmid:16143414

37. Bottaro A, Yasutake Y, Nomura T, Casadio M, Morasso P. Bounded stability of the quiet standing posture: an intermittent control model. Human Movement Science. 2008;27(3):473–495. pmid:18342382

38. Kowalczyk P, Glendinning P, Brown M, Medrano-Cerda G, Dallali H, Shapiro J. Modelling human balance using switched systems with linear feedback control. Journal of The Royal Society Interface. 2012;9(67):234–245.

39. Gawthrop P, Loram I, Lakie M, Gollee H. Intermittent control: a computational theory of human control. Biological Cybernetics. 2011;104(1–2):31–51. pmid:21327829

40. Suzuki Y, Nomura T, Morasso P. Stability of a double inverted pendulum model during human quiet stance with continuous delay feedback control. In: Engineering in Medicine and Biology Society, EMBC, 2011 Annual International Conference of the IEEE. IEEE; 2011. p. 7450–7453.

41. Suzuki Y, Nomura T, Casadio M, Morasso P. Intermittent control with ankle, hip, and mixed strategies during quiet standing: a theoretical proposal based on a double inverted pendulum model. Journal of Theoretical Biology. 2012;310:55–79. pmid:22732276

42. Gawthrop P, Loram I, Lakie M. Predictive feedback in human simulated pendulum balancing. Biological Cybernetics. 2009;101(2):131–146. pmid:19588160

43. Günther M, Grimmer S, Siebert T, Blickhan R. All leg joints contribute to quiet human stance: A mechanical analysis. Journal of Biomechanics. 2009;42(16):2739–2746.

44. Sasagawa S, Shinya M, Nakazawa K. Interjoint dynamic interaction during constrained human quiet standing examined by induced acceleration analysis. Journal of Neurophysiology. 2014;111(2):313–322. pmid:24089399

45. Hsu WL, Scholz JP, Schöner G, Jeka JJ, Kiemel T. Control and estimation of posture during quiet stance depends on multijoint coordination. Journal of Neurophysiology. 2007;97(4):3024–3035. pmid:17314243

46. Horlings CG, van Engelen BG, Allum JH, Bloem BR. A weak balance: the contribution of muscle weakness to postural instability and falls. Nature Clinical Practice Neurology. 2008;4(9):504–515. pmid:18711425

47. Clark AE. Biarticular Muscles Influence Postural Responses: Implications for Treatment of Stiff-Knee Gait. Master’s Thesis, University of Tennessee. 2012;.

48. Zebenay M, Lippi V, Mergener T. Human-like humanoid robot posture control. In: Informatics in Control, Automation and Robotics (ICINCO), 2015 12th International Conference on. IEEE; 2015. p. 304–309.

49. Delp SL, Anderson FC, Arnold AS, Loan P, Habib A, John CT, et al. OpenSim: open-source software to create and analyze dynamic simulations of movement. Biomedical Engineering, IEEE Transactions on. 2007;54(11):1940–1950.

50. Au C. Gait 2392 and 2354 Models. 2013;Available from: http://simtk-confluence.stanford.edu:8080/display/OpenSim/Gait+2392+and+2354+Models.

51. Hicks J. Simulation-based design to prevent ankle injuries. 2014;Available from: http://simtk-confluence.stanford.edu:8080/display/OpenSim/Simulation-Based+Design+to+Prevent+Ankle+Injuries.

52. Millard M, Uchida T, Seth A, Delp SL. Flexing computational muscle: Modeling and simulation of musculotendon dynamics. Journal of Biomechanical Engineering. 2013;135(2):021005. pmid:23445050

53. Loram ID, Maganaris CN, Lakie M. Paradoxical muscle movement in human standing. The Journal of physiology. 2004;556(3):683–689. pmid:15047776

54. Liu MQ, Anderson FC, Schwartz MH, Delp SL. Muscle contributions to support and progression over a range of walking speeds. Journal of Biomechanics. 2008;41(15):3243–3252. pmid:18822415

55. Hamner SR, Delp SL. Muscle contributions to fore-aft and vertical body mass center accelerations over a range of running speeds. Journal of Biomechanics. 2013;46(4):780–787. pmid:23246045

56. Gopalakrishnan A, Modenese L, Phillips AT. A novel computational framework for deducing muscle synergies from experimental joint moments. Frontiers in Computational Neuroscience. 2014;8:153. pmid:25520645

57. Skalshøi O, Iversen CH, Nielsen DB, Jacobsen J, Mechlenburg I, Søballe K, et al. Walking patterns and hip contact forces in patients with hip dysplasia. Gait & Posture. 2015;42(4):529–533.

58. Hunt K, Crossley F. Coefficient of restitution interpreted as damping in vibroimpact. Journal of Applied Mechanics. 1975;42(2):440–445.

59. Sherman MA, Seth A, Delp SL. Simbody: Multibody dynamics for biomedical research. Procedia Iutam. 2011;2:241–261. pmid:25866705

60. Dorn TW, Wang JM, Hicks JL, Delp SL. Predictive simulation generates human adaptations during loaded and inclined walking. PLoS ONE. 2015 04;10(4):e0121407. pmid:25830913

61. Zajac FE. Muscle and tendon: Properties, models, scaling, and application to biomechanics and motor control. Critical Reviews in Biomedical Engineering. 1988;17(4):359–411.

62. Winters JM. An improved muscle-reflex actuator for use in large-scale neuromusculoskeletal models. Annals of Biomedical Engineering. 1995;23(4):359–374. pmid:7486344

63. Jacobs DA. First-order activation dynamics. 2015;Available from: http://simtk-confluence.stanford.edu:8080/display/OpenSim/First-Order+Activation+Dynamics.

64. Fitzpatrick R, Burke D, Gandevia SC. Loop gain of reflexes controlling human standing measured with the use of postural and vestibular disturbances. Journal of Neurophysiology. 1996;76(6):3994–4008. pmid:8985895

65. Gatev P, Thomas S, Kepple T, Hallett M. Feedforward ankle strategy of balance during quiet stance in adults. The Journal of Physiology. 1999;514(3):915–928. pmid:9882761

66. Hof A, Gazendam M, Sinke W. The condition for dynamic stability. Journal of biomechanics. 2005;38(1):1–8. pmid:15519333

67. Pai YC, Patton J. Center of mass velocity-position predictions for balance control. Journal of biomechanics. 1997;30(4):347–354. pmid:9075002

68. Heitmann S, Ferns N, Breakspear M. Muscle co-contraction modulates damping and joint stability in a three-link biomechanical limb. Front Neurorobot. 2011;5(5). pmid:22275897

69. Sousa AS, Silva A, Tavares JMR. Biomechanical and neurophysiological mechanisms related to postural control and efficiency of movement: A review. Somatosensory & motor research. 2012;29(4):131–143.

70. McClearn D. Anatomy of raccoon (Procyon lotor) and coati (Nasua narica and N. nasua) forearm and leg muscles: Relations between fiber length, moment-arm length, and joint-angle excursion. Journal of Morphology. 1985;183(1):87–115. pmid:3973925

71. Panzer VP, Bandinelli S, Hallett M. Biomechanical assessment of quiet standing and changes associated with aging. Archives of Physical Medicine and Rehabilitation. 1995;76(2):151–157. pmid:7848073

72. Hansen N, Müller SD, Koumoutsakos P. Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES). Evolutionary Computation. 2003;11(1):1–18. pmid:12804094

73. The CMA Evolution Strategy. 2015;Available from: https://www.lri.fr/~hansen/cmaesintro.html.

74. Warnica MJ, Weaver TB, Prentice SD, Laing AC. The influence of ankle muscle activation on postural sway during quiet stance. Gait & Posture. 2014;39(4):1115–1121.

75. Maurer C, Peterka RJ. A new interpretation of spontaneous sway measures based on a simple model of human postural control. Journal of Neurophysiology. 2005;93(1):189–200. pmid:15331614