Coordinated alpha and gamma control of muscles and spindles in movement and posture

Mounting evidence suggests that both α and γ motoneurons are active during movement and posture, but how does the central motor system coordinate the α-γ controls in these tasks remains sketchy due to lack of in vivo data. Here a computational model of α-γ control of muscles and spindles was used to investigate α-γ integration and coordination for movement and posture. The model comprised physiologically realistic spinal circuitry, muscles, proprioceptors, and skeletal biomechanics. In the model, we divided the cortical descending commands into static and dynamic sets, where static commands (αs and γs) were for posture maintenance and dynamic commands (αd and γd) were responsible for movement. We matched our model to human reaching movement data by straightforward adjustments of descending commands derived from either minimal-jerk trajectories or human EMGs. The matched movement showed smooth reach-to-hold trajectories qualitatively close to human behaviors, and the reproduced EMGs showed the classic tri-phasic patterns. In particular, the function of γd was to gate the αd command at the propriospinal neurons (PN) such that antagonistic muscles can accelerate or decelerate the limb with proper timing. Independent control of joint position and stiffness could be achieved by adjusting static commands. Deefferentation in the model indicated that accurate static commands of αs and γs are essential to achieve stable terminal posture precisely, and that the γd command is as important as the αd command in controlling antagonistic muscles for desired movements. Deafferentation in the model showed that losing proprioceptive afferents mainly affected the terminal position of movement, similar to the abnormal behaviors observed in human and animals. Our results illustrated that tuning the simple forms of α-γ commands can reproduce a range of human reach-to-hold movements, and it is necessary to coordinate the set of α-γ descending commands for accurate and stable control of movement and posture.


