Modulation of gluteus medius activity reflects the potential of the muscle to meet the mechanical demands during perturbed walking

Mediolateral stability during walking can be controlled by adjustment of foot placement. Reactive activity of gluteus medius (GM) is modulated during the gait cycle. However, the mechanisms behind the modulation are yet unclear. We measured reactive GM activity and kinematics in response to a mediolateral platform translation during different phases of the gait cycle. Forward simulations of perturbed walking were used to evaluate the isolated effect of the perturbation and the GM response on gait stability. We showed that the potential of GM to adjust lateral foot placement and prevent collisions during swing varies during the gait cycle and explains the observed modulation. The observed increase in stance, swing or combined GM activity causes an outward foot placement and therefore compensates for the loss of stability caused by a perturbation early in the gait cycle. GM activity of the swing leg in response to a platform translation late in the gait cycle counteracts foot placement, but prevents collision of the swing foot with the stance leg. This study provides insights in the neuromechanics of reactive control of gait stability and proposes a novel method to distinguish between the effect of perturbation force and reactive muscle activity on gait stability.


Muscle moment arms
Computation speed in the forward simulations was increased by approximating muscle tendon length -joint angle relationships with polynomial functions. Coefficients were estimated by minimzing the difference between the muscle tendon length computed with the polynomial approximation an OpenSim's Muscle Analysis Tool. The order of the polynomial was increased stepwise until the difference between the muscle tendon length and moment arms computed with the polynomial approximation and OpenSim's Muscle Analysis Tool was smaller than 1.5 mm for the kinematics observed during normal and perturbed during walking [6]. Moment arms were computed as the partial derivative of the muscle-tendon length to the selected degree of freedom (equation 2).
where dM j is the moment arm of the muscle for the joint angle q j , L M T is the muscle tendon length. Muscle tendon velocity was computed as the product of the moment arm and the angular velocity (Equation 3) where V M T is the muscle-tendon velocity and N dof the number of degrees of freedom spanned by the muscle.

Ground contact model
Ground contact forces were simulated with a spring-damper based contact model. Contact geometry was modelled by one contact sphere at the heel and three at the metatarsal heads ( Figure  1). Vertical reaction force F y acting on each sphere was computed based on the indentation and indentation velocity of the contact sphere in the ground [4].
where k is the stiffness parameter, r is the radius of the sphere, y is the indentation of the sphere in the ground, c is the damping parameter andẏ is the indentation velocity of the sphere in the ground. Details on the contact parameters are shown in table 2 . Friction forces (F x and F z ) were computed as a function of the vertical ground reaction force and the horizontal velocity of the contact sphere.
where F x and F z is the friction force in respectively anterior-posterior and medio-lateral direction,ẋ andż are respectively the velocity of the contact sphere in the anterior-posterior and medio-lateral direction.

Tracking skeletal dynamics
We created a model-based simulation of unperturbed walking by tracking the experimental data for a gait cycle of unperturbed walking. First, we solved for the joint torques that, when input to the skeleton dynamics, minimized the difference between simulated and measured movement (details below). Next, we solved for the muscle excitations that produced these joint torques by solving the muscle redundancy problem (see section 3). We used the model described in section 1 and replaced the muscle-tendon actuators by joint actuators. Hence, model inputs were the joint torques and the model's state (x) consisted of the generalized coordinates (q) and velocities (q).
Model inputs (joint torques) were optimized to minimize the difference between the measured and simulated marker positions, ground reaction forces and joint torques: where t 0 is the time of left heelstrike and t end the time of the consecutive left heelstrike, p is the vector containing the 3D marker coordinates of all 48 markers on the skeletal model,p is a vector with the corresponding measured marker coordinates, F andF are vectors containing respectively the simulated and measured ground contact forces and moments and nF equals two (ground reaction force under left and right foot).T is a vector containing the inverse dynamic joint torques based on the measured kinematics and ground reaction forces, T is a vector containing the corresponding simulated joint torques. The optimization was subject to skeleton dynamics. We solved this dynamic optimization problem using direct collocation. To improve the numerical condition of the resulting non-linear programming problem, we first reformulated the optimization problem as described below.

Problem formulation
First, we used an implicit formulation of skeleton dynamics [3]. To this aim, we introduced joint accelerations uq as additional controls, simplifying the dynamic equations: Skeleton dynamics were then imposed by a path constraint: Second, we defined the kinematic chain with the foot as a floating base, the position and orientation of the model was described by 6 degrees of freedom between the foot and the ground. By choosing the foot instead of the pelvis as a floating base, we obtained a more direct relation between the model's states (kinematics of the foot) and the ground reaction forces, which improved convergence of the optimization problem. As a consequence, we had to define seperate kinematic chains for the left and right leg, which were coupled at the pelvis [3]. Each chain had a reduced pelvis segment with a mass and inertia that equalled half of the mass and inertia of the original pelvis segment. We introduced forces F p and moments M p acting on the right reduced pelvis as additional controls. Reaction forces and moments were applied on the left reduced pelvis. We then added additional path constraints to impose that the two reduced pelvises had the same kinematics. For both models, skeleton dynamics was evaluated through OpenSim's Inverse Dynamics Tool during the optimization. The residual forces and moments (i.e. external forces and moments) computed with inverse dynamics were constrained to be equal to the ground reaction forces computed with the contact model implemented in Matlab.