Introduction
The physiological system of human motor control is not only highly redundant (Bernstein, 1967;Martin et al., 2009), but also endowed with intricate dual α and γ sensorimotor control (Granit, 1975;Pierrot-Deseilligny and Burke, 2005;Lan and He, 2012;Prochazka and Ellaway, 2012). Our understanding about how movements are organized and how muscles are coordinated in performing different tasks remains incomplete due to lack of in vivo data during behaviors. One of the remaining issues in sensorimotor control is to account for the role of γ motor system in movement and posture. In spite of thorough elucidation of peripheral efferent and afferent innervations of the spindle organ (Matthews, 1964;Boyd, 1980;Hulliger, 1984;Prochazka, 1999), the importance of γ motor system in motor control is still not well understood.
A consistently observed phenomenon during both movement and posture control is α-γ co-activation (Vallbo, 1971;Taylor et al., 2004Taylor et al., , 2006Prochazka and Ellaway, 2012). Direct recording by Taylor et al. (2004Taylor et al. ( , 2006) revealed a co-varying pattern of gamma-static and dynamic firings with joint angle during locomotion. There was plenty evidence of independent control of γ -motoneurons during movement (Prochazka et al., 1985;Dimitriou and Edin, 2008). The co-activation of γ motoneurons with α motor activity was generally viewed to compensate for the unloading effects of muscle contraction to spindle sensitivity. Using a realistic virtual arm (VA) model , Lan and He (2012) suggested a plausible function of γ s fusimotor control to convey centrally encoded joint angle information to the periphery through regulating spindle sensitivity, so that the Ia signaling from spindle afferents is kept faithfully proportional to the joint angle during movement and muscle contraction. This not only explains γ co-activation with α activity, but also supports the hypothesis that the γ s command reinforces the centrally planned kinematics of movement and posture by way of spinal circuits.
Neurophysiological studies have identified separate spinal pathways and circuits of sensorimotor system in details, where the α and γ commands interact with each other to produce sensory and motor outputs (Baldissera et al., 1981;Lemon et al., Abbreviations: α, γ , Alpha, gamma motor nerves and neurons; PN, Propriospinal neuron; EMG, Electromyography; GTO, Golgi Ten Organ; CS-VA, Corticospinal virtual arm; VA, Virtual arm; α s , γ s , Static alpha, gamma command; α d , γ d , Dynamic alpha, gamma command; αMN, γ MN, Alpha, gamma motoneuron; U m , Muscle input in the CS-VA model; U s , Spindle input in the CS-VA model; F m , Muscle force; L ce , Muscle fascicle length; L mt , Musculo-tendon length; Ia, Primaryafferent from spindle; II, Secondary afferent from spindle; Ib, Afferent from GTO; PC, Pectoralis Major Clavicle; DP, Deltoid Posterior; BS, Brachialis; Tlt, Triceps Lateral; Bsh, Biceps Short Head; Tlh, Triceps Long Head; θ el , Elbow angle; θ sh , Shoulder angle; γ -PN, Dynamic gamma inhibition on PN; Ia(−), Ia reciprocal inhibition on αMN; Ia-PN, Ia afferent gain on PN; Ia(+), Stretch reflex; Ib(−), Ib inhibition from GTO; RC, Renshaw Cell; Error P , Error M , mean difference between simulation and experiment of terminal posture, movement; V peak , Peak velocity; 3SD, Three times of standard deviation of static EMG; AG1, First agonistic muscle burst phase; AG2, Second agonistic muscle burst phase; ANT, Antagonistic muscle burst phase; RS, Ratio of static EMG (Bsh/Tlh); RP, Ratio of peak EMG (ANT/AG1); amp AG1 amp ANT , Amplitude of AG1 and ANT phase on α d ; pw AG1, pw ANT , Pulse width of AG1 and ANT phase on α d ; MJT, Minimum jerk trajectory. 2004;Pierrot-Deseilligny and Burke, 2005). In addition to the mono-synaptic cortico-motoneuronal pathway (Lawrence and Kuypers, 1968;Lemon, 1999;Quallo et al., 2012), disynaptic excitatory and inhibitory cortico-motoneuronal pathways via PNs in C3-C4 were found to exist in cats and in macaque monkeys, as well as in humans (Malmgren and Pierrot-Deseilligny, 1988;Gracies et al., 1994;Alstermark et al., 2007;Isa et al., 2007). The PN was shown to play an important role in reaching movement of upper limb (Alstermark et al., 1981;Alstermark and Isa, 2012). In a computational analysis to emulate involuntary oscillatory movements in human upper extremity (Hao et al., 2013), it is postulated that movement signals are transmitted via the disynaptic pathway of the PN network, where the γ -dynamic (γ d ) command is integrated with the α-dynamic (α d ) command to produce pre-motoneuronal outputs. The γ d command encodes kinetic information of joint acceleration (or deceleration), and gates the α d command of double frequency at the PN network to determine the timing of activation for a pair of flexor and extensor muscles during oscillatory movements.
These studies suggest that while the α motor system provides the main drive for muscles, the γ motor system executes a more subtle control for movement dynamics and the maintenance of posture. In this paper, we extend the α-γ model (Hao et al., 2013) to investigate the modular control of voluntary movement and posture, and to demonstrate coordination of a set of α-γ descending commands in the control of movement and posture. Human reach-and-hold movements and muscle EMG activities were recorded and analyzed to guide the specification of the central descending commands. The model behavior was matched to the data of human movement and posture in all subjects. A variety of modifications were introduced in model structures to assess the functional significance of α-γ coordination by "deafferenting" or "deefferenting" the model. Results of this study provided a quantitative evaluation of the relative importance that each proprioceptive afferent and descending α-γ command may have on movement and posture control. Agreement between model predictions and human movement data further corroborated the α-γ coordination as a rudimentary means of sensorimotor control. The results argue that coordinated γ control with α activation is essential for accurate and stable control of movement and posture. Preliminary analysis of this study appeared in a conference proceeding (Li et al., 2014).

Corticospinal Virtual Arm (CS-VA) Model
The module of posture and movement control based on α-γ dual control system was implemented in the computational CS-VA model (Figure 1). This model is based on realistic physiological studies (Cheng et al., 2000;Taylor et al., 2006;Alstermark et al., 2007;Isa et al., 2007) and has been validated to capture the realistic neuromechanical properties of human upper limb in previous work (Song et al., 2008a,b;Lan and He, 2012;He et al., 2013). It consists of four parts, the primary motor and sensory cortex, PN network, spinal cord circuitry and virtual arm (VA) model. The spinal cord network with PNs innervating a pair of antagonistic FIGURE 1 | Corticospinal virtual arm (CS-VA) model. Descending α and γ commands from motor cortex were processed in PN network and spinal cord circuitry, the coordinated output of α and γ motor neurons are transmitted to virtual muscles and spindles, the virtual arm is then activated. PN represents propriospinal interneuron. The subscripts "s" and "d" in α and γ motor signals stand for static and dynamic commands, respectively. U m is muscle input, which is equivalent to human EMG; U s is spindle input; L mt and F m represent muscle tendon length and muscle force computed form virtual arm model; L ce represents fascicle length; Ia and II are sensory feedback of primary and secondary afferents from muscle spindle; and Ib is the feedback of Golgi Tendon Organ (GTO). The virtual arm has two joints including shoulder and elbow in the horizontal plane, with two degrees of freedom, i.e., shoulder flexion/extension and elbow flexion/extension. Three pairs of antagonistic muscles are included in the model: shoulder flexor Pectoralis Major Clavicle (PC), extensor Deltoid Posterior (DP), elbow flexor Brachialis (BS), extensor Triceps Lateral (Tlt), and biarticular flexor Biceps Short Head (Bsh), extensor Triceps Long Head (Tlh). Simulated joint kinematics, especially elbow angle is compared to elbow angle (θ el ) recorded in human experiment. Modified from Hao et al. (2013) with permission. muscles (Figure 2) has been validated in the mechanism study of Parkinsonian tremor (Hao et al., 2013). The VA has two joints including shoulder and elbow in the horizontal plane, two degrees of freedom, namely shoulder flexion/extension and elbow flexion/extension. The six dominating muscles includes two pairs of monoarticular muscles, shoulder flexor Pectoralis Major Clavicle (PC), extensor Deltoid Posterior (DP), elbow flexor Brachialis (BS), extensor Triceps Lateral (Tlt) and one pair of biarticular muscles, flexor Biceps Short Head (Bsh), extensor Triceps Long Head (Tlh). Muscle spindle and Golgi Tendon Organ (GTO) are implemented within the VA to provide proprioceptive feedback. Three parts of the CS-VA model, including PN network, spinal cord circuitry and VA, have been integrated in SIMULINK/MATLAB platform. In this study, biarticular muscles of Bsh and Tlh were mainly used to realized movement and posture.
Within the CS-VA model, central motor commands from motor cortex are transmitted to spinal α and γ motoneurons (αMN, γ MN) in two pathways, mono-synaptic pathway carrying static α and γ motor commands to corresponding motoneurons for posture control, and multi-synaptic pathway transmitting dynamic α and γ motor commands to MNs via PN network for movement control. Muscles and spindles are activated by α and γ motor signals respectively after the regulation of PN and reflex network within spinal cord. The muscular dynamics computed by SIMM model of arm give a real-time muscle tendon length to virtual muscle, at the same time, the dynamics are reflected by Ia and II afferents from muscle spindle, and Ib afferent from GTO. Feedback from proprioceptive afferents participates in signal processing within spinal cord network and cortex simultaneously. In this framework, central inputs of static and dynamic α and γ motor commands are specialized in Section Specification of Central Descending Commands, the input of virtual muscle (U m ) is equivalent to recorded EMG, and the output joint kinematics gives a behavioral indicator to compare with elbow angle (θ el ) in human experiment.
PNs in C3-C4 receive wide descending control from cortico-, rubro-, reticulo and tectospinal tracts, and project inhibition or excitation to downstream αMN , at the same time receive feedback from Ia and cutaneous afferents (Alstermark et al., 1984a,b,c). Based on neurophysiology structures of PN network, αMN, γ MN, and the spinal reflexes, a model of spinal circuitry was implemented within the CS-VA model as is shown in Figure 2 (Hao et al., 2013;Li et al., 2014). The central descending commands in the model were passed down to αMN and γ MN through mono-synaptic and multi-synaptic PN pathways, respectively . There was experimental evidence indicating that PNs mediate motor commands for reaching movement Alstermark and Isa, 2012). Computational analysis (Lan and He, 2012) suggested that γ s conveys kinematic information of joint angle. Therefore, in this model, we assume that the descending commands can be divided into static set (α s and γ s ) for posture and dynamic set (α d and γ d ) for movement. This division is consistent to the evidence of individual posture and movement modules in central motor control system (Kurtzer et al., 2005).
As illustrated in Figure 2, α and γ commands control muscles and spindles through a set of spinal circuits. α d is integrated at the FIGURE 2 | α-γ control model in spinal cord circuitry for a pair of antagonistic muscles. α and γ commands are coordinated by PN network and reflex circuitry within spinal cord. The subscripts "e" and "f" stand for extensor and flexor; MN represents motor neuron pool; U represents muscle input; d e and d f are gains of γ d inhibition on PN(γ -PN); p e and p f are gains of PN to reciprocal inhibition; g e and g f are gains of recurrent inhibition from Renshaw Cell (RC); b e and b f represent inhibition gains on α motor neurons from Ib afferents (Ib(−)); s e and s f are stretch reflex gains from Ia afferents ((Ia(+)); r e and r f are reciprocal inhibition gains on α motor neurons (Ia(−)); a e and a f are Ia afferent gains on PN(Ia-PN). Excitatory synapse at neurons is indicated by a "y" termination, and inhibitory synapse is indicated by a filled dot termination. The function of the feedback from Ia afferents in dashed lines are discussed later.
PN with γ d in a pattern of mirror inhibition (γ d or 1-γ d ). The recording from dynamic γ MN in Taylor et al. (2006) revealed a reciprocal change of γ d with joint angle. The PN also receives excitation from autogenic Ia afferent (Ia-PN) (Alstermark et al., 1984b;Malmgren and Pierrot-Deseilligny, 1988). After the PN, α s and α d commands converge at αMN, and its output is further regulated by reciprocal inhibition (Ia(−)) from antagonistic muscle, recurrent inhibition of Renshaw Cell (RC), autogenic stretch reflexes from Ia (Ia(+)), and Ib afferents (Ib(−)) (Eccles et al., 1957(Eccles et al., , 1960Windhorst, 2007). The output of PN (Y PN ) and αMN (Y αMN ) are given in Equation (1), in which subscript "e" and "f " represent extensor and flexor, is the sum of descending commands to αMN, C(t) represents background activation of αMN pools, τ is the time constant of excitation with the value 0.029 (sec), σ (t) represents the signal dependent noise, which is a Gaussian distributed random signal (Fuglevand et al., 1993;Jones et al., 2002). C(t) was multiplied by σ (t) to yield the noise corrupted background activations (C′(t)) in Equations (1g) and (1h). The output of αMN (Equations 1i and 1j) is the sum of all excitatory and inhibition inputs; here, the values of υ′ and ϕ′ are proportional to Ia and Ib afferent discharge frequencies, respectively; and g, s, r, b are reflex gains of RC, Ia(+), Ia(−), Ib(−). The values of spinal reflex gains are listed in Table 1.

Subjects and Experiments
Seven healthy adults participated in this elbow joint's reach and hold experiment. The human subject study was approved by the Internal Review Board (IRB) of University of Southern California (USC). All of them obtained a brief explanation of this study before the experiment, and signed the informed consent. During the experiment, the subject sat comfortably with the upper arm maintained perpendicular to the trunk in the horizontal plane (Figure 3A), the hand and forearm were secured on the manipulation to make single joint movements around the elbow. All subjects initially held their elbow at 90 • , and performed successive reach and hold movements triggered by the visual cue of LED light, following three blocks. In Block 1, the elbow was extended with the range of 30 degrees, and after Analysis of kinematic and EMG data during posture and movement. During movement, the elbow moves from initial posture to terminal posture, and the movement onset (t 0 ) and offset (t 1 ) are defined as the time at which the velocity change (increases or decreases) is 10% of the peak velocity (V peak ). Initial posture period starts at time (t 0 -1.5) and terminal posture is defined from time (t 0 + 4) after the angle stabilizes. The time window is 1 s. The tri-phasic EMG pattern of extensor Triceps long head (Tlh) and flexor Biceps short head (Bsh with inverse value) is shown in the bottom (low passed at a cut off frequency of 50 Hz with digital Butterworth filter). three successive extensions, the elbow stopped at 0 • , the reversal flexion with the same movement range and holding posture was performed after a 5 s holding period at 0 • . In Block 2, the range of extension and flexion was 45 • , thus one holding period at the mid-posture of 45 • was performed, and in Block 3 the elbow had direct extension and flexion movement between 90 • and 0 • . During the experiments, the subjects moved as fast as them could between postures, and a 5 s holding was required at each posture. Each block was repeated for five trials, and between blocks, the subjects had more than 10 min to rest. During the movement, the elbow angle was recorded and EMGs of biceps short head (Bsh) and triceps long head (Tlh) were collected using bi-polar surface electrodes, the EMG signals were pre-amplified at the gain of 1000, band-pass filtered with cut-off frequencies at 5 Hz and 1000 Hz, and then sampled at 2000 Hz for post-processing.

Analysis of Kinematics and EMG Data
Recorded joint trajectory and EMG were analyzed to provide guidance for simulation. Since Subject 3 couldn't follow the experimental protocol, the data was excluded from result analysis. The joint angles were low-pass filtered with a cutoff frequency of 10 Hz to remove high-frequency noise, and differentiated to obtain velocity. The collected EMG of biceps and triceps were band-pass filtered with a cut-off frequency between 20 and 500 Hz to remove motion artifacts and high-frequency noise, rectified, and then low passed at the cut-off frequency 50 Hz for further analysis. As shown in Figure 3, the time for movement onset (t 0 ) and offset (t 1 ) was defined as the time where velocity increased and decreased to 10% of peak velocity (V peak ) respectively (Atkeson and Hollerbach, 1985). The static angles and EMG during initial posture were then defined as the averaged value from time (t 0 − 1.5) (sec) to time (t 0 − 0.5) (sec) before movement onset. Terminal posture phase started from time (t 1 + 4) (sec), and steady state posture angle was obtained in a window of 1 (sec). The onset of muscle firing was calculated as the time at which EMG amplitude exceeded static EMG for three times of standard deviation (3SD) of it, and maintained higher than it for at least 25 ms (Hodges and Bui, 1996;Takatoku and Fujiwara, 2010), the offset of muscle firing was calculated as the time EMG decreased lower than 3 SD below static EMG. According to this threshold criterion, the duration of muscle firing was obtained, and a low passed EMG at the cut off frequency of 6 Hz (Steele, 1994) was used to detect the amplitude of EMG. Off-line Digital Butterworth filters were used in this study with forward and reverse passes to avoid phase shift. The static and dynamic features of movements and EMGs were used to tune parameters of descending commands for corresponding posture and movement.

Specification of Central Descending Commands
To fit the behaviors of the model in Section Corticospinal Virtual Arm (CS-VA) Model to the experimental data in Section Human Reach and Hold Experiment, the parameters of the model are fixed throughout the simulation. Only parameters of central descending commands are adjusted to capture the realistic feature of human movements and postures, as described in the following sub-sections.
Determining Static Command set (α s , γ s ) for Posture Posture was realized by adjusting static α and γ commands in this model. Earlier results (Lestienne et al., 1981) indicated that terminal posture could be coded as the ratio of activation levels between antagonistic muscles of a joint, thus α s of Bsh and Tlh in the simulation was set according to the relationship between ratio of static EMG (Bsh/Tlh) (RS) and static elbow angle (θ el ) of human experiment (Equation 2). The initial and terminal α s of Bsh and Tlh were then decided by corresponding RS at angle θ el from the experiment. The value of α s was adjusted in the model to keep a low activation on muscles during postures. The transition of α s between initial and terminal postures was assumed as a ramped changing pattern. To ensure a single joint movement in the simulation, the shoulder was fixed through adding a constant high co-activation level on α s of PC and DP.
γ s is related to centrally planned kinematics. In this study, we adopt an experimentally approved criterion for the central planned angle trajectory, the minimal-jerk criterion (Hogan, 1984;Flash and Hogan, 1985). The objective function is given in Equation (3a), and the minimal-jerk trajectory (MJT) is obtained by integration from t 0 and t 1 , which are times of movement onset and end. The outcome is a smooth angle trajectory in Equation (3b), as follows.
Determining Dynamic Command Set (α d , γ d ) for Movement According Taylor's recording (Taylor et al., 2000(Taylor et al., , 2004(Taylor et al., , 2006, γ d efferent had a constant baseline of firing rate during both posture and movement, and its firing rate sensitively changed with onset of muscle lengthening, thus presented a leading phase before movement. In the previous analysis of involuntary oscillatory movements (Hao et al., 2013), it was hypothesis that γ d represented joint acceleration, and was used to steer activations of antagonistic muscles for joint acceleration and deceleration. This reproduced all features of involuntary oscillatory movements such as Parkinsionian tremor. Thus in the present study, we also adopted this hypothesis, and determined the γ d command as the acceleration of the MJT as follows: In which θ el (t) was adopted from Equation (3b). The constant bias was set at 0.5 (normalized) to give a balanced (or mirror) inhibition on antagonistic PNs (see Figure 2). Parameter k was chosen so that the value of γ d was within 0 and 1. Movement was performed through integrating dynamic α and γ commands at the PN network. It has been proposed that, the motor system modulated the pulse amplitude and duration of muscle activations to produce movement with different velocities and ranges Gottlieb et al., 1989), and the pulse strategy has successfully generated scaled movements (Lan, 1997;Lan et al., 2005). Thus, a pair of pulses was utilized on α d to act as agonistic acceleration (AG1) and antagonistic deceleration (ANT) phase respectively. The pulse waveform was given in Equations (6a) and (6b): where amp AG1 and amp ANT were pulse amplitudes of AG1 and ANT (for extension movement, AG1 is on Tlh and ANT is on Bsh respectively), and their values were based on the ratio of peak EMG (ANT/AG1) (RP) obtained in human movement. pw AG1 and pw ANT were the pulse widths for AG1 and ANT, which were set as the durations of experimental EMGs according to the 3SD threshold criterion in Section Analysis of Kinematics and EMG Data. t AG1 andt ANT were the starting times for AG1 phase and ANT phase, t sim was the terminal time of simulation. Therefore, by adjusting the amplitude and width of the bi-phasic pulses, movements with different ranges and velocities could be achieved by the model.

Determining the Stabilization Pulse
A third pulse with a ramped decreasing form was necessary to stabilize the joint after movement. It was implemented in α s of agonist muscle (second agonistic muscle burst, AG2), since tri-phasic patterned EMG is found a common feature in fast reaching movement (Ghez and Martin, 1982;Flanders et al., 1994;Berardelli et al., 1996), and AG2 has been widely accepted as a phase to help stabilize the elbow at terminal posture after movement (Hannaford and Stark, 1985;Takatoku and Fujiwara, 2010). The amplitude and duration of AG2 were adjusted in the simulation to stabilize the joint after movement.

Simulation with Intact Feedforward and Feedback Controls
Normal condition with intact feedforward and feedback control was simulated to match typical trials in human experiment. For each subject, different postures were simulated and compared to the experimental relationship between RS and steady state angle. Stiffness control by varying co-activation level of antagonist muscles at different postures was demonstrated. To determine joint stiffness, a force perturbation was added directly on single joint muscle BS, and the stiffness was calculated as the ratio of changes in elbow torque with respect to joint angle . For each subject, the extension movement from 60 • to 30 • (Movement 1) was specifically simulated and matched to experimental data. More simulations were performed to fit experimental movements of larger range and in opposite direction for Subject 1, Simulation of 11 s was usually performed.
After the system had converged at a steady initial state, a random, signal dependent noise (Jones et al., 2002) was added on U m , the elbow was driven at time 6 (sec). Reflex gains used in this study were set within the range that kept the system stable ( Table 1) . All α and γ commands used in simulation were normalized to values between 0 and 1. Nominal parameter values of spinal circuits were adopted from Hao et al. (2013), and were listed in Table 1.

Simulation with Abnormal Feedforward and Feedback Controls
In this study, we examined the effects of altering model structure on movement and posture by deafferenting and deefferenting the model. Deafferented conditions were modeled by assigning the related gains of Ia, and Ib to zero. Deefferented conditions of the model were obtained by setting one of the descending commands to zero at a time, while keeping others intact. The effect of abnormal ratio (Bsh/Tlh) of α s on terminal postures was also examined by assigning static ratio out of normal ranges of experiments. Note that γ d was always maintained a bias of 0.5 in the analysis of deefferentation study. The abnormal conditions were applied to the model 0.1 (sec) before movement initiation in simulations. The difference between experimental and simulated movements was used to quantify the effects of abnormal control of movement and posture. Angular errors for movements and terminal steady state posture were evaluated as in Equations (7a) and (7b): θ el(P sim (j)) −θ el(P exp (j)) (7a) where Error P and Error M gave the averaged differences in terminal steady state posture and movement between simulated and experiment respectively (Figure 3), W and V were the number of resampled data points during terminal steady state posture and movement phases.

Features in Kinematics and EMGs of Human Posture and Movement
All subjects presented a common pattern in movement and posture at the elbow joint. The joint angle was stabilized after a slight overshoot, and the velocity displayed a bell-shaped profile. Tri-phasic EMG was observed in fast movements, with AG1 firing at accelerating phase, ANT at deceleration phase. The AG2 of agonist muscle appeared during movement offset and lasted until the joint was stabilized at terminal posture. When the elbow was stabilized after movement, the EMG of Tlh and Bsh were maintained at different steady state levels, which varied with elbow angle. An example of the relations between EMG ratio and joint angle during posture and movement was presented in Figure 4. The EMG of Tlh and Bsh were found to vary reciprocally with the elbow angle from 0 • to 110 • . The ratio of static EMG (Bsh/Tlh) (RS) was well fitted against the elbow angle in a linear manner ( Figure 4A, P < 0.01). During movement although  the peak velocity increased with movement range and speed, The relationship between peak velocity and ratio of peak EMG (ANT/AG1) (RP) was almost flat (P > 0.05), as shown in Figure 4B, and the ratio of peak EMG had an average value of 0.67 (±0.29). These features were present in all subjects. A summary of experimental relation between the ratio of static EMGs and elbow angle for all 6 subjects was given in Table 2. This indicates that the central module for posture control purposefully varies the static levels of antagonistic muscle activation in such particular linear manner to compensate for the change in moment arms of muscles at different joint angles. This experimental linear relationship is sufficient to guide the specification of α s descending commands in the model for posture maintenance.

Simulated Posture Control
Posture control was then simulated by tuning α s inputs to Bsh and Tlh for each subject, and by setting γ s according to the static angle θ el . Using the linear relation to guide the specification of α s commands for Bsh and Tlh muscles in the model, a comparable relationship was obtained in simulated postures for all subjects. The results in Table 2 showed that the simulated linear equations for each subject matched to the experimental relationship closely (P < 0.05). Overall an average linear relationship existed in experiment and simulated groups for all six subjects (P < 0.05), and they had a similar slope and intercept. This indicates that posture control can be achieved by tuning α s , as well as setting γ s inputs to relevant antagonistic muscles.
Thus, the posture module of motor system can control joint stiffness while maintaining the same joint posture. Figure 5 demonstrated such a strategy by increasing the co-activation levels of antagonistic muscles for Bsh and Tlh, while keeping the γ s commands at the values corresponding to elbow angles of 30 • , 50 • , 65 • respectively. Figure 5A showed that the simulated posture angles against the experimental linear relation. Since multiple activation levels of muscles can achieve the same elbow posture, joint stiffness can be controlled to desired levels for different motor tasks. We calculated joint stiffness by applying force perturbations on BS and measured the ratio between change of joint torque and angle. Figure 5B showed that joint stiffness increased with the muscle activation, while the same posture was maintained. The results indicate that the central motor system can control joint stiffness and posture independently by tuning the levels of α s commands with programmed γ s control for postures.

Simulated Movement Control
The extension movement from 60 • to 30 • (Movement 1) was simulated for six subjects individually. An extension movement from 90 • to 45 • (Movement 2), a flexion movement from 30 • to 60 • (Movement 3) were simulated for Subject 1. Figures 6A-C depict the descending commands tuned for Movement 1, 2, and 3 respectively. During commands assignment, γ s and γ d commands were calculated from the MJT of target trajectory (Equations 4 and 5). And the initial and terminal postures were maintained by setting the ratio of α s of Bsh and Tlh as experimental values. Thus, movements between postures were simulated with tuning of the pulse height and width of the α d command. The amplitudes and pulse widths were adopted from experimental EMG (Equation 6). The third pulse in the α s of Tlh acted to help stabilize the elbow after fast movement. The posture commands were also specified for the other four muscles in the model to keep the simulation running properly.
The simulated angle, velocity and muscle activation of Movement 1, 2, and 3 were compared with the experimental movements (Figures 7A-C). The initial and terminal angles could finely match the experimental postures, and the three movements displayed the signature bell-shaped velocity profile, and tri-phasic firing pattern of Bsh and Tlh. Results showed that the elbow could be maintained closely at terminal postures, and the errors during movements were relatively small compared to the range of movements.
Furthermore, all simulated movements matched well to those of experimental trajectories and EMGs in all subjects by tuning the α d pulse command. The errors of trajectory fitting for Movement 1 were calculated according to Equation (7) and listed in Table 3. It is demonstrated that the movement module of the central motor system could control movements of a range of angles and velocities by tuning the α d pulse command while coordinating posture commands.

Simulation with Abnormal Afferent and Efferent Controls
Movements under abnormal afferent and efferent controls for Movement 1 were simulated to assess the individual contribution of each descending command to movement and posture. Figure 8 illustrated the various effects for Subject 1. Figure 8A showed that deafferentation after Ia(−) was removed overshot the target of movement with a smaller steady state angle; while removing Ia(+) and Ia-PN undershot the target of movement, note that removing Ia(+) resulted in a larger steady state angle, while Ia-PN had no effect on final posture. When all Ia afferents to αMN were removed in deafferented state, the joint undershot the target of movement, and stabilized gradually to a slightly larger steady state posture.
The results of abnormal efferent commands were shown in Figure 8B. When the ratio of α s was doubled, the movement and steady state posture were slightly affected compared to those in  normal condition. When γ s was removed, the joint undershot the target of movement similar to that of deafferentation. After α d was removed or γ d kept a constant inhibition on PNs, the elbow couldn't make the fast movement as in normal condition. The slow approach to the target posture was brought about by the action of spinal reflexes. Thus, feedforward control of descending commands (α d and γ d ) is absolutely essential for fast movements.
This general pattern of behavior in all subjects is summarized in the phase diagram of errors in Figure 9. As illustrated in Figure 9A, removing Ia(+), Ia(−) and all Ia afferents generally impacted the terminal postures, these results were consistent with the behaviors observed in humans and animals (Polit and Bizzi, 1979;Gentilucci et al., 1994;Gordon et al., 1995). But removing Ia-PN seemed to reduce errors in both movement and posture. Static and dynamic descending commands showed distinct impacts on posture and movement ( Figure 9B). The deviated ratios of α s for Bsh and Tlh from normal values yielded a wider departure from the target posture, which was different for individual subjects. The absence of γ s resulted in an error primarily in posture. The importance of dynamic commands for movement was clearly seen from the large Error M in Figure 9B. This is consistent with the observation in cats after the PNs innervating the upper extremity muscles were removed (Alstermark et al., 1981;Alstermark and Isa, 2012).

Discussion
An important issue in motor control has been to differentiate the role of the γ motor system from that of the α motor system (Stein, 1974;Houk and Rymer, 1981;Bizzi et al., 1991). Although increasing evidence revealed α-γ co-activation during movement and posture control (Vallbo, 1971;Taylor et al., 2006;Prochazka and Ellaway, 2012), the dominant effects of the α motor system tended to diminish the function significance of the γ motor system. Using a computational virtual arm (VA) model , Lan and He (2012) re-interpreted a set of experiment data from Stein et al. (2004), Cordo et al. (2002), Taylor et al. (2006), and suggested that γ s may encode centrally planned information of joint angle and reinforce the planned joint angle though regulating spindle sensitivity. By coupling the VA model with a corticospinal (CS) network of PN in analyzing involuntary oscillatory movements (Hao et al., 2013), we further hypothesized that γ d represents centrally planned joint acceleration. In this paper, the combined CS-VA model was used to implicate the necessity of α-γ coordination in movement and posture control. Only the set of α-γ descending commands were adjusted to fit human movement data. Close match of model behaviors to those observed in human experiments as demonstrated in the results (Figure 7, Table 3) established the validity of the model, therefore, providing a neurophysiologically realistic, multi-scale computational model to evaluate the contribution of various components of descending control in sensorimotor functions. In particular, this model will also be valuable to understand sensorimotor dysfunctions (Hao et al., 2013) and to design novel rehabilitation strategies for motor relearning (Zhuang et al., 2015).
A distinct feature of this computational model is the division of sensorimotor control into movement module and posture module by the dynamic and static descending command sets (Bizzi et al., 2008;Diedrichsen and Classen, 2012). Results of this study indicate that it is possible to coordinate the two sets of descending α-γ commands to achieve accurate control of movement dynamics and stable maintenance of final posture. Methods adopted in this study to determine descending commands are novel, in that these descending commands are specified according to proven rules of central planning for kinematics (Hogan, 1984;Flash and Hogan, 1985) and EMG data collected in human subjects performing reach-and-hold tasks. Assumptions are also adopted from previous analytical studies regarding γ s encoding of joint angle (Lan and He, 2012), and γ d representation of joint acceleration of planned movement trajectory (Taylor et al., 2006;Hao et al., 2013). Results in Figure 5 and Table 2 indicate that α s and α d commands can be tuned based on EMG signals of human data at steady state and during movement. Figure 5 shows that tuning α s based on the linear ratio of antagonistic muscles in Figure 4 is necessary to achieve independent control of joint posture and stiffness, which is an important aspect of regulation of motor functions by the central motor system (Mussa-Ivaldi et al., 1985). For movement control, the α d command is integrated with γ d command at the PN to distribute properly the activation to flexor or extensor acting at the joint, and its pulse amplitude and width can be tuned according to the speed of movement and duration of EMG bursts. We illustrated that adjusting these descending commands can fit reach-and-hold movements for a range of amplitude and direction in different subjects. The ability of the model to fit experiment movements suggests that the computational model captures the neural mechanism of corticospinal computation, as well as the modular nature of organization and coordination of descending α-γ commands by the central motor control system Scheidt and Ghez, 2007;Scheidt et al., 2011;Poston et al., 2013).
The model predicted the contribution of descending α-γ commands to movement and posture by deefferenting the CS-VA model in simulation. With abnormal α s command out of the range of experimental ratio (Figures 8B, 9B), the terminal angle deviated from its targeted position. The large variation shown at the terminal position suggests that accurate γ s command is essential for posture maintenance. The tardy movements under abnormal α d and γ d (Figure 9B) indicated that fast reaching movement must be performed with proper coordination of α d and γ d commands. This is consistent with the observation that cats were unable to carry out skilled reaching movement without PN (Alstermark et al., 1981;Alstermark and Isa, 2012).
Proprioceptive afferents from muscle spindle are important for motor learning (Jeannerod, 1988;Schmidt and Lee, 2011), but are not found indispensible for control of movements, since deafferentation in human patients and animals did not entirely disable their movement execution. Early study in deafferented patients indicated obvious motor dysfunction only with larger errors in terminal position acquisition (Polit and Bizzi, 1979;Gentilucci et al., 1994;Gordon et al., 1995) and lower joint stiffness during posture maintenance (Bizzi et al., 1984), but changes in movement kinematic and muscle firing pattern were not obvious (Taub et al., 1966(Taub et al., , 1975Vaughan et al., 1970;Rothwell et al., 1982;Bizzi et al., 1984;Gordon et al., 1995). This implied that proprioceptive afferent played an important role in posture maintenance and contributed to fine control of movement.
This functional role of proprioceptive afferents is reiterated by deafferenting the CS-VA model (Figures 8A, 9A). It was shown that, when Ia(+) and Ia-PN was removed respectively, the movement slowed down, and the terminal posture shifted, and after Ia(−) was removed, the movement speeded up, and the terminal posture targeted at a lower angle. This confirmed the positive feedback of Ia(+) and Ia-PN, and inhibition role of Ia(−). However, despite the difference of peak velocity, these movements presented similar bell-shaped velocity profile, only settled at different terminal angles. Deafferented simulation showed slowed movement and increased errors in terminal angle. This was similar to the behaviors observed in deafferented primates (Polit and Bizzi, 1979), which showed that the deafferented primate after intensive training was able to control pointing movement with normal like kinematics and muscle activation, but inaccurate positions.

Conclusion
A corticospinal computational model based on the modular organization of movement and posture was validated in this study by fitting the model to experimental human data. Analysis of simulated movement and posture with intact and altered model structures demonstrated that it is necessary to coordinate the set of α-γ descending commands in order to achieve effective control of accurate movement dynamics and stable postures. Results suggest that the central commands of posture module are mediated via a mono-synaptic corticospinal pathway, while those of movement module are transmitted to spinal motoneruons through a multi-synaptic corticospinal pathway involving the propriospinal neurons (PN). The PN network plays the pivotal role to integrate the α d and γ d commands for movement generation. The model is able to capture many essential aspects of motor behaviors, such as independent regulation of joint angle and stiffness and the signature temporal pattern of EMGs, by simply tuning the α s and α d commands. This study suggests a plausible neural computational mechanism for the central motor system to control movement and posture. The model will be useful as a complementary tool to understand neural control of movements, as well as a valuable platform to aid design of novel rehabilitation strategy for motor disabilities.
Author Contributions SL performed model simulation, analyzed experimental data, prepared figures and tables and drafted the manuscript; XH and MH contributed to set up the computational model; JM and CZ contributed in the analysis of human experiment data. CZ also carried out part of simulation work. CN offered intellectual suggestions and edited the manuscript. NL conceived the computational approach, designed the human subject experiment, proposed analytical method and edited the final version of the manuscript.