Summary optimal control problem
In summary, we solved the following dynamic optimization problem through direct collocation: 2. Controls: 3. Dynamic constraints:ẋ = q uq (13) 4. Path constraints: • Input joint torques should equal the inverse dynamics torques computed based on the model's kinematics.
• Residual forces and moments at the foot from OpenSim's Inverse Dynamics Tool F ground−f oot should equal the forces and moments computed based on the contact model.
• position and orientation of the pelvis in the left and right model should be equal where F p and M p are the forces and torques to simulate the action of one model on the other, T ID is a vector containing the inverse dynamics torques, F ground−f oot is a vector containing the residual inverse dynamics forces and moments, F contact is vector containing the ground reaction forces and moments computed with the ground contact model, P P elvisLef t and P P elvisRight are vectors containing the position and orientation of the left and right pelvis in the ground.

Implementation
The direct collocation software GPOPS II was used to formulate the non-linear programming problem, which was solved using IPOPT. Inverse dynamics and the computation of the marker locations was solved using the OpenSim libraries trough a Mex file. Computation time was decreased by distributing the collocation points between multiple cores using OpenMPI. The ground contact forces and moments were computed in Matlab.

Tracking results
The root-mean square error of the tracking simulation of unperturbed walking was 0.03 cm for the marker positions, 29N for the ground reaction forces and 6.7 Nm for the joint torques. The tracking of the ground reaction forces ( Figure 2) and joint torques (Figure 3) is shown for one representative trial.

Tracking muscle dynamics
In a second step, the muscle excitations that generate the inverse dynamic joint torques were computed using dynamic optimization. We modified the dynamic optimization approach proposed by [2] to include tracking of the gluteus medius activity (measured with electromyography). The objective function contains a part that minimizes muscle excitations squared, a second part the minimizes the differences between the measured and simulated excitations of the gluteus medius and a third part that minimizes the torque generated by ideal torque actuators (i.e. reserve actuators according to OpenSim terminology) (Equation 17).
where nMus=92 is the number of muscles in the model, e i is the simulated excitation of muscle i,ê GM is the processed EMG of the gluteus medius, s is an optimization variable and scales the EMG from voltage to muscle excitations, e GM is the simulated excitation of the gluteus medius, T max equals 150Nm and is the reserve actuator of degree of freedom j, e T j is the activation of the reserve actuator of DOF k. Note that the weight w 2 is high to penalize the use of the reserve actuators. The reserve actuators were added to guarantee feasibility in the presence of fast changes in joint torques in the tracking simulation that can not be generated by the muscles.This objective function was minimized subject to dynamic constraints following from muscle activation and contraction dynamics and to path constraints imposing that the model's actuators (muscles, passive torque, reserve actuators) should produce the inverse dynamics joint torques. In this case, the inverse dynamics joint torques were the result of the tracking simulation described in section 2. Reserve actuator torques did not exceed 2 Nm in the simulation. Joint torques resulting from the muscles, passive forces and reserve actuators were imposed to be equal to the joint torques from tracking the skeletal dynamics. The tracking of the gluteus medius activity in one representative trail can be found in Figure 4. The source code of the muscle redundancy solver can be found on the simtk project optcntrlmuscle (https://simtk.org/projects/optcntrlmuscle).

Forward simulation
The initial state and simulated muscle excitations from the tracking simulations were used to reconstruct the unperturbed walking in a forward simulation. The passive response to the perturbation was simulated by imposing the platform translation as an external force in the forward simulation of unperturbed walking. The active response was simulated by imposing the perturbation force and adding the measured muscle response to the perturbation (i.e. the difference between measured activity during unperturbed walking and after perturbation) to the muscle excitations from the tracking simulation of unperturbed walking. Muscle and skeletal 5 Sensitivity analysis

Stiffness contact model
The sensitivity of the simulated active and passive response to contact model properties was evaluated by modifying the stiffness of the contact model in the tracking and forward simulation. The stiffness k was stepwise increased in the tracking simulation of one representative trial (Equation 4). The results from the tracking simulation and the modified contact model were used in the forward integration of the passive and active response. The stiffness had only a small influence on the simulated ground reaction forces in the tracking simulation and had only a limited effect on the simulated stride width and margin of stability in the active and passive response ( Figure 5). Since the contact stiffness had only a small influence on the simulation results, we selected a relatively low stiffness value k = 10 to decrease the stiffness of the differential equations.

Lower bounds on activity
It has been hypothesized that muscle co-contraction is used as a strategy to control balance during walking by increasing the joint stiffness [5]. The influence of this co-contraction strategy on the simulated passive response was evaluated by constraining the minimal muscle activity in the muscle redundancy solver. The minimal muscle activity was increased stepwise in the muscle redundancy solver. The stride width and margin of stability at first heelstrike after perturbation was computed in the passive and active response using forward simulations for the different levels of minimal muscle activity. Stride width decreased with increased muscle activity and a small change in the margin of stability was found in the active and passive response simulations ( Figure  6).

Heelstrike detection
Left and right heel strikes were determined from the ground reaction force. When the subjects were walking with one foot on each belt, a simple vertical force threshold of 30N was used to detect initial contact of the foot with the ground. An additional criteria was needed to detect heelstrikes when the perturbation caused an inward stepping strategy (subjects were therefore walking with both feet on one belt). Peaks in anterior-posterior movement of the combined center of pressure were identified using Matlab's peakfind function and were used to detect heelstrikes (Figure 7). The difference between the center of pressure based and vertical force threshold is on average 0.007s for normal split belt treadmill walking (standard deviation equals 0.002s).  Figure 7: Peaks in anterior-posterior movement of the combined center of pressure (blue line) were used to detect heelstrikes when subjects were walking with both feet on one belt. The black and green dots are the peak in the anterior-posterior movement of the combined center of pressure and represent respectively the left heelstrike (lhs) and right heelstrike (rhs). The red and blue dots represent respectively the left toe-off (lto) and right toe-off (rto).

Marker protocol
Whole-body motion was recorded with an extended plug in gait marker set (figure 8).